如何解决绘制二维高斯椭圆轮廓
我想画一个与水平集(轮廓)对应的椭圆
遵循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])
但是这显然行不通。
解决方法
-
您应该在伽马中对 sx2 和 sy2 进行平方。
-
gamma 应该是平方根。
-
将得到的椭圆乘以 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()
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。