前回に引き続き、時系列モデルについてのINLA(積分段階的ラプラス近似)の解説書。全14章だが、第1,2章だけはメモしながら読むことにした。いよいよ第2章。

Wang, X., Yue, Y.R., Faraway, J. (2023) Dynamic Time Series Models using R-INLA: An Applied Perspective. Chapter 2. A Review of INLA.

2.1 イントロダクション
[略]

2.2 ラプラス近似
　\(n\)個の観察\(\mathbf{y}^n = (y_1, \ldots, y_n)^\top\)と確率モデル\(p(\mathbf{y}|\Theta)\)に基づく、\(p\)次元パラメータ\(\Theta\)の対数尤度関数$$ \ell(\Theta; \mathbf{y}^n) = \sum_{t=1}^n \log p(y_t | \Theta)$$ があるとしよう。以下、\(\ell(\Theta)\)と略記する。
　なんだかしらんが\(w(\Theta), \pi(\Theta)\)という関数があるとしよう。ベイジアンの枠組みでいうと\(\pi(\Theta)\)は事前分布、\(w(\Theta) = u(\Theta) \pi(\Theta)\)で、\(u(\Theta)\)というのはなんらか関心ある関数である。
　で、次の比率について考えよう。$$ \frac{\int w(\Theta) \exp(\ell(\Theta)) d\Theta}{\int \pi(\Theta) \exp(\ell(\Theta)) d\Theta} $$ ベイジアンの枠組みでいうとこれは\(u(\Theta)\)の事後期待値\(E(u(\Theta)|\mathbf{y}^n)\)である。たとえば\(u(\Theta) = \Theta\)なら\(\Theta\)の事後平均だし、\(u(\Theta) = p(\mathbf{y}|\Theta)\) なら\(\mathbf{y}^n\)のもとでの事後予測である。
　簡略に、\(\Lambda(\Theta) = \ell(\Theta) + \log \pi(\Theta)\)と繰り込めば、$$ \frac{u(\Theta) \exp(\Lambda(\Theta)) d \Theta}{\int \exp(\Lambda(\Theta)) d\Theta} $$ とも書ける。
　ここまではいいっすね？

　Lindley(1980)はこの式の漸近的展開を考えた。事後モード\(\hat{\Theta}\)のまわりでテイラー展開すると$$ E(u(\Theta) | \mathbf{y}^n) \approx u(\hat{\Theta}) + \frac{1}{2} \left( \sum_{i,j = 1}^p u_{ij} h_{ij} + \sum_{i,j,k,l=1}^p \Lambda_{ijk} u_i h_{ij} h_{kl} \right) $$ $$ u_i \equiv \left. \frac{\partial u(\Theta)}{\partial \Theta_i} \right|_{\Theta^n = \hat{\Theta}} $$ $$ u_{ij} \equiv \left. \frac{\partial^2 u(\Theta)}{\partial \Theta_i \partial \Theta_j} \right|_{\Theta^n = \hat{\Theta}} $$ $$ \Lambda_{ijk} \equiv \left. \frac{\partial^3 \Lambda(\Theta)}{\partial \Theta_i \partial \Theta_j \partial \Theta_k} \right|_{\Theta^n = \hat{\Theta}} $$ ただし\(h_{ij}\)は、\(\Lambda\)の\(\hat{\Theta}\)におけるヘシアンの逆行列に-1を掛けた行列の要素である。
　[\(\Theta^n\)というのがなんだかよくわからない。おそらく、サイズ\(n\)のデータに基づく推定量という意味だと思うけど、添え字のみにあって式のなかにないのは変だと思う。誤植かな？ まあそれはともかく、Lindleyの近似ってこれまで読んだ奴でも何度も引用されてたけど、3次まで展開するようなややこしい話だったのか、とびっくり]

　これだとややこしすぎるので、Tierney & Kadane (1986)は2次導関数まででどうにかしようと考えた。彼らが提案した近似は、\(\Lambda^*(\Theta) = \log u(\Theta) + \Lambda(\Theta)\)とし、\(\Lambda\)の\(\hat{\Theta}\)におけるヘシアンの逆行列に-1を掛けた行列を\(\Sigma(\hat{\Theta})\)、\(\Lambda^*\)について同様に\(\Sigma^*(\hat{\Theta})\)として$$ E(u(\Theta)|\mathbf{y}^n) \approx \left( \frac{|\Sigma^*(\hat{\Theta})|}{|\Sigma(\hat{\Theta})|} \right) \exp(n(\Lambda^*(\hat{\Theta}) – \Lambda(\hat{\Theta}) )) $$ というもの。[Tierney & Kadane (1986) の2章の話だ]

