2次元正規分布サンプリングと確率等高線の図示

2次元正規分布からデータをサンプルしよう!

 正規分布からサンプリングするにはどうすればよいのでしょうか?1次元の場合はすぐに実装できますが、2次元正規分布の場合は、共分散を考慮しないといけないので、少しサンプル方法がトリッキーになります。{x_1}{x_2}をそれぞれ独立に1次元正規分布からサンプリングしてしまうと上手く行きません。条件付き分布の考え方を導入する必要があります。このあたりの話を参考にしつつ、実装しました。
 まずは、基本的な情報をセッティングします。平均ベクトル、共分散行列の値は『はじめてのパターン認識』47ページを基に作成しました。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab

#info of x1
mu_x1 = 2
sd_x1 = np.sqrt(0.7068966)

#info of x2
mu_x2 = 1
sd_x2 = np.sqrt(1.7931034)

#info of cov
cov = 0.5172414

#lo
lo = cov / (sd_x1 * sd_x2)

ここで、一変数目({x_1})をサンプリングします。ここでは通常の正規乱数を使うだけです。

#samples of x1
n = 200
x1 = np.random.normal(loc=mu_x1,scale=sd_x1,size=n)

次に、{x_1}が与えられた時の{x_2}を考えていきます。ここの考え方は、下記資料が参考になります。
http://www012.upp.so-net.ne.jp/doi/sas/simulation/multi_co/bivariate_normal.pdf

#x2 given x1
x2 = np.ones(n)
for i in range(n):
    mean = mu_x2 + lo * (sd_x1/sd_x2) * (x1[i]-mu_x1)
    stdev = (1-lo**2)*(sd_x2**2)
    x2[i] = np.random.normal(loc=mean,scale=stdev,size=1)

plt.plot(x1,x2,'.')
plt.xlim([-2,6])
plt.ylim([-3,5])
plt.xlabel('x1')
plt.ylabel('x2')
plt.title('Samples from Normal Distribution')

これで、このようなグラフが書けます。2次元正規分布からのサンプリングができました。
f:id:decompose:20160325233855p:plain

確率の等高線を図示しよう!

最後に、確率の等高線も同じグラフに表示させていきます。

x = np.linspace(-7.0,7.0,200)
y = np.linspace(-5.0,5.0,200)
X,Y=np.meshgrid(x,y)
Z = mlab.bivariate_normal(X,Y,sd_x1,sd_x2,mu_x1,mu_x2,cov)
CS = plt.contour(X,Y,Z)
plt.clabel(CS,inline=1,fontsize=10)
plt.show()

これで、このようなグラフが書けます。
f:id:decompose:20160325233859p:plain

参考

この本の第4章(確率モデルと識別関数)の4.2節(確率モデル)の4.2.1(正規分布関数)に基づいています。
www.amazon.co.jp