流れ流れてデータ解析とかで生計を立てていますが、一日に三回くらい、なぜ私はこんなことをしているのか… もっとましな人生があったのではないか… という疑念に囚われる。もっとましな人生ってなに？ よくわかんないけど。
さらに、Stanのコンパイルでイライラするたび、なぜ私はこんな面倒なことをしているのか… これはMCMCじゃなくてなにか別の方法でも解けるのではないか… という疑念に囚われる。別の方法ってなに？ 知らんけど。INLAとか？
そんなこんなで、数年前から何度かINLAを実戦投入しようとし、その度に挫折している。検索すると2017年、R-INLAについての解説を読んでいるようだが、メモを読み返すと内容をあまり理解できていないことが丸分かりである。切ないのう。
Rue, H., Martino, s., Chopin, N. (2009) Approximate Bayesian Inference for Latent Gaussian Models by Using Integrated Nested Laplace Approximations. Journal of Royal Statistical Society, B. 71(2), 319-392.
そういうわけで、先日も仕事の都合で発作的にINLAの勉強をはじめた(仕事のほうは時間切れでそれどころではなくなり、結局はStanで切り抜けた)。その際のメモ。長くなるので2パートくらいにわける。
Gaussianは原則としてガウシアンと訳しているが、GMRFはガウス・マルコフ確率場と訳している。一貫性がないね。
1. イントロダクション
1.1 本論文の目標
本論文は潜在ガウシアンモデルについての近似ベイジアン推論を行う方法について論じる。
潜在ガウシアンモデルというのは、構造化加法予測子を持つベイジアン加法回帰モデルの下位クラスであって、こういうモデルである。反応変数\( y_i \)は指数分布族に属し、その平均\( \mu_i \)は、構造化加法予測子 \( \eta_i \)とリンク関数 \(g(\cdot)\)でつながっている。つまり \( g(\mu_i) = \eta_i\)である。でもって、$$ \eta_i = \alpha + \sum_{j=1}^{n_f} f^{(j)}(u_{ji}) + \sum_{k=1}^{n_\beta} \beta_k z_{ki} + \varepsilon_i$$ \( \{f^{(j)}(\cdot)\} \)はどんな関数でもよい。\(\alpha, \{f^{(j)}(\cdot)\}, \{\beta_k\}, \{\varepsilon_i\}\)にはガウシアン事前分布が与えられている。ハイパーパラメータのベクトル\(\theta\)はガウシアンでなくてよい。
本論文の目的は次の2つである。(1)\(\mathbf{x}\)や\(\theta\)の要素の事後周辺分布の決定論的な近似を正確に速く得る方法を提供する。(2)それらの周辺分布の使い方を示す。
1.2 潜在ガウシアンモデルの用途
[メモは省略するけど、いろんな回帰モデルも、動的モデルも、空間・時空間モデルもカバーできる]
1.3 潜在ガウシアンモデル: 表記と基本的性質
以下ではなにかの条件付き密度を\( \pi(\cdot|\cdot) \)と書きます。
すべての潜在ガウシアン変数、つまり\( \{ \eta_i \}, \alpha, \{f^{(j)}(\cdot)\}, \{\beta_k\} \)のベクトルを\(\mathbf{x}\)とし (話の都合上、\(\{\varepsilon_i\}\)の代わりに\(\{\eta_i\}\)をいれる)、その長さを\(n\)とする。
密度\(\pi(\mathbf{x}|\theta_1)\)はガウシアンで、平均0で、精度行列は\(\mathbf{Q}(\theta_1)\)とする。
\(n_d\)個の観察変数\(\mathbf{y} = \{y_i; i \in \mathcal{I}\}\)が、\(\pi(\mathbf{y}|\mathbf{x}, \theta_2)\)にiidに従うとする。
\( \theta = (\theta^\top_1, \theta^\top_2)^\top \)とし、長さを\(m\)とする。
では、事後分布はどうなるでしょうか。\(\mathbf{Q}(\theta)\)は正則として、$$ \pi(\mathbf{x}, \theta | \mathbf{y}) $$ $$\propto \pi(\theta) \pi(\mathbf{x} | \theta) \prod_{i \in \mathcal{I}} \pi(y_i | x_i, \theta) $$ [ここの\(x_i\)という表記がよくわかんないんですよね。\(i\)は観察個体のインデクスで\(1\)から\(n_d\)まで動くのに、\(\mathbf{x}\)は長さ\(n\)の潜在変数ベクトルじゃないですか。よくわかんないけど、気持ちとしては、\(x_i\)は\(\mathbf{x}\)のなかで\(i\)番目の個体に関係する要素からなるベクトル、つまり\( \eta_i, \alpha, \{f^{(j)}(\cdot)\}, \{\beta_k\}\) からなるベクトルなのではないかと思うのだけれど、そんなら太字で書くべきだと思うんですよね…] $$ \propto \pi(\theta) |\mathbf{Q}(\theta)|^{1/2} \exp \left( -\frac{1}{2} \mathbf{x}^\top \mathbf{Q}(\theta) \mathbf{x} + \sum_{i \in \mathcal{I}} \log(\pi(y_i|x_i, \theta)) \right) $$ [なんでかというと、話を簡単にするために\(\Sigma = \mathbf{Q}^{-1}(\theta)\)として、\(\pi(x|\theta)\)は\(MVN(\mathbf{0}, \Sigma)\)の密度関数、つまり$$ \frac{\exp \left(-\frac{1}{2}\mathbf{x}^\top \Sigma^{-1} \mathbf{x} \right)}{\sqrt{(2\pi)^n |\Sigma|}} $$だからね。分母の\(\sqrt{(2\pi)^n}\)をネグっているわけ]
もし線形制約をかけるなら、\(k \times n\), ランク\(k\)の行列\(\mathbf{A}\)を使って\( \mathbf{Ax} = \mathbf{e} \)と書く。
やりたいのは、事後周辺分布 \(\pi(x_i|\mathbf{y}), \pi(\theta|\mathbf{y}), \pi(\theta_j | \mathbf{y}) \)の近似である。
潜在ガウシアンモデルというのはたいてい次の2つの性質を持っている。本論文でもこれらを仮定する。
- 潜在場\(\mathbf{x}\)は条件付き独立性を持っている。したがってそれはガウス・マルコフ確率場(GMRF)で、その精度行列はスパースである。
- ハイパーパラメータの数\(m\)は小さい。せいぜい6個くらい。
1.4 推論: MCMCによるアプローチ
[メモ省略。計算時間的に無理なことが多いよねって話]
1.5 推論: 決定論的近似
[MC誤差と近似バイアスの比較。パス]
1.6 機械学習における近似手法
[変分ベイズ、expectation-propagationの話。パス。そんなことよりさあ、INLAの話が聞きたいのよ～]
1.7 推論: 新しいアプローチ
いま関心ある事後周辺分布は$$ \pi(x_i|\mathbf{y}) = \int \pi(x_i|\theta, \mathbf{y}) \pi(\theta|\mathbf{y}) d \theta$$ $$ \pi(\theta_j|\mathbf{y}) = \int \pi(\theta|\mathbf{y}) d \theta_{-j}$$ ですね。
我々の新しいアプローチのポイントは、この形式を使ってnested近似[???]をつくるという点である。密度\(\pi(\cdot|\cdot)\)の近似を\(\tilde{\pi}(\cdot|\cdot)\)として、$$ \tilde{\pi}(x_i|\mathbf{y}) = \int \tilde{\pi}(x_i|\theta, \mathbf{y}) \tilde{\pi}(\theta|\mathbf{y}) d\theta$$ $$ \tilde{\pi}(\theta_i|\mathbf{y}) = \int \tilde{\pi} (\theta|\mathbf{y}) d\theta_{-j} $$ つまり近似を数値積分するわけです(だから\(\theta\)の次元が大きいと困る)。
[ここから話の流れを見失い苦戦した。おかげでほぼ逐語訳。まずは \(\pi(\theta | \mathbf{y})\)をどう推定するかという話から]
我々のアプローチは、\(\theta\)の周辺事後分布についての次の近似に基づく。$$ \tilde{\pi}(\theta | \mathbf{y}) \propto \left. \frac{\pi(\mathbf{x}, \theta, \mathbf{y}) }{\tilde{\pi}_G (\mathbf{x} | \theta, \mathbf{y})} \right|_{\mathbf{x} = \mathbf{x}^*(\theta)} $$ \(\tilde{\pi}_G\)というのはガウシアン近似、\(\mathbf{x}^*(\theta)\)とは\(\theta\)のもとでの\(\pi(\mathbf{x}|\theta, \mathbf{y})\)のモードである。
[なにが困るといって、この式が！わからない！！
$$ \tilde{\pi}(\theta | \mathbf{y}) \propto \frac{\pi(\mathbf{x}, \theta, \mathbf{y}) }{\tilde{\pi}_G (\mathbf{x} | \theta, \mathbf{y})}$$ っていうだけならなんとなく分かるんですよ。そこに\(\mathbf{x} = \mathbf{x}^*(\theta)\)を代入するという意味がよくわからない…]
これはTierney & Kadane (1986) が提案した周辺事後分布のラプラス近似と等価である。彼らによれば、近似誤差は相対的で[?]、再基準化後のオーダーは\(O(n_d^{-3/2})\)である。しかし、\(n\)は固定されておらず\(n_d\)に依存するから、ラプラス展開においてふつう援用される標準的漸近性の仮定はここでは検証できない。4章を参照。[わからん。なにいってんだろうか]
\(\tilde{\pi}(\theta | \mathbf{y})\)それ自体は正規性からはほど遠い。従って、\(\pi(\theta | \mathbf{y})\)へのガウシアン近似に基づく近似は、我々の目的のためには正確性が十分でない。このことは、上の式のモードとにおいて評価される「\(\mathbf{x}^*\)のまわりでの等価なガウシアン観察」に基づく近似にもあてはまる。[この段落も全くわからん。なにいってんだろうか…]
我々のアプローチの重要な側面は、\(\tilde{\pi}(\theta | \mathbf{y})\)と\(\tilde{\pi}(x_i|\mathbf{y})\)を「ノンパラメトリックな」やりかたで検討・操作する点である[???]。
Rue & Martino (2007)は上の式に基づき、さまざまな潜在ガウシアンモデルに関して\(\theta\)の事後周辺分布を近似し、\(\tilde{\pi}(\theta|\mathbf{y})\)は正確だと結論した。潜在場の事後周辺分布に関しては、彼らは\(\tilde{\pi}_G(\mathbf{x}|\theta, \mathbf{y})\)から出発し、そこから導出される正規周辺分布を使って\(x_i | \theta, \mathbf{y}\)の密度を近似するという方法を提案している。すなわち、ガウス近似の平均と周辺分散を\(\mu(\theta), \sigma^2(\theta)\)として$$ \tilde{\pi} (x_i | \theta, \mathbf{y}) = N\{x_i; \mu_i(\theta), \sigma^2_i(\theta)\} $$ この近似の誤差は大きい。で、この近似を数値積分して、潜在場に関する関心ある周辺分布の近似を得る。$$ \hat{\pi}(x_i | \mathbf{y}) = \sum_k \tilde{\pi}(x_i | \theta_k, \mathbf{y}) \tilde{\pi}(\theta_k | \mathbf{y}) \Delta_k $$ この近似は正確である。
このアプローチの問題は、\(x_i\)の近似が正確でないことに気づくのが難しく、近似を正確にすることもできないという点である。近似の誤差をコントロールして積分点を適応的・自動的に決めるということもできない。
本論文では、これらの問題点を解決し、潜在ガウシアンモデルの完全に自動的な近似推論の方法を提案する。名付けて積分段階的ラプラス近似(INLA)。
1.8 本論文の構成
[パス]
2. 前提
以下では\(\mathbf{x}\)から\(i\)番目の要素を抜いたのを\(\mathbf{x}_{-i}\)と書く。\(\Gamma(a, b)\)(平均は\(a/b\))の点\(\tau\)における密度を\(\Gamma(\tau; a, b)\)と書く。
2.1 ガウス・マルコフ確率場
マルコフ性を持つガウシアン確率変数\(\mathbf{x} = (x_i, \ldots, x_n)\)をGMRFという。マルコフ性とは、なんらかの\(i \neq j\)において\(x_i\)と\(x_j\)が\(\mathbf{x}_{-ij}\)のもとで独立だということである
[ここでちょっと混乱したんだけど、Wikipediaによれば、マルコフ性というのは3種類あるんだってさ。(1)ペアワイズマルコフ性。隣接しない任意の\(i \neq j\)について、他のすべての変数を与えられたときに条件付き独立になる。(2)局所マルコフ性。変数\(i\)に直接つながっている変数が条件付けられたとき、\(i\)が任意の変数\(j\)と条件付き独立になる。(3)大域マルコフ性。変数\(i\)から変数\(j\)に至るすべての経路が変数集合\(S\)を通るとき、\(S\)で条件付ければ、\(i\)と\(j\)が条件付き独立になる。(3), (2), (1)の順に強い仮定なんだけど、なんだかよくわかんないけど同時分布が狭義正測度であれば、(1)(2)(3)は同値なのだそうだ。ここでは(1)のマルコフ性の話をしているんでしょうね]
マルコフ特性は、精度行列(=共分散行列の逆行列)\(Q\)においては簡単に表現できる。\(x_i\)と\(x_j\)が\(\mathbf{x}_{-ij}\)のもとで条件付き独立であるとき、その時に限り、\(Q_{ij} = 0\)である。いま、\(x\)の条件付き独立性を無向グラフ\(G\)で表すと[原文では筆記体]、\(\mathbf{x}\)は\(G\)に関してGRMFであるという。
[あれ？ もし\(G\)が完全無向グラフだったらGMRFにはならなくない？ だってどの\(x_i, x_j\)も、\(\mathbf{x}_{-ij}\)のもとで条件付き独立にならないから]
\(x\)の密度はこうなる。平均を\(\mu\)として、$$ \pi(\mathbf{x}) = (2\pi)^{-n/2} |\mathbf{Q}|^{1/2} \exp \left( -\frac{1}{2} (\mathbf{x}-\mu)^\top \mathbf{Q} (\mathbf{x}-\mu) \right) $$ [びびっちゃだめ、深呼吸。単にMVNの密度関数を書いているだけだ]
\(\mathbf{Q}\)はたいていスパースである。ってことはコレスキー分解\(\mathbf{Q} = \mathbf{L} \mathbf{L}^\top\)は速い。\(\mathbf{L}\)もまたスパースである。ノードの順番を並び替えて、\(\mathbf{L}\)のなかの非ゼロの要素をさらに減らすことができる。[…中略…]
\(\mathbf{Q}\)が含まれる式をコレスキー分解をつかって解くこともできる。たとえば、\(\mathbf{Q}\mathbf{x} = \mathbf{b}\)を解くとき、まず\(\mathbf{L} \mathbf{v} = \mathbf{b}\)を解いて、次に\(\mathbf{L}^\top \mathbf{x} = \mathbf{v}\)を解く[1本目に2本目を代入すると\(\mathbf{LL}^\top = \mathbf{b}\)だから]。
もし\(\mathbf{z} \sim N(\mathbf{0}, \mathbf{I})\)なら、\(\mathbf{L}^\top \mathbf{x} = \mathbf{z}\)の解は精度行列\(\mathbf{Q}\)を持つ。これはGRMFからランダムなサンプルを得る一般的な方法である。[恥ずかしながら、ここ、わからない…]
周辺分散も効率的に計算できる。上述のように、\(\mathbf{L}^\top \mathbf{x} = \mathbf{z}, \mathbf{z} \sim N(\mathbf{0}, \mathbf{I})\) の解\(x\)は精度行列\(Q\)を持つ。で、$$ L_{ii} x_i = z_i – \sum_{k=1+i}^n L_{ki} x_{k} $$ である。両側に\(x_j, j \geq i\)を掛けて期待値をとると、$$ \Sigma_{ij} = \frac{\delta_{ij}^2}{L_{ii}} – \frac{1}{L_{ii}} \sum^n_{k=i+1} L_{ki} \Sigma_{kj}$$ ただし\(\Sigma = \mathbf{Q}^{-1}\)は共分散行列、\(\delta_{ij}\)は\(i=j\)のときだけ1になるダミー変数である。[正直にいおう、ここも写経である。まあ\(x\)の共分散行列が求まりますってことね]
GMRFに線形制約がかかっている場合は以下の戦略が使える。[…中略…]
2.2 ガウシアン近似
[きっと肝になる箇所なので、ほぼ逐語訳]
我々のアプローチは、次の形式の密度に対するガウシアン近似に基づいている。$$ \pi(\mathbf{x}) \propto \exp \left( -\frac{1}{2} \mathbf{x}^\top \mathbf{Q} \mathbf{x} + \sum_i g_i (x_i) \right) $$ \(g_i(x_i)\)は我々のセッティングでは\(\log(\pi(y_i | x_i, \theta))\)である。
[さあ深呼吸。なにをゆうとるのかというと… きっと、事後密度は事前密度と尤度の積に比例するってことでしょうね。事前分布は正規だから、事前密度は精度行列を\(\mathbf{Q}\)として\( \propto \exp \left( -\frac{1}{2} (\mathbf{x} – \mu)^\top \mathbf{Q} (\mathbf{x} – \mu) \right) \)ですよね。あ、\(\mu = 0\)と置いているってことか！ 観察事例\(i\)の個体尤度は\(\pi(y_i|x_i, \theta)\)で、対数尤度はその対数の和\(\sum_i g_i(x_i)\)だ。その指数を事前密度に掛けている、つまり指数記号のなかで足しているわけだ。なるほどね！ 完全に理解したよ！]
ガウシアン近似\(\tilde{\pi}_G(\mathbf{x})\)は、matching the modal configration and the curvature at the modeによって得られる。モードはNewton-Raphson法で反復的に計算する。この方法はスコアリング・アルゴリズムとか、その変異体であるフィッシャー・スコアリング・アルゴリズムとして知られている。
いま、最初のguessを\(\mu^{(0)}\)としよう[この\(\mu^{(0)}\)は太字、つまりベクトル]。\(g_i(x_i)\)を\(\mu_i^{(0)}\)のまわりで2次まで展開する。$$ g_i(x_i) \approx g_i(\mu_i^{(1)}) + b_i x_i – \frac{1}{2} c_i x^2_i $$ \(\{b_i\}, \{c_i\}\)は\(\mu_{(0)}\)に依存する。
[ここがよくわかんない。テイラー展開のこと？ だったら第2項は\( g’_i(\mu_i^{(0)}) (x_i – \mu_i^{(0)}) \), 第3項は\(\frac{1}{2} g^{\prime\prime}_i(\mu_i^{(0)}) (x_i – \mu_i^{(0)})^2 \)になるはずだよね。あっ、もしかすると、第2項と第3項をバラして整理すれば \(b_i x_i – \frac{1}{2} c_i x^2_i\)という形に変形できます、ってこと？ きっとそうだね！ 完全に理解した！]
こうして、精度行列\(\mathbf{Q}+\mathrm{diag}(\mathbf{c})\)が得られる。\(\{\mathbf{Q} + \mathrm{diag}(\mathbf{c})\} \mu^{(1)} = \mathbf{b} \)の解としてモードが得られる。
[正直にいおう。ここがわからん。なにいってんだかさっぱりわからん。事前密度の精度行列が\(\mathbf{Q}\)だとわかっているとして、事後密度\(\pi(\mathbf{x})\)のモードを\(\mu^{(0)}\)に仮置きしたら、事後密度の精度行列はちょっと大きくなる、というのはわかる。それが\(\mathbf{Q} + \mathrm{diag}(\mathbf{c})\) だ、ってこと？ そこがわからん。で、もし事後密度の精度行列を\(\mathbf{Q} + \mathrm{diag}(\mathbf{c})\)に仮置きしたら、こんどはモードは\(\{\mathbf{Q} + \mathrm{diag}(\mathbf{c})\} \mu^{(1)} = \mathbf{b} \)を満たすような\(\mu^{(1)}\)になる、ってこと？ その理屈もわからん。うぐぐぐぐ。なにもわからない…]
これを収束するまで繰り返す。収束した平均を\(\mathbf{x}^*\), 精度行列を\(\mathbf{Q}^* = \mathbf{Q} + \mathrm{diag} (\mathbf{c}^*)\)としよう。 もし線形制約があったら[…中略…]。
上の\(\pi(\mathbf{x})\)の式に出てくるnon-quadraticな項は [\(\sum_i g_i(x_i)\)のことかな]、\(x_i\)の関数であって、たとえば\(x_i\)と\(x_j\)の関数ではない。だから、ガウシアン近似における精度行列は、\(\mathbf{Q} + \mathrm{diag}(\mathbf{c})\)の形式となる。GRMFのマルコフ性が保存されるので、計算しやすい。
[うん、雰囲気としてわかる。上の式を事前分布から事後分布への更新としてみると、個体の尤度を掛けただけなので、個体ごとの事後分布の分散は小さくなるけど共分散は小さくならないよ、つまり精度行列でいうと対角成分だけが更新されるよ、ってことね？ 俺は雰囲気だけで読み進めている…]
モードとモードにおけるcurvatureをマッチングすることによって得られる分布に関して上の式をガウシアン近似させる方法の改善については、以下をみよ。Rue(2001 JRSS), Rue & Held(2005 書籍), Kuss & Rasmussen(2005 J.Mach.Learn.Res.)。ここではこの問題については追及しない。
——–
いよいよINLAの登場だが、疲れたのでこの辺までにしておこう。続きはまた後日。
[2023/09/15追記: 2章のメモを書き直した。自分がどこを理解していないかがよくわかった。進歩だ！]
[2023/09/22追記: 2章のメモをさらに追記。]