pythonで多変量正規分布をグラフにしてイメージを掴む

世の中の色々な現象を表現する正規分布。
その多変量版である多変量正規分布についてpythonで色々やってみました。
jupyter notebookでやると良いかもです。

※この記事は私が運営するブログの下記記事を書き直したものです。
ふぐの子「pythonで多変量正規分布」

目次

• 多変量正規分布とは？
• 2次元のケース
• 3次元のケース
• まとめ

多変量正規分布とは？

n次元の多変量正規分布は、以下のように表されます。
$$ f(\vec{x}) = \frac{1}{(2\pi)^{n/2}\sqrt{|\Sigma|}}\exp(\frac{1}{2}(\vec{x}-\vec{\mu} )^T\Sigma^{-1}(\vec{x}-\vec{\mu} )) $$

1次元の正規分布は平均と分散の2つのパラメータで表現できるが、
多変量の場合、これが式の中で出てきている、次のベクトルと行列に置き換えられる。

 \begin{equation}\vec{\mu}=\begin{pmatrix}\mu_1 \\ \mu_2 \\ \vdots \\ \mu_n \\  \end{pmatrix} , \ \ \ \Sigma =  \begin{pmatrix} \sigma_{1}^2 & \cdots & \sigma_{1i} & \cdots & \sigma_{1n}\\ \vdots & \ddots & & & \vdots \\ \sigma_{i1} & & \sigma_{i}^2 & & \sigma_{in} \\ \vdots & & & \ddots & \vdots \\ \sigma_{n1} & \cdots & \sigma_{ni} & \cdots & \sigma_{n}^2 \end{pmatrix} \end{equation}

ここで、$ \mu_{i}, \ \ \ \sigma_{i}^2 $ は $i$ 番目の変数の平均、分散で、
$ \sigma_{ij}=\sigma_{ji} \ \ (i\neq j) $ は $i$ 番目と $j$ 番目の変数間の共分散。
この非対角項の成分があると、その変数間に相関があるようなイメージになります。

それでは次に、pythonで多変量席分布のグラフをプロットして、どのようなものなのかイメージをつかんでみます。

■2次元のケース

from numpy.random import *
from matplotlib import pyplot as plt
from pylab import rcParams
import numpy as np
rcParams['figure.figsize'] = 20,20

# 2次元正規分布version

# 平均
mu1=10
mu2=20
# 分散
sigma1=30
sigma2=30
sigma12=20
sigma1_2=90
sigma2_2=10
sigma12_2=30

# 平均ベクトルと共分散行列
mu = [mu1, mu2]
sigma = [[sigma1, sigma12], [sigma12, sigma2]]
sigma_2 = [[sigma1_2, sigma12], [sigma12, sigma2_2]]
sigma_3 = [[sigma1, sigma12_2], [sigma12_2, sigma2]]

# 2次元正規乱数
values = multivariate_normal(mu, sigma, 5000)
values_2 = multivariate_normal(mu, sigma_2, 5000)
values_3 = multivariate_normal(mu, sigma_3, 5000)

# グラフをプロット
fig = plt.figure()
ax0 = fig.add_subplot(2,2,1)
ax0.plot(values[:,0], values[:,1], 'o', c='b',label="2d", alpha=0.1)
ax0.set_title("2d normalplot")
ax0.set_xlim(-50, 50)
ax0.set_ylim(-50, 50)
ax0.grid()
# パラメータをいじってみる
ax = fig.add_subplot(2,2,2)
ax.plot(values[:,0], values[:,1], 'o', c='b',label="2d", alpha=0.1)
ax.plot(values_2[:,0], values_2[:,1], 'o', c='r',label="2d", alpha=0.1)
ax.plot(values_3[:,0], values_3[:,1], 'o', c='g',label="2d", alpha=0.1)
ax.set_title("2d normalplot")
ax.set_xlim(-50, 50)
ax.set_ylim(-50, 50)
ax.grid()
fig.show()

