まずは2次元の正規分布を描画する。簡便のために対角成分以外は0とする。
import numpy as np from scipy.stats import multivariate_normal import matplotlib.pyplot as plt np.random.seed(2) mu_x,mu_y = 5,10 sigma_x,sigma_y = 1,5 cov = np.array([[1, 0], [0, 1]]) # 共分散行列 n_samples = 100000 samples = np.random.multivariate_normal(np.array([0, 0]), cov, n_samples) x = samples.T[0]*sigma_x+mu_x y = samples.T[1]*sigma_y+mu_y fig = plt.figure(figsize=(6, 6)) grid = plt.GridSpec(4, 4, hspace=0.2, wspace=0.2) main_ax = fig.add_subplot(grid[1:4, 0:3]) main_ax.scatter(x, y, alpha=0.5, s=1) xlim = main_ax.get_xlim() ylim = main_ax.get_ylim() x_hist = fig.add_subplot(grid[0, 0:3], sharex=main_ax) x_hist.hist(x, bins=40, color='tab:gray', range=xlim) y_hist = fig.add_subplot(grid[1:4, 3], sharey=main_ax) y_hist.hist(y, bins=40, orientation='horizontal', color='tab:gray', range=ylim) main_ax.set_xlabel('X') main_ax.set_ylabel('Y') plt.show()
この正規分布における、90%信頼区間はどう引けばよいだろうか?
対角成分のみの2次元正規分布の中心からのXY成分の標準偏差で規格化した距離の平方がカイ二乗分布に従うことより、自由度2のカイ二乗分布のp値が0.1になるときのカイ二乗の平方根をとればよいと分かる。つまり sigma * chi2.isf(0.1,2) ** 0.5 である。99% C.L.であれば、sigma * chi2.isf(0.01,2) ** 0.5 である。
main_ax.scatter(x, y, alpha=0.5, s=1) の後に以下のコードを挿入しよう。matplotlibの楕円描画では、長軸と短軸の全長をそれぞれ指定するため、標準偏差の2倍を与える。
mu_x0,mu_y0 = mu_x,mu_y # 初期値を記録 sigma_x0,sigma_y0 = sigma_x,sigma_y # 初期値を記録 from scipy.stats import chi2 sigma_x = sigma_x0 * chi2.isf(0.1,2)**0.5 # 自由度2、p値が0.1 sigma_y = sigma_y0 * chi2.isf(0.1,2)**0.5 # 自由度2、p値が0.1 import matplotlib.patches as patches ellipse = patches.Ellipse((mu_x, mu_y), sigma_x*2, sigma_y*2, edgecolor='r', facecolor='none',ls="-",alpha=0.5) main_ax.add_patch(ellipse)
それっぽくできた。本当にこの楕円の中に90%のイベントが入っているか確かめてみよう。楕円内にあるかどうかを判定するコードを書く。
先程のsigma_y = sigma_y0 * chi2.isf(0.1,2)**0.5 # 自由度2、p値が0.1 の行の直後に、以下のコードを挿入しよう
def is_in_ellipse(x, y, mu_x, mu_y, sigma_x, sigma_y): return ((x - mu_x)**2 / sigma_x**2 + (y - mu_y)**2 / sigma_y**2) <= 1 hitx = x[is_in_ellipse(x,y,mu_x,mu_y,sigma_x, sigma_y)] hity = y[is_in_ellipse(x,y,mu_x,mu_y,sigma_x, sigma_y)] main_ax.scatter(hitx, hity, alpha=0.5, s=1, color="tab:green") main_ax.text(mu_x,mu_y,f"{len(hitx)/n_samples:.2%}",ha="center",va="center")
統計上のゆらぎの範囲で、確かに90%がこの楕円 (図上では真円に近いが軸を見ると確かに楕円) の中に入っていることが分かる。
この楕円を、任意の形、任意の場所に持っていったときの確率はどう求めればよいか。おそらく2つ方法、モンテカルロ法と積分法がある。
モンテカルロ法は簡単で、先程の関数を使い回すことで求めることができる。
sigma_x sigma_y に代入した部分を下記のように書き換えよう。
# sigma_x = sigma_x0 * chi2.isf(0.1,2)**0.5 # 自由度2、p値が0.1 # sigma_y = sigma_y0 * chi2.isf(0.1,2)**0.5 # 自由度2、p値が0.1 mu_x,mu_y=3,5 sigma_x,sigma_y=1,2
この楕円の中には約2.23%のイベントがあることが分かった。
次に、積分法で求めてみる。2次元の正規分布の平均値を0、標準偏差を1に規格化してしまっていたので、与える平均値と標準偏差も同じく規格化して与える。
from scipy.integrate import dblquad def pdf(x, y): return multivariate_normal([0, 0], cov).pdf([x, y]) norm_mu_x = (mu_x-mu_x0)/sigma_x0 norm_mu_y = (mu_y-mu_y0)/sigma_y0 norm_sigma_x = sigma_x/sigma_x0 norm_sigma_y = sigma_y/sigma_y0 P, error = dblquad( lambda x, y: pdf(x, y), norm_mu_x-norm_sigma_x, norm_mu_x+norm_sigma_x, lambda x: norm_mu_y - norm_sigma_y * np.sqrt(1 - ((x-norm_mu_x) / norm_sigma_x)**2), # y の下限 lambda x: norm_mu_y + norm_sigma_y * np.sqrt(1 - ((x-norm_mu_x) / norm_sigma_x)**2) # y の上限 ) print(f"楕円領域内の確率: {P:.3%}±{error:.3%}")
楕円領域内の確率: 2.213%±0.000%
と出力されるだろう。先程の 約2.23%と誤差の範囲で一致している。積分法は誤差も算出してくれるので、特に確率が低い領域での計算で役に立つだろう。
積分法と、描画を全て統合したコード
import numpy as np from scipy.stats import multivariate_normal import matplotlib.pyplot as plt np.random.seed(2) mu_x,mu_y = 5,10 sigma_x,sigma_y = 1,5 cov = np.array([[1, 0], [0, 1]]) # 共分散行列 n_samples = 100000 samples = np.random.multivariate_normal(np.array([0, 0]), cov, n_samples) x = samples.T[0]*sigma_x+mu_x y = samples.T[1]*sigma_y+mu_y fig = plt.figure(figsize=(6, 6)) grid = plt.GridSpec(4, 4, hspace=0.2, wspace=0.2) main_ax = fig.add_subplot(grid[1:4, 0:3]) main_ax.scatter(x, y, alpha=0.5, s=1) mu_x0,mu_y0 = mu_x,mu_y # 初期値を記録 sigma_x0,sigma_y0 = sigma_x,sigma_y # 初期値を記録 from scipy.stats import chi2 # sigma_x = sigma_x0 * chi2.isf(0.1,2)**0.5 # 自由度2、p値が0.1 # sigma_y = sigma_y0 * chi2.isf(0.1,2)**0.5 # 自由度2、p値が0.1 mu_x,mu_y=3,5 sigma_x,sigma_y=1,2 norm_mu_x = (mu_x-mu_x0)/sigma_x0 norm_mu_y = (mu_y-mu_y0)/sigma_y0 norm_sigma_x = sigma_x/sigma_x0 norm_sigma_y = sigma_y/sigma_y0 x_min, x_max = norm_mu_x-norm_sigma_x, norm_mu_x+norm_sigma_x y_lower = lambda x: norm_mu_y - norm_sigma_y * np.sqrt(1 - ((x - norm_mu_x) / norm_sigma_x) ** 2) # y の下限 y_upper = lambda x: norm_mu_y + norm_sigma_y * np.sqrt(1 - ((x - norm_mu_x) / norm_sigma_x) ** 2) # y の上限 def pdf(x, y): return multivariate_normal([0, 0], cov).pdf([x, y]) from scipy.integrate import dblquad P, error = dblquad( lambda x, y: pdf(x, y), x_min, x_max,y_lower,y_upper) # 下限と上限のy値を計算 x_values = np.linspace(x_min,x_max,1000) y_values_lower = y_lower(x_values) y_values_upper = y_upper(x_values) main_ax.fill_between(x_values*sigma_x0+mu_x0, y_values_lower*sigma_y0+mu_y0, y_values_upper*sigma_y0+mu_y0, color='tab:green', alpha=0.7) main_ax.text(mu_x,mu_y,f"{P:.2%}",ha="center",va="center") import matplotlib.patches as patches ellipse = patches.Ellipse((mu_x, mu_y), sigma_x*2, sigma_y*2, edgecolor='r', facecolor='none',ls="-",alpha=0.5) main_ax.add_patch(ellipse) xlim = main_ax.get_xlim() ylim = main_ax.get_ylim() x_hist = fig.add_subplot(grid[0, 0:3], sharex=main_ax) x_hist.hist(x, bins=40, color='tab:gray', range=xlim) y_hist = fig.add_subplot(grid[1:4, 3], sharey=main_ax) y_hist.hist(y, bins=40, orientation='horizontal', color='tab:gray', range=ylim) main_ax.set_xlabel('X') main_ax.set_ylabel('Y') plt.show()
