多変量正規分布をPythonでplotして理解する

#はじめに
統計を勉強していた際に出てきた「多変量正規分布」のイメージを掴むためpythonでplotしてみました。今回は可視化して際にわかりやすいよう$n$数を2にして二次元正規分布をplotしています。

##参考
多変量正規分布の理解とそのplotを行うに当たって下記を参考にさせていただきました。

- [多変量正規分布を理解する](https://math-note.com/multivariate-normal-distribution/)
- [pythonで多変量正規分布をグラフにしてイメージを掴む](https://qiita.com/fugunoko/items/g)
- [二次元ガウス分布をplotしてみる](https://openbook4.me/projects/231/sections/1538)

#多変量正規分布の概要
$n$変数の多変量正規分布は下記のように表されます。

```math

f(\vec{x}) = \frac{1}{\sqrt{(2\pi)^n |\sum|}}exp \left \{-\frac{1}{2}{}^t (\vec{x}-\vec{\mu}) {\sum}^{-1} (\vec{x}-\vec{\mu}) \right \}

```
変数が$n$個あるためデータを$n$次元ベクトル表記で表します。さらに平均値$\mu$も変数の数だけ存在するため同様にベクトル表記で表します。

```math

{ \begin{equation}\vec{x}=\begin{pmatrix}x_1 \\ x_2 \\ \vdots \\ x_n \\  \end{pmatrix}, \vec{\mu}=\begin{pmatrix}\mu_1 \\ \mu_2 \\ \vdots \\ \mu_n \\  \end{pmatrix}   \end{equation}
}

```
ある一つの要素$x_{i}$が確率変数$X_{i}$のデータを表し、平均値$\mu_i$が確率変数$X_{i}$の平均値を表します。
続いて分散についてですが、多変量となる場合は各データの分布のみならずデータ間の相関を考慮する必要があるため**分散共分散行列$\sum$**を用います。

```math

{ \begin{equation}\ \ \ \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}
}

```
$\sigma^2_i$は$i$番目の変数の分散で、$\sigma_{ij} =\sigma_{ji} (i≠j)$は$i$番目と$j$番目の変数間の共分散になっています。
そして$n$が$2$の時の二次元正規分布は下記のように表されます。

```math
N_2 \left ( \begin{pmatrix}  \mu_x \\  \mu_y \\  \end{pmatrix} , \begin{pmatrix}  \sigma_{x}^2 & \sigma_{xy}\\  \sigma_{xy} & \sigma_{y}^2\\  \end{pmatrix} \right  )

```

それでは二次元正規分布をplotしてみたいと思います。

#二次元正規分布のplot
二次元正規分布をplotするスクリプトが下記になります。
まずは2変数とも標準正規分布に従い、互いに独立な場合で出力してみます。

```python
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm

#関数に投入するデータを作成
x = y = np.arange(-20, 20, 0.5)
X, Y = np.meshgrid(x, y)

z = np.c_[X.ravel(),Y.ravel()]

#二次元正規分布の確率密度を返す関数
def gaussian(x):
    #分散共分散行列の行列式
    det = np.linalg.det(sigma)
    print(det)
    #分散共分散行列の逆行列
    inv = np.linalg.inv(sigma)
    n = x.ndim
    print(inv)
    return np.exp(-np.diag((x - mu)@inv@(x - mu).T)/2.0) / (np.sqrt((2 * np.pi) ** n * det))

#2変数の平均値を指定
mu = np.array([0,0])
#2変数の分散共分散行列を指定
sigma = np.array([[1,0],[0,1]])

Z = gaussian(z)
shape = X.shape
Z = Z.reshape(shape)

#二次元正規分布をplot
fig = plt.figure(figsize = (15, 15))
ax = fig.add_subplot(111, projection='3d')
    
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm)
plt.show()

```
出力結果は下記です。2変数とも正規分布のため偏りのない尖ったグラフになります。

![多変量正規分布1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/395230/290abae2-093e-b90d-9c2d-c6d143d144b1.png)

それでは異なる形のグラフもplotします。
2変数の分布が下記となる場合の二次元正規分布をplotしてみましょう。

```python
#2変数の平均値を指定
mu = np.array([3,1])
#2変数の分散共分散行列を指定
sigma = np.array([[10,5],[5,10]])
```
下記は先ほどのplotと同様です。

```python

Z = gaussian(z)
shape = X.shape
Z = Z.reshape(shape)

#二次元正規分布をplot
fig = plt.figure(figsize = (15, 15))
ax = fig.add_subplot(111, projection='3d')
    
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm)
plt.show()

```

出力結果は下記になります。今回は互いに相関している分布をplotしたため、少し斜めに歪んだ形になっていることがわかります。

![ダウンロード (1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/395230/a19771c3-9b0e-6359-a9f4-c231ab5df5fc.png)

数式上ではわかりにくかったことも、こうやって可視化するとイメージが掴みやすいですね。

#Next
統計を勉強していると中々数式だけでイメージが掴めないことも多いので、pythonで自分で書いてみたりplotして可視化するのを積極的にやっていきたいと思います。