Demonstration of plotting an error ellipse (probability ellipse) in Python.

plot_error_ellipse.py

from matplotlib import pyplot as plt
from matplotlib.patches import Ellipse
import numpy as np

fig, ax = plt.subplots(1)
ax.set_aspect('equal')

m = np.array([8, 12])
S = np.array([[16, np.sqrt(78)], [np.sqrt(78), 9]])
r = np.random.multivariate_normal(m, S, size=100) # random points to plot
plt.scatter(r[:, 0], r[:, 1], color="C0", s=2, zorder=1)

su2 = ((S[0, 0] + S[1, 1]) + np.sqrt((S[0, 0] - S[1, 1])**2 + 4*S[0, 1]**2)) / 2
sv2 = ((S[0, 0] + S[1, 1]) - np.sqrt((S[0, 0] - S[1, 1])**2 + 4*S[0, 1]**2)) / 2
slope1 = (su2 - S[0, 0]) / S[0, 1]
slope2 = (sv2 - S[0, 0]) / S[0, 1]
theta = np.arctan(slope1) * 180.0 / np.pi

# P = 0.500
ax.add_patch(Ellipse(m, 1.177*np.sqrt(su2)*2, 1.177*np.sqrt(sv2)*2, theta, edgecolor="C1", fill=False, zorder=0))
# P = 0.900
ax.add_patch(Ellipse(m, 2.146*np.sqrt(su2)*2, 2.146*np.sqrt(sv2)*2, theta, edgecolor="C2", fill=False, zorder=0))
# P = 0.950
ax.add_patch(Ellipse(m, 2.448*np.sqrt(su2)*2, 2.448*np.sqrt(sv2)*2, theta, edgecolor="C3", fill=False, zorder=0))

# Plot basis vector e1
intercept = m[1] - slope1 * m[0]
sx = m[0] - np.sqrt(su2)*3
sy = sx*slope1 + intercept
ex = m[0] + np.sqrt(su2)*3
ey = ex*slope1 + intercept
plt.quiver(sx, sy, ex-sx, ey-sy, angles='xy', scale_units='xy', scale=1, alpha=0.3, zorder=-1)

# Plot basis vector e2
intercept = m[1] - slope2 * m[0]
sy = m[1] - np.sqrt(sv2)*3
sx = (sy - intercept) / slope2
ey = m[1] + np.sqrt(sv2)*3
ex = (ey - intercept) / slope2
plt.quiver(sx, sy, ex-sx, ey-sy, angles='xy', scale_units='xy', scale=1, alpha=0.3, zorder=-1)

plt.xlim(-20, 30)
plt.ylim(-10, 25)
plt.xlabel("x")
plt.ylabel("y")
plt.savefig("ellipse.png")