3資産のポートフォリオを組むとき、私たちは過去の収益率から $3\times 3$ の共分散行列を推定して、リスクとリターンのトレードオフを計算します。ところが標本数が少ないと、推定された共分散行列はノイズだらけで、最適化が「過去のノイズに過剰反応した極端なポートフォリオ」を吐き出してしまう——これは資産運用の現場でよく知られた失敗です。同じことは、共分散構造を持つガウス過程のハイパーパラメータ推定や、多変量回帰の誤差分散の推定、トピックモデルの潜在変数モデリングにも起こります。
こうした「共分散行列の推定が不安定になる」問題への王道の処方箋が、共分散行列に事前分布を置くベイズ推定です。1次元の分散 $\sigma^2$ に対して逆ガンマ分布が共役事前分布になることはよく知られていますが、これを多次元の共分散行列 $\bm{\Sigma}$ に拡張したのが Wishart分布 とその逆数版である 逆Wishart分布 です。前者は精度行列（共分散行列の逆）の事前分布として、後者は共分散行列そのものの事前分布として、ベイズ多変量解析の中心に居座っています。
Wishart/逆Wishartを理解すると、世界が広がります。ポートフォリオ最適化での収縮推定（James-Stein型の改善）、ガウス過程の共分散ハイパーパラメータ事後推論、LDAのガウス版（潜在トピックを多変量ガウスで表現）、さらにはマルチタスク学習でのタスク間共分散構造の推定まで、共分散がからむあらゆる場所で顔を出します。
本記事の内容
- 「共分散行列に事前分布を置く」とはどういうことかの直感
- Wishart分布 $W_p(\bm{V}, n)$ の定義、密度関数、自由度 $n$ と尺度行列 $\bm{V}$ の意味
- 逆Wishart分布 $W_p^{-1}(\bm{\Psi}, \nu)$ との関係（精度⇔共分散の双対）
- 多変量正規モデルでの共役性と、事後ハイパーパラメータ更新式の省略なしの導出
- Normal-Inverse-Wishart 事前分布（平均 $\bm{\mu}$ と共分散 $\bm{\Sigma}$ を同時にベイズ化）
scipy.stats.wishart/invwishartによる標本生成、事後更新シミュレーション、ポートフォリオ応用までのPython実装
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
直感: 共分散行列に事前分布を置くとはどういうことか
「分散」に事前分布を置く話は、1次元のベイズ推定ですでに馴染みがあるはずです。観測 $x_i \sim N(\mu, \sigma^2)$ で $\sigma^2$ を推定したいとき、$\sigma^2$ そのものに逆ガンマ分布 $\mathrm{IG}(\alpha, \beta)$ を被せます。$\sigma^2$ は正の実数なので、台が $(0, \infty)$ の分布が必要で、その中で正規モデルと相性が良い（共役になる）のが逆ガンマでした。
では多次元、つまり $p$ 次元の観測 $\bm{x}_i \sim N(\bm{\mu}, \bm{\Sigma})$ で $\bm{\Sigma}$ を推定したいときはどうでしょうか。今度の $\bm{\Sigma}$ は単なる「正の数」ではなく、$p \times p$ の対称正定値行列です。事前分布を置きたいなら、その台は「$p \times p$ の対称正定値行列全体の空間」でなければなりません。
ここで思い出して欲しいのが、ガンマ分布の構成です。$Z_1, \dots, Z_n$ が独立に $N(0,1)$ なら、$\sum Z_i^2 \sim \chi^2_n$、これがガンマ分布の特殊形でした。$p$ 次元に拡張して $\bm{z}_1, \dots, \bm{z}_n$ が独立に $N(\bm{0}, \bm{V})$ なら、
$$ \bm{S} = \sum_{i=1}^n \bm{z}_i \bm{z}_i^\top $$
が定義する $p \times p$ 行列が従う分布——これがまさに Wishart 分布 です。$\bm{z}_i \bm{z}_i^\top$ は外積で対称半正定値、十分なサンプルを足せば正定値になります。「ガウスベクトルの外積和の分布」と覚えれば、Wishartは怖くありません。
逆Wishartは、Wishart従う行列の逆行列が従う分布です。直感的には：
- Wishart: 精度行列 $\bm{\Lambda} = \bm{\Sigma}^{-1}$ の事前分布として自然
- 逆Wishart: 共分散行列 $\bm{\Sigma}$ そのものの事前分布として自然
1次元の対応で言えば、ガンマ ↔ Wishart、逆ガンマ ↔ 逆Wishart、という双対です。分散 $\sigma^2$ に逆ガンマを置くのを多次元化したのが、$\bm{\Sigma}$ に逆Wishartを置くこと、と理解できます。
イメージが掴めたところで、Wishart分布の数学的な定義に踏み込んでいきましょう。
Wishart分布の定義と性質
構成的定義
まず構成的定義から入ります。$\bm{V}$ を $p \times p$ の正定値対称行列、$n$ を $n \geq p$ の整数とします。$\bm{z}_1, \dots, \bm{z}_n$ を独立同分布で $\bm{z}_i \sim N_p(\bm{0}, \bm{V})$ とするとき、
$$ \bm{S} = \sum_{i=1}^n \bm{z}_i \bm{z}_i^\top $$
の従う分布を Wishart 分布 と呼び、$\bm{S} \sim W_p(\bm{V}, n)$ と書きます。$n$ を自由度、$\bm{V}$ を尺度行列と呼びます。$p=1$ で $\bm{V}=v$ が単なる正数なら、$\bm{S} = \sum z_i^2$ は $v \chi^2_n$ 分布、つまりガンマ分布の特殊形に一致します。Wishartは「カイ二乗の多次元版」というわけです。
定義から、$\bm{S}$ は $n$ 個の対称半正定値行列の和なので、$n \geq p$ ならほぼ確実に正定値になります。台は $\{\bm{S} : \bm{S}\ \text{は}\ p\times p\ \text{正定値対称}\}$ で、これは $p(p+1)/2$ 次元の凸錐（正定値錐）です。
密度関数
Wishart 分布 $W_p(\bm{V}, n)$ の密度関数は、正定値対称行列 $\bm{S}$ に対して
$$ f(\bm{S}) = \frac{1}{2^{np/2} |\bm{V}|^{n/2} \Gamma_p(n/2)} |\bm{S}|^{(n-p-1)/2} \exp\!\left(-\frac{1}{2}\mathrm{tr}(\bm{V}^{-1}\bm{S})\right) $$
で与えられます。ここで $\Gamma_p(\cdot)$ は多変量ガンマ関数：
$$ \Gamma_p(a) = \pi^{p(p-1)/4} \prod_{j=1}^{p} \Gamma\!\left(a – \frac{j-1}{2}\right) $$
です。$p=1$ なら通常のガンマ関数に一致し、$\Gamma_p$ は次元 $p$ が大きいほど積で「補正」を加える役割を果たします。
密度の中心構造を読み解くと、$|\bm{S}|^{(n-p-1)/2}$ が「行列のサイズの効果」、$\exp(-\mathrm{tr}(\bm{V}^{-1}\bm{S})/2)$ が「$\bm{V}$ から離れすぎないよう減衰させる項」です。$p=1$ で $\bm{V} = v$ なら $\mathrm{tr}(\bm{V}^{-1}\bm{S}) = s/v$ となり、ガンマ分布 $\mathrm{Gamma}(n/2, 1/(2v))$ の密度と完全に一致します。
平均と分散
$\bm{S} \sim W_p(\bm{V}, n)$ の期待値は
$$ E[\bm{S}] = n\bm{V} $$
です。導出は構成的定義から直接：$E[\bm{z}_i\bm{z}_i^\top] = \bm{V}$ なので、$E[\bm{S}] = \sum_{i=1}^n E[\bm{z}_i\bm{z}_i^\top] = n\bm{V}$。シンプルです。自由度 $n$ を大きくすると、$\bm{S}$ は $n\bm{V}$ の周りに集中していき、$n \to \infty$ で確定論的に $n\bm{V}$ に近づきます。
分散構造は要素ごとに
$$ \mathrm{Var}(S_{ij}) = n(V_{ij}^2 + V_{ii}V_{jj}) $$
となります。$p=1$ なら $\mathrm{Var}(S) = n(v^2 + v^2) = 2nv^2$ で、$\chi^2_n$ 分布の分散 $2n$ をスケール $v^2$ で補正した値と一致します。
加法性と回転不変性
Wishart には便利な性質がいくつかあります。
加法性: $\bm{S}_1 \sim W_p(\bm{V}, n_1)$、$\bm{S}_2 \sim W_p(\bm{V}, n_2)$ が独立なら、$\bm{S}_1 + \bm{S}_2 \sim W_p(\bm{V}, n_1 + n_2)$。これも構成的定義から明らかで、$n_1$ 個と $n_2$ 個の $\bm{z}_i\bm{z}_i^\top$ を足すだけです。
線形変換: $\bm{A}$ を $q \times p$ の行列としたとき、$\bm{S} \sim W_p(\bm{V}, n)$ なら $\bm{A}\bm{S}\bm{A}^\top \sim W_q(\bm{A}\bm{V}\bm{A}^\top, n)$。これは、$\bm{z}_i \sim N(\bm{0}, \bm{V})$ から $\bm{A}\bm{z}_i \sim N(\bm{0}, \bm{A}\bm{V}\bm{A}^\top)$ への変換から従います。
これらの性質は、後で事後分布を計算する際に「個別の和を組み合わせる」ことの正当性を担保してくれます。
ここまでで Wishart の輪郭が見えました。次に、これと表裏一体の関係にある逆Wishart分布を導入し、共分散行列の事前分布として使う準備を整えましょう。
逆Wishart分布と関係性
定義
逆Wishart 分布 は、Wishartに従う行列の逆行列の分布として定義されます。$\bm{\Lambda} \sim W_p(\bm{V}, \nu)$ なら、$\bm{\Sigma} = \bm{\Lambda}^{-1}$ は逆Wishart分布に従い、これを $\bm{\Sigma} \sim W_p^{-1}(\bm{\Psi}, \nu)$ と書きます。ここで尺度行列のパラメータは $\bm{\Psi} = \bm{V}^{-1}$ と取り替えるのが慣習です。つまり：
$$ \bm{\Sigma} \sim W_p^{-1}(\bm{\Psi}, \nu) \iff \bm{\Sigma}^{-1} \sim W_p(\bm{\Psi}^{-1}, \nu) $$
この双対関係が記号上のすべてです。1次元での $\sigma^2 \sim \mathrm{IG}(\alpha, \beta) \iff 1/\sigma^2 \sim \mathrm{Gamma}(\alpha, \beta)$ の多次元版と思えば自然です。
密度関数
変数変換 $\bm{\Sigma} = \bm{\Lambda}^{-1}$ のヤコビアンを計算すると（行列値変数なので注意が必要、$|d\bm{\Lambda}/d\bm{\Sigma}| = |\bm{\Sigma}|^{-(p+1)}$）、逆Wishart の密度は
$$ f(\bm{\Sigma}) = \frac{|\bm{\Psi}|^{\nu/2}}{2^{\nu p/2}\Gamma_p(\nu/2)} |\bm{\Sigma}|^{-(\nu+p+1)/2} \exp\!\left(-\frac{1}{2}\mathrm{tr}(\bm{\Psi}\bm{\Sigma}^{-1})\right) $$
となります。$p=1$ なら $|\bm{\Sigma}|^{-(\nu+p+1)/2} = \sigma^{-(\nu+2)}$、$\mathrm{tr}(\bm{\Psi}\bm{\Sigma}^{-1}) = \psi/\sigma^2$ で、逆ガンマ分布 $\mathrm{IG}(\nu/2, \psi/2)$ の密度に一致します。多変量化した形が逆Wishartと理解できます。
期待値
$\bm{\Sigma} \sim W_p^{-1}(\bm{\Psi}, \nu)$ の期待値は
$$ E[\bm{\Sigma}] = \frac{\bm{\Psi}}{\nu – p – 1}, \quad \nu > p + 1 $$
です。1次元なら $E[\sigma^2] = \psi/(\nu – 2)$ で、これも逆ガンマの期待値 $\beta/(\alpha-1)$ で $\alpha = \nu/2$, $\beta = \psi/2$ とした値に一致します。$\nu – p – 1$ という奇妙な分母は、$p$ 次元化に伴う「自由度の調整」と思ってください。
期待値の存在条件 $\nu > p + 1$ は、事前分布の自由度を選ぶ際の重要な指針になります。「無情報に近い事前分布が欲しい」と思って $\nu$ を小さくしすぎると、期待値すら定義されない暴れた分布になってしまうので注意が必要です。実用上は $\nu = p + 2$（期待値が存在する最小に近い値）程度を「弱情報事前分布」として選ぶことが多いです。
モード
$\bm{\Sigma} \sim W_p^{-1}(\bm{\Psi}, \nu)$ のモード（密度の最大点）は
$$ \mathrm{Mode}[\bm{\Sigma}] = \frac{\bm{\Psi}}{\nu + p + 1} $$
です。期待値と分母が異なる（$\nu – p – 1$ vs $\nu + p + 1$）ことから、逆Wishart は右に裾を引いた非対称分布だとわかります。事前分布として使う際、「事前のベスト推測」を $\bm{\Psi}/(\nu+p+1)$ に合わせるか $\bm{\Psi}/(\nu-p-1)$ に合わせるかは設計判断です。
ここで一度、逆Wishart を「分散の事前分布」として使う動機を再確認しておきましょう。次のセクションで、多変量正規モデルに逆Wishart事前を置いたときに、事後分布が再び逆Wishartになる——という共役性を導出します。
多変量正規モデルでの共役性
設定
ベイズの王道の問題設定です。観測 $\bm{x}_1, \dots, \bm{x}_N$ を独立同分布で
$$ \bm{x}_i \sim N_p(\bm{\mu}, \bm{\Sigma}) $$
から得たとします。本セクションでは $\bm{\mu}$ は既知とし、共分散行列 $\bm{\Sigma}$ のみを推定対象とします（$\bm{\mu}$ も未知の場合は次節の Normal-Inverse-Wishart で扱います）。事前分布を
$$ \bm{\Sigma} \sim W_p^{-1}(\bm{\Psi}_0, \nu_0) $$
と置きます。$\nu_0$ は「事前知識の強さ」、$\bm{\Psi}_0$ は「事前に想定する分散構造」を表すハイパーパラメータです。
尤度の整理
まず尤度を整理します。多変量正規の密度は
$$ p(\bm{x}_i | \bm{\Sigma}) = (2\pi)^{-p/2} |\bm{\Sigma}|^{-1/2} \exp\!\left(-\frac{1}{2}(\bm{x}_i – \bm{\mu})^\top \bm{\Sigma}^{-1}(\bm{x}_i – \bm{\mu})\right) $$
です。$N$ サンプル分の同時尤度は
$$ p(\bm{X}|\bm{\Sigma}) = (2\pi)^{-Np/2} |\bm{\Sigma}|^{-N/2} \exp\!\left(-\frac{1}{2}\sum_{i=1}^N (\bm{x}_i – \bm{\mu})^\top \bm{\Sigma}^{-1}(\bm{x}_i – \bm{\mu})\right) $$
となります。指数部の二次形式の和は、トレースの性質 $\bm{a}^\top \bm{M} \bm{a} = \mathrm{tr}(\bm{M}\bm{a}\bm{a}^\top)$ を使うと
$$ \sum_{i=1}^N (\bm{x}_i – \bm{\mu})^\top \bm{\Sigma}^{-1}(\bm{x}_i – \bm{\mu}) = \mathrm{tr}\!\left(\bm{\Sigma}^{-1}\sum_{i=1}^N (\bm{x}_i – \bm{\mu})(\bm{x}_i – \bm{\mu})^\top\right) = \mathrm{tr}(\bm{\Sigma}^{-1}\bm{S}) $$
と書けます。ここで $\bm{S} = \sum_{i=1}^N (\bm{x}_i – \bm{\mu})(\bm{x}_i – \bm{\mu})^\top$ は散布行列（scatter matrix）です。トレースに押し込んだことで、尤度は $\bm{\Sigma}$ について
$$ p(\bm{X}|\bm{\Sigma}) \propto |\bm{\Sigma}|^{-N/2} \exp\!\left(-\frac{1}{2}\mathrm{tr}(\bm{\Sigma}^{-1}\bm{S})\right) $$
という形にまとまります。$\bm{\Sigma}$ に依存しない $(2\pi)^{-Np/2}$ は比例定数として落としました。
事後分布の導出
ベイズの定理より、事後分布は事前 $\times$ 尤度に比例します。
$$ p(\bm{\Sigma}|\bm{X}) \propto p(\bm{\Sigma}) \cdot p(\bm{X}|\bm{\Sigma}) $$
事前と尤度の形をそれぞれ書き下すと：
$$ p(\bm{\Sigma}) \propto |\bm{\Sigma}|^{-(\nu_0 + p + 1)/2} \exp\!\left(-\frac{1}{2}\mathrm{tr}(\bm{\Psi}_0 \bm{\Sigma}^{-1})\right) $$
$$ p(\bm{X}|\bm{\Sigma}) \propto |\bm{\Sigma}|^{-N/2} \exp\!\left(-\frac{1}{2}\mathrm{tr}(\bm{\Sigma}^{-1}\bm{S})\right) $$
これらを掛け合わせます。指数部はトレースの線形性 $\mathrm{tr}(\bm{A}) + \mathrm{tr}(\bm{B}) = \mathrm{tr}(\bm{A}+\bm{B})$ を使って合算します。
$$ p(\bm{\Sigma}|\bm{X}) \propto |\bm{\Sigma}|^{-(\nu_0 + p + 1)/2 – N/2} \exp\!\left(-\frac{1}{2}\mathrm{tr}\big((\bm{\Psi}_0 + \bm{S})\bm{\Sigma}^{-1}\big)\right) $$
指数の係数を整理すると：
$$ -\frac{\nu_0 + p + 1 + N}{2} = -\frac{(\nu_0 + N) + p + 1}{2} $$
これは逆Wishart密度 $W_p^{-1}(\bm{\Psi}, \nu)$ の冪指数 $-(\nu + p + 1)/2$ の形と完全に一致しています。比較すると：
$$ \nu_{\mathrm{post}} = \nu_0 + N, \quad \bm{\Psi}_{\mathrm{post}} = \bm{\Psi}_0 + \bm{S} $$
つまり事後分布は
$$ \boxed{\bm{\Sigma} | \bm{X} \sim W_p^{-1}(\bm{\Psi}_0 + \bm{S},\ \nu_0 + N)} $$
となります。事前分布も事後分布も同じ「逆Wishart」族で表現できる——これが共役性の正体です。
更新式の解釈
更新式 $\nu_{\mathrm{post}} = \nu_0 + N$、$\bm{\Psi}_{\mathrm{post}} = \bm{\Psi}_0 + \bm{S}$ を解釈してみます。
- 自由度の更新: 事前の「擬似観測数」$\nu_0$ に実観測数 $N$ を足す。事前情報を「$\nu_0$ サンプル分の知識」とみなしている。
- 尺度行列の更新: 事前の散布行列 $\bm{\Psi}_0$ に観測の散布行列 $\bm{S}$ を足す。事前と観測のエビデンスを単純に重ね合わせる形。
事後期待値は
$$ E[\bm{\Sigma}|\bm{X}] = \frac{\bm{\Psi}_0 + \bm{S}}{\nu_0 + N – p – 1} $$
となり、$N \to \infty$ で
$$ E[\bm{\Sigma}|\bm{X}] \to \frac{\bm{S}}{N} = \hat{\bm{\Sigma}}_{\mathrm{MLE}} $$
すなわち最尤推定量に収束します。これは大標本では事前の影響が消えるという、ベイズ推定の標準的な性質です。$N$ が小さい領域では、事前 $\bm{\Psi}_0$ が「サンプルからの推定」を引き寄せる正則化として機能します。
ここまでは $\bm{\mu}$ を既知としてきましたが、現実の問題では $\bm{\mu}$ も $\bm{\Sigma}$ も同時に未知です。両者に同時に事前分布を置く枠組みが、次に見る Normal-Inverse-Wishart 事前分布です。
Normal-Inverse-Wishart事前分布
階層的構造
$\bm{\mu}$ と $\bm{\Sigma}$ の両方を未知としたとき、共役な同時事前分布を作るには、次の階層的構造を採用します。
$$ \bm{\Sigma} \sim W_p^{-1}(\bm{\Psi}_0, \nu_0), \quad \bm{\mu} | \bm{\Sigma} \sim N_p\!\left(\bm{\mu}_0, \frac{\bm{\Sigma}}{\kappa_0}\right) $$
これを Normal-Inverse-Wishart (NIW) 分布 と呼び、$(\bm{\mu}, \bm{\Sigma}) \sim \mathrm{NIW}(\bm{\mu}_0, \kappa_0, \bm{\Psi}_0, \nu_0)$ と書きます。ハイパーパラメータは4つ：
- $\bm{\mu}_0$: $\bm{\mu}$ の事前期待値
- $\kappa_0$: $\bm{\mu}$ の事前精度の強さ（事前の「擬似観測数」）
- $\bm{\Psi}_0$: $\bm{\Sigma}$ の事前尺度行列
- $\nu_0$: $\bm{\Sigma}$ の事前自由度
$\bm{\mu}$ の事前共分散を $\bm{\Sigma}/\kappa_0$ という形にする点が肝心です。$\bm{\Sigma}$ が大きいときは $\bm{\mu}$ の不確実性も大きい——という自然なスケーリングを表現しており、これが共役性を生みます。
NIW の密度
同時密度は
$$ p(\bm{\mu}, \bm{\Sigma}) = N_p(\bm{\mu} | \bm{\mu}_0, \bm{\Sigma}/\kappa_0) \cdot W_p^{-1}(\bm{\Sigma} | \bm{\Psi}_0, \nu_0) $$
で、$\bm{\Sigma}$ について整理すると
$$ p(\bm{\mu}, \bm{\Sigma}) \propto |\bm{\Sigma}|^{-(\nu_0 + p + 2)/2} \exp\!\left(-\frac{1}{2}\mathrm{tr}(\bm{\Psi}_0\bm{\Sigma}^{-1}) – \frac{\kappa_0}{2}(\bm{\mu}-\bm{\mu}_0)^\top \bm{\Sigma}^{-1}(\bm{\mu}-\bm{\mu}_0)\right) $$
の形になります。$|\bm{\Sigma}|^{-1/2}$ がガウス分布から、$|\bm{\Sigma}|^{-(\nu_0+p+1)/2}$ が逆Wishart から出てきて、合計で $|\bm{\Sigma}|^{-(\nu_0+p+2)/2}$ となるところがポイントです。
NIW の事後更新
$N$ サンプル $\bm{x}_1, \dots, \bm{x}_N$ を観測したときの事後分布を計算しましょう。標本平均と標本散布を
$$ \bar{\bm{x}} = \frac{1}{N}\sum_{i=1}^N \bm{x}_i, \quad \bm{S} = \sum_{i=1}^N (\bm{x}_i – \bar{\bm{x}})(\bm{x}_i – \bar{\bm{x}})^\top $$
と定義します（散布行列は $\bar{\bm{x}}$ 中心であることに注意——$\bm{\mu}$ ではなく標本平均）。
導出は少し長いですが、要点は二次形式の和を「平均項」と「標本散布項」に分解することです。
$$ \sum_{i=1}^N (\bm{x}_i – \bm{\mu})(\bm{x}_i – \bm{\mu})^\top = \bm{S} + N(\bar{\bm{x}} – \bm{\mu})(\bar{\bm{x}} – \bm{\mu})^\top $$
この恒等式は、各 $\bm{x}_i – \bm{\mu}$ を $(\bm{x}_i – \bar{\bm{x}}) + (\bar{\bm{x}} – \bm{\mu})$ と分解して外積をとり、交差項が消える（$\sum (\bm{x}_i – \bar{\bm{x}}) = \bm{0}$）ことから従います。
これを尤度と事前に押し込んで二次形式を平方完成すると、事後ハイパーパラメータは：
$$ \kappa_N = \kappa_0 + N $$
$$ \bm{\mu}_N = \frac{\kappa_0 \bm{\mu}_0 + N \bar{\bm{x}}}{\kappa_0 + N} $$
$$ \nu_N = \nu_0 + N $$
$$ \bm{\Psi}_N = \bm{\Psi}_0 + \bm{S} + \frac{\kappa_0 N}{\kappa_0 + N}(\bar{\bm{x}} – \bm{\mu}_0)(\bar{\bm{x}} – \bm{\mu}_0)^\top $$
となり、事後分布は同じ NIW 族で
$$ (\bm{\mu}, \bm{\Sigma}) | \bm{X} \sim \mathrm{NIW}(\bm{\mu}_N, \kappa_N, \bm{\Psi}_N, \nu_N) $$
と書けます。
更新式の解釈
各更新式を読み解きます。
- $\bm{\mu}_N$: 事前 $\bm{\mu}_0$ と標本平均 $\bar{\bm{x}}$ の、事前擬似観測数 $\kappa_0$ と実観測数 $N$ で重み付けした加重平均。$N$ が大きくなるほど標本平均に近づく。
- $\kappa_N = \kappa_0 + N$: $\bm{\mu}$ の事前精度に観測数を足す（観測すればするほど $\bm{\mu}$ の不確実性は減る）。
- $\nu_N = \nu_0 + N$: $\bm{\Sigma}$ の自由度の更新。共役性節と同じ。
- $\bm{\Psi}_N$: 事前 $\bm{\Psi}_0$、標本散布 $\bm{S}$、そして事前平均と標本平均の乖離項 $\frac{\kappa_0 N}{\kappa_0 + N}(\bar{\bm{x}} – \bm{\mu}_0)(\bar{\bm{x}} – \bm{\mu}_0)^\top$ の合計。最後の項が重要で、「事前の平均から標本平均が外れている度合いを、共分散の不確実性に取り込む」役割を果たします。
特に最後の補正項は、Normal-Inverse-Wishart の妙味です。「事前で $\bm{\mu}$ をかなり信じていた（$\kappa_0$ が大きい）のに、観測がそこから外れていた」場合、共分散の事後尺度が大きくなり、「予想外なことが起きたから分散も大きいかも」と自動的に修正されます。
ここまで純粋な理論を見てきました。Wishart/逆Wishart/NIW は実用上どこで使われるのか、応用面を次に見ていきます——その前に、まずはPythonで実際に標本を生成して感覚を掴みましょう。
Python実装 — 標本生成と事後推論
Wishart と逆Wishart の標本生成
scipy.stats には wishart と invwishart が用意されています。まずは標本を生成して、平均と分散が理論値と合うかを確認します。
import numpy as np
from scipy.stats import wishart, invwishart
np.random.seed(0)
# パラメータ設定: 3x3 の例
p = 3
nu = 10 # 自由度
V = np.array([[2.0, 0.5, 0.3],
[0.5, 1.0, 0.2],
[0.3, 0.2, 1.5]]) # 尺度行列（正定値対称）
# Wishart 標本生成（1000サンプル）
W_samples = wishart.rvs(df=nu, scale=V, size=1000, random_state=0)
print("Wishart サンプル形状:", W_samples.shape) # (1000, 3, 3)
# 経験平均と理論平均の比較
W_mean_emp = W_samples.mean(axis=0)
W_mean_th = nu * V
print("経験平均:\n", np.round(W_mean_emp, 3))
print("理論平均 nu*V:\n", np.round(W_mean_th, 3))
実行すると、経験平均と理論平均 $n\bm{V}$ がほぼ一致することが確認できます。たとえば $(1,1)$ 要素は理論値 $10 \times 2.0 = 20.0$ に対し経験値も $20$ 前後の値を返します。これは Wishart 分布の構成的定義「ガウスベクトルの外積和」が正しくシミュレーションされていることの裏付けです。サンプル数を増やせば一致はさらに精緻になります。
続いて逆Wishart のサンプリングを行い、期待値 $\bm{\Psi}/(\nu-p-1)$ を確認します。
import numpy as np
from scipy.stats import invwishart
np.random.seed(1)
p = 3
nu = 10 # 自由度（nu > p+1 = 4 で期待値が存在）
Psi = np.array([[2.0, 0.5, 0.3],
[0.5, 1.0, 0.2],
[0.3, 0.2, 1.5]])
# 逆Wishart サンプル
IW_samples = invwishart.rvs(df=nu, scale=Psi, size=2000, random_state=1)
IW_mean_emp = IW_samples.mean(axis=0)
IW_mean_th = Psi / (nu - p - 1)
print("逆Wishart 経験平均:\n", np.round(IW_mean_emp, 4))
print("理論平均 Psi/(nu-p-1):\n", np.round(IW_mean_th, 4))
経験平均と理論値 $\bm{\Psi}/(\nu-p-1) = \bm{\Psi}/6$ がほぼ一致します。$(1,1)$ 要素なら理論 $2.0/6 \approx 0.333$ に対し経験値も同程度です。逆Wishart は右に裾を引いた分布なので、有限サンプル数では「外れ値」によって経験平均が理論より大きく振れることがあります。サンプルを増やすか、ロバストな中央値を見るのも一手です。
固有値分布の可視化
Wishart の振る舞いを直感的につかむには、生成された行列の固有値の分布を見るのが効果的です。固有値は行列の「主軸方向のスケール」を表します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import wishart, invwishart
np.random.seed(2)
p = 3
nu_list = [5, 10, 30, 100] # 自由度を変えながら比較
V = np.eye(p) # 単位行列を尺度に
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
for nu in nu_list:
# Wishart
W_samples = wishart.rvs(df=nu, scale=V, size=1000, random_state=nu)
eigs_W = np.linalg.eigvalsh(W_samples).flatten()
axes[0].hist(eigs_W, bins=60, alpha=0.5, density=True, label=f'nu={nu}')
# 逆Wishart
IW_samples = invwishart.rvs(df=nu+p+2, scale=V, size=1000, random_state=nu)
eigs_IW = np.linalg.eigvalsh(IW_samples).flatten()
axes[1].hist(eigs_IW, bins=60, alpha=0.5, density=True,
range=(0, 2), label=f'nu={nu+p+2}')
axes[0].set_xlabel('Eigenvalue')
axes[0].set_ylabel('Density')
axes[0].set_title('Wishart eigenvalue distribution (V=I)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[1].set_xlabel('Eigenvalue')
axes[1].set_ylabel('Density')
axes[1].set_title('Inverse-Wishart eigenvalue distribution (Psi=I)')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('wishart_eigenvalues.png', dpi=150, bbox_inches='tight')
plt.show()
このヒストグラムから、自由度 $\nu$ の役割が一目でわかります。Wishart の場合、$\nu$ が小さいほど固有値分布は広がり（極端な大きさの主軸ができやすい）、$\nu$ が大きいと $\nu$ 付近に鋭く集中します（行列がほぼ $\nu \bm{I}$ に近づく）。逆Wishart は逆に、$\nu$ が小さいと固有値の裾が右に長く伸び（共分散の不確実性が大）、$\nu$ が大きいと固有値が小さい値に集中します。実用上、「事前情報をどれくらい強く信じるか」を $\nu$ の選択で表現できる、ということが視覚的に確認できます。
共分散行列推定: MLE vs ベイズ
理論的な利点を確認するため、共分散行列の最尤推定（MLE）とベイズ推定（事後平均）を比較します。サンプル数が少ない領域でベイズが有利になるはずです。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal, invwishart
np.random.seed(3)
# 真の共分散行列（4x4）
p = 4
Sigma_true = np.array([[2.0, 0.8, 0.3, 0.1],
[0.8, 1.5, 0.4, 0.2],
[0.3, 0.4, 1.0, 0.5],
[0.1, 0.2, 0.5, 1.2]])
mu_true = np.zeros(p)
# 事前分布のハイパーパラメータ（弱情報事前: 事前平均 = 単位行列）
nu_0 = p + 2 # 期待値が存在する最小に近い
Psi_0 = (nu_0 - p - 1) * np.eye(p) # 事前期待値が I になる設定
def frobenius_err(A, B):
return np.linalg.norm(A - B, 'fro')
N_list = [10, 20, 50, 100, 200, 500, 1000]
err_mle, err_bayes = [], []
n_trials = 200 # 各 N で複数回試行して平均
for N in N_list:
e_mle, e_bayes = [], []
for trial in range(n_trials):
X = multivariate_normal.rvs(mean=mu_true, cov=Sigma_true,
size=N, random_state=trial)
# MLE: 標本共分散
Sigma_mle = np.cov(X.T, bias=True)
# ベイズ事後平均: (Psi_0 + S) / (nu_0 + N - p - 1)
S = (X - X.mean(axis=0)).T @ (X - X.mean(axis=0))
Sigma_bayes = (Psi_0 + S) / (nu_0 + N - p - 1)
e_mle.append(frobenius_err(Sigma_mle, Sigma_true))
e_bayes.append(frobenius_err(Sigma_bayes, Sigma_true))
err_mle.append(np.mean(e_mle))
err_bayes.append(np.mean(e_bayes))
plt.figure(figsize=(9, 5))
plt.semilogx(N_list, err_mle, 'o-', label='MLE')
plt.semilogx(N_list, err_bayes, 's-', label='Bayesian (posterior mean)')
plt.xlabel('Sample size N')
plt.ylabel('Frobenius error to true Sigma')
plt.title('Covariance estimation: MLE vs Bayesian')
plt.legend()
plt.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.savefig('mle_vs_bayes.png', dpi=150, bbox_inches='tight')
plt.show()
for N, em, eb in zip(N_list, err_mle, err_bayes):
print(f"N={N:4d}: MLE err={em:.3f}, Bayes err={eb:.3f}, improvement={(em-eb)/em*100:+.1f}%")
このグラフから、サンプル数 $N$ が小さい領域（$N=10$〜$50$）でベイズ推定が顕著にMLEより低い誤差を達成することが読み取れます。$N=10$ では MLE の誤差が大きい一方、ベイズは事前分布の「単位行列」に引き寄せられて誤差が抑えられます。サンプルが増えるほど両者は接近し、$N=1000$ ではほぼ同じになります。これは事前情報の影響が大標本で消えるという理論的予想と完全に整合しており、共分散の事前分布が小標本での正則化として機能していることを定量的に確認できます。
Normal-Inverse-Wishart の事後更新シミュレーション
NIW で $\bm{\mu}$ と $\bm{\Sigma}$ を同時推定する例です。事後分布から繰り返しサンプリングし、真値を「信用区間」が覆っているかを確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal, invwishart
np.random.seed(4)
# 真値
p = 2
mu_true = np.array([1.0, -0.5])
Sigma_true = np.array([[1.5, 0.6],
[0.6, 1.0]])
# 観測データ
N = 30
X = multivariate_normal.rvs(mean=mu_true, cov=Sigma_true, size=N, random_state=4)
# NIW 事前
mu_0 = np.zeros(p)
kappa_0 = 1.0
nu_0 = p + 2
Psi_0 = np.eye(p)
# 事後ハイパーパラメータ
x_bar = X.mean(axis=0)
S = (X - x_bar).T @ (X - x_bar)
kappa_N = kappa_0 + N
mu_N = (kappa_0 * mu_0 + N * x_bar) / kappa_N
nu_N = nu_0 + N
Psi_N = Psi_0 + S + (kappa_0 * N / kappa_N) * np.outer(x_bar - mu_0, x_bar - mu_0)
print(f"事後 mu_N = {mu_N}, 真値 = {mu_true}")
print(f"事後 Sigma の期待値:\n{Psi_N/(nu_N-p-1)}")
print(f"真の Sigma:\n{Sigma_true}")
# 事後分布からのサンプリング: Sigma をサンプル → mu をサンプル
n_post = 2000
mu_samples = np.zeros((n_post, p))
Sigma_samples = np.zeros((n_post, p, p))
for i in range(n_post):
Sigma_i = invwishart.rvs(df=nu_N, scale=Psi_N, random_state=i)
mu_i = multivariate_normal.rvs(mean=mu_N, cov=Sigma_i/kappa_N, random_state=i+10000)
Sigma_samples[i] = Sigma_i
mu_samples[i] = mu_i
# 可視化
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# mu の散布図と真値
axes[0].scatter(mu_samples[:, 0], mu_samples[:, 1], alpha=0.2, s=10, label='posterior samples')
axes[0].plot(mu_true[0], mu_true[1], 'r*', markersize=20, label='true mu')
axes[0].plot(x_bar[0], x_bar[1], 'g^', markersize=12, label='sample mean')
axes[0].set_xlabel('mu_1')
axes[0].set_ylabel('mu_2')
axes[0].set_title('Posterior of mu')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# Sigma の対角成分のヒストグラム
axes[1].hist(Sigma_samples[:, 0, 0], bins=40, alpha=0.5, label='Sigma_11')
axes[1].hist(Sigma_samples[:, 1, 1], bins=40, alpha=0.5, label='Sigma_22')
axes[1].axvline(Sigma_true[0, 0], color='red', ls='--', label='true Sigma_11')
axes[1].axvline(Sigma_true[1, 1], color='blue', ls='--', label='true Sigma_22')
axes[1].set_xlabel('Variance')
axes[1].set_ylabel('Frequency')
axes[1].set_title('Posterior of Sigma (diagonal)')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('niw_posterior.png', dpi=150, bbox_inches='tight')
plt.show()
このプロットから、NIW 事後分布の振る舞いが明瞭に見えます。$\bm{\mu}$ の事後サンプル（散布図）は真値（赤い星）を中心に分布し、標本平均（緑の三角）と真値の間に挟まれる位置に事後平均が来ます。これは事前 $\bm{\mu}_0 = \bm{0}$ への引き寄せが効いている証拠です。共分散の対角成分も真値の周りに分布し、$N=30$ という中程度のサンプル数では事後分布の幅は無視できない大きさです——これが不確実性の定量化であり、MLE の点推定では得られない情報です。
理論と実装で基礎を固めたので、最後に実用的な応用例として、金融ポートフォリオ最適化を見ていきましょう。
応用: ポートフォリオ最適化と他の応用先
Markowitz 最適化での共分散推定の問題
金融ポートフォリオ最適化（Markowitzの平均-分散最適化）では、$n$ 個の資産の収益率 $\bm{r} \in \mathbb{R}^n$ に対して、共分散行列 $\bm{\Sigma}$ から「リスク最小ポートフォリオ重み」が次で与えられます。
$$ \bm{w}^* = \frac{\bm{\Sigma}^{-1} \bm{1}}{\bm{1}^\top \bm{\Sigma}^{-1} \bm{1}} $$
問題は、$\bm{\Sigma}$ を標本共分散（MLE）で推定すると、特に資産数 $n$ が観測期間に対して大きい場合（典型: 50資産 × 12ヶ月）、$\bm{\Sigma}$ がほぼ特異になり、$\bm{\Sigma}^{-1}$ が暴走することです。結果として「ある資産に+200%、別の資産に-180%」みたいな現実離れした重みが出てしまいます。
ここで逆Wishart 事前分布によるベイズ推定が威力を発揮します。事前 $\bm{\Sigma} \sim W_p^{-1}(\bm{\Psi}_0, \nu_0)$ を「対角行列に近い形」で置けば、事後平均は標本共分散と対角行列の重み付き平均になり、特異性が解消されます。これは実は古典的な Ledoit-Wolf 収縮推定 とほぼ同じ形になり、Wishart事前は「収縮推定のベイズ的解釈」を与えていると言えます。
実際にシミュレーションで確認しましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
np.random.seed(5)
# 10資産の真の共分散行列（適度な相関構造）
n_assets = 10
A = np.random.randn(n_assets, n_assets)
Sigma_true = A @ A.T / n_assets + np.eye(n_assets) * 0.1
mu_returns = np.random.randn(n_assets) * 0.01
def min_var_weights(Sigma):
"""最小分散ポートフォリオの重み"""
inv_Sigma = np.linalg.inv(Sigma)
ones = np.ones(Sigma.shape[0])
return inv_Sigma @ ones / (ones @ inv_Sigma @ ones)
def portfolio_variance(w, Sigma):
return w @ Sigma @ w
# 観測サンプル数を変えてMLEとベイズの結果を比較
N_list = [15, 20, 30, 50, 100, 300]
mle_realized_var, bayes_realized_var = [], []
# 事前分布: 対角を仮定（事前期待値が対角行列になる設定）
nu_0 = n_assets + 5
Psi_0 = (nu_0 - n_assets - 1) * np.eye(n_assets) # 期待値が I になる
n_trials = 100
for N in N_list:
mle_vars, bayes_vars = [], []
for trial in range(n_trials):
X = multivariate_normal.rvs(mean=mu_returns, cov=Sigma_true,
size=N, random_state=trial * 7)
Sigma_mle = np.cov(X.T, bias=True)
S_scatter = (X - X.mean(axis=0)).T @ (X - X.mean(axis=0))
Sigma_bayes = (Psi_0 + S_scatter) / (nu_0 + N - n_assets - 1)
try:
w_mle = min_var_weights(Sigma_mle)
w_bayes = min_var_weights(Sigma_bayes)
# 真の Sigma に対する実現分散を評価
mle_vars.append(portfolio_variance(w_mle, Sigma_true))
bayes_vars.append(portfolio_variance(w_bayes, Sigma_true))
except np.linalg.LinAlgError:
pass
mle_realized_var.append(np.mean(mle_vars))
bayes_realized_var.append(np.mean(bayes_vars))
# 真の最適ポートフォリオの分散（下限）
true_optimal_var = portfolio_variance(min_var_weights(Sigma_true), Sigma_true)
plt.figure(figsize=(9, 5))
plt.semilogx(N_list, mle_realized_var, 'o-', label='MLE-based weights')
plt.semilogx(N_list, bayes_realized_var, 's-', label='Bayesian weights')
plt.axhline(true_optimal_var, color='red', ls='--', label='True optimal')
plt.xlabel('Sample size N (observation periods)')
plt.ylabel('Realized portfolio variance')
plt.title(f'Min-var portfolio: MLE vs Bayes ({n_assets} assets)')
plt.legend()
plt.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.savefig('portfolio_optimization.png', dpi=150, bbox_inches='tight')
plt.show()
for N, em, eb in zip(N_list, mle_realized_var, bayes_realized_var):
print(f"N={N:3d}: MLE realized var={em:.4f}, Bayes={eb:.4f}, "
f"true optimal={true_optimal_var:.4f}")
このグラフから、$N$ が資産数 $n=10$ に近い領域（$N=15$〜$30$）で MLE ベースのポートフォリオが真の最適より大きく劣化する（実現分散が膨らむ）一方、ベイズ事前で正則化したポートフォリオは真の最適に近い性能を保つことが見て取れます。サンプル数が増えれば両者は収束しますが、現実の資産運用では「使えるヒストリーが短い」ことが多いため、共分散行列のベイズ推定は実用的価値が大きいわけです。
その他の応用先
Wishart/逆Wishart は金融以外にも広く使われています。
ガウス過程回帰のハイパーパラメータ: ガウス過程ではカーネル関数のスケールやノイズ分散を推定しますが、複数の長さスケールが組み合わさる場合、それらをまとめて多変量ガウス的に扱い、逆Wishart で事前分布を置く方法が研究されています。
潜在ディリクレ配分 (LDA) のガウス変種: 標準的なLDAではトピックが離散分布で表現されますが、トピックを連続的な「埋め込み空間の点」とし、多変量ガウスでモデル化する変種（Gaussian-LDA）では、各トピックの共分散行列に逆Wishart事前を置きます。
多変量ベイズ回帰: $\bm{Y} = \bm{X}\bm{B} + \bm{E}$ という多変量回帰モデル（$\bm{E}$ の行が多変量正規）では、誤差共分散 $\bm{\Sigma}$ に逆Wishart 事前を置くと、係数行列 $\bm{B}$ と共分散の同時事後分布が解析的に書けます。
マルチタスク学習: $K$ 個のタスクのパラメータ間に共通の事前分布を置く際、タスク間共分散行列に逆Wishart 事前を置いて「タスクの類似性」をデータから推定するアプローチがあります。
カルマンフィルタの適応化: プロセス雑音や観測雑音の共分散行列を未知としたとき、それらに逆Wishart事前を置いてオンラインベイズ更新する適応カルマンフィルタが提案されています。
いずれも「共分散行列が未知のパラメータとして現れ、その不確実性を真面目に扱いたい」という共通の動機から、Wishart 族が出てくる構図です。
まとめ
本記事では、多変量正規分布の共分散行列に対するベイズ事前分布として、Wishart 分布と逆Wishart 分布を導入し、共役性に基づく事後更新式の導出と実装、そして応用までを解説しました。
- 構成的定義: $\bm{z}_i \sim N_p(\bm{0}, \bm{V})$ の外積和 $\sum \bm{z}_i \bm{z}_i^\top$ が Wishart $W_p(\bm{V}, n)$。逆Wishart はその逆行列の分布。
- 基本性質: $E[W_p(\bm{V}, n)] = n\bm{V}$、$E[W_p^{-1}(\bm{\Psi}, \nu)] = \bm{\Psi}/(\nu-p-1)$。1次元化するとそれぞれカイ二乗・逆ガンマに帰着。
- 共役性: $\bm{\mu}$ 既知の多変量正規モデル $\bm{x}_i \sim N(\bm{\mu}, \bm{\Sigma})$ に $\bm{\Sigma} \sim W_p^{-1}(\bm{\Psi}_0, \nu_0)$ を置くと、事後分布は $W_p^{-1}(\bm{\Psi}_0 + \bm{S}, \nu_0 + N)$。
- Normal-Inverse-Wishart: $\bm{\mu}$ と $\bm{\Sigma}$ を同時にベイズ化する共役事前。事後更新は標本平均と事前平均の加重平均、および「乖離項」を含む散布行列の更新。
- 応用: 小標本下の共分散推定の正則化、ポートフォリオ最適化の収縮推定、ガウス過程ハイパーパラメータ、Gaussian-LDA、多変量ベイズ回帰、マルチタスク学習など。
数値実験では、サンプル数が少ない領域で逆Wishart事前による正則化がMLEを明確に上回ること、ポートフォリオ最適化で「資産数に対してサンプル数が小さい」現実的な状況で実現分散を抑えられることを確認しました。共分散行列のベイズ推定は、教科書的な理論であると同時に、現場の実問題で確かな成果を出す実用技術です。
次のステップとして、以下の記事も参考にしてください。