2.2.1 単純化ラプラス近似
　関心あるのが期待値\(E(u(\Theta)|\mathbf{y}^n)\)じゃなくて、事後周辺分布\(\pi(\Theta_i|\mathbf{y}^n)\)である場合について考えよう。積分の比率の形で書くと以下となる。\(\Theta_{(-i)} = (\Theta_j | j=1,\ldots,n; j \neq i)\)として$$ \frac{\int \pi(\Theta) \exp(\ell(\Theta)) d\Theta_{(-i)}}{\int \pi(\Theta) \exp(\ell(\Theta)) d\Theta} $$

　Tierney & Kadane (1986)はこのラプラス近似も示している。ある固定された\(\Theta_i\)について、\(\pi(\Theta_i, \Theta_{(-i)}) \exp(\ell(\Theta_i, \Theta_{(-i}))\)のmaximaを \(\hat{\Theta}_{(-i)}(\Theta_i)\)とし、そこで評価したヘシアンの逆行列に-1を掛けたのを\(\Sigma^*(\Theta_i)\)と書く。先ほどと同様に、事後モード\(\hat{\Theta}\)で評価した\(\Theta\)の対数事後分布のヘシアンの逆行列に-1を掛けた奴を\(\Sigma(\hat{\Theta})\)と書く。$$ \pi(\Theta_i | \mathbf{y}^n) \approx \left( \frac{|\Sigma^*(\Theta_i)|}{2\pi n |\Sigma(\hat{\Theta}^n)|} \right)^{1/2} \frac{ \pi(\Theta_i, \hat{\Theta}_{(-i)}(\Theta_i) ) \exp(\ell(\Theta_i, \hat{\Theta}_{(-1)}(\hat{\Theta})}{\pi(\hat{\Theta}^n) \exp(\ell(\hat{\Theta}^n))} $$ [Tierney & Kadane (1986) の4章に出てきた式。ここでも\(\hat{\Theta}^n\)ってのが出てくるけど、上添字はあまり気にしなくてよさそうだな]

　ここからはRue, Martino, Chopin(2009)の話。
　事後周辺分布は同時分布\(\pi(\Theta_i, \mathbf{y}^n)\) と比例している。つまり、$$ \pi(\Theta_i | \mathbf{y}^n) \propto \frac{ \pi(\Theta, \mathbf{y}^n)}{\pi(\Theta_{(-i)} | \Theta_i, \mathbf{y}^n)} = \pi(\Theta_i, \mathbf{y}^n)$$ とも書ける。
　そこでRueらはこう考えた。まず、分母を正規分布で近似する。$$ \pi_G(\Theta_{(-i)} | \Theta_i, \mathbf{y}^n) = N(\hat{\Theta}_{(-i)}(\Theta_i) | \Theta_i, \mathbf{y}^n) $$ で、事後周辺分布はこうやって近似する。$$ \pi_{LA}(\Theta_i, \mathbf{y}^n) \propto \left. \frac{\pi(\Theta, \mathbf{y}^n)}{\pi_G(\Theta_{(-i)} | \Theta_i, \mathbf{y}^n)} \right|_{\Theta_{(-i)} = \hat{\Theta}_{(-i)}(\Theta_i)} $$　なお、この近似はTierney & Kadane (1986)によるラプラス近似と等価である。なぜなら、書き換えると $$ \pi_{LA}(\Theta_i | \mathbf{y}^n) \propto \pi(\Theta_i, \hat{\Theta}_{(-i)}(\Theta_i), \mathbf{y}^n) |\Sigma^*(\Theta_i)|^{1/2} $$ となるからである。
　この方法だと\(i\)ごとにヘシアンを評価しないといけないので、\(n\)が大きいときには計算が大変である。

　そこで、Rueらはより単純化した近似手法(単純化ラプラス近似; SLA)を考えた。
　完全事後分布\(\pi(\Theta|\mathbf{y}^n)\)は正規分布で近似される。そこで、分母の\(\pi_G(\Theta_{(-i)} | \Theta_i, \mathbf{y}^n)\) を、\(\Theta_{(-i)}\)の下での\(\Theta_{(-i)}\)の条件つき平均で評価する。$$ \pi_{SLA}(\Theta_i | \mathbf{y}^n) \propto \pi(\Theta_i, \hat{\Theta}_{(-i)}(\Theta_i), \mathbf{y}^n)$$

　[やれやれ… 個別の式についてはいまだ理解できないところがあるのだが、Tierney & KadaneからINLAに至る話の流れがようやくわかったよ…]

2.3 時系列のINLA構造
　オリジナルのINLAでは、潜在ベクトル\(\mathbf{x}^n\)はガウシアン・マルコフ確率場(GRMF)で定義されていた。GRMFは次の3レベルの階層で表現できる。上から、データモデル、GRMF事前分布、ハイパーパラメータ。$$ \mathbf{y}^n | \mathbf{x}_t, \Theta \sim \prod_t p(y_t | \mathbf{x}_t, \theta) $$ $$ \mathbf{x}^n | \Theta \sim N(\mathbf{0}, \Sigma(\theta)) $$ $$ \theta \sim \pi(\theta)$$ ただし、\(x_\ell\)と\(x_m\)は、\(\mathbf{x}^n\)からその二つを抜いた奴の下で条件付き独立である。2本目の共分散行列は\(\theta\)の関数なのでそういう書き方になっている。

　さて、時系列のローカル・レベル・モデルってあるじゃないですか。$$ y_t = x_t + v_t; \ \ v_t \sim N(0, \sigma^2_v) $$ $$ x_t = x_{t-1} + w_t; \ \ w_t \sim N(0, \sigma^2_w)$$ というモデル。これも上記の3階層で表せる。ハイパーパラメータは\(\theta = (\sigma^2_v, \sigma^2_w)^\top\)。潜在変数はガウシアンである。ハイパーパラメータはガウシアンでなくてよい。
　[上のモデルを拡張して] データの分布は指数分布族に属していればなんでもよいとしよう。

　[ああ、最初から時系列に絞って勉強すればよかった。GRMFと個別のモデルの関係がわかりやすい…]

2.3.1 INLAのステップ
　INLAはこの3階層をつかって、事後周辺分布を次のように定義する。$$ \pi(x_i | \mathbf{y}^n) = \int \pi(x_i | \theta, \mathbf{y}^n) \pi(\theta | \mathbf{y}^n) d \theta $$ $$ \pi(\theta_j | \mathbf{y}^n) = \int \pi(\theta \ \mathbf{y}^n) d\theta_{(-j)} $$ この右辺の\(\pi(x_i|\theta, \mathbf{y}^n), \pi(\theta|\mathbf{y}^n), \pi(\theta | \mathbf{y}^n) \)のところにその近似を入れる。

　INLAアプローチは次の3ステップからなる。

1. ハイパーパラメータの事後周辺分布 \(\pi(\theta | \mathbf{y}^n)\) をラプラス近似する。
   1. まず、完全条件付き分布 \( \pi(\mathbf{x}^n | \theta, \mathbf{y}^n)\)を、MVN \(\hat{\pi}_G (\mathbf{x} | \theta, \mathbf{y}^n)\)で近似する。
   2. 次に、ハイパーパラメータの事後密度を近似する。ここでラプラス近似が登場する。[2.2.1 の\(\pi_{LA}(\Theta_i|\mathbf{y}^n)\)の式であろう]
      　\(\pi(\mathbf{x}^n | \theta, \mathbf{y}^n)\)のモード[あれ？その近似であるMVNのモード？ 混乱してきた]を\(\hat{x}^n\)として、$$ \tilde{\pi} (\theta | \mathbf{y}^n) \propto \left. \frac{ \pi(\mathbf{x}^n, \theta, \mathbf{y}^n) }{ \tilde{\pi}_G(\mathbf{x}^n | \theta, \mathbf{y}^n)} \right|_{\mathbf{x}^n = \hat{\mathbf{x}}^n} $$ これはNewton-Raphsonアルゴリズムで求まる。
2. 潜在ガウシアン変数の条件付き分布 \(\pi(x_i | \theta, \mathbf{y}^n)\) を近似する。
   方法は3つある。(a)ガウシアン近似、(b)単純化ラプラス近似、(c)ラプラス近似。[説明は一切ない。そうか、ここはあまりこだわらなくていいのか]
3. 潜在ガウシアン変数の事後分布 \(\pi(x_i | \theta)\)を近似する。
   数値積分を使って、ハイパーパラメータを積分消去する。$$ \tilde{\pi}(x_i | \mathbf{y}^n) = \sum_j \tilde{\pi}(x_i | \theta_i, \mathbf{y}^n) \tilde{\pi}(\theta_j | \mathbf{y}) \Delta_j $$

2.4 INLAにおける予測
　1期先予測は次の式で得られる。$$ p(y_{n+1} | \mathbf{y}^n) = \int p(y_{n+1} | x_{n+1}, \theta, \mathbf{y}^n) \pi(x_{n+1}, \theta | \mathbf{y}^n) d x_{n+1} d\theta $$ [あれ？ 積分記号は2つ必要だよね。こういう略記もあるのかしらん？]
　\(x_t\)の下で\(Y_t\)は条件付き独立だから、こう書き換えられる。$$ p(y_{n+1} | \mathbf{y}^n) = \int p(y_{n+1} | x_{n+1}, \theta) \pi(x_{n+1}, \theta | \mathbf{y}^n) d x_{n+1} d\theta $$　
　右辺の二つ目、1期先の\(x\)の事後分布は下式で得られる。$$ \pi(x_{n+1}, \theta | \mathbf{y}^n) = \pi(x_{n+1} | \theta, \mathbf{y}^n) \pi(\theta | \mathbf{y}^n) $$ […中略…]

　INLAは\(x_t\)のGMRF構造を使ったベイジアン分析をやるのであって、予測分布をシミュレートするために\(x_t\)のマルコフ進化を活用するわけではないことに注意。

2.5 INLAにおける周辺尤度
　Rue et al.(2009)いわく、周辺尤度は$$ m(\mathbf{y}^n | M_k) \approx \int \frac{p(\mathbf{y}^n, \theta, \mathbf{x}^n | M_k)}{\tilde{\pi}_G (\mathbf{x}^n | \mathbf{y}, \theta, M_k)} d\theta $$ 右辺は\(\mathbf{x}^n = \hat{\mathbf{x}}^n(M_k)\)において評価する。

　モデル\(M_k\)の交差妥当化周辺尤度も近似できる。
　観察\(y_t\)のconditional predictive ordinate (CPO)は $$ p(y_t | \mathbf{y}_{(-t)}, M_k) = \int p(y_t | \mathbf{y}_{(-t)}, \theta, M_k) \pi(\theta | \mathbf{y}_{(-t)}, M_k) d\theta $$ 1個目は次の式で評価できる。$$ p(y_t | \mathbf{y}_{(-t)}, \theta, M_k) = 1 / \int \frac{\pi(x_t | \mathbf{y}^n, \theta, M_k)}{p(y_t | x_t, \theta, M_k)} d x_t$$ \(\pi(x_t | \mathbf{y}^n, \theta, M_k)\)はINLAで手に入る。上の積分は数値積分で評価する。つまり $$ \tilde{p}(y_t | \mathbf{y}_{(-t)}, M_k) = 1/ \sum_j \frac{\tilde{\pi}(\theta_j|\mathbf{y}^n, M_k)}{\tilde{p}(y_t | \mathbf{y}_{(-t)}, \theta_j, M_k)} \Delta_j$$ この推定量は\(\tilde{p}(y_t | \mathbf{y}_{(-t)}, \theta_j, M_k)\)の調和平均になっている。

2.6 R-INLAパッケージの基礎
　R-INLAパッケージのメインの関数は inla()である。formula()を受け取ってデータに当てはめる。lm(), glm(), arima()みたいなものである。

　モデルのあてはめは以下のステップからなる。

1. 線形予測子をformulaオブジェクトとして表現する。
   • 固定効果の指定はlm(), glm()と同様。係数はvague priorを持つ(自分で指定することもできる)。実は固定効果にはlinear型とclinear型があるんだけど、その話はあとで。
   • ランダム効果はf()で指定し、formulaに”+”でつなぐ。係数はふつう正規事前分布をもつ。

   たとえばこんな感じ: formula <- y + z1 + z_2 + f(x, model = "iid) + f(id, model = "iid)

2. inla()をコールする。引数として以下がある。
   • formula.
   • family. 正規、ポワソン、二項、などなどいろいろ指定できる。
   • data.
   • control. モデル選択基準とか、リンク関数とか、パラメータの事前分布の指定とか。

予測したかったら、predict()を呼ぶのではなくて、反応にNAを入れる。

　なお、R-INLAは周辺分布を扱う関数inla.tmarginal(), inla.zmarginal()を提供している。[...後略...]

-------
　まだ理解できていないところも多いんだけど、原典のRue et al.(2009)よりもはるかにわかりやすかった。こっちを先に読めばよかった... と思うんだけど、もしこれを先に読んでたら、やっぱり「わからん...」と落ち込んでいただろうな。
　ここから先は、実際にR-INLAを使い倒さないとわかんないだろうな。INLAの理屈についてはこのくらいにしておこう。