绘制二维高斯椭圆轮廓

问题描述

假设我有一个带有 pdf 的二维高斯

我想画一个与水平集（轮廓）对应的椭圆

遵循here我知道我可以用它的特征分解替换精度矩阵来获得

其中伽马是

然后要找到椭圆上点的坐标，我必须这样做

我尝试绘制此图，但它不起作用。

绘制轮廓

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()

上面的代码显示了轮廓

绘制椭圆

# 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])

但是这显然行不通。

解决方法

1. 您应该在伽马中对 sx2 和 sy2 进行平方。

2. gamma 应该是平方根。

3. 将得到的椭圆乘以 P^-1 得到原始坐标系中的点。这是在链接的帖子中提到的。您必须转换回原始坐标系。我实际上不知道如何编码，或者它是否真的有效，所以我把编码留给你。

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()

我的编码尝试：

gamma = math.log(1 / (4*(np.pi**2)*(sx2**2)*(sy2**2)*(1 - rho**2)*(target.pdf(xy)**2)))
eigenvalues,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),:])[1]

# Plot
plt.scatter(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)))
orig_coord1=np.zeros((100,2))
for i in range(len(uv)):
    orig_coord1[i,uv1[i,:])[0]
    orig_coord1[i,:])[1]
plt.scatter(orig_coord1[:,orig_coord1[:,1])

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

plt.show()

有时绘图不起作用，您会收到错误无效的 sqrt，但是当它起作用时，它看起来很好。