多変量正規分布（多次元ガウス分布）は、機械学習・統計学・信号処理において最も重要な確率分布の一つです。その特に強力な性質として、「多変量正規分布を変数のグループに分割したとき、周辺分布も条件付き分布も正規分布になる」というものがあります。
この性質は、ガウス過程回帰、カルマンフィルタ、ベイズ線形回帰など、多くのアルゴリズムの理論的基盤を提供しています。
本記事では、ブロック行列の逆行列とシューア補行列を使って、周辺分布と条件付き分布の公式を完全に導出します。
本記事の内容
- 多変量正規分布の分割表現
- 周辺分布が正規分布になることの導出
- 条件付き分布の導出（シューア補行列を使った完全な証明）
- 条件付き平均 $\bm{\mu}_{a|b}$ と条件付き共分散 $\bm{\Sigma}_{a|b}$ の公式
- ガウス過程回帰への接続
- Pythonで2変量正規分布の可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
多変量正規分布の分割表現
分割の設定
$n$ 次元の確率変数ベクトル $\bm{x}$ を、$n_a$ 次元の $\bm{x}_a$ と $n_b$ 次元の $\bm{x}_b$（$n_a + n_b = n$）に分割します。
$$ \bm{x} = \begin{pmatrix} \bm{x}_a \\ \bm{x}_b \end{pmatrix} $$
対応して、平均ベクトルと共分散行列も分割します。
$$ \bm{\mu} = \begin{pmatrix} \bm{\mu}_a \\ \bm{\mu}_b \end{pmatrix}, \quad \bm{\Sigma} = \begin{pmatrix} \bm{\Sigma}_{aa} & \bm{\Sigma}_{ab} \\ \bm{\Sigma}_{ba} & \bm{\Sigma}_{bb} \end{pmatrix} $$
ここで、
- $\bm{\mu}_a \in \mathbb{R}^{n_a}$, $\bm{\mu}_b \in \mathbb{R}^{n_b}$: 各グループの平均ベクトル
- $\bm{\Sigma}_{aa} \in \mathbb{R}^{n_a \times n_a}$: $\bm{x}_a$ の共分散行列
- $\bm{\Sigma}_{bb} \in \mathbb{R}^{n_b \times n_b}$: $\bm{x}_b$ の共分散行列
- $\bm{\Sigma}_{ab} = \bm{\Sigma}_{ba}^T \in \mathbb{R}^{n_a \times n_b}$: $\bm{x}_a$ と $\bm{x}_b$ の交差共分散行列
多変量正規分布の確率密度関数
$$ p(\bm{x}) = \frac{1}{(2\pi)^{n/2} |\bm{\Sigma}|^{1/2}} \exp\left( -\frac{1}{2} (\bm{x} – \bm{\mu})^T \bm{\Sigma}^{-1} (\bm{x} – \bm{\mu}) \right) $$
シューア補行列
条件付き分布の導出で重要な役割を果たすのがシューア補行列（Schur complement）です。
定義
ブロック行列 $\bm{M} = \begin{pmatrix} \bm{A} & \bm{B} \\ \bm{C} & \bm{D} \end{pmatrix}$ に対して、
- $\bm{D}$ に関する $\bm{M}$ のシューア補行列: $\bm{M}/\bm{D} = \bm{A} – \bm{B}\bm{D}^{-1}\bm{C}$
- $\bm{A}$ に関する $\bm{M}$ のシューア補行列: $\bm{M}/\bm{A} = \bm{D} – \bm{C}\bm{A}^{-1}\bm{B}$
ブロック行列の逆行列
$\bm{\Sigma}$ の逆行列（精度行列）$\bm{\Lambda} = \bm{\Sigma}^{-1}$ をブロック形式で求めます。
$$ \bm{\Lambda} = \bm{\Sigma}^{-1} = \begin{pmatrix} \bm{\Sigma}_{aa} & \bm{\Sigma}_{ab} \\ \bm{\Sigma}_{ba} & \bm{\Sigma}_{bb} \end{pmatrix}^{-1} = \begin{pmatrix} \bm{\Lambda}_{aa} & \bm{\Lambda}_{ab} \\ \bm{\Lambda}_{ba} & \bm{\Lambda}_{bb} \end{pmatrix} $$
ブロック行列の逆行列の公式を導出します。$\bm{\Sigma}\bm{\Lambda} = \bm{I}$ より、
$$ \begin{pmatrix} \bm{\Sigma}_{aa} & \bm{\Sigma}_{ab} \\ \bm{\Sigma}_{ba} & \bm{\Sigma}_{bb} \end{pmatrix} \begin{pmatrix} \bm{\Lambda}_{aa} & \bm{\Lambda}_{ab} \\ \bm{\Lambda}_{ba} & \bm{\Lambda}_{bb} \end{pmatrix} = \begin{pmatrix} \bm{I}_{n_a} & \bm{0} \\ \bm{0} & \bm{I}_{n_b} \end{pmatrix} $$
これを4つのブロック方程式に分解します。
$$ \begin{align} \bm{\Sigma}_{aa}\bm{\Lambda}_{aa} + \bm{\Sigma}_{ab}\bm{\Lambda}_{ba} &= \bm{I}_{n_a} \tag{1} \\ \bm{\Sigma}_{aa}\bm{\Lambda}_{ab} + \bm{\Sigma}_{ab}\bm{\Lambda}_{bb} &= \bm{0} \tag{2} \\ \bm{\Sigma}_{ba}\bm{\Lambda}_{aa} + \bm{\Sigma}_{bb}\bm{\Lambda}_{ba} &= \bm{0} \tag{3} \\ \bm{\Sigma}_{ba}\bm{\Lambda}_{ab} + \bm{\Sigma}_{bb}\bm{\Lambda}_{bb} &= \bm{I}_{n_b} \tag{4} \end{align} $$
(3) より $\bm{\Lambda}_{ba} = -\bm{\Sigma}_{bb}^{-1}\bm{\Sigma}_{ba}\bm{\Lambda}_{aa}$ を (1) に代入:
$$ \begin{align} \bm{\Sigma}_{aa}\bm{\Lambda}_{aa} – \bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}\bm{\Sigma}_{ba}\bm{\Lambda}_{aa} &= \bm{I}_{n_a} \\ (\bm{\Sigma}_{aa} – \bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}\bm{\Sigma}_{ba})\bm{\Lambda}_{aa} &= \bm{I}_{n_a} \end{align} $$
したがって、
$$ \boxed{\bm{\Lambda}_{aa} = (\bm{\Sigma}_{aa} – \bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}\bm{\Sigma}_{ba})^{-1} = (\bm{\Sigma}/\bm{\Sigma}_{bb})^{-1}} $$
同様に (2) より $\bm{\Lambda}_{ab} = -\bm{\Sigma}_{aa}^{-1}\bm{\Sigma}_{ab}\bm{\Lambda}_{bb}$ を (4) に代入:
$$ \begin{align} -\bm{\Sigma}_{ba}\bm{\Sigma}_{aa}^{-1}\bm{\Sigma}_{ab}\bm{\Lambda}_{bb} + \bm{\Sigma}_{bb}\bm{\Lambda}_{bb} &= \bm{I}_{n_b} \\ (\bm{\Sigma}_{bb} – \bm{\Sigma}_{ba}\bm{\Sigma}_{aa}^{-1}\bm{\Sigma}_{ab})\bm{\Lambda}_{bb} &= \bm{I}_{n_b} \end{align} $$
$$ \boxed{\bm{\Lambda}_{bb} = (\bm{\Sigma}_{bb} – \bm{\Sigma}_{ba}\bm{\Sigma}_{aa}^{-1}\bm{\Sigma}_{ab})^{-1} = (\bm{\Sigma}/\bm{\Sigma}_{aa})^{-1}} $$
非対角ブロックは:
$$ \bm{\Lambda}_{ab} = -\bm{\Lambda}_{aa}\bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}, \quad \bm{\Lambda}_{ba} = -\bm{\Lambda}_{bb}\bm{\Sigma}_{ba}\bm{\Sigma}_{aa}^{-1} $$
行列式の公式
ブロック行列の行列式は、シューア補行列を用いて
$$ |\bm{\Sigma}| = |\bm{\Sigma}_{bb}| \cdot |\bm{\Sigma}_{aa} – \bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}\bm{\Sigma}_{ba}| $$
と書けます。この公式は後の導出で使います。
周辺分布の導出
定理
$\bm{x} = (\bm{x}_a^T, \bm{x}_b^T)^T \sim \mathcal{N}(\bm{\mu}, \bm{\Sigma})$ のとき、
$$ \boxed{\bm{x}_a \sim \mathcal{N}(\bm{\mu}_a, \bm{\Sigma}_{aa})} $$
$$ \boxed{\bm{x}_b \sim \mathcal{N}(\bm{\mu}_b, \bm{\Sigma}_{bb})} $$
証明
周辺分布 $p(\bm{x}_a)$ は、同時分布 $p(\bm{x}_a, \bm{x}_b)$ を $\bm{x}_b$ について積分することで得られます。
$$ p(\bm{x}_a) = \int p(\bm{x}_a, \bm{x}_b) \, d\bm{x}_b $$
指数部分の二次形式を $\bm{x}_b$ について平方完成します。$\bm{y}_a = \bm{x}_a – \bm{\mu}_a$, $\bm{y}_b = \bm{x}_b – \bm{\mu}_b$ と変数変換すると、
$$ Q = (\bm{x} – \bm{\mu})^T \bm{\Sigma}^{-1} (\bm{x} – \bm{\mu}) = \begin{pmatrix} \bm{y}_a \\ \bm{y}_b \end{pmatrix}^T \begin{pmatrix} \bm{\Lambda}_{aa} & \bm{\Lambda}_{ab} \\ \bm{\Lambda}_{ba} & \bm{\Lambda}_{bb} \end{pmatrix} \begin{pmatrix} \bm{y}_a \\ \bm{y}_b \end{pmatrix} $$
展開すると、
$$ Q = \bm{y}_a^T \bm{\Lambda}_{aa} \bm{y}_a + 2\bm{y}_a^T \bm{\Lambda}_{ab} \bm{y}_b + \bm{y}_b^T \bm{\Lambda}_{bb} \bm{y}_b $$
$\bm{y}_b$ について平方完成します。
$$ \begin{align} Q &= \bm{y}_b^T \bm{\Lambda}_{bb} \bm{y}_b + 2\bm{y}_a^T \bm{\Lambda}_{ab} \bm{y}_b + \bm{y}_a^T \bm{\Lambda}_{aa} \bm{y}_a \\ &= (\bm{y}_b + \bm{\Lambda}_{bb}^{-1}\bm{\Lambda}_{ba}\bm{y}_a)^T \bm{\Lambda}_{bb} (\bm{y}_b + \bm{\Lambda}_{bb}^{-1}\bm{\Lambda}_{ba}\bm{y}_a) \\ &\quad + \bm{y}_a^T (\bm{\Lambda}_{aa} – \bm{\Lambda}_{ab}\bm{\Lambda}_{bb}^{-1}\bm{\Lambda}_{ba}) \bm{y}_a \end{align} $$
ここで $\bm{\Lambda}_{ab} = \bm{\Lambda}_{ba}^T$ を使いました。平方完成のステップを詳しく示します。
$\bm{y}_b$ の二次形式 $\bm{y}_b^T \bm{\Lambda}_{bb} \bm{y}_b + 2\bm{y}_a^T \bm{\Lambda}_{ab} \bm{y}_b$ を完成させます。
$$ \begin{align} &\bm{y}_b^T \bm{\Lambda}_{bb} \bm{y}_b + 2\bm{y}_a^T \bm{\Lambda}_{ab} \bm{y}_b \\ &= \bm{y}_b^T \bm{\Lambda}_{bb} \bm{y}_b + \bm{y}_b^T \bm{\Lambda}_{ba}\bm{y}_a + \bm{y}_a^T \bm{\Lambda}_{ab} \bm{y}_b \\ &= (\bm{y}_b + \bm{\Lambda}_{bb}^{-1}\bm{\Lambda}_{ba}\bm{y}_a)^T \bm{\Lambda}_{bb} (\bm{y}_b + \bm{\Lambda}_{bb}^{-1}\bm{\Lambda}_{ba}\bm{y}_a) – \bm{y}_a^T \bm{\Lambda}_{ab}\bm{\Lambda}_{bb}^{-1}\bm{\Lambda}_{ba}\bm{y}_a \end{align} $$
したがって、
$$ Q = \underbrace{(\bm{y}_b + \bm{\Lambda}_{bb}^{-1}\bm{\Lambda}_{ba}\bm{y}_a)^T \bm{\Lambda}_{bb} (\bm{y}_b + \bm{\Lambda}_{bb}^{-1}\bm{\Lambda}_{ba}\bm{y}_a)}_{\text{(I): } \bm{y}_b \text{に依存}} + \underbrace{\bm{y}_a^T (\bm{\Lambda}_{aa} – \bm{\Lambda}_{ab}\bm{\Lambda}_{bb}^{-1}\bm{\Lambda}_{ba}) \bm{y}_a}_{\text{(II): } \bm{y}_a \text{のみに依存}} $$
$\bm{x}_b$ について積分すると、(I) の部分がガウス積分になります。$\bm{z} = \bm{y}_b + \bm{\Lambda}_{bb}^{-1}\bm{\Lambda}_{ba}\bm{y}_a$ と置換すると（ヤコビアンは $1$）、
$$ \int \exp\left(-\frac{1}{2}\bm{z}^T\bm{\Lambda}_{bb}\bm{z}\right) d\bm{z} = (2\pi)^{n_b/2} |\bm{\Lambda}_{bb}|^{-1/2} $$
したがって、
$$ p(\bm{x}_a) = \frac{(2\pi)^{n_b/2} |\bm{\Lambda}_{bb}|^{-1/2}}{(2\pi)^{n/2} |\bm{\Sigma}|^{1/2}} \exp\left(-\frac{1}{2}\bm{y}_a^T (\bm{\Lambda}_{aa} – \bm{\Lambda}_{ab}\bm{\Lambda}_{bb}^{-1}\bm{\Lambda}_{ba}) \bm{y}_a\right) $$
正規化定数を整理します。$n = n_a + n_b$ なので $(2\pi)^{n_b/2} / (2\pi)^{n/2} = (2\pi)^{-n_a/2}$ です。
ブロック行列の逆行列から $\bm{\Lambda}_{bb} = (\bm{\Sigma}_{bb} – \bm{\Sigma}_{ba}\bm{\Sigma}_{aa}^{-1}\bm{\Sigma}_{ab})^{-1}$ であり、行列式の公式より
$$ |\bm{\Sigma}| = |\bm{\Sigma}_{aa}| \cdot |\bm{\Sigma}_{bb} – \bm{\Sigma}_{ba}\bm{\Sigma}_{aa}^{-1}\bm{\Sigma}_{ab}| = |\bm{\Sigma}_{aa}| \cdot |\bm{\Lambda}_{bb}|^{-1} $$
したがって $|\bm{\Lambda}_{bb}|^{-1/2} / |\bm{\Sigma}|^{1/2} = 1 / |\bm{\Sigma}_{aa}|^{1/2}$ です。
また、$\bm{\Lambda}_{aa} – \bm{\Lambda}_{ab}\bm{\Lambda}_{bb}^{-1}\bm{\Lambda}_{ba} = \bm{\Sigma}_{aa}^{-1}$ であることを示します。
ブロック行列の逆行列の公式から、$\bm{\Lambda}_{aa} = (\bm{\Sigma}_{aa} – \bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}\bm{\Sigma}_{ba})^{-1}$ であり、$\bm{\Lambda}_{ab} = -\bm{\Lambda}_{aa}\bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}$, $\bm{\Lambda}_{ba} = -\bm{\Sigma}_{bb}^{-1}\bm{\Sigma}_{ba}\bm{\Lambda}_{aa}$ です。
$$ \begin{align} \bm{\Lambda}_{ab}\bm{\Lambda}_{bb}^{-1}\bm{\Lambda}_{ba} &= (-\bm{\Lambda}_{aa}\bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1})\bm{\Lambda}_{bb}^{-1}(-\bm{\Sigma}_{bb}^{-1}\bm{\Sigma}_{ba}\bm{\Lambda}_{aa}) \\ &= \bm{\Lambda}_{aa}\bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}(\bm{\Sigma}_{bb} – \bm{\Sigma}_{ba}\bm{\Sigma}_{aa}^{-1}\bm{\Sigma}_{ab})\bm{\Sigma}_{bb}^{-1}\bm{\Sigma}_{ba}\bm{\Lambda}_{aa} \end{align} $$
ここで $\bm{\Lambda}_{bb}^{-1} = \bm{\Sigma}_{bb} – \bm{\Sigma}_{ba}\bm{\Sigma}_{aa}^{-1}\bm{\Sigma}_{ab}$ を代入しました。これを展開すると長くなるため、別のアプローチを使います。
$(\bm{\Sigma}^{-1})$ 全体の逆行列の $(a,a)$ ブロックのシューア補行列を考えます。
$$ (\bm{\Lambda}/\bm{\Lambda}_{bb})^{-1} = (\bm{\Lambda}_{aa} – \bm{\Lambda}_{ab}\bm{\Lambda}_{bb}^{-1}\bm{\Lambda}_{ba})^{-1} $$
$\bm{\Lambda} = \bm{\Sigma}^{-1}$ なので、$(\bm{\Lambda}/\bm{\Lambda}_{bb})^{-1}$ は $\bm{\Sigma}$ の $(a,a)$ ブロック、すなわち $\bm{\Sigma}_{aa}$ に等しくなります。
これはブロック行列の逆行列の一般的な性質です。$\bm{M}^{-1}$ のブロック $(1,1)$ のシューア補行列は、$\bm{M}$ のブロック $(1,1)$ そのものです。
$$ \bm{\Lambda}_{aa} – \bm{\Lambda}_{ab}\bm{\Lambda}_{bb}^{-1}\bm{\Lambda}_{ba} = \bm{\Sigma}_{aa}^{-1} $$
以上をまとめると、
$$ p(\bm{x}_a) = \frac{1}{(2\pi)^{n_a/2}|\bm{\Sigma}_{aa}|^{1/2}} \exp\left(-\frac{1}{2}(\bm{x}_a – \bm{\mu}_a)^T \bm{\Sigma}_{aa}^{-1} (\bm{x}_a – \bm{\mu}_a)\right) $$
これは $\mathcal{N}(\bm{\mu}_a, \bm{\Sigma}_{aa})$ の密度関数です。 $\square$
条件付き分布の導出
定理
$$ \boxed{\bm{x}_a \mid \bm{x}_b \sim \mathcal{N}(\bm{\mu}_{a|b}, \bm{\Sigma}_{a|b})} $$
ここで、
$$ \boxed{\bm{\mu}_{a|b} = \bm{\mu}_a + \bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}(\bm{x}_b – \bm{\mu}_b)} $$
$$ \boxed{\bm{\Sigma}_{a|b} = \bm{\Sigma}_{aa} – \bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}\bm{\Sigma}_{ba}} $$
証明
条件付き分布は $p(\bm{x}_a \mid \bm{x}_b) = p(\bm{x}_a, \bm{x}_b) / p(\bm{x}_b)$ で定義されます。$p(\bm{x}_b)$ は $\bm{x}_a$ に依存しないので、$\bm{x}_a$ に関する関数形を見れば十分です。
同時分布の指数部分の二次形式を再び考えます。先ほどの平方完成の結果を使います。
$$ Q = (\bm{y}_b + \bm{\Lambda}_{bb}^{-1}\bm{\Lambda}_{ba}\bm{y}_a)^T \bm{\Lambda}_{bb} (\bm{y}_b + \bm{\Lambda}_{bb}^{-1}\bm{\Lambda}_{ba}\bm{y}_a) + \bm{y}_a^T \bm{\Sigma}_{aa}^{-1} \bm{y}_a $$
ただし、ここでは $\bm{x}_b$ を固定して $\bm{x}_a$ について平方完成したいので、別の分解を行います。
$Q$ を $\bm{y}_a$ について整理します。
$$ Q = \bm{y}_a^T \bm{\Lambda}_{aa} \bm{y}_a + 2\bm{y}_a^T \bm{\Lambda}_{ab}\bm{y}_b + \bm{y}_b^T \bm{\Lambda}_{bb} \bm{y}_b $$
$\bm{y}_a$ について平方完成します。
$$ \begin{align} Q &= (\bm{y}_a + \bm{\Lambda}_{aa}^{-1}\bm{\Lambda}_{ab}\bm{y}_b)^T \bm{\Lambda}_{aa} (\bm{y}_a + \bm{\Lambda}_{aa}^{-1}\bm{\Lambda}_{ab}\bm{y}_b) \\ &\quad + \bm{y}_b^T (\bm{\Lambda}_{bb} – \bm{\Lambda}_{ba}\bm{\Lambda}_{aa}^{-1}\bm{\Lambda}_{ab}) \bm{y}_b \end{align} $$
平方完成のステップを詳しく示します。
$$ \begin{align} &\bm{y}_a^T \bm{\Lambda}_{aa} \bm{y}_a + 2\bm{y}_a^T \bm{\Lambda}_{ab}\bm{y}_b \\ &= \bm{y}_a^T \bm{\Lambda}_{aa} \bm{y}_a + \bm{y}_a^T \bm{\Lambda}_{ab}\bm{y}_b + \bm{y}_b^T \bm{\Lambda}_{ba}\bm{y}_a \\ &= (\bm{y}_a + \bm{\Lambda}_{aa}^{-1}\bm{\Lambda}_{ab}\bm{y}_b)^T \bm{\Lambda}_{aa} (\bm{y}_a + \bm{\Lambda}_{aa}^{-1}\bm{\Lambda}_{ab}\bm{y}_b) – \bm{y}_b^T\bm{\Lambda}_{ba}\bm{\Lambda}_{aa}^{-1}\bm{\Lambda}_{ab}\bm{y}_b \end{align} $$
条件付き分布では第2項（$\bm{y}_b$ のみに依存する項）は正規化定数に吸収されるので、$\bm{x}_a$ に依存する部分だけを取り出すと、
$$ p(\bm{x}_a \mid \bm{x}_b) \propto \exp\left(-\frac{1}{2}(\bm{y}_a + \bm{\Lambda}_{aa}^{-1}\bm{\Lambda}_{ab}\bm{y}_b)^T \bm{\Lambda}_{aa} (\bm{y}_a + \bm{\Lambda}_{aa}^{-1}\bm{\Lambda}_{ab}\bm{y}_b)\right) $$
これは $\bm{x}_a$ に関する正規分布の形をしています。
条件付き共分散:
精度行列が $\bm{\Lambda}_{aa}$ なので、共分散行列は
$$ \bm{\Sigma}_{a|b} = \bm{\Lambda}_{aa}^{-1} = \bm{\Sigma}_{aa} – \bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}\bm{\Sigma}_{ba} $$
ここで $\bm{\Lambda}_{aa} = (\bm{\Sigma}_{aa} – \bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}\bm{\Sigma}_{ba})^{-1}$ を使いました。
条件付き平均:
$\bm{y}_a = \bm{x}_a – \bm{\mu}_a$, $\bm{y}_b = \bm{x}_b – \bm{\mu}_b$ を代入すると、中心は $\bm{y}_a = -\bm{\Lambda}_{aa}^{-1}\bm{\Lambda}_{ab}\bm{y}_b$ で、つまり
$$ \begin{align} \bm{x}_a – \bm{\mu}_a &= -\bm{\Lambda}_{aa}^{-1}\bm{\Lambda}_{ab}(\bm{x}_b – \bm{\mu}_b) \end{align} $$
$\bm{\Lambda}_{aa}^{-1}\bm{\Lambda}_{ab}$ を $\bm{\Sigma}$ の成分で表します。
$$ \bm{\Lambda}_{ab} = -\bm{\Lambda}_{aa}\bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1} $$
なので、
$$ -\bm{\Lambda}_{aa}^{-1}\bm{\Lambda}_{ab} = \bm{\Lambda}_{aa}^{-1}\bm{\Lambda}_{aa}\bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1} = \bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1} $$
したがって、
$$ \bm{\mu}_{a|b} = \bm{\mu}_a + \bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}(\bm{x}_b – \bm{\mu}_b) \qquad \square $$
結果のまとめと解釈
条件付き分布の公式をまとめます。
$$ \bm{x}_a \mid \bm{x}_b \sim \mathcal{N}(\bm{\mu}_{a|b}, \bm{\Sigma}_{a|b}) $$
条件付き平均:
$$ \bm{\mu}_{a|b} = \bm{\mu}_a + \bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}(\bm{x}_b – \bm{\mu}_b) $$
- $\bm{\mu}_a$: 事前の平均
- $\bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}$: 回帰係数行列
- $\bm{x}_b – \bm{\mu}_b$: 観測の偏差
- 条件付き平均は $\bm{x}_b$ の線形関数
条件付き共分散:
$$ \bm{\Sigma}_{a|b} = \bm{\Sigma}_{aa} – \bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}\bm{\Sigma}_{ba} $$
- $\bm{\Sigma}_{aa}$: 事前の共分散
- $\bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}\bm{\Sigma}_{ba}$: $\bm{x}_b$ の観測によって減少する不確実性
- 条件付き共分散は $\bm{x}_b$ の値に依存しない（正規分布の重要な性質）
具体例: 2変量正規分布
2変量正規分布 $\bm{x} = (x_1, x_2)^T$ で具体的に計算します。
$$ \bm{\mu} = \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \quad \bm{\Sigma} = \begin{pmatrix} 1 & 0.8 \\ 0.8 & 2 \end{pmatrix} $$
$x_2 = 2$ が観測されたときの $x_1$ の条件付き分布を求めます。
条件付き平均:
$$ \mu_{1|2} = 0 + 0.8 \cdot 2^{-1} \cdot (2 – 1) = 0.4 $$
条件付き分散:
$$ \sigma_{1|2}^2 = 1 – 0.8 \cdot 2^{-1} \cdot 0.8 = 1 – 0.32 = 0.68 $$
したがって $x_1 \mid x_2 = 2 \sim \mathcal{N}(0.4, 0.68)$ です。
ガウス過程回帰への接続
ガウス過程回帰の設定
ガウス過程回帰（Gaussian Process Regression, GPR）は、本記事の条件付き分布の公式をそのまま適用したものです。
$n$ 個の訓練データ $(\bm{X}, \bm{y})$ と $n_*$ 個のテストデータ $\bm{X}_*$ が与えられたとき、
$$ \begin{pmatrix} \bm{y} \\ \bm{f}_* \end{pmatrix} \sim \mathcal{N}\left( \begin{pmatrix} \bm{0} \\ \bm{0} \end{pmatrix}, \begin{pmatrix} \bm{K} + \sigma_n^2 \bm{I} & \bm{K}_* \\ \bm{K}_*^T & \bm{K}_{**} \end{pmatrix} \right) $$
ここで $\bm{K}$, $\bm{K}_*$, $\bm{K}_{**}$ はカーネル行列です。
条件付き分布の公式を適用すると、
$$ \bm{f}_* \mid \bm{y} \sim \mathcal{N}(\bm{\mu}_*, \bm{\Sigma}_*) $$
$$ \bm{\mu}_* = \bm{K}_*^T (\bm{K} + \sigma_n^2 \bm{I})^{-1} \bm{y} $$
$$ \bm{\Sigma}_* = \bm{K}_{**} – \bm{K}_*^T (\bm{K} + \sigma_n^2 \bm{I})^{-1} \bm{K}_* $$
これはまさに先ほどの $\bm{\mu}_{a|b}$ と $\bm{\Sigma}_{a|b}$ の公式そのものです。
Pythonでの可視化
2変量正規分布の周辺分布と条件付き分布
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal, norm
# パラメータの設定
mu = np.array([0, 1])
Sigma = np.array([[1.0, 0.8],
[0.8, 2.0]])
# 条件付き分布の計算
def conditional_normal(mu, Sigma, x_b, idx_a=0, idx_b=1):
"""
多変量正規分布の条件付き分布を計算
x_a | x_b の平均と分散を返す
"""
mu_a = mu[idx_a]
mu_b = mu[idx_b]
sigma_aa = Sigma[idx_a, idx_a]
sigma_bb = Sigma[idx_b, idx_b]
sigma_ab = Sigma[idx_a, idx_b]
# 条件付き平均
mu_cond = mu_a + sigma_ab / sigma_bb * (x_b - mu_b)
# 条件付き分散
sigma_cond = sigma_aa - sigma_ab**2 / sigma_bb
return mu_cond, sigma_cond
# メッシュグリッドの作成
x1 = np.linspace(-3, 3, 200)
x2 = np.linspace(-3, 5, 200)
X1, X2 = np.meshgrid(x1, x2)
pos = np.dstack((X1, X2))
# 同時分布の密度
rv = multivariate_normal(mu, Sigma)
Z = rv.pdf(pos)
# 条件付き分布の例（x2 = 2）
x2_cond = 2.0
mu_cond, sigma_cond = conditional_normal(mu, Sigma, x2_cond)
print(f"条件付き分布: x1 | x2={x2_cond}")
print(f" 平均: {mu_cond:.4f}")
print(f" 分散: {sigma_cond:.4f}")
print(f" 標準偏差: {np.sqrt(sigma_cond):.4f}")
# 可視化
fig = plt.figure(figsize=(14, 10))
# メインの等高線プロット
ax_main = fig.add_axes([0.15, 0.15, 0.55, 0.55])
ax_main.contourf(X1, X2, Z, levels=20, cmap='Blues', alpha=0.7)
ax_main.contour(X1, X2, Z, levels=10, colors='steelblue', linewidths=0.5)
# 条件付けの線
ax_main.axhline(x2_cond, color='red', linewidth=2, linestyle='--',
label=f'$x_2 = {x2_cond}$')
# 条件付き分布の平均を点で表示
ax_main.plot(mu_cond, x2_cond, 'ro', markersize=8)
ax_main.set_xlabel('$x_1$', fontsize=13)
ax_main.set_ylabel('$x_2$', fontsize=13)
ax_main.set_title('Bivariate normal distribution', fontsize=14)
ax_main.legend(fontsize=11, loc='lower right')
# 上側: x1の周辺分布
ax_top = fig.add_axes([0.15, 0.72, 0.55, 0.2])
# 周辺分布 p(x1)
marginal_x1 = norm(mu[0], np.sqrt(Sigma[0, 0]))
ax_top.plot(x1, marginal_x1.pdf(x1), 'b-', linewidth=2,
label=f'$p(x_1)$: $\\mathcal{{N}}({mu[0]}, {Sigma[0,0]})$')
# 条件付き分布 p(x1|x2=2)
cond_dist = norm(mu_cond, np.sqrt(sigma_cond))
ax_top.plot(x1, cond_dist.pdf(x1), 'r-', linewidth=2,
label=f'$p(x_1|x_2={x2_cond})$: $\\mathcal{{N}}({mu_cond:.2f}, {sigma_cond:.2f})$')
ax_top.fill_between(x1, cond_dist.pdf(x1), alpha=0.2, color='red')
ax_top.set_xlim(-3, 3)
ax_top.set_ylabel('Density', fontsize=11)
ax_top.legend(fontsize=10)
ax_top.set_title('Marginal and conditional distribution of $x_1$', fontsize=12)
# 右側: x2の周辺分布
ax_right = fig.add_axes([0.72, 0.15, 0.2, 0.55])
marginal_x2 = norm(mu[1], np.sqrt(Sigma[1, 1]))
ax_right.plot(marginal_x2.pdf(x2), x2, 'b-', linewidth=2,
label=f'$p(x_2)$')
ax_right.axhline(x2_cond, color='red', linewidth=1.5, linestyle='--')
ax_right.set_ylim(-3, 5)
ax_right.set_xlabel('Density', fontsize=11)
ax_right.legend(fontsize=10)
ax_right.set_title('Marginal $p(x_2)$', fontsize=12)
plt.show()
条件付き分布の変化
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
# パラメータの設定
mu = np.array([0, 1])
Sigma = np.array([[1.0, 0.8],
[0.8, 2.0]])
def conditional_normal_1d(mu, Sigma, x_b):
"""x1 | x2 = x_b の条件付き分布"""
mu_cond = mu[0] + Sigma[0, 1] / Sigma[1, 1] * (x_b - mu[1])
sigma_cond = Sigma[0, 0] - Sigma[0, 1]**2 / Sigma[1, 1]
return mu_cond, sigma_cond
# さまざまなx2の値での条件付き分布
x2_values = [-1, 0, 1, 2, 3]
x1_range = np.linspace(-4, 4, 300)
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# 左: 条件付き分布の密度関数
ax = axes[0]
colors = plt.cm.viridis(np.linspace(0, 1, len(x2_values)))
# 周辺分布をグレーで表示
marginal = norm(mu[0], np.sqrt(Sigma[0, 0]))
ax.plot(x1_range, marginal.pdf(x1_range), 'k--', linewidth=2,
alpha=0.5, label='Marginal $p(x_1)$')
for x2_val, color in zip(x2_values, colors):
mu_c, sigma_c = conditional_normal_1d(mu, Sigma, x2_val)
cond = norm(mu_c, np.sqrt(sigma_c))
ax.plot(x1_range, cond.pdf(x1_range), color=color, linewidth=2,
label=f'$x_2 = {x2_val}$: $\\mu = {mu_c:.2f}$')
ax.set_xlabel('$x_1$', fontsize=12)
ax.set_ylabel('Density', fontsize=12)
ax.set_title('Conditional distribution $p(x_1 | x_2)$', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
# 右: 相関係数の影響
ax = axes[1]
rho_values = [0.0, 0.3, 0.6, 0.9]
x2_fixed = 2.0
sigma1, sigma2 = 1.0, np.sqrt(2.0)
for rho in rho_values:
Sigma_rho = np.array([[sigma1**2, rho * sigma1 * sigma2],
[rho * sigma1 * sigma2, sigma2**2]])
mu_c = mu[0] + Sigma_rho[0, 1] / Sigma_rho[1, 1] * (x2_fixed - mu[1])
sigma_c = Sigma_rho[0, 0] - Sigma_rho[0, 1]**2 / Sigma_rho[1, 1]
cond = norm(mu_c, np.sqrt(sigma_c))
ax.plot(x1_range, cond.pdf(x1_range), linewidth=2,
label=f'$\\rho = {rho}$: $\\mu = {mu_c:.2f}$, $\\sigma^2 = {sigma_c:.2f}$')
marginal = norm(mu[0], sigma1)
ax.plot(x1_range, marginal.pdf(x1_range), 'k--', linewidth=2,
alpha=0.5, label='Marginal ($\\rho = 0$)')
ax.set_xlabel('$x_1$', fontsize=12)
ax.set_ylabel('Density', fontsize=12)
ax.set_title(f'Effect of correlation $\\rho$ ($x_2 = {x2_fixed}$)', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
3次元での可視化（等高線 + 条件付きスライス）
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal, norm
from mpl_toolkits.mplot3d import Axes3D
# パラメータの設定
mu = np.array([0, 1])
Sigma = np.array([[1.0, 0.8],
[0.8, 2.0]])
# メッシュグリッド
x1 = np.linspace(-3, 3, 150)
x2 = np.linspace(-3, 5, 150)
X1, X2 = np.meshgrid(x1, x2)
pos = np.dstack((X1, X2))
rv = multivariate_normal(mu, Sigma)
Z = rv.pdf(pos)
# 3Dプロット
fig = plt.figure(figsize=(14, 6))
# 左: 3D表面プロット
ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(X1, X2, Z, cmap='viridis', alpha=0.7,
edgecolor='none')
# 条件付き分布のスライス (x2=2)
x2_cond = 2.0
mu_c = mu[0] + Sigma[0, 1] / Sigma[1, 1] * (x2_cond - mu[1])
sigma_c = Sigma[0, 0] - Sigma[0, 1]**2 / Sigma[1, 1]
cond_pdf = norm(mu_c, np.sqrt(sigma_c)).pdf(x1)
# p(x1, x2=2) = p(x1|x2=2) * p(x2=2) のスケーリング
p_x2 = norm(mu[1], np.sqrt(Sigma[1, 1])).pdf(x2_cond)
joint_slice = cond_pdf * p_x2
ax1.plot(x1, np.full_like(x1, x2_cond), joint_slice,
color='red', linewidth=3, label=f'Slice at $x_2={x2_cond}$')
ax1.set_xlabel('$x_1$', fontsize=11)
ax1.set_ylabel('$x_2$', fontsize=11)
ax1.set_zlabel('$p(x_1, x_2)$', fontsize=11)
ax1.set_title('Joint density with conditional slice', fontsize=13)
ax1.view_init(elev=25, azim=-60)
# 右: 周辺分布と条件付き分布のまとめ
ax2 = fig.add_subplot(122)
# さまざまなx2の値で条件付き平均と95%信頼区間を表示
x2_range = np.linspace(-2, 4, 100)
mu_conds = mu[0] + Sigma[0, 1] / Sigma[1, 1] * (x2_range - mu[1])
sigma_cond = Sigma[0, 0] - Sigma[0, 1]**2 / Sigma[1, 1]
std_cond = np.sqrt(sigma_cond)
ax2.plot(x2_range, mu_conds, 'r-', linewidth=2,
label='Conditional mean $\\mu_{1|2}$')
ax2.fill_between(x2_range, mu_conds - 1.96 * std_cond,
mu_conds + 1.96 * std_cond,
alpha=0.2, color='red', label='95% credible interval')
ax2.fill_between(x2_range, mu_conds - std_cond,
mu_conds + std_cond,
alpha=0.3, color='red', label='68% credible interval')
# サンプル点を表示
np.random.seed(42)
samples = np.random.multivariate_normal(mu, Sigma, 200)
ax2.scatter(samples[:, 1], samples[:, 0], alpha=0.3, s=15, color='steelblue',
label='Samples')
ax2.axhline(mu[0], color='blue', linewidth=1, linestyle=':', alpha=0.5)
ax2.set_xlabel('$x_2$', fontsize=12)
ax2.set_ylabel('$x_1$', fontsize=12)
ax2.set_title('Conditional mean and credible interval', fontsize=13)
ax2.legend(fontsize=10, loc='upper left')
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
多変量の場合の数値検証
import numpy as np
from scipy.stats import multivariate_normal
# 4次元の多変量正規分布
np.random.seed(42)
n = 4
n_a = 2 # x_a は最初の2次元
n_b = 2 # x_b は残りの2次元
# 正定値対称行列を生成
A = np.random.randn(n, n)
Sigma = A @ A.T + np.eye(n)
mu = np.array([1.0, 2.0, -1.0, 0.5])
print("=== パラメータ ===")
print(f"mu = {mu}")
print(f"Sigma =\n{Sigma}")
# ブロックに分割
Sigma_aa = Sigma[:n_a, :n_a]
Sigma_ab = Sigma[:n_a, n_a:]
Sigma_ba = Sigma[n_a:, :n_a]
Sigma_bb = Sigma[n_a:, n_a:]
mu_a = mu[:n_a]
mu_b = mu[n_a:]
# 条件付き分布の解析解
x_b_obs = np.array([0.0, 1.0]) # 観測値
mu_cond = mu_a + Sigma_ab @ np.linalg.solve(Sigma_bb, x_b_obs - mu_b)
Sigma_cond = Sigma_aa - Sigma_ab @ np.linalg.solve(Sigma_bb, Sigma_ba)
print(f"\n=== 条件付き分布（解析解）===")
print(f"x_b の観測値: {x_b_obs}")
print(f"条件付き平均: mu_a|b = {mu_cond}")
print(f"条件付き共分散:\n{Sigma_cond}")
# モンテカルロ法による検証
n_samples = 1000000
samples = np.random.multivariate_normal(mu, Sigma, n_samples)
x_a_samples = samples[:, :n_a]
x_b_samples = samples[:, n_a:]
# x_b が x_b_obs に近いサンプルを抽出
tolerance = 0.1
mask = np.all(np.abs(x_b_samples - x_b_obs) < tolerance, axis=1)
x_a_conditional = x_a_samples[mask]
print(f"\n=== モンテカルロ法による検証 ===")
print(f"条件を満たすサンプル数: {mask.sum()} / {n_samples}")
print(f"条件付き平均（MC）: {x_a_conditional.mean(axis=0)}")
print(f"条件付き平均（解析）: {mu_cond}")
print(f"条件付き共分散（MC）:\n{np.cov(x_a_conditional.T)}")
print(f"条件付き共分散（解析）:\n{Sigma_cond}")
# 周辺分布の検証
print(f"\n=== 周辺分布の検証 ===")
print(f"周辺平均（MC）: {x_a_samples.mean(axis=0)}")
print(f"周辺平均（解析）: {mu_a}")
print(f"周辺共分散（MC）:\n{np.cov(x_a_samples.T)}")
print(f"周辺共分散（解析）:\n{Sigma_aa}")
まとめ
本記事では、多変量正規分布の周辺分布と条件付き分布を厳密に導出しました。
- 周辺分布: $\bm{x}_a \sim \mathcal{N}(\bm{\mu}_a, \bm{\Sigma}_{aa})$（$\bm{x}_b$ を積分消去すると、ブロック共分散がそのまま残る）
- 条件付き分布: $\bm{x}_a \mid \bm{x}_b \sim \mathcal{N}(\bm{\mu}_{a|b}, \bm{\Sigma}_{a|b})$
- 条件付き平均: $\bm{\mu}_{a|b} = \bm{\mu}_a + \bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}(\bm{x}_b – \bm{\mu}_b)$
- 条件付き共分散: $\bm{\Sigma}_{a|b} = \bm{\Sigma}_{aa} – \bm{\Sigma}_{ab}\bm{\Sigma}_{bb}^{-1}\bm{\Sigma}_{ba}$
- シューア補行列: ブロック行列の逆行列を求める上で不可欠な道具
- ガウス過程回帰: 条件付き分布の公式がそのまま予測分布の公式になる
条件付き平均が $\bm{x}_b$ の線形関数であること、条件付き共分散が $\bm{x}_b$ の値に依存しないことは、正規分布の最も重要な性質です。この性質がカルマンフィルタ、ベイズ線形回帰、ガウス過程など、多くのアルゴリズムを解析的に扱えるものにしています。
次のステップとして、以下の記事も参考にしてください。