『Graphics』绘制二维数据集的置信椭圆

2020-12-31

字数统计: 1.8k | 阅读时长≈ 7 分钟

如何绘制二维正态分布数据的误差椭圆，又名置信椭圆。

在二维或者更高维的空间里，数据的聚类往往需要添加一个“置信区间”。仿照一维空间的数据，置信区间往往相对于点估计而来的，在统计学中，一个概率样本的置信区间（Confidence interval）是对这个样本的某个总体参数的区间估计。置信区间展现的是这个参数的真实值有一定概率落在测量结果的周围的程度，其给出的是被测量参数的测量值的可信程度，一般我们用95%置信区间来表示。
那么二维空间的数据也是如此，二维空间的置信区间往往利用置信椭圆来描述，展现真实值的可信程度，一般我们用95%置信椭圆来表示。
为什么是置信椭圆呢，往往是因为我们的而数据分布形似椭圆

置信椭圆绘制原理

一、数据中心在原点且轴对齐型

先看最简单的，轴对齐型的，比方说椭圆长轴平行于x轴

假设我有两列随机变量 x 和 y
对于该二维数据，x和y，我们不妨计算下x和y的协方差矩阵：

由于我们的数据，形似椭圆，并且椭圆长轴平行于x轴，那么协方差矩阵中x轴方向的方差为8.4213，y轴方向的方差为0.9387。由于x方向与y方向正交，所以x和y这两个随机变量的协方差为0，也可以理解为相关性为0。

由于我们的数据点分布形似椭圆，所以我们定义置信区间就把它定义为椭圆，我们知道椭圆的标准方程为：
$$
\left(\frac{x}{a}\right)^{2}+\left(\frac{y}{b}\right)^{2}=1
$$
那么我们不妨将椭圆的半长轴a定义为σx，半短轴b定义为σy，那么有：
$$
\left(\frac{x}{\sigma_{x}}\right)^{2}+\left(\frac{y}{\sigma_{y}}\right)^{2}=s
$$
σ_x与σ_y分别代表x和y的标准差，s与置信度有关。

显然，该椭圆的数据中心为(0,0)。
其中s定义椭圆的规模，可以是任意的数（例如，s=1）。现在的问题是如何选择s，使得所得到的椭圆规模代表我们所选择的置信水平（例如，95％的置信水平对应于s =5.991）。(不太明白这个椭圆的规模是什么意思)
这个s=5.991可以通过卡方分布（随机变量平方和）进行计算（假设随机变量x和y服从高斯分布），并且自由度为2.
$$
\left(\frac{x}{\sigma_{x}}\right)^{2}+\left(\frac{y}{\sigma_{y}}\right)^{2}=s
$$

联想到卡方分布，上式左侧可以看作为随机变量的平方和，通过查询卡方表，有：
$$
P(s<5.991)=1-0.05=0.95
$$
我们知道当 s = 5.991 时，P(s < 5.991) = 0.95，也就是说，我们要寻找一个在椭圆上的点，使得 s = 5.991，那么将 s = 5.991 带入椭圆方程：
$$
\left(\frac{x}{\sigma_{x}}\right)^{2}+\left(\frac{y}{\sigma_{y}}\right)^{2}<5.991
$$
这个等式表明数据点小于5.991的概率为95%（在椭圆内表示数据点小于5.991）。
那么，s=5.991即为95%置信区间，表示有一组数据，其x方向方差为σx，y方向方差为σy时，如果s=5.991，那么95%的数据点在该椭圆内。
那么置信椭圆的长轴2a：
$$
2 \sigma_{x} \sqrt{5.991}
$$
置信椭圆的短轴2b：
$$
2 \sigma_{y} \sqrt{5.991}
$$

二、数据中心在原点且任意方向的置信椭圆

在数据是相关的情况下，例如存在协方差，所产生的误差椭圆不会是轴对齐的。在这种情况下，如果我们暂时定义一个新的坐标系，使得所述椭圆变为轴对齐，然后旋转所产生的椭圆，那么上面的结论依然有效。

换句话说，之前我们计算平行于x轴和y轴的方差，现在我们计算平行于置信椭圆长轴和短轴的方差，需要计算的方差方向由图1粉红和绿色箭头表示出来。

首先计算随机变量x和y的协方差矩阵，根据协方差矩阵计算该矩阵的特征向量（上图绿色为椭圆短轴方向的特征向量；粉红色为椭圆长轴方向的特征向量），我们定义椭圆长轴方向的特征向量为v1，椭圆短轴方向的特征向量为v2；而特征向量大小为特征值，我们定义椭圆长轴方向的特征值为λ1，椭圆短轴方向的特征值为λ2。
此时椭圆的长轴为：
$$
2 \sqrt{5.991\lambda_1}
$$
椭圆的短轴为：
$$
2 \sqrt{5.991\lambda_2}
$$
并且定义椭圆长轴与x轴正方向夹角为α：(为了获得椭圆的方向，我们简单地计算最大特征向量的角度)

$$ \alpha=\arctan \frac{\mathbf{v}_{1}(y)}{\mathbf{v}_{1}(x)} $$ 那么，我们就可以根据角度来对置信椭圆进行“旋转”。

三、任意位置置信椭圆

有了上面的介绍，那么任意位置的椭圆无非是满足数据中心不在原点，且椭圆长轴与x轴正方向夹角为任意角α。并以此来建立椭圆方程，画出置信椭圆（可见数据中心化可以简便很多计算）

python 绘制置信椭圆

# -*- encoding: utf-8 -*-
"""
@software: PyCharm
@file : 112.py
@time : 2020/12/31 
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
import matplotlib.transforms as transforms
import pandas as pd


def confidence_ellipse(x, y, ax, n_std=3.0, facecolor='none', **kwargs):
    if x.size != y.size:
        raise ValueError("x and y must be the same size")

    cov = np.cov(x, y)
    lambda_, v = np.linalg.eig(cov)  # 计算特征值和特征向量
    angle = np.rad2deg(np.arccos(v[0, 0]))  # 计算旋转角度

    ell_radius_x = np.sqrt(n_std*lambda_[0])
    ell_radius_y = np.sqrt(n_std*lambda_[1])

    # 默认椭圆中心为原点
    mean_x,mean_y = np.mean(x), np.mean(y)
    ellipse = Ellipse((mean_x,mean_y),
                      width=ell_radius_x * 2,
                      height=ell_radius_y * 2,
                      facecolor=facecolor,
                      angle=angle,
                      **kwargs)

    return ax.add_patch(ellipse)

if __name__ == '__main__':
    data = pd.read_excel('data.xlsx', index_col=None)
    x = data['x']
    y = data['y']

    # 卡方值， 95％的置信水平对应于s =5.991
    ss = 5.991

    fig, ax_kwargs = plt.subplots(figsize=(6, 6))

    # ax_kwargs.axvline(c='grey', lw=1)
    # ax_kwargs.axhline(c='grey', lw=1)


    confidence_ellipse(x, y, ax_kwargs,n_std=ss,
                       alpha=0.9,  edgecolor='red', zorder=0)
    ax_kwargs.scatter(x, y, s=2)
    ax_kwargs.scatter(np.mean(x),np.mean(y), c='red', s=9)
    plt.show()