环境
Ubuntu 22.04
Python 3.10.13
numpy 1.26.3
matplotlib 3.7.2
核心代码
def groupcir(x, y, level=0.95, df=2, **kwargs):
'''
置信椭圆:描述多元正态分布的置信区间。
椭圆的长轴和短轴表示方差,而椭圆的旋转角度则表示协方差。
本函数计算了p={level},自由度df={df}的置信椭圆
:type x: array like
:param x: x postion of samples
:type y: array like
:param y: y postion of samples
:type level: float
:param level: confidence level
:type df: int
:param df: degrees of freedom
:param kwargs:pass to mpl.patches.Ellipse
'''
critical_value = chi2.ppf(level, df=df)
x0 = np.mean(x) # 椭圆中心x坐标
y0 = np.mean(y) # 椭圆中心y坐标
cov = np.cov(x, y) # 样本协方差矩阵
lambda_, v = np.linalg.eig(cov) # 求特征值和特征向量
a = np.sqrt(lambda_[0]) * critical_value # 长轴
b = np.sqrt(lambda_[1]) * critical_value # 短轴
theta = np.arctan2(v[1, 0], v[0, 0]) * 180 / np.pi # 旋转角度
ellipse = mpl.patches.Ellipse( (x0, y0), a, b, angle=theta, **kwargs)
return ellipse
整体代码如下
import matplotlib as mpl
from scipy.stats import chi2
from matplotlib.lines import Line2D
import numpy as np
import matplotlib.pyplot as plt
def groupcir(x, y, level=0.95, df=2, **kwargs):
'''
置信椭圆:描述多元正态分布的置信区间。
椭圆的长轴和短轴表示方差,而椭圆的旋转角度则表示协方差。
本函数计算了p={level},自由度df={df}的置信椭圆
:type x: array like
:param x: x postion of samples
:type y: array like
:param y: y postion of samples
:type level: float
:param level: confidence level
:type df: int
:param df: degrees of freedom
:param kwargs:pass to mpl.patches.Ellipse
'''
critical_value = chi2.ppf(level, df=df)
x0 = np.mean(x) # 椭圆中心x坐标
y0 = np.mean(y) # 椭圆中心y坐标
cov = np.cov(x, y) # 样本协方差矩阵
lambda_, v = np.linalg.eig(cov) # 求特征值和特征向量
a = np.sqrt(lambda_[0]) * critical_value # 长轴
b = np.sqrt(lambda_[1]) * critical_value # 短轴
theta = np.arctan2(v[1, 0], v[0, 0]) * 180 / np.pi # 旋转角度
ellipse = mpl.patches.Ellipse( (x0, y0), a, b, angle=theta, **kwargs)
return ellipse
if __name__ == '__main__':
data = np.random.normal(2, 1, (10000, 2))
distances = np.abs(data[:, 1] - data[:, 0])
median_distance = np.median(distances)
y = np.where(distances > median_distance, 1, 0)
fig=plt.figure()
ax = fig.add_subplot(111)
data_0 = data[y==0]
data_1 = data[y==1]
plt.scatter(data_0[:, 0], data_0[:, 1], color='#90C549', s=5, alpha=1)
plt.scatter(data_1[:, 0], data_1[:, 1], color="#C04000", s=5,alpha=1)
ellipse = groupcir(data_0[:, 0], data_0[:, 1], level=0.8, fill=False, linewidth=2,linestyle='--',edgecolor='#268388')
ax.add_patch(ellipse)
ellipse = groupcir(data_1[:, 0],data_1[:, 1], level=0.8, fill=False, linewidth=2,linestyle='--',edgecolor="#800000")
ax.add_patch(ellipse)
plt.show()