python - 绘制二维高斯的椭圆轮廓
问题描述
假设我有一个带有 pdf 的 2D 高斯
我想绘制一个对应于水平集(轮廓)的椭圆
在这里我知道我可以用它的特征分解替换精度矩阵来获得 伽马在哪里
然后要找到椭圆上点的坐标,我必须这样做
我尝试绘制此图,但它不起作用。
绘制等高线
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])
然而,这显然是行不通的。
解决方案
您应该在 gamma 中平方 sx2 和 sy2。
伽玛应该是平方根。
将得到的椭圆乘以 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, 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()
推荐阅读
- sequence - Teradata 根据前几行创建日期范围窗口
- node.js - 我正在尝试登录迪士尼商店。我尝试了许多不同的选择器,但它们似乎都不起作用。有什么帮助或提示吗?
- java - Spring Boot:延迟获取不适用于 getById()
- sql - 多次使用同一列时,SQL Server 无法将文本转换为时间
- python - 姜戈 | 在同一个模板上多次使用同一个表单
- reactjs - 在 Emotion 11/Next js 10 应用程序中启用 css 道具的正确方法是什么
- amazon-web-services - 如何将图像从 React Native 上传到 AWS S3?错误 403
- python - 未找到带有参数 '('',)' 的 'admin-color-product-map-edit' 的反向操作。尝试了 1 种模式:['admin1/colorProductMap_edit/(?P
[0-9]+)$'] - java - > 无法隔离参数 com.android.build.gradle
- json - 是否可以使用 curl 设置 json 抓取,使用当前日期抓取并在脚本中更改接下来 5 天的日期?