How to get bivariate normal probability distribution with specified standard deviation and angle in python

Question

I want to get a bivariate normal probability distribution(hereafter referred to as BNPD) from ellipse parameters. As shown in the illustraion below, the length of the major axis of the ellipse will be taken as 3 times the standard deviation of BNPD, the same to the minor axis, and the rotation angle of ellipse is taken as the rotation angle of BNPD.

This is the python test code, not success yet:


import numpy as np
import cv2 as cv

# initialize ellispe parameters
size = 500
ctx = 280
cty = 220
width_radius = 40
height_radius = 120
theta = 45 

# draw a reference line and ellipse
bg = np.ones((size,size, 3))
cv.line(bg,(0,int(size/2)),(size,int(size/2)),(0,0,0),1)
cv.line(bg,(int(size/2),0),(int(size/2),size),(0,0,0),1)
cv.ellipse(bg,(ctx,cty),(width_radius,height_radius),theta,0,360,(0,0,255),2)
cv.imshow("result", bg)
cv.waitKey(0)

# convert ellipse parameters to bivariate normal probability distribution parameters
mean_x = ctx - size//2
mean_y = cty - size//2
sigma_x = width_radius/(3*np.cos(theta)) # I doubt here
sigma_y = height_radius/(3*np.sin(theta))

gaussian_x = lambda x: np.exp(-(x-mean_x)**2/(2*(sigma_x**2)))/(np.sqrt(2*np.pi) * sigma_x)
gaussian_y = lambda y :np.exp(-(y-mean_y)**2/(2*(sigma_y**2)))/(np.sqrt(2*np.pi) * sigma_y)
X = np.linspace(-size/2, size/2, size)
Y = np.linspace(-size/2, size/2, size)

result_x = gaussian_x(X).reshape(-1,1)
result_y = gaussian_y(Y).reshape(1,-1)
result_xy = np.dot(result_x, result_y)

# visualize bivariate normal probability distribution
import matplotlib.pyplot as plt
fig, axes = plt.subplots()
axes.imshow(result_xy)
plt.show()

Run above code, it shows: As you can see, the center and angle of ellipse is not correspond to the mean and angle of BNPD.

-------------------------update-----------------------------

As @Matt Pitkin suggest, I try out the scipy.stats.multivariate_normal with deduced covariance matrix, but from the plotted graph, it seems the rotation angle is slightly smaller than 45°. The code is:

phi = theta
semiMajorAxis = width_radius
semiMinorAxis = height_radius

varX1 = semiMajorAxis**2 * np.cos(phi)**2 + semiMinorAxis**2 * np.sin(phi)**2
varX2 = semiMajorAxis**2 * np.sin(phi)**2 + semiMinorAxis**2 * np.cos(phi)**2
cov12 = (semiMajorAxis**2 - semiMinorAxis**2) * np.sin(phi) * np.cos(phi) 
mean = [mean_x, mean_y]
covariance = [[varX1, cov12], 
              [cov12,varX2]]
from scipy.stats import multivariate_normal
X, Y = np.meshgrid(X, Y)
xy = np.stack((X, Y), -1)
fxy = multivariate_normal.pdf(xy, mean=mean, cov=covariance)
fig, axes = plt.subplots()
axes.imshow(fxy)
plt.show()

Are you happy to use the [`scipy.stats.multivariate_normal`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.multivariate_normal.html#scipy.stats.multivariate_normal) function? If so, I'd suggest looking at https://stackoverflow.com/questions/41807958/convert-position-confidence-ellipse-to-covariance-matrix to calculate your covariance matrix and then input that into `multivariate_normal` as in the example shown in the docs. — Matt Pitkin, Feb 22 '23 at 12:48
@MattPitkin Hi, it is OK to use `scipy.stats.multivariate_normal` function, I updated my question, there is still some problem. — Wade Wang, Feb 22 '23 at 16:10

Wade Wang · Answer 1 · 2023-02-23T06:08:45.243

Thanks's to @Matt Pitkin's guidance in the comment, the final solution is :

phi =  theta*np.pi/180
semiMajorAxis = width_radius
semiMinorAxis = height_radius

varX1 = semiMajorAxis**2 * np.cos(phi)**2 + semiMinorAxis**2 * np.sin(phi)**2
varX2 = semiMajorAxis**2 * np.sin(phi)**2 + semiMinorAxis**2 * np.cos(phi)**2
cov12 = (semiMajorAxis**2 - semiMinorAxis**2) * np.sin(phi) * np.cos(phi) 
mean = [mean_x, mean_y]
covariance = [[varX1, cov12], 
              [cov12,varX2]]
from scipy.stats import multivariate_normal
X, Y = np.meshgrid(X, Y)
xy = np.stack((X, Y), -1)
fxy = multivariate_normal.pdf(xy, mean=mean, cov=covariance)
fig, axes = plt.subplots()
axes.imshow(fxy)
plt.show()

This is the newest plot result:

How to get bivariate normal probability distribution with specified standard deviation and angle in python

1 Answers1