使用python绘制置信椭圆

Fr0_o

于 2024-03-20 17:21:55 发布

阅读量227

点赞数 3

文章标签： python 开发语言

本文链接：https://blog.csdn.net/Fr0_o/article/details/136882621

版权

本文介绍了一个在Ubuntu22.04环境中，使用Python3.10.13版本的numpy1.26.3和matplotlib3.7.2库实现的置信椭圆绘制函数，用于描述多元正态分布的置信区间。函数通过计算样本协方差和特征值生成椭圆图形。

摘要由CSDN通过智能技术生成

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