「多次元正規分布の式ってどうなってるの？」の疑問に丁寧に答える

# 多次元正規分布の式、複雑すぎ問題
統計学の基本とも言える**正規分布**。

この正規分布は、下図のように**多次元の変数でも定義することができます**（**多次元正規分布**）
<img width="400" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/b4373333-12e3-9bc5-a944-21541477d4e7.png">

多次元正規分布は多くの機械学習アルゴリズム（[主成分分析](https://quest-research.co.jp/research/principal-component-analysis)、[ナイーブベイズ分類器](https://qiita.com/aflc/items/13fe52243c35d3b678b0)、[マハラノビスタグチ法](https://www.toshiba-dme.co.jp/dme/coretec/hin/hin_mt.htm)、[混合ガウスモデル](https://datachemeng.com/gaussianmixturemodel/)等）で活用される確率分布であり、特に**教師なし学習**においては**最重要理論の一つ**と言える存在であり、データ分析において避けては通れない理論と言えます。

一方で以下の**多次元正規分布の式**(**確率密度関数**)は多くの書籍で天下り的に与えられており、

```math
N(\boldsymbol{x} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma})=\frac{1}{(2\pi)^\frac{n}{2}}\frac{1}{|\boldsymbol{\Sigma}|^\frac{1}{2}} \exp \Bigl(-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^T \boldsymbol{\Sigma}^{-1} (\boldsymbol{x}-\boldsymbol{\mu}) \Bigr)
```

<br>

**「複雑すぎて分からん！ちゃんと解説してくれ！」**
<img width=320 src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/8aeb0491-e864-420b-c3db-1c3f5e84daa3.png">
とプンプンになってしまった方も多いかと思います。私もその一人です。

というわけで、精神の平静？を取り戻すべく**上記の式の意味**について調べたので、**なるべく丁寧に解説**したいと思います。

# 本文に進む前に
本記事の内容を理解するためには、**「同時分布」「変数の独立性」「共分散」などに関する知識が必要**です。
これらについては**以下の記事**で詳細に解説しましたので、まずは**こちらをご一読頂ければ**と思います。

https://qiita.com/c60evaporator/items/93295d3890c9e20e521d

# 多次元正規分布の解説
多次元正規分布について、具体的な解説を進めていきます。

2次元の正規分布 → 3次元以上の正規分布

の順で進めます。

## 2次元正規分布
まずは理解しやすいよう2次元で解説します。

### 独立変数の2次元正規分布
2次元の正規分布をシンプルに考えると、以下のように**1つ目の確率変数X1の正規分布と2つ目の確率変数X2の正規分布を掛け合わせた同時分布** $p(x_1,x_2)$が思い浮かびます。

```math
p_1(x_1)=\frac{1}{\sqrt{2\pi\sigma_1^2}}\exp\Bigl(-\frac{1}{2\sigma_1^2}(x_1-\mu_1)^2\Bigr)\\
p_2(x_2)=\frac{1}{\sqrt{2\pi\sigma_2^2}}\exp\Bigl(-\frac{1}{2\sigma_2^2}(x_2-\mu_2)^2\Bigr)
```

```math
p(x_1,x_2)=p_1(x_1)p_2(x_2)
```

上記の分布に関してPythonで等高線を作成し、可視化してみます（$\mu_1=1, \mu_2=1, \sigma_1=1, \sigma_2=2$とする）

```python:独立2変数の正規分布
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# パラメータ
MU1 = 1
MU2 = 1
SIGMA1 = 2
SIGMA2 = 1
# (x1,x2)格子点を作成
x1_grid = np.linspace(-6, 6, 100)
x2_grid = np.linspace(-6, 6, 100)
X1, X2 = np.meshgrid(x1_grid, x2_grid)

# 同時分布関数
def p(x1, x2):
    p1 = np.exp(-np.square(x1-MU1)/(2*SIGMA1**2))/np.sqrt(2*np.pi*SIGMA1)
    p2 = np.exp(-np.square(x2-MU2)/(2*SIGMA2**2))/np.sqrt(2*np.pi*SIGMA2)
    p = p1*p2
    return p

# 同時分布をプロットするメソッド
def plot_joint_distribution(x1, x2):
    # 同時分布を計算
    P = p(X1, X2)
    # 曲面をプロット
    fig = plt.figure(figsize = (8, 8))
    ax = fig.add_subplot(111, projection="3d")
    ax.set_xlabel("x1", size = 16)  # x1軸
    ax.set_ylabel("x2", size = 16)  # x2軸
    ax.set_zlabel("p(x1,x2)", size = 16)  # p軸
    ax.plot_surface(X1, X2, P, cmap = "YlGn_r")
    plt.show()
    # 等高線をプロット
    fig, axes = plt.subplots(1, 1, figsize=(6, 6))
    plt.contour(X1, X2, P, cmap="YlGn_r")
    ax.set_xlabel("x1", size = 16)  # x1軸
    ax.set_ylabel("x2", size = 16)  # x2軸
    plt.show()

# 同時分布プロット実行
plot_joint_distribution(X1, X2)
```
**・3Dグラフ**
![norm3d_indpendent.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/a5810136-b9ec-2520-d9a4-9a17b16e9cce.png)
**・等高線**
![norm_indpendent.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/cfe689fc-1819-e321-ab6d-646c7adb95a6.png)

この分布は$p(x_1,x_2)=p_1(x_1)p_2(x_2)$という式から、変数X1とX2は独立であることが分かります。
よって上の等高線のように、**この式では**どこで切っても断面形状が同じ分布しか再現できず、冒頭のような**各変数が相関を持った斜めの分布**が再現できません。

斜めの確率分布を再現するにはどうすれば良いでしょうか？

### 2次元正規分布と変数変換
斜めの確率分布を再現するには、**行列による変数変換**が有効です。

変換前の変数を(y1, y2)、変換後の変数を(x1, x2)とすると、変換行列

```math
\boldsymbol{A}
=
\begin{pmatrix}
a & b \\
c & d
\end{pmatrix}
```
を用いて

```math
\left(\begin{array}{c}
x_1 \\
x_2 \\
\end{array}\right)
=
\begin{pmatrix}
a & b \\
c & d
\end{pmatrix}
\left(\begin{array}{c}
y_1 \\
y_2 \\
\end{array}\right)
```
と表すことで、行列による変数変換が実現できます。

「行列」と聞いただけで頭が痛くなってきた方も多いかと思いますが、そんな時は**変換要素を分解**して、**グラフで変換過程を可視化**すると、理解が進むかと思います。

#### ・QR分解と変換要素の分解
全てのnxn行列**A**は、直交行列**Q**と上三角行列**R**を使って以下のように分解できます。

```math
\boldsymbol{A}=\boldsymbol{Q}\boldsymbol{R}
```
<img width=540 src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/93457960-88ff-e7e9-7699-ef8bdbdedcac.png">

この分解を[QR分解](https://ja.wikipedia.org/wiki/QR%E5%88%86%E8%A7%A3)と呼び、[Aが正則であればRの対角成分が全て正となって](http://nalab.mind.meiji.ac.jp/~mk/labo/text/eigenvalues/node32.html)扱いやすくなります。

2x2行列においては任意のa,b,c,dに対して以下のように表せます（以後、Aが非正則な場合は片方の変数が消えて2次元正規分布が表現できないため、考慮しないこととします）

```math
\begin{pmatrix}
a & b \\
c & d
\end{pmatrix}
=
\boldsymbol{Q}
\begin{pmatrix}
m_1 & m_3 \\
0 & m_2
\end{pmatrix}
```

#### ・直交行列Qの分解
2x2行列では直交行列Qは
左右反転行列**Mir**→回転行列**Rot**

```math
\boldsymbol{Q}=\boldsymbol{Rot}\cdot\boldsymbol{Mir}=
\begin{pmatrix}
\cos{\theta} & -\sin{\theta} \\
\sin{\theta} & \cos{\theta}
\end{pmatrix}
\begin{pmatrix}
1 & 0 \\
0 & (-1)^k
\end{pmatrix}
```
で表されます（[参考1](http://web.econ.keio.ac.jp/staff/tose/cours/2016/slin/STlin106Chapter07_20160630a.pdf)、[参考2](http://www.eds.it-hiroshima.ac.jp/koyama/19enshuu/FPkouki/11_P.pdf)）
（kは0か1を取り、0ならば反転なし、1ならば反転ありを表します）

式を見ると分かりますが、反転がある場合

```math
\boldsymbol{Q}=\begin{pmatrix}\cos{\theta} & \sin{\theta} \\ \sin{\theta} & -\cos{\theta} \end{pmatrix}
```

となり、反転がない場合

```math
\boldsymbol{Q}=\begin{pmatrix}\cos{\theta} & -\sin{\theta} \\ \sin{\theta} & \cos{\theta} \end{pmatrix}
```

となるため、対角成分が一致しているかで反転の有無を見分ける事ができます (一致していなければ反転あり)

#### ・上三角行列Rの分解
2x2行列では上三角行列Rは、以下のようにスキュー行列**Skw**と拡大行列**Mag**で表されます。
(tanφの値域が-∞〜∞なので、m1tanφは任意のm3を表現できます)

```math
\boldsymbol{R}=
\begin{pmatrix}
m_1 & m_3 \\
0 & m_2
\end{pmatrix}
=
\begin{pmatrix}
m_1 & m_1\tan\phi \\
0 & m_2
\end{pmatrix}
=
\begin{pmatrix}
m_1 & 0 \\
0 & m_2
\end{pmatrix}
\begin{pmatrix}
1 & \tan\phi \\
0 & 1
\end{pmatrix}
=\boldsymbol{Mag}\cdot\boldsymbol{Skw}
```

<br>
よって任意の2x2正則行列Aは、

```math
\boldsymbol{A}=
\boldsymbol{Rot}\cdot\boldsymbol{Mir}\cdot\boldsymbol{Mag}\cdot\boldsymbol{Skw}
```
と「スキュー」→「拡大」→（「反転」）→「回転」の組み合わせで表せる事が分かります。

### 2次元正規分布における変数変換の可視化
2次元正規分布においては、下図のように独立2変数の標準正規分布（平均が全て0, 分散が全て1）をスタートとして
**「1.標準正規分布」→「2.スキュー」→「3.拡大」→（「4.反転」）→「5.回転」**
の順の変換を、前述のように2x2行列**A**で表現できます。

上記変換に「**6.平行移動**」を加えた下図のような変換を実施することで、2次元の正規分布を定義することができます。

![正規分布変数変換.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/12286718-1f1a-e6b1-d1bf-4cbd133f57a5.png)

図を見れば明らかなように、行列**A**による**変数変換で斜めの正規分布が再現**できていますね！

### 各変換要素の解説
画像処理における座標変換（[参考](https://imagingsolution.blog.fc2.com/blog-entry-284.html)）をイメージすると分かりやすいですが、各変換要素の内容と行列について解説します。

#### 1. 2変数の標準正規分布
2変数における標準正規分布は、[先ほどの独立2変数の正規分布](https://qiita.com/c60evaporator/items/d53053358105b0117f2c#%E7%8B%AC%E7%AB%8B%E5%A4%89%E6%95%B0%E3%81%AE2%E6%AC%A1%E5%85%83%E6%AD%A3%E8%A6%8F%E5%88%86%E5%B8%83)において、μ=0、σ=1を代入したものとなります。
この標準正規分布の関数を$q(y_1,y_2)$とすると、以下のように表せます。

```math
q_1(y_1)=\frac{1}{\sqrt{2\pi}}\exp(-\frac{1}{2}y_1^2)\\
q_2(y_2)=\frac{1}{\sqrt{2\pi}}\exp(-\frac{1}{2}y_2^2)
```

```math
q(y_1,y_2)=q_1(y_1)q_2(y_2)=\frac{1}{2\pi}\exp\bigl(-\frac{1}{2}(y_1^2+y_2^2)\bigr)
```
標準正規分布は各変数の周辺分布の掛け算で表されるため、確率変数**y1, y2は独立**であることが分かります。

また$y_1^2+y_2^2$の部分は、ベクトル**y**=(y1, y2)同士の内積、すなわち

```math
\langle\boldsymbol{y},\boldsymbol{y}\rangle
```
となるため、**2変数標準正規分布はベクトルで以下のようにも表記**できます。

```math
q(\boldsymbol{y})=\frac{1}{2\pi}\exp\bigl(-\frac{1}{2}\langle\boldsymbol{y},\boldsymbol{y}\rangle\bigr)
```

#### 2. スキュー(せん断)
「スキュー」は馴染みのない言葉ですが、以下のように長方形を平行四辺形に変化させる変換を表します。
<img width="400" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/e3099837-8423-2103-d5b4-31100813a60b.png">
**スキューを表す行列**は、以下の式で表されます（φはスキュー角度）

```math
\boldsymbol{Skw}=\begin{pmatrix}
1 & \tan{\phi} \\
0 & 1
\end{pmatrix}
```

#### 3. 拡大
**拡大・縮小を表す行列**は、以下の式で表されます

```math
\boldsymbol{Mag}=\begin{pmatrix}
m_1 & 0 \\
0 & m_2
\end{pmatrix}
```

$m_1, m_2$は、

```math
m_1: y_1方向の拡大倍率\\
m_2: y_2方向の拡大倍率
```

を表しています。

#### 4. 反転
**反転を表す行列**は、以下の式で表されます
（kは0か1を取り、0ならば反転なし、1ならば反転ありを表します）

```math
\boldsymbol{Mir}=\begin{pmatrix}
1 & 0 \\
0 & (-1)^k
\end{pmatrix}
```

#### 5. 回転
**回転を表す行列**は、以下の式で表されます（θは回転角度）

```math
\boldsymbol{Rot}=\begin{pmatrix}
\cos{\theta} & -\sin{\theta} \\
\sin{\theta} & \cos{\theta}
\end{pmatrix}
```

#### 6. 平行移動
平行移動は行列でも表現できますが、ベクトルで表現した方がシンプルです。
（[前述のアフィン変換](https://imagingsolution.blog.fc2.com/blog-entry-284.html)のように行列のサイズを1つ大きくして平行移動を行列で表現可能）

平行移動前の変数を(x1',x2')とすると、平行移動後の変数(x1, x2)は以下のように表されます。

```math
\left(\begin{array}{c}
x_1 \\
x_2 \\
\end{array}\right)
=
\left(\begin{array}{c}
x_1' \\
x_2' \\
\end{array}\right)
+
\left(\begin{array}{c}
\mu_1 \\
\mu_2 \\
\end{array}\right)
```

以下のように、$\boldsymbol{\mu}=(\mu_{1}, \mu_{2})$は変換後の座標軸(x1, x2)における平均値を表しています

```math
\mu_{1}: x_1方向の平均値\\
\mu_{2}: x_2方向の平均値
```

この式は以下のように一般化して表現可能です

```math
\boldsymbol{x}=\boldsymbol{x'}+\boldsymbol{\mu}
```

#### ・変数変換の一般式
前述の「2.スキュー」→「3.拡大」→（「4.反転」）→「5.回転」→「6.平行移動」の変換に対して、

```math
変換前の確率変数\boldsymbol{y}=\left(\begin{array}{c}
y_1 \\
y_2 \\
\end{array}\right)
```

```math
変換後の確率変数\boldsymbol{x}=\left(\begin{array}{c}
x_1 \\
x_2 \\
\end{array}\right)
```
の関係は、以下の式で表せます（反転がない場合は行列**Mir**は不要）

```math
\boldsymbol{x}=(\boldsymbol{Rot}\cdot\boldsymbol{Mir}\cdot\boldsymbol{Mag}\cdot\boldsymbol{Skw})\boldsymbol{y}+\boldsymbol{\mu}\\
=\boldsymbol{A}\boldsymbol{y}+\boldsymbol{\mu}
```
Aの各成分を詳細に表記すると、以下のようになります。

```math
\boldsymbol{A}=
\begin{pmatrix}
\cos{\theta} & -\sin{\theta} \\
\sin{\theta} & \cos{\theta}
\end{pmatrix}
\begin{pmatrix}
1 & 0 \\
0 & (-1)^k
\end{pmatrix}
\begin{pmatrix}
m_1 & 0 \\
0 & m_2
\end{pmatrix}
\begin{pmatrix}
1 & \tan{\phi} \\
0 & 1
\end{pmatrix}
```

```math
=\begin{pmatrix}
m_1\cos{\theta} & m_1\cos{\theta}\tan{\phi}-(-1)^km_2\sin{\theta} \\
m_1\sin{\theta} & m_1\sin{\theta}\tan{\phi}+(-1)^km_2\cos{\theta}
\end{pmatrix}
```

以上で、行列Aによって座標がどのように変化するか可視化できたので、**以降は簡単のため**

```math
a=m_1\cos{\theta} \\
b=m_1\cos{\theta}\tan{\phi}-(-1)^km_2\sin{\theta} \\
c=m_1\sin{\theta} \\
d=m_1\sin{\theta}\tan{\phi}+(-1)^km_2\cos{\theta}
```
とおき、

```math
\boldsymbol{A}=
\begin{pmatrix}
a & b \\
c & d
\end{pmatrix}
```
**と表すこととします**。

### 【参考】特異値分解による変数変換の可視化

上ではQR分解により変数変換を可視化しましたが。別の行列分解法である**特異値分解**によっても可視化することができます (むしろQR変換よりも一般的かもしれません。[参考](https://qiita.com/osadataku/items/cb778b49799e718c68b8))

[特異値分解](https://nisshingeppo.com/ai/singular-value-decomposition/)を用いると、行列**A** は以下のように分解する事ができます。

```math
\boldsymbol{A} = \boldsymbol{U}\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{V}^T\\
{\left\{\begin{array}{ll}
\boldsymbol{U}: \boldsymbol{A}\boldsymbol{A}^Tの固有ベクトルを列ベクトルとして並べた行列 (直交行列となる)\\
\boldsymbol{\Gamma}: \boldsymbol{A}\boldsymbol{A}^Tの固有値を並べた対角行列\\
\boldsymbol{V}: \boldsymbol{A}^T\boldsymbol{A}の固有ベクトルを列ベクトルとして並べた行列 (直交行列となる)
\end{array}\right.
}
```

特に今回のようにAが正方かつ正則な行列の場合、$\boldsymbol{\Gamma}^{\frac{1}{2}}$は$\boldsymbol{A}\boldsymbol{A}^T$の[固有値の平方根を並べた対角行列となります](https://ja.wikipedia.org/wiki/%E7%89%B9%E7%95%B0%E5%80%A4%E5%88%86%E8%A7%A3)

[先ほどのQR分解と同様]()に、直交行列UおよびVは左右反転行列**Mir**→回転行列**Rot**に分解する事ができ (反転は適用されないこともある)、対角行列

```math
\boldsymbol{\Gamma}^{\frac{1}{2}}=\begin{pmatrix} \Sigma_1 & 0 \\ 0 & \Sigma_2 \end{pmatrix}
```

は拡大行列**Mag**とみなすことができます。

よってこれらを特異値分解の式$\boldsymbol{A} = \boldsymbol{U}\boldsymbol{\Gamma}^{\frac{1}{2}}\boldsymbol{V}^T$ に代入すると、

```math
\boldsymbol{A}=
\boldsymbol{Rot}\cdot\boldsymbol{Mir}\cdot\boldsymbol{Mag}\cdot\boldsymbol{Mir}\cdot\boldsymbol{Rot}
```

と「回転」→（「反転」）→「拡大」→（「反転」）→「回転」の組み合わせで表せる事が分かります。
($\boldsymbol{V}$ は転置をとっているので**Mir**と**Rot**を逆向きにしています)

さらに、最初の2つの「回転」→「反転」の操作では、円形の標準正規分布をスタートとすると形状が変化しないため、これらは無視する事ができます。

<img width=240 src=https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/9aa639fe-591d-fc71-6395-67fb9496353c.png>

また「拡大」では円形の線対称は崩れないので、4個目の「反転」も無視できます。

<img width=240 src=https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/70b95bf6-8e3f-2c0c-2a73-81008af0f554.png>

よってこれらを

~~「回転」→（「反転」）→~~「拡大」→（~~「反転」~~）→「回転」

のように除いた

**「拡大」→「回転」**

の2操作のみで表せる事が分かります
(特異値分解後の行列**Σ**が「拡大」を、**U**が「回転」を表します)

上記変換に「**平行移動**」を加えた下図のような変換を実施することで、QR分解と同様に2次元の正規分布を定義することができます。

<img width=480 src=https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/610167/8e84eac4-1ee4-e23b-5c84-1615fe3a6af9.png>

QR分解よりもシンプルな操作で表せる事が分かりました。

### 非独立変数の2次元正規分布
前節のように、変換前の標準正規分布で表せる座標系を**y**=(y1, y2)、変換後の座標系を**x**=(x1, x2)とすると、**x**は以下の式で表せます。

```math
\boldsymbol{x}=\boldsymbol{A}\boldsymbol{y}+\boldsymbol{\mu}
```
この式を変形して**y**を**x**の式で表すと、

```math
\boldsymbol{A}\boldsymbol{y}=\boldsymbol{x}-\boldsymbol{\mu}
```

```math
\boldsymbol{y}=\boldsymbol{A}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})
```
となります。

これを[前述の**y**の標準正規分布のベクトル式に代入すると](https://qiita.com/c60evaporator/items/d53053358105b0117f2c#1-2%E5%A4%89%E6%95%B0%E3%81%AE%E6%A8%99%E6%BA%96%E6%AD%A3%E8%A6%8F%E5%88%86%E5%B8%83)（<**y**,**y**>は**y**同士の内積を表します）

```math
q(\boldsymbol{y})=\frac{1}{2\pi}\exp\bigl(-\frac{1}{2}\langle\boldsymbol{y},\boldsymbol{y}\rangle\bigr)\\
=\frac{1}{2\pi}\exp\bigl(-\frac{1}{2}\langle\boldsymbol{A}^{-1}(\boldsymbol{x}-\boldsymbol{\mu}),\boldsymbol{A}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\rangle\bigr)\\
```
[転置と内積の関係](https://mathlandscape.com/transposed-matrix/)$\langle\boldsymbol{A}\boldsymbol{v},\boldsymbol{w}\rangle=\langle\boldsymbol{v},\boldsymbol{A}^T\boldsymbol{w}\rangle$を適用して

```math
=\frac{1}{2\pi}\exp\bigl(-\frac{1}{2}\langle(\boldsymbol{x}-\boldsymbol{\mu}),(\boldsymbol{A}^{-1})^T\boldsymbol{A}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\rangle\bigr)\\
```
[転置と逆行列の関係](https://mathlandscape.com/transposed-matrix/)$(\boldsymbol{A}^{-1})^T=(\boldsymbol{A}^T)^{-1}$を適用して

```math
=\frac{1}{2\pi}\exp\bigl(-\frac{1}{2}\langle(\boldsymbol{x}-\boldsymbol{\mu}),(\boldsymbol{A}^T)^{-1}\boldsymbol{A}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\rangle\bigr)\\
```
[逆行列同士の積の関係](https://risalc.info/src/inverse-matrix.html#prod)$\boldsymbol{A}^{-1}\boldsymbol{B}^{-1}=(\boldsymbol{B}\boldsymbol{A})^{-1}$を適用して

```math
=\frac{1}{2\pi}\exp\bigl(-\frac{1}{2}\langle(\boldsymbol{x}-\boldsymbol{\mu}),(\boldsymbol{A}\boldsymbol{A}^T)^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\rangle\bigr)\\
```
内積の表記法を$\langle\boldsymbol{v},\boldsymbol{w}\rangleから\boldsymbol{v}^T\boldsymbol{w}$（転置したベクトルとの積）に変更して

```math
=\frac{1}{2\pi}\exp\bigl(-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^T(\boldsymbol{A}\boldsymbol{A}^T)^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\bigr)\cdots④
```

#### ・変数変換後の積分が1となるよう規格化
ここで、変数変換後も確率密度関数の「積分値=1」の条件を満たすよう**規格化**する必要があります。
すなわち、上記④=`q(y1, y2)`に対して、変数変換後の確率密度関数`p(x1, x2)`が

```math
p(x_1, x_2)=Cq(y_1, y_2)
```
を満たすような規格化定数Cを求めます。

`q(y1, y2)`を`p(x1, x2)`に変数変換したとき、[積分値はヤコビ行列（正方行列で変数変換する場合、その行列そのもの）を掛ける必要がある](https://www.momoyama-usagi.com/entry/math-analysis24)、すなわち$dx_1dx_2=\bigl|\boldsymbol{A}\bigr|dy_1dy_2$となります。

変数変換後の確率密度関数`p(x1, x2)`の全定義域での積分に、上記ヤコビ行列による変換を適用すると

```math
\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}p(x_1,x_2)dx_1dx_2\\
=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}Cq(y_1,y_2)\bigl|\boldsymbol{A}\bigr|dy_1dy_2\\
=C\bigl|\boldsymbol{A}\bigr|\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}q(y_1,y_2)dy_1dy_2
```
ここで、`q(y1, y2)`の全定義域での積分$\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}q(y_1,y_2)dy_1dy_2=1$より

```math
=C\bigl|\boldsymbol{A}\bigr|
```
確率密度関数`p(x1, x2)`の全定義域での積分は1となるので、上式=1となる必要があるため、

```math
C\bigl|\boldsymbol{A}\bigr|=1\\
\therefore C=\frac{1}{\bigl|\boldsymbol{A}\bigr|}
```
これを④に適用して、

```math
p(\boldsymbol{x})=p(x_1, x_2)=\frac{1}{\bigl|\boldsymbol{A}\bigr|}q(y_1, y_2)\\
=\frac{1}{2\pi}\frac{1}{\bigl|\boldsymbol{A}\bigr|}\exp\bigl(-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^T(\boldsymbol{A}\boldsymbol{A}^T)^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\bigr)\cdots⑤
```
※上記では計算で規格化定数を求めましたが、[こちらの記事](https://www.geisya.or.jp/~mwm48961/electro/multi_integral5.htm)のように、変数変換により**行列式で表される平行四辺形分だけ面積が引き伸ばされる**ため、積分値(体積)=1を維持するためには引き伸ばした分高さを縮めるため、p(**x**)を**行列式で割る必要**がある、と理解すると直感的に分かりやすいかと思います。

#### ・行列式の変換
天下り的ですが、$\boldsymbol{A}\boldsymbol{A}^T=\boldsymbol{\Sigma}$とおきます。
このとき、

```math
\bigl|\boldsymbol{\Sigma}\bigr|=\bigl|\boldsymbol{A}\boldsymbol{A}^T\bigr|
```
[積の行列式](https://w3e.kanazawa-it.ac.jp/math/category/gyouretu/senkeidaisu/henkan-tex.cgi?target=/math/category/gyouretu/senkeidaisu/seki_no_gyouretusiki.html)$\bigl|\boldsymbol{A}\boldsymbol{B}\bigr|=\bigl|\boldsymbol{A}\bigr|\bigl|\boldsymbol{B}\bigr|$より

```math
=\bigl|\boldsymbol{A}\bigr|\bigl|\boldsymbol{A}^T\bigr|
```
[転置行列の行列式](https://w3e.kanazawa-it.ac.jp/math/category/gyouretu/senkeidaisu/gyouretusiki_no_seisitu/henkan-tex.cgi?target=/math/category/gyouretu/senkeidaisu/gyouretusiki_no_seisitu/gyouretusiki_no_tenti.html)$\bigl|\boldsymbol{A}^T\bigr|=\bigl|\boldsymbol{A}\bigr|$より

```math
=\bigl|\boldsymbol{A}\bigr|\bigl|\boldsymbol{A}\bigr|\\
=\bigl|\boldsymbol{A}\bigr|^2\\
\therefore \bigl|\boldsymbol{A}\bigr|=\bigl|\boldsymbol{\Sigma}\bigr|^\frac{1}{2}
```
となります。これを⑤に代入して

```math
p(\boldsymbol{x})=\frac{1}{2\pi}\frac{1}{|\boldsymbol{\Sigma}|^\frac{1}{2}}\exp\bigl(-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\bigr)
```

これは[冒頭の多次元正規分布の式](https://qiita.com/c60evaporator/items/d53053358105b0117f2c#%E5%A4%9A%E6%AC%A1%E5%85%83%E6%AD%A3%E8%A6%8F%E5%88%86%E5%B8%83%E3%81%AE%E5%BC%8F%E8%A4%87%E9%9B%91%E3%81%99%E3%81%8E%E5%95%8F%E9%A1%8C)のn=2の場合（nは変数の数）と等しくなることが分かります。
あとは天下り的に定義した$\boldsymbol{\Sigma}=\boldsymbol{A}\boldsymbol{A}^T$が分散共分散行列と等しくなることを証明すれば完了です。

なお、ここまでの式変形に使用した行列の公式は全て2次元に限らない正則行列に対して成立するため、
上式は**3変数以上の正規分布についても成立**します(nの値が変わる事にご注意下さい)

### 変数変換行列と分散共分散行列
前節で定義した$\boldsymbol{\Sigma}=\boldsymbol{A}\boldsymbol{A}^T$ですが、

これが分散共分散行列

```math
\boldsymbol{\Sigma}=
\begin{pmatrix}
\sigma_1 & \sigma_{12} \\
\sigma_{12} & \sigma_2
\end{pmatrix}
```

```math
{\left\{\begin{array}{ll}
\sigma_1 \cdots x_1の分散 \\
\sigma_2 \cdots x_2の分散 \\
\sigma_{12} \cdots x_1とx_2の共分散\\
\end{array}\right.
}
```
と等しいことを証明します

#### ・前節の行列Σの各成分の値
まずは**Σ**の具体的な値を求めてみます。

```math
\boldsymbol{\Sigma}=\boldsymbol{A}\boldsymbol{A}^T\\
=
\begin{pmatrix}
a & b \\
c & d
\end{pmatrix}
\begin{pmatrix}
a & c \\
b & d
\end{pmatrix}\\
=
\begin{pmatrix}
a^2+b^2 & ac+bd \\
ac+bd & c^2+d^2
\end{pmatrix}
```
この各成分が、分散共分散行列の各成分と等しいことを示します。

#### ・分散共分散行列の各成分の値
[前述の変数変換の式](https://qiita.com/c60evaporator/items/d53053358105b0117f2c#2%E6%AC%A1%E5%85%83%E6%AD%A3%E8%A6%8F%E5%88%86%E5%B8%83%E3%81%A8%E5%A4%89%E6%95%B0%E5%A4%89%E6%8F%9B)

```math
\left(\begin{array}{c}
x_1 \\
x_2 \\
\end{array}\right)
=
\begin{pmatrix}
a & b \\
c & d
\end{pmatrix}
\left(\begin{array}{c}
y_1 \\
y_2 \\
\end{array}\right)
```
から、

```math
{\left\{\begin{array}{ll}
x_1=ay_1+by_2 \\
x_2=cy_1+dy_2 \\
\end{array}\right.
}
```
となることを踏まえ、変換後の変数x1,x2に対する分散、共分散を求めます。

##### x1の分散
まずx1の分散$\sigma_1$を求めると

```math
\sigma_1=V[x_1]=V[ay_1+by_2]
```
[和の分散の公式](https://qiita.com/c60evaporator/items/93295d3890c9e20e521d#%E5%92%8C%E3%81%AE%E5%88%86%E6%95%A3)より

```math
=V[ay_1]+V[by_2]+\rm{Cov}[ay_1,by_2]
```
定数倍の[分散](https://qiita.com/c60evaporator/items/93295d3890c9e20e521d#%E5%AE%9A%E6%95%B0%E5%80%8D%E3%81%AE%E5%88%86%E6%95%A3)および[共分散](https://qiita.com/c60evaporator/items/93295d3890c9e20e521d#%E5%AE%9A%E6%95%B0%E5%80%8D%E3%81%AE%E5%85%B1%E5%88%86%E6%95%A3)の公式より

```math
=a^2V[y_1]+b^2V[y_2]+ab\rm{Cov}[y_1,y_2]
```
[変数y1とy2は独立](https://qiita.com/c60evaporator/items/d53053358105b0117f2c#1-2%E5%A4%89%E6%95%B0%E3%81%AE%E6%A8%99%E6%BA%96%E6%AD%A3%E8%A6%8F%E5%88%86%E5%B8%83)なので共分散は0となり

```math
=a^2V[y_1]+b^2V[y_2]
```
変数y1とy2の確率分布q(y1,y2)は標準正規分布なので、分散は1となり

```math
=a^2+b^2
```
よって$V[x_1]=a^2+b^2$となることが分かります

##### x2の分散
x2の分散$\sigma_1$は、x1のときと同様に

```math
\sigma_2=V[x_2]=V[cy_1+dy_2]\\
=V[cy_1]+V[dy_2]+\rm{Cov}[cy_1,dy_2]\\
=c^2V[y_1]+d^2V[y_2]+cd\rm{Cov}[y_1,y_2]\\
=c^2V[y_1]+d^2V[y_2]\\
=c^2+d^2
```
よって$V[x_2]=c^2+d^2$となることが分かります

##### x1とx2の共分散
x2の分散$\sigma_1$は、x1のときと同様に

```math
\sigma_{12}=\rm{Cov}[x_1,x_2]\\
=E\bigl[(x_1-E[x_1])(x_2-E[x_2]) \bigr]
```
[共分散の計算法の公式](https://qiita.com/c60evaporator/items/93295d3890c9e20e521d#%E5%85%B1%E5%88%86%E6%95%A3%E3%81%AE%E8%A8%88%E7%AE%97%E6%B3%95)より

```math
=E[x_1x_2]-E[x_1]E[x_2]\\
=E[(ay_1+by_2)(cy_1+dy_2)]-E[ay_1+by_2]E[cy_1+dy_2]\\
=acE[y_1^2]+(ad+bc)E[y_1y_2]+bdE[y_2^2]-(aE[y_1]+bE[y_2])(cE[y_1]+dE[y_2])
```
ここで、確率分布q(y1,y2)は標準正規分布なので、y1およびy2の平均$E[y_1]=E[y_2]=0$となるので

```math
=acE[y_1^2]+(ad+bc)E[y_1y_2]+bdE[y_2^2]
```
変数y1とy2は独立なので、$E[y_1y_2]=E[y_1]E[y_2]=0$が成り立つので

```math
=acE[y_1^2]+bdE[y_2^2]
```
また、[分散の公式](https://qiita.com/c60evaporator/items/93295d3890c9e20e521d#%E5%88%86%E6%95%A3)$V[f(x)]=E[f(x)^2]-E[f(x)]^2より、E[y_1^2]=V[y_1]+E[y_1]^2=V[y_1]=1$から（y2についても同様）

```math
=ac\cdot1+bd\cdot1\\
=ac+bd
```

よって

```math
{\left\{\begin{array}{ll}
\sigma_1=a^2+b^2 \\
\sigma_2=c^2+d^2 \\
\sigma_{12}=ac+bd \\
\end{array}\right.
}
```
となり、前節の行列Σと各成分が等しい、すなわち**前節の行列**$\boldsymbol{\Sigma}=\boldsymbol{A}\boldsymbol{A}^T$**は分散共分散行列である**ことが分かりました。

また、共分散の式$\sigma_{12}=ac+bd$を見ると分かりますが、共分散は負になることもあります（正負のイメージは[こちら](https://qiita.com/c60evaporator/items/93295d3890c9e20e521d#%E5%85%B1%E5%88%86%E6%95%A3%E3%81%AE%E7%9B%B4%E6%84%9F%E7%9A%84%E3%81%AA%E3%82%A4%E3%83%A1%E3%83%BC%E3%82%B8)参照）

※ なお、[前述の](https://qiita.com/c60evaporator/items/d53053358105b0117f2c#1-2%E5%A4%89%E6%95%B0%E3%81%AE%E6%A8%99%E6%BA%96%E6%AD%A3%E8%A6%8F%E5%88%86%E5%B8%83)標準正規分布の分散共分散行列**Σ**は、単位行列**I**となります。

### 3次元以上の場合
QR分解を始めとして、ここまで使用した行列の公式は全て2次元に限らない正則行列に対して成立するため、
分散共分散行列$\boldsymbol{\Sigma}=\boldsymbol{A}\boldsymbol{A}^T$は、**3次元以上の正規分布に関しても成立**します。

3次元以上においても直交行列は回転と反転を、上三角行列は拡大とスキューを表すので、行列Aはこれら（＋平行移動）を組み合わせた変数変換行列となります。

このとき、分散共分散行列$\boldsymbol{\Sigma}$は各変数の分散`V`および共分散`Cov`を用いて

```math
\boldsymbol{\Sigma}=
\left(
\begin{array}{ccc}
 {\rm V}[x_1] & {\rm Cov}[x_1,x_2] & \ldots &  {\rm Cov}[x_1,x_n]\\
 {\rm Cov}[x_2,x_1] & {\rm V}[x_2] & \ldots & {\rm Cov}[x_2,x_n] \\
\vdots & \vdots & \ddots & \vdots \\
{\rm Cov}[x_n,x_1] & {\rm Cov}[x_n,x_2] & \ldots & {\rm V}[x_n]
\end{array}
\right)
```
のように求められます。

また行列$\boldsymbol{A}$は、n次元標準正規分布

```math
q_1(y_1)=\frac{1}{\sqrt{2\pi}}\exp(-\frac{1}{2}y_1^2)\\
q_2(y_2)=\frac{1}{\sqrt{2\pi}}\exp(-\frac{1}{2}y_2^2)\\
:
q_n(y_n)=\frac{1}{\sqrt{2\pi}}\exp(-\frac{1}{2}y_n^2)
```

```math
q(y_1,y_2 \cdots,y_n)=q_1(y_1)q_2(y_2)\cdots q_n(y_n)=\frac{1}{(2\pi)^{ \frac{n}{2}} }\exp\bigl(-\frac{1}{2}(y_1^2+y_2^2\cdots +y_n^2)\bigr)
```
で表せる独立な変数**y**=(y1,y2..,yn)を、実際に使用する変数**x**=(x1,x2..,xn)に変換するための変換行列です。

<br>

上記分散共分散行列**Σ**および平均値ベクトル**μ**を用いて、多次元正規分布$p(\boldsymbol{x} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma})$は以下のように表されます（冒頭の式と同じ）

```math
p(\boldsymbol{x} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma})=\frac{1}{(2\pi)^\frac{n}{2}}\frac{1}{|\boldsymbol{\Sigma}|^\frac{1}{2}} \exp \Bigl(-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^T \boldsymbol{\Sigma}^{-1} (\boldsymbol{x}-\boldsymbol{\mu}) \Bigr)
```



# まとめ
**・多次元正規分布は正規分布を2次元以上の確率変数に拡張したもの**
**・多次元正規分布は、独立変数による標準正規分布を行列Aで変数変換することで表現できる**
**・変数変換行列Aは、スキュー、拡大、反転、回転、平行移動の組み合わせで表される**
**・分散共分散行列Σは、変数変換行列Aを使用して$\boldsymbol{\Sigma}=\boldsymbol{A}\boldsymbol{A}^T$と表せる**
**・多次元正規分布の確率密度関数は、分散共分散行列Σおよび平均ベクトルμから以下のように表せる**

```math
p(\boldsymbol{x})=\frac{1}{(2\pi)^\frac{n}{2}}\frac{1}{|\boldsymbol{\Sigma}|^\frac{1}{2}} \exp \Bigl(-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^T \boldsymbol{\Sigma}^{-1} (\boldsymbol{x}-\boldsymbol{\mu}) \Bigr)
```

# 出典
**・統計学入門　東京大学出版会**

http://www.utp.or.jp/book/b300857.html

**・英語版Wikipedia(他次元正規分布)**

https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Bivariate_case