【python】95%信頼楕円/確率楕円を描画する

　「ライブラリあるやろｗ」と思ったら、なくて顔面蒼白になった。

　しょうがないから調べて実装した。

理論的なもの

　ちゃんと数式を書いて説明する気概がないので、アバウトに説明する。

　適当な二次元正規分布のデータがあるとする。PCAと同じ要領で分散共分散行列を対角化する。

　対角化した行列の対角成分（＝固有値）は、データを軸同士の相関がないような空間に変換して（＝要するにぐるっと回して）あげたときの変換先の軸上における分散である。

　分散がわかれば、一次元のときの信頼区間的なものがわかる。それを決めるためのデータの中心位置からの距離にはマハラノビス距離という名前が付いている。そして、これのニ乗は $\chi^2$ 分布になるので、けっきょく累積確率だけ決めれば適当に定まる、ということがわかる*1。相関がなくなるように軸を作ったので、各軸ごとで求めてやれば（その軸上での）信頼楕円の径がわかる。

　この段階で楕円の幅と高さがわかっていることになるので、あとは適当に回転角を計算すると、

楕円の中心位置（単純に全データの平均で良い）
幅と高さ
回転角

　
　がわかることになり、楕円が描ける。

実装

　実装したものをそのまま貼っておきます。説明はしないので、参考にしたければしてください。

　一応入力データを色々変えてみておかしくなさそうなことは確認していますが、ミスがないとは言い切れません。ご注意ください。見つけたら指摘して頂けると嬉しいです。

# coding: UTF-8

import numpy as np
from scipy.stats import chi2

import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse

class ConfidenceEllipse:
    def __init__(self, data, p=0.95):
        self.data = data
        self.p = p

        self.means = np.mean(data, axis=0)
        self.cov = np.cov(data[:,0], data[:,1])

        lambdas, vecs = np.linalg.eigh(self.cov)
        order = lambdas.argsort()[::-1]
        lambdas, vecs = lambdas[order], vecs[:,order]

        c = np.sqrt(chi2.ppf(self.p, 2))
        self.w, self.h = 2 * c * np.sqrt(lambdas)
        self.theta = np.degrees(np.arctan(
            ((lambdas[0] - lambdas[1])/self.cov[0,1])))
        
    def get_params(self):
        return self.means, self.w, self.h, self.theta

    def get_patch(self, line_color="black", face_color="none", alpha=0):
        el = Ellipse(xy=self.means,
                     width=self.w, height=self.h,
                     angle=self.theta, color=line_color, alpha=alpha)
        el.set_facecolor(face_color)
        return el

def gen_data():
    return np.random.multivariate_normal([3,3], [[0.3,-0.2],[-0.2,1]], size=100)
def main():
    data = gen_data()

    fig = plt.figure()
    ax = fig.add_subplot(1,1,1)
    ax.scatter(data[:,0], data[:,1], color="b", marker=".", s=3)

    el = ConfidenceEllipse(data, p=0.95)
    ax.add_artist(el.get_patch(face_color="blue", alpha=0.5))
    plt.savefig("img.png")

if __name__ == "__main__":
    main()

　こんな感じの絵が出てきます。

f:id:hayataka2049:20180214235058p:plain

　正しいはずだけど、詳しくチェックはしていないので、自己責任でご利用ください。

参考文献

*1:てきとーすぎる説明だな

静かなる名辞

pythonとプログラミングのこと

【python】95%信頼楕円/確率楕円を描画する

理論的なもの

実装

参考文献