多次元確率分布　多次元正規分布

確率変数$${y_1, \cdots, y_d}$$がそれぞれ独立な正規分布に従う時、$${{\bf y}=(y_1, \cdots, y_d)^T}$$と記せば、この$${d}$$次元の確率変数$${{\bf y}}$$の確率密度関数は、
$${g({\bf y})=\displaystyle{\Pi^d_{j=1}\frac{1}{\sqrt{2\pi}} exp(-\frac{y_j^2}{2})=\frac{1}{(2\pi)^{d/2}} e^{-\frac{1}{2}{\bf y \cdot y}}}}$$
で与えられ、それぞれが正規分布に従うことから、
$${E[{\bf y}]=0}$$であり、分散共分散行列は$${d\times d}$$の単位行列を用いて
$${V[{\bf y}]= {\bf I}}$$
である。
　$${d\times d}$$の正則行列$${{\bf K}}$$と$${d}$$次元のベクトル$${{\bf \mu}}$$を使って、$${{\bf y}}$$を
$${{\bf x}={\bf y K}+{\bf \mu}}$$と変換した$${{\bf x}}$$の確率密度関数は、$${\bf y}$$の確率密度関数$${g({\bf y})}$$を用いて、
$${f({\bf x})=\displaystyle{\frac{1}{|det{\bf K}|} g({\bf x})}}$$
と与えられる。
$${{\bf y}={\bf K}^{-1}({\bf x} -{\bf \mu}) , \ {\bf y}^T=({\bf x} -{\bf \mu})^T({\bf K}^{-1})^T}$$を用いて、
$${f({\bf x})=\displaystyle{\frac{1}{|det{\bf K}|} \frac{1}{(2\pi)^{d/2}} exp\Big(-\frac{1}{2}({\bf x-\mu})^T ({\bf K}^{-1})^T{\bf K}^{-1}({\bf x- \mu})\Big) }}$$
ここで、$${({\bf K K}^{T})^{-1}=\Sigma^{-1}}$$とおけば、$${(det{\bf K})^2=det(\Sigma)}$$となり、
$${f({\bf x})=\displaystyle{ \frac{1}{\sqrt{det{\bf \Sigma}}} \frac{1}{(2\pi)^{d/2}} exp\Big(-\frac{1}{2}({\bf x-\mu})^T {\bf \Sigma}^{-1}({\bf x- \mu})\Big)}}$$で与えられる。

$${\bf x}$$の期待値と分散は、
$${E[{\bf x}]={\bf K}E[{\bf y}]+\bf \mu= \mu}$$
$${V[{\bf x}]=V[{\bf K y} + {\bf \mu}]={\bf K} V[{\bf y}]{\bf K}^{T}={\bf\Sigma}}$$
となる。
よって、$${\bf x}$$は、期待値が$${\bf \mu}$$で分散共分散行列$${\bf \Sigma}$$の多変量正規分布となる。

多変量正規分布に従う乱数のpythonでの実装は、numpy.randomのmultivariate_normalで行うことができる。以下のコードでは、平均がそれぞれ$${1,2}$$、共分散が
$${\begin{pmatrix}6 & 3 \\3 & 5 \\\end{pmatrix}}$$
とした２変量の正規分布に従う乱数を$${1000000}$$生成し、$${1000000}$$行$${2}$$列の行列に入れている。
それぞれの列の平均と、分散とをとれば、ほぼ指定の平均と共分散が確かめられている。

import numpy as np
from numpy import random as rand
import matplotlib.pyplot as plt
rand.seed(0)
randnum=1_000_000
mean = (1, 2)
cov = np.array([[6, 3], [3, 5]])
x=rand.multivariate_normal(mean, cov, randnum)
print(x.shape)
print(x[:,0].mean(),x[:,1].mean(),np.cov(x.T))
plt.plot(x[:, 0], x[:, 1], '.', alpha=0.5)
plt.axis('equal')
plt.grid()
plt.show();