仕事の都合で、分散共分散行列に情報的な逆ウィシャート事前分布を与えるとき、そのパラメータをどう解釈したらいいのかがわからずに困っていた。いろいろ探した結果、なんと灯台もと暗し、Muthen一家の若頭Asparouhovさんが、哀れなMplusユーザ向けに親切な解説を書いてくれていた。
事前分布として逆ウィシャート分布を使うことが自体がモダンじゃない、というご意見もあろうかと思うが、そんな頭の良い君たちにいいたい。その良く回る口を閉じ、おうちに帰ってStanで遊んでな、と。真の男が愛するMplusの世界では、宇宙を統べる方程式$$ y = \nu + \Lambda \eta + K x + \varepsilon, \ \ \varepsilon \sim MVN(0, \Theta) $$ $$ \eta = \alpha + B \eta + \Gamma x + \zeta, \ \ \zeta \sim MVN(0, \Psi)$$の\(\Psi\)の事前分布は逆ウィシャート分布一択なのだ。
おさらい
えーっと、いきなり兄貴の意見を聴くのは畏れ多いので、逆ウィシャート分布についてちょっと復習しておこう。以下よりメモ。
Rossi, Allenby, McCulloch (2005) “Bayesian Statistics and Marketing“, Wiley. 2.8.3 Bayesian Inference for Covariance Metrices.
まずは\(m\)次元の多変量正規分布について考えよう。話を簡単にするために平均を無視し、共分散行列を\(\Sigma\)として、ベクトル\(y\)の同時密度は $$ f(y) = \frac{1}{(\sqrt{2 \pi})^m \sqrt{|\Sigma|}} \exp\left(-\frac{1}{2} y^\top \Sigma^{-1} y \right)$$ である[ここはwikipediaをみながら書きました]。ってことは、$$ p(y_1, \ldots, y_n | \Sigma) \propto \prod_{i=1}^n |\Sigma|^{-\frac{1}{2}} \exp \left(-\frac{1}{2} y^\top_i \Sigma^{-1} y_i \right) $$ ですわね。
右辺を書き換えるぜ。定数を外に出し、\(\prod\)を\(\exp\)のなかに繰り込んで $$ = |\Sigma|^{-\frac{n}{2}} \exp \left( -\frac{1}{2} \sum_i^n y_i^\top \Sigma^{-1} y_i \right) $$ この式の \(y_i^\top \Sigma^{-1} y_i\)というところは、結局は一本の多項式になるわけで、だからそこにトレース記号をかぶせても変わらない。さて、トレースの内側では積の順番を変えられるから$$ \mathrm{tr} (y^\top_i \Sigma^{-1} y_i) = \mathrm{tr}(y_i y^\top_i \Sigma^{-1}) $$ 清書します。\(S = \sum y_i y^\top_i\)とする。簡略のため \(\mathrm{etr}(\cdot) = \exp(\mathrm{tr}(\cdot)) \)とする。$$ p(y_1, \ldots, y_n | \Sigma) \propto |\Sigma|^{-\frac{n}{2}} \mathrm{etr} \left( -\frac{1}{2} S \Sigma^{-1} \right) $$
それではご紹介しましょう。逆ウィシャート分布\(IW(\upsilon_0, V_0)\)です！
$$ p(\Sigma | \upsilon_0, V_0) \propto |\Sigma|^{-\frac{\upsilon_0 + m + 1}{2}} \mathrm{etr} \left( – \frac{1}{2} V_0 \Sigma^{-1} \right)$$ つまり、\(V_0\)は分散共分散行列\(S\)に相当し、\(\upsilon_0 + m + 1\)が標本サイズ \(n\)に相当するわけね。
密度関数は超面倒くさい式になるので書かないけど、以下の性質をもつ。
- \(\upsilon_0 > m\) のときに積分可能。
- \(\upsilon_0 \geq m + 2\)のとき、\(E(\Sigma) = \frac{V_0}{\upsilon – m – 1} \)となる。
- \(\Sigma\)の対角に沿って正方ブロックを切り出してその周辺分布をみると、\(V_0\)からその部分を切り出したのと同じ逆ウィシャート分布に従う。
- \(V_0\)が分布の位置、\(\upsilon_0\)が分布の広がりを表している。ただし、特に\(\upsilon_0\)が小さいとき、強く歪んだ分布になる。逆カイ二乗分布の行列版だと思うとよい。
これは\(\Sigma\)の自然共役事前分布になっております。$$ \Sigma \sim IW(\upsilon_0, V_0) $$ としておくと、事後分布は $$ \Sigma | Y \sim IW(\upsilon_0 + n, V_0 + S) $$ となる。
自然共役事前分布の常として、「なにかの事前分布とどこかよそのデータから得られた事後分布」だと捉えることもできる。Zellnerの事前分布を使ったと考えれば、\(\upsilon_0\)は有効標本サイズで\(V_0\)は二乗和とクロス積の行列[\(S\)ってことね]と解釈できる。つまり\(V_0 = \upsilon \hat{\Sigma}\)と解釈できる。[…]
Asparouhov兄貴からお言葉を賜ります
以下よりメモする。原文の記号に従っているので、上記と記号が異なるし、前半と後半でも記号の意味が異なることに注意。
Asparouhov, T., Muthen, B. (2021) Bayesian Analysis of Latent Variable Models using Mplus. Mplus Technical Appendix. 10. Appdendix A. Priors.
まずはガンマ分布\(IG(\alpha, \beta)\)について考えよう。
この分布は\((0, \infty)\)上に密度を持ち、密度関数は$$ \frac{\beta^\alpha}{\Gamma(\alpha)} x^{-\alpha-1} \exp(-\frac{\beta}{X}) $$ である。平均は\(\alpha > 1\)のときに \(\frac{\beta}{\alpha-1}\), そうでないときは無限となる。分散は\(\alpha > 2\)のときに\(\frac{\beta^2}{(\alpha-1)^2 (\alpha-2)}\), そうでないときに無限となる。
ある分散パラメータに対して情報事前分布として逆ガンマ分布を与えたい場合は、その平均を\(m\), 分散を\(\upsilon\)としたいなら、$$ m = \frac{\beta}{\alpha -1}, \ \upsilon = \frac{\beta^2}{(\alpha-1)^2 (\alpha-2)} $$ とすればいいわけだ。幸いこれは簡単に解けて、$$ \alpha = 2 + \frac{m^2}{\upsilon}, \beta = m + \frac{m^3}{\upsilon} $$ となるね。\(\alpha\)が小さくなりすぎてしまう場合は、かわりに分布の最頻値が\(\frac{\beta}{\alpha + 1}\)であることを思い出すとよいだろう。
パラメータ\(\beta\)はただのスケールパラメータで、\(X \sim IG(\alpha, 1)\)ならば\(\beta X \sim IG(\alpha, \beta)\)である。
ここまではいいね？
では、今度は逆ウィシャート分布\(IW(\Psi, m)\)について考える。
ある期待値を持つ情報事前分布を指定したい場合、以下を思い出すと良い。\(IW(\Psi, m)\)の平均は\(m > p+1\)のときのみ存在し \( \frac{\Psi}{m-p-1}\)である。\(m \leq p+1\)の場合については、かわりに、最頻値が \( \frac{\Psi}{m+p+1}\)であることを思い出すと良い。
分散(つまり情報性のレベル)はパラメータ\(m\)によってコントロールされる。\(m\)が大きいほど情報的になる。
- 対角要素の周辺分布について考えると、\(i\)番目の対角要素の周辺分布は\(IG(\frac{m-p+1}{2}, \frac{\psi_{ii}}{2})\)となる。その平均は\(m > p+1\)のときに\(\frac{\psi_{jj}}{m-p-1}\), その分散は\(m > p + 3\)のときに\(\frac{2 \psi^2_{jj}}{(m-p-1)^2(m-p-3)}\)である。だから、ある分散を持つ情報事前分布を設定したい場合は、その分散から\(m\)を逆算し、それに\(m-p-1\)を掛けて\(\Psi\)を求めればよい。
- いっぽう、非対角要素の周辺分布は閉形式では書けない[自分向けに追記しておくが、\(\Psi\)が対角行列なら書ける。下記を参照]。その平均は\(m > p+1\)のときに\(\frac{\psi_{jj}}{m-p-1}\)であり[対角要素とおなじ]、その分散は\(m > p + 3\)のときに$$ \frac{(m-p+1)\psi^2_{ij} + (m-p-1) \psi_{ii} \psi_{jj}}{(m-p)(m-p-1)^2 (m-p-3)} $$ となる。
このように、逆ウィシャート事前分布の情報性のレベルはリジッドに決まる。行列のなかのあるパラメータの情報性が他の全パラメータの情報性を決める。場面によっては柔軟性が足りないかもしれない。
諸君のために特に便利な事前分布を教えてあげよう。
- \(D\)を対角行列として\(IW(D, p+1)\)。こうすると、すべての相関の周辺分布が\((-1, 1)\)上で一様となり、分散の周辺分布は\(IG(1, \frac{d_{ii}}{2})\)となる。その最頻値は\( \frac{d_{ii}}{4} \)だから、お望みの最頻値が\( \frac{d_{ii}}{4} \)になるように\(D\)を決めれば良い。
- もっと一般的にいうと、\(IW(D, m)\)。こうすると、すべての相関の周辺分布が\((-1, 1)\)上でベータ分布 \(B(\frac{m-p+1}{2}, \frac{m-p+1}{2})\)となる[←これが知りたかったんだよーー！ありがとう兄貴！]。\(m \geq p\)ならばその平均は0, 分散は\(\frac{1}{m-p+2}\)である。
云々。