Drawing Ellipse Contour of 2D Gaussian

Question

Suppose I have a 2D Gaussian with pdf

I want to draw an ellipse corresponding to the level-set (contour)

Following here I know that I can replace the precision matrix with its eigendecomposition to obtain where gamma is

Then to find coordinates of the points on the ellipse I would have to do

I tried plotting this but it is not working.

Plotting the Contours

from scipy.stats import multivariate_normal
import numpy as np
from numpy.linalg import eigh
import math
import matplotlib.pyplot as plt

# Target distribution
sx2 = 1.0
sy2 = 2.0
rho = 0.6
Sigma = np.array([[sx2, rho*math.sqrt(sx2)*math.sqrt(sy2)], [rho*math.sqrt(sx2)*math.sqrt(sy2), sy2]])
target = multivariate_normal(mean=np.zeros(2), cov=Sigma)
# Two different contours
xy = target.rvs()
xy2 = target.rvs() 
# Values where to plot the density
x, y = np.mgrid[-2:2:0.1, -2:2:0.1]
zz = target.pdf(np.dstack((x, y)))
fig, ax = plt.subplots()
ax.contour(x,y, zz, levels=np.sort([target.pdf(xy), target.pdf(xy2)]))
ax.set_aspect("equal")
plt.show()

The code above shows the contour

Plotting the Ellipse

# Find gamma and perform eigendecomposition
gamma = math.log(1 / (4*(np.pi**2)*sx2*sy2*(1 - rho**2)*(target.pdf(xy)**2)))
eigenvalues, P = eigh(np.linalg.inv(Sigma))
# Compute u and v as per link using thetas from 0 to 2pi
thetas = np.linspace(0, 2*np.pi, 100)
uv = (gamma / np.sqrt(eigenvalues)) * np.hstack((np.cos(thetas).reshape(-1,1), np.sin(thetas).reshape(-1, 1)))
# Plot
plt.scatter(uv[:, 0], uv[:, 1])

However this clearly doesn't work.

Answer 1

1. You should square sx2 and sy2 in gamma.

2. gamma should be square rooted.

3. Multiply the resulting ellipse by P^-1 to get points in the original coordinate system. That's mentioned in the linked post. You have to convert back to the original coordinate system. I don't know actually how to code this, or if it actually works, so I leave the coding to you.

gamma = math.log(1 / (4*(np.pi**2)*(sx2**2)*(sy2**2)*(1 - rho**2)*(target.pdf(xy)**2)))
eigenvalues, P = eigh(np.linalg.inv(Sigma))
# Compute u and v as per link using thetas from 0 to 2pi
thetas = np.linspace(0, 2*np.pi, 100)
uv = (np.sqrt(gamma) / np.sqrt(eigenvalues)) * np.hstack((np.cos(thetas).reshape(-1,1), np.sin(thetas).reshape(-1, 1)))


orig_coord=np.linalg.inv(P) * uv #I don't how to code this in python

plt.scatter(orig_coord[:,0], orig_coord[:,1])
plt.show()

My attempt at coding it:

gamma = math.log(1 / (4*(np.pi**2)*(sx2**2)*(sy2**2)*(1 - rho**2)*(target.pdf(xy)**2)))
eigenvalues, P = eigh(np.linalg.inv(Sigma))
# Compute u and v as per link using thetas from 0 to 2pi
thetas = np.linspace(0, 2*np.pi, 100)
uv = (np.sqrt(gamma) / np.sqrt(eigenvalues)) * np.hstack((np.cos(thetas).reshape(-1,1), np.sin(thetas).reshape(-1, 1)))


orig_coord=np.zeros((100,2))
for i in range(len(uv)):
    orig_coord[i,0]=np.matmul(np.linalg.inv(P), uv[i,:])[0]
    orig_coord[i,1]=np.matmul(np.linalg.inv(P), uv[i,:])[1]

# Plot
plt.scatter(orig_coord[:, 0], orig_coord[:, 1])

gamma1 = math.log(1 / (4*(np.pi**2)*(sx2**2)*(sy2**2)*(1 - rho**2)*(target.pdf(xy2)**2)))
uv1 = (np.sqrt(gamma1) / np.sqrt(eigenvalues)) * np.hstack((np.cos(thetas).reshape(-1,1), np.sin(thetas).reshape(-1, 1)))
orig_coord1=np.zeros((100,2))
for i in range(len(uv)):
    orig_coord1[i,0]=np.matmul(np.linalg.inv(P), uv1[i,:])[0]
    orig_coord1[i,1]=np.matmul(np.linalg.inv(P), uv1[i,:])[1]
plt.scatter(orig_coord1[:, 0], orig_coord1[:, 1])

plt.axis([-2,2,-2,2])

plt.show()

Sometimes the plots don't work and you get the error invalid sqrt, but when it works it looks fine.