左のグラフを見ると、平均 $(x_1, x_2) = (10, 20)$ の周りで分散に従って分布しているのが分かる。
変数 $x_1 x_2$ 間の共分散も sigma12$ = 20$ に設定しているため、相関関係を持つようなプロットになった。

右のグラフではさらにパラメータを変動させてみた。
赤色は分散を横長に大きくしている。
緑色は変数 $x_1 x_2$ 間の共分散を sigma12_2$ = 20$ にしていて、ちょうど直線に乗るようなグラフになっている。

相関係数を設定値とグラフの実際の値で比較してみると以下のようになり
同様値に確かになっている。
特に緑色は相関係数が1となるので、直線に乗るようなグラフになった。

# 相関係数の理論計算値
corr1_theory =sigma12 / np.sqrt(sigma1 * sigma2)
corr2_theory =sigma12 / np.sqrt(sigma1_2 * sigma2_2)
corr3_theory =sigma12_2 / np.sqrt(sigma1 * sigma2)
print(corr1_theory,corr2_theory,corr3_theory)

# 相関係数の実測値
corr1 = np.corrcoef(values[:,0],values[:,1])[0,1]
corr2 = np.corrcoef(values_2[:,0],values_2[:,1])[0,1]
corr3 = np.corrcoef(values_3[:,0],values_3[:,1])[0,1]
print(corr1,corr2,corr3)
0.6666666666666666 0.6666666666666666 1.0
0.6670848283319503 0.6527793760878543 0.9999999999999999

■3次元のケース

3次元でも同様にグラフをプロットした。

%matplotlib inline
from numpy.random import *
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from pylab import rcParams
rcParams['figure.figsize'] = 20,20

# 3次元正規分布version

# 平均
mu = [0, 0, 0]

# 分散
sigma1 = 40
sigma2 = 40
sigma3 = 40
sigma12 = 20
sigma23 = 20
sigma31 = 20
sigma1_2 = 90
sigma2_2 = 10
sigma3_2 = 30
sigma12_2 = 40
sigma23_2 = 40
sigma31_2 = 40

sigma =  [[sigma1, sigma12, sigma31], [sigma12, sigma2, sigma23], [sigma31, sigma23, sigma3]]
sigma_2 = [[sigma1_2, sigma12, sigma31], [sigma12, sigma2_2, sigma23], [sigma31, sigma23, sigma3_2]]
sigma_3 = [[sigma1, sigma12_2, sigma31_2], [sigma12_2, sigma2, sigma23_2], [sigma31_2, sigma23_2, sigma3]]

# 3次元正規乱数
values1 = multivariate_normal(mu, sigma, 5000)
values2 = multivariate_normal(mu, sigma_2, 5000)
values3 = multivariate_normal(mu, sigma_3, 5000)

# グラフをプロット
fig = plt.figure()
ax0 = fig.add_subplot(221, projection='3d')
ax0.set_xlim([-50,50])
ax0.set_ylim([-50,50])
ax0.set_zlim([-50,50])
ax0.scatter(values1[:,0], values1[:,1], values1[:,2], s=10, alpha=0.1, c='b')
ax0.legend()
# パラメータをいじってみる
ax1 = fig.add_subplot(222, projection='3d')
ax1.set_xlim([-50,50])
ax1.set_ylim([-50,50])
ax1.set_zlim([-50,50])
ax1.scatter(values1[:,0], values1[:,1], values1[:,2], s=10, alpha=0.1, c='b')
ax1.scatter(values2[:,0], values2[:,1], values2[:,2], s=10, alpha=0.1, c='r')
ax1.scatter(values3[:,0], values3[:,1], values3[:,2], s=10, alpha=0.1, c='g')
ax1.legend()
plt.show()

まとめ

多変量正規分布について概念を学習。pythonを用いて、2次元と3次元のグラフをプロットしてイメージを掴むことができた。