Tierney, L, Kadane, J.B. (1986) Accurate Approximations for Posterior Moments and Marginal Densities. Journal of the American Statistical Association, 82(393), 82-86.
都合により調べ物をしていて(後日のためにメモしておくとRue, Martino, & Chopin(2009), INLAについての論文)、途中で話についていけなくなったので、引用を遡って読んでみた奴。
こんなの読むなんて柄じゃないんですけどね。数学できないんですけどね。いったいどういう罰ゲームなのか。
1. イントロダクション
パラメータ空間上になめらかな正値関数\(g\)があるとします。対数尤度関数を\(LL\), 事前密度を\(\pi\)として、\(g(\theta)\)の事後平均はこうですよね。
[どこにも断りがないけど、サイズ\(n\)のデータ\(X^{(n)}\)のもとでパラメータ\(\theta\)をベイズ推論するなら、という話であろう。原文では対数尤度関数\( \mathscr{L} \)としているが、わかりにくいので勝手に書き換えた] $$ E_n[g] = E[g(\theta) | x^{(n)}] = \frac{\int g(\theta) \exp(LL(\theta)) \pi(\theta) d\theta}{\int \exp(LL(\theta)) \pi(\theta) d\theta} $$ [学力のない私のためにメモしておくと、分子の積分と分母の積分は意味が異なると思う。まず、事後密度関数 $$ p(\theta) = \frac{\exp(LL(\theta)) \pi(\theta)}{\int \exp(LL(\theta)) \pi(\theta) d\theta} $$ がある。で、\(\theta\)を変換する関数\(g(\cdot)\)がある。では、\(g(\theta)\)の期待値がどうなるかというと、$$ E[g(\theta)] = \int g(\theta)p(\theta) d\theta$$ だ、という話の流れである。そうですよね？]
上の式の分母の積分をどう近似するか。よくある方法は、ラプラスの方法をつかって、事後モードを平均、その点での対数事後密度の二次導関数の逆数かける-1 を分散にした正規分布を積分するという方法で、Lindley(1980), Misteller & Wallace (1964)がそうしている。
[なにをいっているのかというと、こういうことであろう。分母で積分しようとしている関数は、尤度関数と事前分布の積\(f(x) = \exp(LL(\theta)) \pi(\theta)\)である。これ自体は密度関数ではないけれど、これを積分1に正規化した密度関数を、\(N(\theta_0, \frac{1}{A})\)で近似したいなあという話である。\(\theta_0\)は密度関数のモードとする(どうやって見つけるのかは別にして)。\(A\)はどう決めるかというと、\(g(\theta) = \log f(x) \)として\(A = -g^{\prime\prime}(\theta_0)\)とする。なぜかというと、\(g(\theta)\)を\(\theta_0\)のまわりでテイラー展開して$$ g(\theta) \approx g(\theta_0) + \frac{g^{\prime\prime}(\theta_0)}{2} (\theta – \theta_0)^2$$ 両辺の指数をとると $$ f(\theta) \approx f(\theta_0) \exp \left(\frac{g^{\prime\prime}(\theta_0)}{2} (\theta – \theta_0)^2 \right) $$ いっぽう正規分布の密度関数は\( \propto \exp \left( -\frac{1}{2\sigma^2} (x – \mu)^2 \right) \)だから、これは平均\(\theta_0\), 分散\(-\frac{1}{g^{\prime\prime}(\theta_0)}\)の正規分布の密度関数のかたちになっているよねという話である。あってますでしょうか]
本論文では分子の積分もうまく近似する方法について考える。
直観的に説明するとこういう方法である。もし\(g\)がゼロから離れた事後中央値のまわりに集まっているなら、分母も分子も似たような形をしているので、両方に同じようにラプラスの方法を適用すれば、同じような誤差をとるので、分母と分子の比においては誤差はキャンセルされるだろう。
2. 事後平均と分散
\(L = \log \pi + LL/n, L^* = \log g + \log \pi + LL/n\)として、上の式を書き換える。$$ E_n[g] = \frac{\int \exp(nL^*) d\theta}{\int \exp(nL) d\theta} $$ [恥ずかしながら、もうここからわからない。$$ E_n[g] = \frac{\int g(\theta) \exp(LL(\theta)) \pi(\theta) d\theta}{\int \exp(LL(\theta)) \pi(\theta) d\theta} $$ $$ = \frac{\int \exp \{ \log(g(\theta)) + LL(\theta) + \log(\pi(\theta)) \} d\theta}{\int \exp \{ LL(\theta) + \log(\pi(\theta)) \} d\theta} $$ となるはずではないか。ひょっとすると、\(L = (\log \pi + LL)/n, L^* = (\log g + \log \pi + LL)/n\)の誤植だろうか？]
分母の近似を考えよう。\(L\)のモードを\(\hat{\theta}\)とし、\(\sigma^2 = -\frac{1}{L^{\prime\prime}(\hat{\theta}) } \)として、$$ \int \exp(nL(\theta)) d\theta \approx \int \exp \left( nL(\hat{\theta}) – \frac{n(\theta – \hat{\theta})^2}{2\sigma^2} \right) d\theta = \sqrt{2\pi} \sigma n^{-1/2} \exp(n L(\hat\theta)) $$ [あああ！速すぎる！ まず\(L(\theta)\)をテイラー展開して $$ L(\theta) \approx L(\hat{\theta}) + \frac{L^{\prime\prime}(\hat{\theta})}{2} (\theta – \hat{\theta})^2 = L(\hat{\theta}) – \frac{(\theta – \hat{\theta})^2}{2\sigma^2} $$ 元の式の\(L(\theta)\)をこれで近似して $$ \int \exp \left( nL(\hat{\theta}) – \frac{n(\theta – \hat{\theta})^2}{2\sigma^2} \right) d\theta $$ $$ = \int \exp \left( nL(\hat{\theta}) \right) \exp \left(- \frac{n(\theta – \hat{\theta})^2}{2\sigma^2} \right) d\theta$$ あ！ひとつめの\(\exp\)は\(\theta\)の関数になってない！ ふたつめの\(\exp\)の積分は、私の手には負えないけれど、よくみると分散\(\sigma^2/n\)の正規分布の密度関数と同じ形になっているから、積分は分散\(\sigma^2/n\)の正規分布の正規化項\(\sqrt{2 \pi (\sigma^2/n)}\)である。よって$$ = \exp(nL(\hat{\theta})) \sqrt{2 \pi} \sigma n^{-1/2} $$ となるわけだ]
同様に分子についても考える。\(L^*\)のモードを\(\hat{\theta}^*\)とし、\( \sigma^{*2} = -\frac{1}{L^{*\prime\prime}}(\hat{\theta}^*) \)として、$$ \int \exp(nL^*) d\theta \approx \sqrt{2\pi} \sigma^* n^{-1/2} \exp(nL^* (\hat{\theta}^*)) $$
比をとって、$$ \hat{E}_n[g] = \frac{\sigma^*}{\sigma} \exp\left( n(L^*(\hat{\theta}^*) – L(\hat{\theta})) \right) $$ 以下ではこれをラプラス近似と呼ぶ。
多パラメータにも拡張できる。
\(\hat{\theta}^*\)における\(L^*\)のヘッセ行列の逆行列に-1を掛けたのを\(S^*\)とし[原文では\(\Sigma\)に縦棒が刺さっている記号]、\(\hat{\theta}\)における\(L\)のヘッセ行列の逆行列に-1を掛けたのを\(S\)として、$$ \hat{E}_n[g] = \left( \frac{ \mathrm{det} S^* }{ \mathrm{det} S } \right)^{1/2} \exp \left( n(L^*(\hat{\theta}^*) – L(\hat{\theta}) \right)$$
[速すぎてむかつくけど、これはまあ信じましょう]
このラプラス近似の漸近的正確性は… [中略]
というわけで、\(g\)の事後平均についての近似が得られた。事後分散については、$$ \hat{V}_n[g] = \hat{E}_n[g^2] – \hat{E}_n[g]^2 $$ と近似すれば良い。その正確性は…[中略]
実務的な観点からは、\(n\)が小さいとき、またすごく大きいとき、事後分散の近似には要注意である…[中略]
最後に、\(g\)は正値だという仮定について。これは、分子と分母で形が似ているために必要な仮定である… [中略]
3. 予測分布
[いまちょっと関心ないのでメモ省略]
4. 周辺事後密度
今度は多パラメータの場合の、個々のパラメータの周辺事後密度を近似する方法について考えよう。
[mathjaxではギリシア文字を太字にしても見た目が変わらないので、以下、\(\theta\)の太字の代わりに\(\Theta\)を使う]
\(\Theta = (\theta_1, \ldots, \theta_p) \)について、その事後モード、つまり尤度を最大化する値を\(\hat{\Theta}\)とする。\(LL + \log \pi\)の\(\hat{\Theta}\)におけるヘッセ行列の逆行列に-1を掛けた\(p \times p\)行列を\(S\)とする。
\(\Theta\)を最初のパラメータ\(\theta_1\)とそれ以外のパラメータ\(\Theta_2\)にわける。\(\theta_1\)を所与として、関数\(h(\cdot) = \pi(\theta, \cdot) \exp(LL(\theta_1, \cdot))\)を最大化するベクトルを\(\hat{\Theta}^*_2 = \hat{\Theta}^*_2(\theta_1)\)とする。\(\log h(\cdot) / n\)のヘッセ行列の逆行列に-1を掛けた\((p-1) \times (p-1)\)行列を\(S^*\)とする。
ラプラス近似を使って$$ \pi_1(\theta | x{(n)}) = \pi_{n,1}(\theta_1) = \frac{ \int \pi(\theta_1, \Theta_2) \exp(LL(\theta_1, \Theta_2)) d\theta }{ \int \pi(\Theta) \exp(LL(\Theta)) d\Theta} $$ だから、周辺密度の近似として以下が得られる。$$ \hat{\pi}_{n,1}(\theta_1) = \left( \frac{\mathrm{det} S^*(\theta_1)}{2\pi n \ \mathrm{det} S} \right)^{1/2} \frac{\pi(\theta_1, \hat{\Theta}^*_2) \exp(LL(\theta_1, \Theta_2^*)) }{ \pi(\hat{\Theta}) \exp(LL(\hat{\Theta})) } $$
[以下、この近似の正確性についての説明が半ページ続くが、力尽きて読み飛ばした]
5. 応用例
[略]
6. 結論
[略]
—————-
途中で力尽きて投げやりになってしまったが、どういう話なのかということはなんとなーく分かったので、まあ、これでよしとしよう…
2023/09/17: 1節と、2節の単変量のパートのメモを大幅に追記した。