Triad sou.

ラプラス近似 (1)

積分ラプラス近似についてです。
だいぶ昔に [1] を引用して非線形混合効果モデルについて書いたことがあったのですが、積み残しを消化 (詳細を理解) するためにラプラス近似について調べてみました。

ラプラス近似を見てなんじゃこりゃと思ったのが [1] p.5 の (6) 式の下の式の変換で、ここをちゃんと理解できればいいなぁと思っていました。
そこで論文を検索してみたところ [2] が見つかり、この論文には一般化線形混合効果モデルのラプラス近似について書いてありました。

一般化線形混合効果モデルで用いられるラプラス近似で近似したいのは周辺尤度関数を求めるための積分で、被積分関数である全尤度関数が
\[
\exp\{h(\mathbf b)\}
\] と指数関数でくるまれた形になっている状態を考えています。
$\mathbf b\in \mathbb{R}^q$ は変量効果に対応する確率変数で、$h(\mathbf b)$ はスカラー関数です (ベクトル関数ではない)。

関数 $h(\mathbf b)$ を変量効果の推定量 $\hat{\mathbf b}$ の周りで 2 次まで Taylor 展開すると
\[
h(\mathbf b)\simeq h(\hat{\mathbf b})+h^{\prime}(\hat{\mathbf b})(\mathbf b-\hat{\mathbf b})+\frac{1}{2}(\mathbf b - \hat{\mathbf b})^t h^{\prime\prime}(\hat{\mathbf b})(\mathbf b - \hat{\mathbf b})
\] となります。
ここで、$\hat{\mathbf b}$ は関数 $h(\mathbf b)$ を最大化する解 (最尤推定値) なので、$h^{\prime}(\hat{\mathbf b})=0$ となることが分かります。
したがって、
\[
h(\mathbf b)\simeq h(\hat{\mathbf b})+\frac{1}{2}(\mathbf b - \hat{\mathbf b})^t h^{\prime\prime}(\hat{\mathbf b})(\mathbf b - \hat{\mathbf b}),
\] という近似式が得られます。

上で求めた Taylor 展開による近似を使って、周辺尤度を求める積分
\[
\int_{\mathbf b\in \mathbb{R}^q}\exp\{h(\mathbf b)\}\mathrm{d}\mathbf b,
\] の近似を求めたものがラプラス近似です。
ということで、Taylor 展開した近似式を代入して整理すると、
\[
\begin{align*}
\int_{\mathbf b\in \mathbb{R}^q}\exp\{h(\mathbf b)\}d\mathbf b & \simeq
\int_{\mathbf b\in \mathbb{R}^q} \exp\left\{ h(\hat{\mathbf b})+\frac{1}{2}(\mathbf b - \hat{\mathbf b})^t h^{\prime\prime}(\hat{\mathbf b})(\mathbf b - \hat{\mathbf b})\right\} \mathrm{d}\mathbf b \\ &=
\exp\{h(\hat{\mathbf b})\} \times \int_{\mathbf b\in \mathbb{R}^q} \exp\left\{ \frac{1}{2}(\mathbf b - \hat{\mathbf b})^t h^{\prime\prime}(\hat{\mathbf b})(\mathbf b - \hat{\mathbf b})\right\} \mathrm{d}\mathbf b
\end{align*}
\] と変形できます。
上式右辺第 2 項の積分は、うまく変形すると多変量正規分布, $N(\hat{\mathbf b}, -h^{\prime\prime}(\hat{\mathbf b})^{-1})$, の密度関数について全領域を積分したもの ($=1$) になるため、
\[
\begin{align*} (2\pi)^{q/2}&|-h^{\prime\prime}(\hat{\mathbf b})|^{-1/2} \times \\ &\quad~
\int_{\mathbf b\in \mathbb{R}^q}\frac{1}{(2\pi)^{q/2}|-h^{\prime\prime}(\hat{\mathbf b})|^{-1/2}}\exp\left\{ -\frac{1}{2}(\mathbf b - \hat{\mathbf b})^t (-h^{\prime\prime}(\hat{\mathbf b})^{-1})^{-1} (\mathbf b - \hat{\mathbf b})\right\} \mathrm{d}\mathbf b \\ &=
(2\pi)^{q/2}|-h^{\prime\prime}(\hat{\mathbf b})|^{-1/2} \times 1
\end{align*}
\] となります。
したがって、前述の第 2 項の積分を変形するとスケーリングの項だけが得られます。

以上から、周辺尤度関数のラプラス近似は
\[
\int_{\mathbf b\in \mathbb{R}^q}\exp\{h(\mathbf b)\}d\mathbf b\simeq (2\pi)^{q/2}|-h^{\prime\prime}(\hat{\mathbf b})|^{-1/2}\exp\{h(\hat{\mathbf b})\}
\] となります。
ちなみに [2] では、3 次以降まで Taylor 展開した場合についても評価しています。

CiteseerX 版の [1] も追ってみるとこの展開をしていますが、Taylor 展開で近似した 2 次の項にくっついてるはずの 1 / 2 が消えているような気がします、今度 print 版をちゃんと読んでみようと思います。
非線形混合効果モデルだと (でも?)、$h^{\prime\prime}(\hat{\mathbf b})$ は特殊な形をしているため、2 階偏微分しなくてはいけない部分が無視できる様な小さな値になり、計算量を 1 階偏微分を求める程度に減らすことが出来るようです。
こうして求められた周辺尤度の対数を元に最適化を行いパラメータの推定値を計算するようです。