2008年のリーマンショックでも、2020年のコロナショックでも、株価の暴落は「100年に1度」と評されました。けれども実際には、こうした暴落は数十年に何度も繰り返し起きています。原因の一つは、私たちがリターンの分布を多変量正規分布でモデル化しがちなことにあります。正規分布の裾は指数関数の二乗で減衰するため、平均から離れた値は事実上「ほぼ起こらない」扱いです。ところが現実の市場リターンは、もっと太い裾——いわゆるファットテール——を持ち、極端事象が正規分布の予言よりはるかに頻繁に起こるのです。
この「裾の重さ」を一つの自由度パラメータで連続的に制御できる多変量モデルが、本記事で扱う多変量t分布（multivariate t-distribution）です。多変量t分布は、平均ベクトル $\bm{\mu}$ と共分散構造 $\Sigma$ に加え、自由度 $\nu$ という第3のパラメータを持ち、$\nu \to \infty$ で多変量正規分布に一致し、$\nu$ が小さいほど裾が重くなります。これにより、外れ値の影響を吸収するロバストな統計モデルが自然に構成できます。
応用先は驚くほど幅広いです。一つは金融リスク管理——多変量t分布でポートフォリオ収益をモデル化することで、正規分布では過小評価されるテールリスク（VaR、CVaR）を現実的な水準で見積もれます。もう一つはロバストクラスタリング——通常のガウス混合モデル（GMM）は外れ値1つで平均と分散が大きく引きずられますが、各成分をt分布に置き換えた Mixture of t（MoT） は外れ値を「自然に小さな重みで扱う」効果を持ちます。さらに、外れ値検出、ロバスト線形回帰、ベイズモデリングの事前分布など、多変量t分布は現代統計学のいたるところで顔を出します。
本記事の内容
- 多変量正規分布が捉えきれない外れ値・極端事象の直感的理解
- 多変量t分布の確率密度関数の定義と、ガンマ関数を用いた導出
- ガウス・スケールミクスチャとしての表現と、自由度 $\nu$ の幾何学的意味
- 平均・共分散の存在条件、Mahalanobis距離とF分布の関係
- EMアルゴリズムによる最尤推定（$\bm{\mu}, \Sigma, \nu$ の同時推定）
- Pythonによるロバスト混合モデル（MoT）のスクラッチ実装とGMMとの比較
- 金融リスク（VaR/CVaR）への応用と、正規分布との差異の数値検証
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- t分布の定義と性質をわかりやすく解説
- 多変量正規分布とは？定義から性質まで徹底解説
- 多変量ガウス分布の周辺化と条件付け
- 条件付き多変量ガウス分布の導出
- ガウス混合モデルとEMアルゴリズム
- ガウス混合モデルによるクラスタリング
直感 — 多変量ガウスでは捕まらない外れ値
多変量正規分布 $N(\bm{\mu}, \Sigma)$ の密度は、Mahalanobis距離 $d^2 = (\bm{x}-\bm{\mu})^\top \Sigma^{-1}(\bm{x}-\bm{\mu})$ に対して $\exp(-d^2/2)$ で減衰します。指数関数のさらに二次で押さえつけるわけですから、$d$ が大きい点はほぼゼロ確率です。たとえば1次元標準正規分布で $|x| > 4$ となる確率はおよそ $6 \times 10^{-5}$、$|x| > 6$ になると $2 \times 10^{-9}$ で「天文学的にまれ」と言ってよい水準です。
ところが、S&P500の日次対数リターンを実際に集めてみると、$|x| > 4\sigma$ の日は数年に1度は出現し、$|x| > 6\sigma$ の日も歴史を遡れば何度かあります。これを正規分布で説明しようとすると、「この10年に何度も100万年に1度の事象が起きた」と言わざるを得ません。直感的には、現実の市場が起こせる「動き」は、正規分布が想定するよりずっと派手なのです。
この観察は1次元だけでなく、ベクトル値の現象でも成り立ちます。複数の銘柄や複数のセンサ計測を同時に扱う多変量データでは、平常時には密度の中心に固まるが、ショック時には複数次元が同時に大きく動く——という構造がよく見られます。多変量正規分布の楕円等密度線は、こうした「同時の極端事象」をほとんど確率ゼロと扱ってしまいます。
t分布はこの問題に対する古典的な処方箋です。1次元のStudentのt分布が標本平均の検定で正規分布の代わりに使われるのと同じ理由で、多変量t分布は「サンプル数が少ない・外れ値がある・分散が不確実」といった状況での自然な拡張になります。直感的なイメージとしては、
- 正規分布の楕円密度を、外側へ滑らかに「裾を引き伸ばした」もの
- 自由度 $\nu$ を回転ノブのように回すと、$\nu$ 小：裾が太い／$\nu$ 大：裾が細い、と連続的に変化
という連続族です。$\nu \to \infty$ で正規分布に戻ることを後ほど見ます。
このイメージを数式に落とし込みましょう。次セクションで、多変量t分布の密度関数を定義し、なぜその形になるのかをガンマ関数の構造から見ていきます。
多変量t分布の定義と導出
密度関数の定義
$p$ 次元の確率ベクトル $\bm{x} \in \mathbb{R}^p$ が、平均ベクトル $\bm{\mu} \in \mathbb{R}^p$、正定値スケール行列 $\Sigma \in \mathbb{R}^{p \times p}$、自由度 $\nu > 0$ の多変量t分布 $\mathcal{T}_p(\bm{\mu}, \Sigma, \nu)$ に従うとは、密度関数が次で与えられることを言います。
$$ p(\bm{x} \mid \bm{\mu}, \Sigma, \nu) = \frac{\Gamma\!\left(\frac{\nu+p}{2}\right)}{\Gamma\!\left(\frac{\nu}{2}\right) (\nu \pi)^{p/2} |\Sigma|^{1/2}} \left[1 + \frac{1}{\nu}(\bm{x}-\bm{\mu})^\top \Sigma^{-1}(\bm{x}-\bm{\mu})\right]^{-(\nu+p)/2} $$
この式は一見複雑ですが、構造を分解すれば自然に読み解けます。
- $|\Sigma|^{1/2}$ と $(\bm{x}-\bm{\mu})^\top\Sigma^{-1}(\bm{x}-\bm{\mu})$（Mahalanobis距離の二乗 $d^2$）の組は、多変量正規分布と全く同じです。スケール行列 $\Sigma$ がベクトル方向ごとの広がりを担います。
- 鍵は最後の $[1 + d^2/\nu]^{-(\nu+p)/2}$ です。多変量正規分布の $\exp(-d^2/2)$ を、多項式（べき乗）型の減衰 に置き換えた形になっています。
- 先頭のガンマ関数の比は規格化定数で、密度の積分が1になるよう調整しています。
ベキ乗減衰こそが「重い裾」の正体です。$d^2$ が大きくなると、正規分布なら指数関数で猛烈に潰されますが、t分布では多項式でゆるやかにしか減衰しません。だから極端値の確率がずっと残ります。
何を「t分布」と呼んでいるのか
1次元のStudentのt分布は、独立な標準正規 $Z$ とカイ二乗 $V \sim \chi^2_\nu$ から構成される $T = Z / \sqrt{V/\nu}$ の分布です。多変量t分布は、これを直接ベクトル化したものとして次のように構成できます。
$$ \bm{T} = \bm{\mu} + \frac{\bm{Z}}{\sqrt{V/\nu}}, \quad \bm{Z} \sim N_p(\bm{0}, \Sigma), \quad V \sim \chi^2_\nu, \quad \bm{Z} \perp V $$
つまり、多変量正規ベクトル $\bm{Z}$ を、独立なカイ二乗の平方根 $\sqrt{V/\nu}$ で割って平均をずらしたものが多変量t分布です。$\sqrt{V/\nu}$ は確率変数なので、ベクトルのスケールがランダムに揺らぐ正規分布——という解釈ができます。
密度の導出（変数変換による）
上の構成から密度を計算してみましょう。同時密度を求めて $V$ を周辺化する戦略です。簡単のため $\bm{\mu} = \bm{0}$ とします。$\bm{Z} \sim N_p(\bm{0}, \Sigma)$ と $V \sim \chi^2_\nu$ の同時密度は独立性から積で書けます。
$$ p(\bm{z}, v) = \frac{1}{(2\pi)^{p/2}|\Sigma|^{1/2}} \exp\!\left(-\tfrac{1}{2}\bm{z}^\top \Sigma^{-1} \bm{z}\right) \cdot \frac{v^{\nu/2 – 1} e^{-v/2}}{2^{\nu/2} \Gamma(\nu/2)} $$
変数変換 $\bm{x} = \bm{z}/\sqrt{v/\nu}$、$v = v$ を行います。$\bm{z} = \sqrt{v/\nu}\, \bm{x}$ なので、Jacobianは $|\partial \bm{z}/\partial \bm{x}| = (v/\nu)^{p/2}$ です。代入すると、
$$ p(\bm{x}, v) = \frac{(v/\nu)^{p/2}}{(2\pi)^{p/2}|\Sigma|^{1/2}} \exp\!\left(-\tfrac{v}{2\nu}\bm{x}^\top \Sigma^{-1} \bm{x}\right) \cdot \frac{v^{\nu/2-1} e^{-v/2}}{2^{\nu/2} \Gamma(\nu/2)} $$
ここで Mahalanobis距離の二乗 $d^2 = \bm{x}^\top \Sigma^{-1} \bm{x}$ と置き、$v$ をまとめて整理します。
$$ p(\bm{x}, v) = \frac{1}{(2\pi)^{p/2} |\Sigma|^{1/2} 2^{\nu/2} \Gamma(\nu/2) \nu^{p/2}} \cdot v^{(\nu+p)/2 – 1} \exp\!\left(-\tfrac{v}{2}\left[1 + \tfrac{d^2}{\nu}\right]\right) $$
$v$ について積分すれば $\bm{x}$ の周辺密度が得られます。ガンマ関数の積分公式 $\int_0^\infty v^{\alpha-1} e^{-\beta v} dv = \Gamma(\alpha)/\beta^\alpha$ を、$\alpha = (\nu+p)/2$、$\beta = \tfrac{1}{2}(1 + d^2/\nu)$ として適用します。
$$ \int_0^\infty v^{(\nu+p)/2 – 1} e^{-v[1 + d^2/\nu]/2}\, dv = \frac{\Gamma\!\big((\nu+p)/2\big)}{\big(\tfrac{1}{2}[1 + d^2/\nu]\big)^{(\nu+p)/2}} = 2^{(\nu+p)/2}\,\Gamma\!\big((\nu+p)/2\big) \left[1 + \tfrac{d^2}{\nu}\right]^{-(\nu+p)/2} $$
この結果を代入し、$2^{(\nu+p)/2}/(2^{\nu/2} \cdot (2\pi)^{p/2}) = 1/\pi^{p/2}$ で $2$ のべきを整理すると、冒頭の密度式
$$ p(\bm{x}) = \frac{\Gamma((\nu+p)/2)}{\Gamma(\nu/2) (\nu \pi)^{p/2} |\Sigma|^{1/2}} \left[1 + \frac{d^2}{\nu}\right]^{-(\nu+p)/2} $$
がきれいに得られます。
平均と共分散の存在条件
t分布は裾が重いので、モーメントが必ずしも有限ではありません。詳細な計算は省きますが、結論だけまとめておくと:
- 平均 $E[\bm{X}] = \bm{\mu}$ が存在するのは $\nu > 1$ のとき。$\nu \leq 1$ では平均自体が定義されない（コーシー分布が $\nu = 1$ の特殊例）。
- 共分散行列 は $\nu > 2$ のとき存在し、$\mathrm{Cov}[\bm{X}] = \dfrac{\nu}{\nu – 2}\, \Sigma$ となる。$\Sigma$ そのものは共分散ではなくスケール行列であって、$\nu$ が小さいほど真の共分散が膨らむ点に注意。
つまり、$\Sigma$ を「正規分布の共分散」と同じ感覚で読み替えてはいけない。これがt分布を扱うときの最初の落とし穴です。データから $\Sigma$ を推定する場合、サンプル共分散 $\hat{S}$ とは $\hat{S} \approx \tfrac{\nu}{\nu-2} \hat{\Sigma}$ の関係になります。
ここまでで密度の形と基本性質が見えました。しかし、上の導出をもう一度眺めると、「正規分布 $\bm{Z}$ をランダムスケール $\sqrt{V/\nu}$ で割る」という構成法がより深い構造を持っていることに気付きます。実はt分布は、無限個の正規分布の混合として表現できる——これがガウス・スケールミクスチャ表現で、後のEMアルゴリズム導出の鍵になります。
ガウス・スケールミクスチャ表現
表現の主張
多変量t分布 $\mathcal{T}_p(\bm{\mu}, \Sigma, \nu)$ は、次の階層構造として書けます。
$$ \begin{aligned} \bm{X} \mid u &\sim N_p\!\left(\bm{\mu},\; \Sigma/u\right) \\ u &\sim \mathrm{Gamma}(\nu/2,\; \nu/2) \end{aligned} $$
ここで $\mathrm{Gamma}(\alpha, \beta)$ は形状 $\alpha$・rate $\beta$ のガンマ分布（密度 $\propto u^{\alpha-1} e^{-\beta u}$）です。$u$ を潜在スケール変数として、それで条件付けた正規分布の共分散が $\Sigma/u$ になる——という構造です。
直感的には、「観測ごとに正規分布の幅がランダムに変わる」という生成過程です。たまたま $u$ が小さい観測ではスケールが膨らんで外れ値らしい点が出やすく、$u$ が大きい観測ではスケールが縮まって中央付近に集まります。この「スケールが揺らぐ正規分布の重ね合わせ」が、多変量t分布なのです。
導出
ガンマ分布 $u \sim \mathrm{Gamma}(\nu/2, \nu/2)$ の密度は
$$ p(u) = \frac{(\nu/2)^{\nu/2}}{\Gamma(\nu/2)} u^{\nu/2 – 1} e^{-\nu u / 2} $$
条件付き密度 $p(\bm{x} \mid u) = N_p(\bm{x}; \bm{\mu}, \Sigma/u)$ は
$$ p(\bm{x} \mid u) = \frac{u^{p/2}}{(2\pi)^{p/2} |\Sigma|^{1/2}} \exp\!\left(-\tfrac{u}{2} d^2\right) $$
ここで $d^2 = (\bm{x}-\bm{\mu})^\top \Sigma^{-1}(\bm{x}-\bm{\mu})$ です。$\Sigma/u$ の行列式が $u^{-p} |\Sigma|$ になり、逆行列が $u \Sigma^{-1}$ になることを使いました。
周辺密度は $u$ について積分して
$$ p(\bm{x}) = \int_0^\infty p(\bm{x} \mid u) p(u)\, du = \frac{(\nu/2)^{\nu/2}}{(2\pi)^{p/2} |\Sigma|^{1/2} \Gamma(\nu/2)} \int_0^\infty u^{(\nu+p)/2 – 1} \exp\!\left(-\tfrac{u}{2}[\nu + d^2]\right) du $$
ガンマ積分 $\int_0^\infty u^{\alpha-1} e^{-\beta u} du = \Gamma(\alpha)/\beta^\alpha$ を $\alpha = (\nu+p)/2$、$\beta = (\nu+d^2)/2$ で適用すると
$$ \int_0^\infty u^{(\nu+p)/2-1} e^{-u(\nu+d^2)/2} du = \Gamma\!\big((\nu+p)/2\big) \left(\tfrac{\nu+d^2}{2}\right)^{-(\nu+p)/2} $$
代入して整理すると、多変量t分布の密度がそのまま現れます。前節のカイ二乗による導出と本質的に同じ計算ですが、こちらの方が「混合分布」としての解釈が前面に出ます。
潜在変数 $u$ の事後分布
ガウス・スケールミクスチャ表現のうれしいところは、潜在変数 $u$ の事後分布が解析的に求まることです。観測 $\bm{x}$ を与えたときの $u$ の事後密度は、ベイズの定理より $p(u \mid \bm{x}) \propto p(\bm{x} \mid u) p(u)$ で、上の積分被積分関数そのものです。整理すると
$$ u \mid \bm{x} \sim \mathrm{Gamma}\!\left(\frac{\nu + p}{2},\; \frac{\nu + d^2}{2}\right) $$
ガンマ分布の平均が $\alpha/\beta$ であることを使うと、
$$ E[u \mid \bm{x}] = \frac{\nu + p}{\nu + d^2} $$
この期待値の式は、後のEMアルゴリズムで観測ごとの重みとして使います。$d^2$ が大きい（中心から遠い）観測ほど $E[u\mid\bm{x}]$ が小さくなり、推定への寄与が抑えられる——これがt分布のロバスト性の正体です。
$\nu \to \infty$ で正規分布に戻る
ガンマ分布 $\mathrm{Gamma}(\nu/2, \nu/2)$ は平均1・分散 $2/\nu$ ですから、$\nu \to \infty$ で $u$ は $1$ に確率収束します。つまり「ランダムスケール $u$」のばらつきがなくなり、$\bm{X} \mid u=1 \sim N_p(\bm{\mu}, \Sigma)$ に収束します。実際、密度式でも
$$ \lim_{\nu \to \infty} \left[1 + \frac{d^2}{\nu}\right]^{-(\nu+p)/2} = e^{-d^2/2} $$
となり（$\lim (1+x/n)^n = e^x$ の応用）、規格化定数のガンマ関数の比も $(2\pi)^{p/2}$ に向かい、密度全体が多変量正規分布の形に収束することが確認できます。
ガウス・スケールミクスチャ表現により、t分布が「正規分布の自然な拡張」であることがはっきりしました。次に、自由度 $\nu$ が密度の裾をどう変えるかを、もう少し定量的に見ていきましょう。
自由度νと裾の制御
裾の減衰の速さ
$d^2$ が十分大きいとき、密度の主要項は $[1 + d^2/\nu]^{-(\nu+p)/2} \approx (d^2/\nu)^{-(\nu+p)/2}$ となります。すなわち
$$ p(\bm{x}) \sim C \cdot d^{-(\nu+p)} \quad (\text{as } d \to \infty) $$
密度がベキ的に減衰する——これが「重い裾」の数学的内容です。指数関数 $e^{-d^2/2}$ よりも遥かに緩やかな減衰で、極端値の確率が残ります。減衰の指数 $\nu + p$ は、$\nu$ と次元 $p$ で決まります。
- $\nu$ が小さいほどベキ指数が小さく、裾が太い。
- $\nu \to \infty$ で指数関数減衰（正規分布）に滑らかに移行する。
- 次元 $p$ が高いほど指数 $\nu + p$ が大きくなり、相対的に裾が締まる。
モーメントが消える境界
裾が重すぎると、モーメント（平均・分散・尖度など）が存在しなくなります。前にも触れた通り:
- $\nu \leq 1$：平均すら定義されない（コーシー分布的）
- $1 < \nu \leq 2$：平均は存在するが、分散が無限大
- $2 < \nu \leq 4$：分散は存在するが、尖度が無限大
- $\nu > 4$：尖度（4次モーメント）まで存在し、$\mathrm{Kurt}[X] = 3 + \dfrac{6}{\nu – 4}$（1次元の場合）
金融データではよく $\nu \approx 4$〜$6$ 程度が推定されると言われ、これは「分散はあるが尖度はギリギリ」という状況に対応します。実は$\nu = \infty$ の正規分布の超過尖度 0 と比べると、$\nu = 5$ なら超過尖度 $6$、$\nu = 7$ でも超過尖度 $2$ と、裾の太さの違いは数値で見ても大きいことがわかります。
Mahalanobis距離の分布
多変量正規分布では Mahalanobis距離の二乗 $D^2 = (\bm{X}-\bm{\mu})^\top \Sigma^{-1}(\bm{X}-\bm{\mu})$ がカイ二乗分布 $\chi^2_p$ に従いました。多変量t分布では、この量がスケーリングされたF分布になります。
$$ \frac{D^2}{p} \sim F_{p,\nu} $$
ここで $F_{p,\nu}$ は自由度 $(p, \nu)$ のF分布です。これは、t分布の構成 $\bm{X} = \bm{\mu} + \bm{Z}/\sqrt{V/\nu}$ から、$D^2 = \bm{Z}^\top \Sigma^{-1} \bm{Z} \cdot \nu/V$ が $\chi^2_p \cdot \nu / \chi^2_\nu$ の比率になることから来ています。これはF分布の定義そのものです。
この性質は 外れ値検出 で有用です。正規分布なら $D^2$ が $\chi^2_p$ の上側分位点を超えたら外れ値、と判定しましたが、多変量t分布の枠組みでは F分布の分位点に置き換わります。同じ閾値でもt分布の方が緩く判定するため、外れ値だらけのデータでも誤検出が抑えられます。
ここまでで多変量t分布の確率的性質が一通り見えました。次は実用上もっとも重要な問題——データから $\bm{\mu}, \Sigma, \nu$ をどう推定するか——に進みます。ガウス・スケールミクスチャ表現が、EMアルゴリズムによる解析的なM-stepを可能にします。
EMアルゴリズムでの最尤推定
完全データ尤度の準備
$N$ 個の観測 $\{\bm{x}_n\}_{n=1}^N$ が i.i.d. に $\mathcal{T}_p(\bm{\mu}, \Sigma, \nu)$ に従うとします。各観測には潜在スケール変数 $u_n$ が対応していると考え、ガウス・スケールミクスチャ表現
$$ \bm{x}_n \mid u_n \sim N_p(\bm{\mu}, \Sigma/u_n), \quad u_n \sim \mathrm{Gamma}(\nu/2, \nu/2) $$
を完全データモデルとします。観測データだけの対数尤度 $\log p(\bm{x}_n \mid \bm{\mu}, \Sigma, \nu)$ は $\bm{\mu}, \Sigma$ について閉形式の最大化ができない（多項式項が分母に入る）ため、潜在変数 $u_n$ を導入してEMで攻めるのが定石です。
完全データ対数尤度 $\log p(\bm{x}_n, u_n \mid \bm{\mu}, \Sigma, \nu) = \log p(\bm{x}_n \mid u_n, \bm{\mu}, \Sigma) + \log p(u_n \mid \nu)$ を書き下すと
$$ \log p(\bm{x}_n, u_n \mid \theta) = -\frac{p}{2}\log(2\pi) – \frac{1}{2}\log|\Sigma| + \frac{p}{2}\log u_n – \frac{u_n}{2} d_n^2 + \frac{\nu}{2}\log\frac{\nu}{2} – \log\Gamma(\nu/2) + \left(\tfrac{\nu}{2} – 1\right)\log u_n – \frac{\nu u_n}{2} $$
ここで $d_n^2 = (\bm{x}_n – \bm{\mu})^\top \Sigma^{-1}(\bm{x}_n – \bm{\mu})$ です。
E-step：潜在変数の事後
現在の推定値 $\theta^{(t)} = (\bm{\mu}^{(t)}, \Sigma^{(t)}, \nu^{(t)})$ のもとで、$u_n$ の事後分布は前節で求めた通り
$$ u_n \mid \bm{x}_n \sim \mathrm{Gamma}\!\left(\frac{\nu^{(t)} + p}{2},\; \frac{\nu^{(t)} + d_n^{2,(t)}}{2}\right) $$
E-stepで必要な十分統計量は次の2つです。
$$ w_n^{(t)} := E[u_n \mid \bm{x}_n, \theta^{(t)}] = \frac{\nu^{(t)} + p}{\nu^{(t)} + d_n^{2,(t)}} $$
$$ E[\log u_n \mid \bm{x}_n, \theta^{(t)}] = \psi\!\left(\frac{\nu^{(t)} + p}{2}\right) – \log\!\left(\frac{\nu^{(t)} + d_n^{2,(t)}}{2}\right) $$
ここで $\psi(\cdot)$ はディガンマ関数です。$w_n^{(t)}$ こそが「観測ごとの重み」で、後で見るように外れ値（$d_n^2$ が大）は自動的に小さな重みになります。
M-step（$\bm{\mu}, \Sigma$）：重み付き正規分布の最尤推定
完全データ対数尤度の期待値を $Q(\theta \mid \theta^{(t)}) = \sum_n E[\log p(\bm{x}_n, u_n \mid \theta)]$ とすると、$\bm{\mu}, \Sigma$ に関する項は $E[u_n] = w_n^{(t)}$ で置き換わって
$$ Q(\bm{\mu}, \Sigma) = -\frac{N}{2}\log|\Sigma| – \frac{1}{2}\sum_n w_n^{(t)} (\bm{x}_n – \bm{\mu})^\top \Sigma^{-1} (\bm{x}_n – \bm{\mu}) + \text{const} $$
これは重み付き多変量正規分布の最尤推定そのものです。$\bm{\mu}$ に関して偏微分してゼロとおき、
$$ \bm{\mu}^{(t+1)} = \frac{\sum_n w_n^{(t)} \bm{x}_n}{\sum_n w_n^{(t)}} $$
続いて $\Sigma$ に関して偏微分してゼロとおくと、
$$ \Sigma^{(t+1)} = \frac{1}{N} \sum_n w_n^{(t)} (\bm{x}_n – \bm{\mu}^{(t+1)})(\bm{x}_n – \bm{\mu}^{(t+1)})^\top $$
正規分布のサンプル平均・サンプル共分散が、重み $w_n$ 付きの加重平均・加重共分散に置き換わっただけです。実装はサンプル統計量の修正版で、極めて素直に書けます。
M-step（$\nu$）：1次元数値解
自由度 $\nu$ に関する項は
$$ Q(\nu) = \sum_n \left[\frac{\nu}{2}\log\frac{\nu}{2} – \log\Gamma(\nu/2) + \frac{\nu}{2}\big(E[\log u_n] – E[u_n]\big)\right] $$
これは $\nu$ について凸ではない超越方程式になり、閉形式の解はありません。実用的には、$\partial Q / \partial \nu = 0$ をブレント法・二分法・ニュートン法などの1次元数値求解で解きます。具体的には
$$ -\psi(\nu/2) + \log(\nu/2) + 1 + \frac{1}{N}\sum_n \big(E[\log u_n] – E[u_n]\big) = 0 $$
を $\nu > 0$ の範囲で解きます。$N$ が大きければよく振る舞う方程式で、scipyの brentq などで安定に解けます。
アルゴリズム全体
以上をまとめると、EMアルゴリズムは次のループになります。
- 初期化：$\bm{\mu}^{(0)}, \Sigma^{(0)}$ をサンプル平均・サンプル共分散、$\nu^{(0)}$ を適当な値（例：10）に。
- E-step：各 $n$ で $d_n^{2,(t)}$ を計算し、$w_n^{(t)} = (\nu^{(t)} + p)/(\nu^{(t)} + d_n^{2,(t)})$ を更新。
- M-step：上の式で $\bm{\mu}^{(t+1)}, \Sigma^{(t+1)}, \nu^{(t+1)}$ を更新。
- 対数尤度の変化が閾値以下になるまで反復。
EMの一般理論より、対数尤度は単調非減少で、局所最適値に収束します。ロバスト統計の文脈では、この $w_n$ がIRLS（重み付き最小二乗の反復解法）における重みと同じ役割を果たし、外れ値の影響を自動で削っていきます。
理論はここまで揃いました。次セクションで、実際にこのEMをPythonで実装し、外れ値に対してどう振る舞うかを GMMと比較しながら確認します。
Python実装 — ロバスト混合モデル
単一の多変量t分布のEM推定
まず混合モデルに行く前に、単一の多変量t分布のパラメータ推定を実装し、外れ値への耐性を確認します。
import numpy as np
from scipy.special import digamma, gammaln
from scipy.optimize import brentq
def mahalanobis_sq(X, mu, Sigma_inv):
"""各行ベクトルの Mahalanobis距離の二乗を計算"""
diff = X - mu
return np.einsum("ni,ij,nj->n", diff, Sigma_inv, diff)
def fit_mvt_em(X, max_iter=100, tol=1e-6, nu_init=10.0):
"""単一多変量t分布のEM推定。
Returns: mu, Sigma, nu, log_likelihood_history
"""
N, p = X.shape
# 初期化：サンプル平均・サンプル共分散
mu = X.mean(axis=0)
Sigma = np.cov(X, rowvar=False) + 1e-6 * np.eye(p)
nu = float(nu_init)
ll_history = []
for it in range(max_iter):
# E-step：観測ごとの重み w_n と E[log u_n]
Sigma_inv = np.linalg.inv(Sigma)
d2 = mahalanobis_sq(X, mu, Sigma_inv)
w = (nu + p) / (nu + d2)
E_log_u = digamma((nu + p) / 2) - np.log((nu + d2) / 2)
# M-step：mu, Sigma
mu_new = (w[:, None] * X).sum(axis=0) / w.sum()
diff = X - mu_new
Sigma_new = (w[:, None, None] * diff[:, :, None] * diff[:, None, :]).sum(axis=0) / N
# M-step：nu（1次元方程式を Brent法で解く）
mean_term = (E_log_u - w).mean()
def eq(nu_val):
return -digamma(nu_val / 2) + np.log(nu_val / 2) + 1 + mean_term
try:
nu_new = brentq(eq, 1e-3, 200.0)
except ValueError:
nu_new = nu # 解が見つからない場合は更新しない
# 対数尤度（観測データ）
Sigma_new_inv = np.linalg.inv(Sigma_new)
d2_new = mahalanobis_sq(X, mu_new, Sigma_new_inv)
log_det = np.linalg.slogdet(Sigma_new)[1]
ll = (gammaln((nu_new + p) / 2) - gammaln(nu_new / 2)
- 0.5 * p * np.log(nu_new * np.pi) - 0.5 * log_det
- 0.5 * (nu_new + p) * np.log1p(d2_new / nu_new))
ll_total = ll.sum()
ll_history.append(ll_total)
# 収束判定
if it > 0 and abs(ll_history[-1] - ll_history[-2]) < tol:
mu, Sigma, nu = mu_new, Sigma_new, nu_new
break
mu, Sigma, nu = mu_new, Sigma_new, nu_new
return mu, Sigma, nu, ll_history
このコードの肝は2つです。一つは E-step で計算する重み w = (nu + p) / (nu + d2) で、これが Mahalanobis距離の大きい点（外れ値候補）に対して自動的に小さな値を与えます。もう一つは M-step が「重み付きサンプル平均・サンプル共分散」になっている点で、ガウス分布のMLEを最小限の修正で流用できる構造です。$\nu$ の更新だけは1次元数値求解が必要ですが、brentq で安定に解けます。
動作確認：正規分布 vs t分布 vs 外れ値混入データ
次に、合成データでEMが正しく動くか確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_t
# 真のパラメータ
np.random.seed(42)
mu_true = np.array([1.0, -1.0])
Sigma_true = np.array([[2.0, 1.0], [1.0, 1.5]])
nu_true = 4.0
N = 500
# 多変量t分布からサンプリング
rv = multivariate_t(loc=mu_true, shape=Sigma_true, df=nu_true, seed=42)
X = rv.rvs(size=N)
# EM推定
mu_hat, Sigma_hat, nu_hat, ll_hist = fit_mvt_em(X)
print(f"真の mu : {mu_true}")
print(f"推定 mu : {np.round(mu_hat, 3)}")
print(f"真の Sigma :\n{Sigma_true}")
print(f"推定 Sigma :\n{np.round(Sigma_hat, 3)}")
print(f"真の nu : {nu_true}")
print(f"推定 nu : {nu_hat:.3f}")
print(f"対数尤度の収束: {ll_hist[0]:.2f} -> {ll_hist[-1]:.2f} ({len(ll_hist)} iters)")
実行すると、$\bm{\mu}$、$\Sigma$、$\nu$ がいずれも真の値の近く（典型的に $\bm{\mu}$ で誤差 0.05〜0.1、$\Sigma$ で 0.1〜0.3、$\nu$ で $\pm 1$ 程度）に収束し、対数尤度が単調増加することが確認できます。$N = 500$ サンプルでこの精度なら実用レベルです。$\nu$ の推定誤差が他より大きいのは、自由度パラメータが裾の形からしか情報を得られず、推定が本質的に難しいためです。
多変量正規 vs 多変量t の裾の可視化
両分布の密度を等高線で重ね描きし、裾の太さの違いを見ます。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal, multivariate_t
mu = np.zeros(2)
Sigma = np.array([[1.0, 0.6], [0.6, 1.0]])
xx, yy = np.meshgrid(np.linspace(-6, 6, 200), np.linspace(-6, 6, 200))
grid = np.stack([xx.ravel(), yy.ravel()], axis=1)
mvn = multivariate_normal(mu, Sigma)
fig, axes = plt.subplots(1, 3, figsize=(15, 5), sharex=True, sharey=True)
# 等高線レベル（対数スケールで広く取って裾を強調）
levels = np.logspace(-6, -1, 10)
for ax, nu in zip(axes, [3, 10, 100]):
mvt = multivariate_t(mu, Sigma, df=nu)
pdf_t = mvt.pdf(grid).reshape(xx.shape)
pdf_n = mvn.pdf(grid).reshape(xx.shape)
ax.contour(xx, yy, pdf_t, levels=levels, colors='red', linewidths=1.5)
ax.contour(xx, yy, pdf_n, levels=levels, colors='blue', linewidths=1.0, linestyles='--')
ax.set_title(f'MVT df={nu} (red) vs MVN (blue dashed)')
ax.set_xlabel('x1')
ax.set_ylabel('x2')
ax.grid(alpha=0.3)
ax.set_aspect('equal')
plt.tight_layout()
plt.savefig('mvt_vs_mvn_contours.png', dpi=150, bbox_inches='tight')
plt.show()
このプロットから3つの重要な特徴が読み取れます。第一に、$\nu = 3$ の多変量t分布（赤実線）の外側等高線は、同じスケールの正規分布（青破線）よりも明らかに大きく外側に広がっています。低密度領域での確率密度がはるかに大きいことを意味し、これがファットテールの正体です。第二に、$\nu = 10$ では赤と青の差はかなり縮まりますが、最外側ではまだt分布の方が太い。第三に、$\nu = 100$ では赤と青がほぼ重なり、$\nu \to \infty$ で正規分布に収束する性質が視覚的に確認できます。中心部の高密度領域では3つともほぼ重なり、「中央は変わらず、裾だけが太くなる」というt分布の役割が明瞭です。
ロバスト混合モデル（Mixture of t）
単一t分布を $K$ 成分の混合に拡張し、外れ値に強いクラスタリングを実装します。GMMのEMにスケール変数 $u_{nk}$ を追加するだけで、構造は素直に拡張できます。
import numpy as np
from scipy.special import digamma, gammaln
from scipy.optimize import brentq
def fit_mot_em(X, K, max_iter=100, tol=1e-5, nu_init=10.0, seed=0):
"""Mixture of multivariate t-distributions（MoT）の EM 推定。
各成分 k は pi_k, mu_k, Sigma_k, nu_k を持つ。
"""
rng = np.random.default_rng(seed)
N, p = X.shape
# 初期化：ランダム選択
idx = rng.choice(N, K, replace=False)
mus = X[idx].copy()
Sigmas = np.stack([np.cov(X, rowvar=False) + 1e-6 * np.eye(p) for _ in range(K)])
nus = np.full(K, float(nu_init))
pis = np.full(K, 1.0 / K)
def comp_log_pdf(X, mu, Sigma, nu):
p_dim = X.shape[1]
Sigma_inv = np.linalg.inv(Sigma)
diff = X - mu
d2 = np.einsum("ni,ij,nj->n", diff, Sigma_inv, diff)
log_det = np.linalg.slogdet(Sigma)[1]
return (gammaln((nu + p_dim) / 2) - gammaln(nu / 2)
- 0.5 * p_dim * np.log(nu * np.pi) - 0.5 * log_det
- 0.5 * (nu + p_dim) * np.log1p(d2 / nu)), d2
ll_hist = []
for it in range(max_iter):
# E-step：成分責任 r_nk と スケール期待値 u_nk
log_pdfs = np.zeros((N, K))
d2_all = np.zeros((N, K))
for k in range(K):
log_pdfs[:, k], d2_all[:, k] = comp_log_pdf(X, mus[k], Sigmas[k], nus[k])
log_weighted = log_pdfs + np.log(pis)
log_norm = np.logaddexp.reduce(log_weighted, axis=1)
log_r = log_weighted - log_norm[:, None]
r = np.exp(log_r) # 責任 (N, K)
u = (nus + p) / (nus + d2_all) # スケール期待 (N, K)
ru = r * u # 有効重み
# M-step：pi, mu, Sigma, nu を成分ごとに更新
Nk = r.sum(axis=0)
pis = Nk / N
for k in range(K):
mus[k] = (ru[:, k:k+1] * X).sum(axis=0) / ru[:, k].sum()
diff = X - mus[k]
Sigmas[k] = (ru[:, k:k+1, None] * diff[:, :, None] * diff[:, None, :]).sum(axis=0) / Nk[k]
Sigmas[k] += 1e-6 * np.eye(p) # 正則化
# nu_k の更新
E_log_u_k = digamma((nus[k] + p) / 2) - np.log((nus[k] + d2_all[:, k]) / 2)
mean_term = (r[:, k] * (E_log_u_k - u[:, k])).sum() / Nk[k]
def eq(nu_val):
return -digamma(nu_val / 2) + np.log(nu_val / 2) + 1 + mean_term
try:
nus[k] = brentq(eq, 1e-3, 200.0)
except ValueError:
pass
ll = log_norm.sum()
ll_hist.append(ll)
if it > 0 and abs(ll_hist[-1] - ll_hist[-2]) < tol:
break
return {"pis": pis, "mus": mus, "Sigmas": Sigmas, "nus": nus,
"responsibilities": r, "log_likelihood": ll_hist}
このMoTのEMはGMMのEMと骨格はほぼ同じです。違いは、各観測 $\bm{x}_n$ が成分 $k$ に対して「責任」$r_{nk}$ と「スケール期待値」$u_{nk}$ の2つの重みを持つ点で、M-stepの平均・共分散は両者の積 $r_{nk} u_{nk}$ で重み付けされます。外れ値は $u_{nk}$ を介して自動的に小さな寄与に抑えられ、GMMでは平均が引きずられる問題が緩和されます。
GMM vs MoT：外れ値混入下のクラスタリング比較
実際に、2つのクラスタに10%の外れ値を加えたデータで両モデルを比較します。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.mixture import GaussianMixture
# データ生成：2クラスタ + 外れ値
np.random.seed(1)
n_inlier = 200
mu1 = np.array([0, 0])
mu2 = np.array([5, 5])
Sigma1 = np.array([[1.0, 0.5], [0.5, 1.0]])
Sigma2 = np.array([[1.5, -0.3], [-0.3, 1.0]])
X1 = np.random.multivariate_normal(mu1, Sigma1, n_inlier)
X2 = np.random.multivariate_normal(mu2, Sigma2, n_inlier)
# 外れ値：広い一様分布から
n_out = 40
X_out = np.random.uniform(-10, 15, size=(n_out, 2))
X = np.vstack([X1, X2, X_out])
# GMM (sklearn)
gmm = GaussianMixture(n_components=2, random_state=0, n_init=5).fit(X)
labels_gmm = gmm.predict(X)
# MoT (自作)
result = fit_mot_em(X, K=2, seed=0)
labels_mot = result["responsibilities"].argmax(axis=1)
# 推定された中心を比較
print("GMM means:")
print(np.round(gmm.means_, 2))
print("MoT means:")
print(np.round(result["mus"], 2))
print("MoT estimated nu:", np.round(result["nus"], 2))
print("True means: [0,0] and [5,5]")
# プロット
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
for ax, labels, means, title in [
(axes[0], labels_gmm, gmm.means_, 'GMM'),
(axes[1], labels_mot, result["mus"], 'MoT (Mixture of t)'),
]:
ax.scatter(X[:, 0], X[:, 1], c=labels, cmap='coolwarm', s=15, alpha=0.6)
ax.scatter(means[:, 0], means[:, 1], c='black', marker='X', s=200, label='center')
ax.scatter([mu1[0], mu2[0]], [mu1[1], mu2[1]], c='gold', marker='*',
s=300, edgecolors='black', label='true center')
ax.set_title(title)
ax.legend()
ax.grid(alpha=0.3)
ax.set_xlim(-12, 17)
ax.set_ylim(-12, 17)
plt.tight_layout()
plt.savefig('gmm_vs_mot.png', dpi=150, bbox_inches='tight')
plt.show()
この比較から、MoTの優位性がはっきり読み取れます。GMMの推定中心（黒のX印）は、外れ値の影響で真の中心（金の星）から目に見える距離だけずれる傾向があります。特に、共分散行列も外れ値を取り込もうとして大きく膨らみ、結果として2クラスタの境界が曖昧になりがちです。一方、MoTの中心はほぼ真の値に張り付きます。これはEM内部で外れ値が小さな $u_{nk}$ を獲得し、平均・共分散の計算でほとんど無視されるからです。推定された $\nu$ も $5$ 未満の小さな値になり、データが「重い裾を持つ」とアルゴリズム自身が認識していることがわかります。
理論と実装を一通り通したので、最後に金融リスクへの応用を見て本記事を締めくくります。
応用 — 金融リスク・異常検知
VaR と CVaR：尾の上のリスク指標
ポートフォリオの損失 $L$ に対し、Value-at-Risk（VaR） とは「確率 $\alpha$ で損失がこれ以下に収まる」という閾値です。
$$ \mathrm{VaR}_\alpha(L) = \inf\{l : P(L \leq l) \geq \alpha\} $$
$\alpha = 0.99$ なら「99%の日は損失がこの値以下」、すなわち「1%の最悪日に超える損失額」を意味します。Conditional VaR（CVaR、別名Expected Shortfall） はその外側の条件付き期待値で、
$$ \mathrm{CVaR}_\alpha(L) = E[L \mid L \geq \mathrm{VaR}_\alpha] $$
VaRが「閾値」、CVaRが「閾値を超えた場合の平均損失」です。バーゼル規制の文脈でも、CVaRの方が裾の情報を捉えると評価され、近年は標準指標になっています。
正規 vs t分布：VaRの差
リターンを多変量正規でモデル化するか多変量tでモデル化するかで、VaR・CVaRの値は劇的に変わります。これを実データ風シミュレーションで確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal, multivariate_t, norm, t as student_t
# 3資産ポートフォリオを想定
np.random.seed(7)
p = 3
mu_r = np.array([0.0005, 0.0008, 0.0003]) # 日次平均リターン
# 共分散行列（年率 25% 程度のボラ、相関 0.4）
sigma_daily = 0.25 / np.sqrt(252)
corr = np.array([[1.0, 0.4, 0.3],
[0.4, 1.0, 0.5],
[0.3, 0.5, 1.0]])
Sigma_r = sigma_daily**2 * corr
# 真のデータ生成は t 分布 (nu=5) — 現実的な金融データ
true_nu = 5.0
# t分布の共分散は Sigma * nu/(nu-2) なので、Sigma を逆スケール
Sigma_t_scale = Sigma_r * (true_nu - 2) / true_nu
n_days = 2000
rv_t = multivariate_t(loc=mu_r, shape=Sigma_t_scale, df=true_nu, seed=7)
returns = rv_t.rvs(n_days)
# 等ウェイトポートフォリオの損失
weights = np.array([1/3, 1/3, 1/3])
pf_returns = returns @ weights
pf_losses = -pf_returns
# モデル1：正規分布でフィット
mu_n = pf_returns.mean()
sig_n = pf_returns.std(ddof=1)
# モデル2：1次元t分布でフィット
nu_fit, loc_fit, scale_fit = student_t.fit(pf_returns)
# VaR / CVaR を解析的に計算（99%）
alpha = 0.99
# 正規分布の VaR / CVaR
var_n = -mu_n + sig_n * norm.ppf(alpha)
cvar_n = -mu_n + sig_n * norm.pdf(norm.ppf(alpha)) / (1 - alpha)
# t分布の VaR / CVaR（t分布は対称なので符号を反転）
t_q = student_t.ppf(alpha, df=nu_fit)
var_t = -loc_fit + scale_fit * t_q
pdf_at_q = student_t.pdf(t_q, df=nu_fit)
cvar_t = -loc_fit + scale_fit * pdf_at_q / (1 - alpha) * (nu_fit + t_q**2) / (nu_fit - 1)
# 経験的 VaR / CVaR（参考）
emp_var = np.quantile(pf_losses, alpha)
emp_cvar = pf_losses[pf_losses >= emp_var].mean()
print(f"推定された t分布の自由度 nu: {nu_fit:.2f}")
print(f"\n--- 99% VaR / CVaR ---")
print(f"正規モデル: VaR={var_n*100:.3f}% CVaR={cvar_n*100:.3f}%")
print(f"t分布モデル: VaR={var_t*100:.3f}% CVaR={cvar_t*100:.3f}%")
print(f"経験(実測): VaR={emp_var*100:.3f}% CVaR={emp_cvar*100:.3f}%")
この出力から、3点が読み取れます。第一に、推定された自由度 $\nu$ は真の値 5 に近い値（典型的に4.5〜5.5程度）になり、t分布が裾の太さを正しく拾っていることがわかります。第二に、99% VaRでも正規モデルとt分布モデルで数値が違いますが、その差はCVaRで決定的に拡大します——t分布のCVaRは正規のそれより2割〜5割大きい値になり、これが規制資本やリスク予算の必要量に直接効きます。第三に、t分布のVaR・CVaRは経験値（実測）にずっと近く、正規モデルは裾を過小評価していることが定量的に確認できます。
損失分布の可視化
最後に、損失のヒストグラムに2つのモデルの密度を重ねて見ます。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, t as student_t
fig, ax = plt.subplots(figsize=(10, 6))
ax.hist(pf_losses, bins=80, density=True, alpha=0.5, color='gray', label='Empirical losses')
x_grid = np.linspace(pf_losses.min(), pf_losses.max(), 500)
ax.plot(x_grid, norm.pdf(x_grid, loc=-mu_n, scale=sig_n),
'b-', lw=2, label='Normal fit')
ax.plot(x_grid, student_t.pdf(x_grid, df=nu_fit, loc=-loc_fit, scale=scale_fit),
'r-', lw=2, label=f't fit (df={nu_fit:.2f})')
ax.axvline(var_n, color='blue', ls='--', alpha=0.7, label=f'Normal VaR99={var_n*100:.2f}%')
ax.axvline(var_t, color='red', ls='--', alpha=0.7, label=f't VaR99={var_t*100:.2f}%')
ax.axvline(emp_var, color='black', ls=':', alpha=0.7, label=f'Empirical VaR99={emp_var*100:.2f}%')
ax.set_xlabel('Daily portfolio loss')
ax.set_ylabel('Density')
ax.set_title('Loss distribution: Normal vs t fit (true df=5)')
ax.legend()
ax.set_yscale('log') # 裾を見やすく
ax.set_ylim(1e-3, ax.get_ylim()[1])
ax.grid(alpha=0.3, which='both')
plt.tight_layout()
plt.savefig('var_loss_distribution.png', dpi=150, bbox_inches='tight')
plt.show()
ヒストグラム（灰色）とt分布フィット（赤）は裾までほぼ重なる一方、正規分布フィット（青）はヒストグラムの裾よりも明らかに早く落ち込みます。対数縦軸で見ると差は一目瞭然で、損失 2% を超えるあたりから正規分布の密度はt分布の半分以下になり、「正規モデルでは数百日に1度しか起こらないはずの大損失が、実際には数十日に1度起きている」という現象が定量的に確認できます。VaRの縦線も、青（正規）が経験値（黒点線）より内側にあり、リスクを過小評価していることが視覚的にわかります。
Mahalanobis距離による異常検知
最後に、多変量t分布の枠組みでの異常検知も簡単に触れておきます。前述したように、$D^2/p \sim F_{p,\nu}$ なので、F分布の上側分位点を閾値として「裾の重さを許容した外れ値検出」ができます。
import numpy as np
from scipy.stats import f as fdist, multivariate_t
# 2次元データを生成（内側はt分布、明確な外れ値を混ぜる）
np.random.seed(3)
p = 2
rv = multivariate_t(loc=[0, 0], shape=[[1, 0.5], [0.5, 1]], df=4, seed=3)
X_inlier = rv.rvs(300)
X_outlier = np.array([[8, -7], [-9, 8], [10, 10]])
X_all = np.vstack([X_inlier, X_outlier])
# 多変量t分布でフィット
mu_h, Sig_h, nu_h, _ = fit_mvt_em(X_all)
Sig_inv = np.linalg.inv(Sig_h)
d2 = mahalanobis_sq(X_all, mu_h, Sig_inv)
# 閾値：F分布の 99% 分位点
thr = p * fdist.ppf(0.99, p, nu_h)
mask = d2 > thr
print(f"推定 nu: {nu_h:.2f}")
print(f"閾値 D^2 (99% F分布): {thr:.2f}")
print(f"検出された外れ値数: {mask.sum()} / {len(X_all)}")
print(f"明確な外れ値 (最後の3つ) の D^2: {np.round(d2[-3:], 2)}")
実行すると、明確に外れ値として注入した3点はいずれも閾値を大きく超え、検出されます。重要なのは、内側のt分布から生成された点で「軽くはみ出している」ものは F分布の閾値内に収まり、誤検出が少ない点です。同じデータを正規分布の $\chi^2_p$ 閾値で判定すると、内側点の多くも「外れ値」と誤判定されてしまうでしょう。裾の重い母集団からのデータでは、裾の重い参照分布で判定するのが筋——これが多変量t分布の異常検知における価値です。
まとめ
本記事では、多変量t分布を理論・実装・応用の3つの観点から解説しました。
- 裾の重さを連続的に制御するパラメータ $\nu$：密度がベキ的に減衰し、$\nu \to \infty$ で多変量正規分布に滑らかに収束する。$\nu$ が小さいほど外れ値の確率が大きく、ファットテールが自然に表現される。
- ガウス・スケールミクスチャ表現：$\bm{X} \mid u \sim N(\bm{\mu}, \Sigma/u)$、$u \sim \mathrm{Gamma}(\nu/2, \nu/2)$ という階層構造として書ける。観測ごとに正規分布の幅がランダムに揺らぐ生成過程として理解でき、潜在変数 $u$ の事後期待値 $E[u \mid \bm{x}] = (\nu+p)/(\nu+d^2)$ がそのまま観測の重みになる。
- 平均・共分散の存在条件：$\nu > 1$ で平均、$\nu > 2$ で共分散 $\mathrm{Cov}=\nu/(\nu-2)\cdot\Sigma$ が存在。$\Sigma$ は共分散ではなくスケール行列であることに注意。Mahalanobis距離の二乗は $D^2/p \sim F_{p,\nu}$。
- EMアルゴリズム：完全データ尤度をスケール変数 $u_n$ で書き、E-stepで重み $w_n = (\nu+p)/(\nu+d_n^2)$ と $E[\log u_n]$ を計算、M-stepで $\bm{\mu}, \Sigma$ は加重サンプル平均・共分散、$\nu$ は1次元数値求解で更新する。
- ロバスト混合モデル（MoT）：GMMのEMに $u_{nk}$ を追加するだけで、外れ値に頑強なクラスタリングが実現する。実験では真の中心がほぼ正確に推定され、GMMでは引きずられる問題が解消されることを確認した。
- 金融リスク・異常検知への応用：VaR・CVaRがt分布モデルで現実的な値になり、正規モデルでは過小評価される。F分布閾値による外れ値検出は、裾の重い母集団でも誤検出を抑えられる。
多変量t分布は、「正規分布が便利すぎて、現実が見えなくなる」というアナリストの罠への解毒剤と言えます。中央付近のモデル化能力を保ったまま、裾だけを現実に合わせて太くできる——この性質が、ロバスト統計・金融工学・機械学習の多くの場面で支持される理由です。
次のステップとして、以下の記事も参考にしてください。