ここしばらく、時系列モデリングの解説書を頼りに、INLA(積分段階的ラプラス近似)の理屈を理解しようと試みていた。目当ての2章は読み終えたので、この本はこれでおしまいにするつもりにするつもりだったのだが、なんとなく3章も目を通してしまった。

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

せっかくなので、各節のごくごく概要についてメモしておく。ちなみに本章は長いのだが、チャートが多く、実質的な内容はみかけほど多くない。

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

3.1 ランダム・ウォーク＋ノイズモデル
またの名をローカルレベルモデル。$$ 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.2.1 R-INLA モデルformula
　こんな感じになる。
　 formula.rw1 <- y.rw1 ~ f(id.w, model = "rw1", constr = FALSE) -1
　R-INLAでは、デフォルトでは項(ここでは\(x_t\))の合計が0に制約されるので、f()のconstrでこれを解除している。

3.2.2 モデル実行
　こんな感じ。
　model.rw1 <- inla(formula.rw1, family = "gaussian", data = rw1.dat, control.predictor = list(compute = TRUE))
　出力の見方は... [略]

3.2.3 ハイパーパラメータの事前分布の指定
　R-INLAでは、事前分布はパラメータの「内的表現」を通じて割り当てられる。たとえば、精度は対数尺度で表現される。観察のノイズ分散\(\sigma^2_v\)の事前分布は\(\sigma^2_v\)じゃなくて\(\theta_v = \log(1/\sigma^2_v\)に対して付与される。状態のノイズの分散も同様である。
　\(\theta_v, \theta_w\)のデフォルトの事前分布は log gamma(shape=1, inverse scale = 5e-05) となる。

3.2.4 ハイパーパラメータの事後分布
　上記の理由で、たとえば\(\sigma^2_v\)の事後分布が欲しかったら変換が必要になる。inla.tmarginal()を使う。こんな感じ。
　sigma2.v.dist <- inla.tmarginal(fun = function(x) exp(-x), marginal = model.rw1$internal.marginals.hyperpar$'Log precision for the Gaussian observations')
　[ひえー。リスト要素の指定がめんどくさそう...]

　事後分布の期待値を取りたかったらinla.emarginal(), 分位点を取りたかったらinla.qmarginal(), スプライン平滑化したかったらinla.smarginal(), サマリーを取りたかったらinla.zmarginal()を使え。[コード例を見ていると、変換しつつ期待値とかを取れるようだ]

3.2.5 潜在状態と反応に対してフィットした値
　\(\pi(x_t | \mathbf{y}^n)\)のサマリーはモデルオブジェクトにもう入っている。ついでにKLダイバージェンスも入っている。これは単純化ラプラス近似と標準正規分布の間の距離を表す。
　\(y_t\)の予測分布の平均ももう入っている。

3.2.6 DLMにおけるフィルタリングとスムージング
　[inla()は全時点一気に推定しちゃう、つまりスムージングしちゃうので、もしフィルタリングしたかったら時点ごとにinla()を繰り返し呼べ、という話。略。そういうニーズってあるのかしらん...]

3.3 AR(1)レベル＋ノイズモデル
　こういうモデル。$$ y_t = \alpha + x_t + v_t; \ \ v_t \sim N(0, \sigma^2_v)$$ $$ x_t = \phi x_{t-1} + w_t; \ \ w_t \sim N(0, \sigma^2_w), \ \ t = 2, \ldots, n$$ $$ x_1 \sim N(0, \sigma^2_w / (1-\phi^2))$$ この\(\phi\)のことをR-INLAでは\(\rho\)と呼ぶので、そこんとこ注意するように。

　[ハイパーパラメータの事前分布についての説明など。略]

　formulaでf(..., model = "ar1", ...)と指定してinla()を呼ぶ。うっかり\(\alpha\)を落とさないようにね。

3.4 高次ARラグを持つDLM
　AR(p)レベル＋ノイズモデルの場合...
　[formulaでf(..., model = "ar", order = 2, ...)と指定するそうである。めんどくさいので読み飛ばす]

3.5 ドリフトつきランダムウォーク＋ノイズモデル
　ドリフトをいれて、
　$$ y_t = x_t + v_t; \ \ v_t \sim N(0, \sigma^2_v)$$ $$ x_t = \delta + x_{t-1} + w_t; \ \ w_t \sim N(0, \sigma^2_w)$$ このモデルはこういう風に書き換えて考えること。$$ y_t = \delta t + z_t + v_t; \ \ v_t \sim N(0, \sigma^2_v)$$ $$ z_t = z_{t-1} + w_t; \ \ w_t \sim N(0, \sigma^2_w)$$ formulaはこんな感じid.deltaはid.wと同じく時点番号。
　formula.rw.drift.1 <- y.rw1.drift ~ f(id.w, model = "rw1", constr = FALSE) - 1 + id.delta
　あるいはこういう風に書くという手もある。
　formula.rw.drift.2 <- y.rw.drift ~ f(id.w, model = "rw1", constr = FALSE) - 1 + f(id.delta, model = "linear", mean.linear = 0, prec.linear = 0.001)

　[結果の取り出し方について。パス]

3.6 二次多項式モデル
　ドリフトが時変するモデルを考えよう。これはふつう二次多項式モデルと呼ばれる。$$ y_t = x_t + v_t$$ $$ x_t = x_{t-1} + \delta_{t-1} + w_{t,1} $$ $$ \delta_t = \delta_{t-1} + w_{t,2}$$
　状態ベクトル\(\mathbf{x}^*_t = (x_t, \delta_t)^\top\)を考えて書き直す。$$ y_t = \left( \begin{array}{c} 1 & 0 \end{array} \right) \mathbf{x}^*_t + v_t; \ \ v_t \sim N(0, V)$$ $$ \mathbf{x}^*_t = \Phi \mathbf{x}^*_{t-1} + w^*_t $$ $$ \Phi = \left( \begin{array}{cc} 1 & 1 \\ 0 & 1 \end{array} \right) $$ $$ \mathbf{w}^*_t \sim N \left( \mathbf{0}, \left( \begin{array}{cc} \sigma^2_{w1} & 0 \\ 0 & \sigma^2_{w2} \end{array} \right) \right) $$

さて、元の式に照らすと、$$ 0 = x_t - x_{t-1} - \delta_{t-1} - w_{t,1}$$ $$ 0 = \delta_t - \delta_{t-1} - w_{t,2}$$ である。そこで、目的変数を\(y_1, \ldots, y_n\)と捉えるのではなくて、次の3つの行列を縦に積んだものだと考える。

• \(n\)行3列の行列で、1列目に\(y_t\)がはいっている。ほかはNA。
• \(n-1\)行3列の行列で、2列目に0が、他にはNAが入っている。
• \(n-1\)行3列の行列で、3列目に0が、他にはNAが入っている。

　モデルのあてはめにあたって、R-INLAは列ごとに異なる尤度を考える。

• 1列目のデータ\(y_t\)に対しては、状態ベクトル\(\mathbf{x}^*_t\)の下で、未知の精度\(V^{-1}\)を持つ正規分布に従うと仮定される。
• 偽のデータである2列目と3列目に対しては、\(\mathbf{x}^*, \mathbf{x}^*_{t-1}, \mathbf{w}^*_t\)の下で、既知の分散0を持つ決定的変数だと仮定される。R-INLAはこれを固定された高い精度を持つ正規分布として表現する。

　\(\mathbf{x}_t\)の事前分布はこうなる。[...略...]

　[以下、この定式化に基づくコードの紹介。あまりにトリッキー過ぎてついていけないのでパス]

3.7 状態・観察の予測
　R-INLAでは学習データの末尾にNAをつけることで予測できる。[...後略...]

3.8 モデル比較
3.8.1 in-sampleの比較
　次の3つが計算できる。inla(..., control.compute = dic(dic = TRUE))という風に指定する。

• DIC. パラメータの事後平均を\(\hat{\theta}\), 対数尤度を\(L(\Theta; \mathbf{y}^n)\)とし、\( D = -2 \log L(\Theta; \mathbf{y}^n)\)とし、\(\bar{D} = E_{\Theta | \mathbf{y}^n} (D)\)として、\(\bar{D} + (\bar{D} - D(\hat{\Theta}))\).
• conditional predictive ordinate(CPO) \(p(y_t|\mathbf{y}_{(-1)})\)
• probability integral transform(PIT) \(p(y^{new}_t \leq y_t | \mathbf{y}_{(-t)})\)

CPOとPITはモデル適合のleave-one-out型の指標である。なお、result$cpo$failureというのがついてくる。これはCPOを求めるときに、INLAが置いている暗黙的な仮定が破られたかどうかのベクトルで、すべて0になっていてほしい。そうでないときは再計算するべし。PITにも同じのがある。
　それから、周辺尤度の対数も出力に入っているよ。

3.8.2 out-of-sampleの比較
　MAPE, MAEを使うがよろしい。コード例は14章をみよ。

3.9 デフォルトでない事前分布の指定
　[具体的な指定方法の話なので省略するけど、デフォルトの事前分布のパラメータを変える、ないし事前分布そのもの指定する方法のほかに、table priors(密度を表で示す)、expression priors(密度をformulaで示す)、という方法もあるのだそうな。
　また、普通の事前分布を指定する代わりに、penalized compexity事前分布というのを指定することもできる由。詳細はinla.doc("pc.prec")をみよとのこと]

3.10 潜在効果・ハイパーパラメータの事後分布からのサンプリング
[パスするけど、コード例がついている。これは助かることがあるかも]

3.11 未知の観察に対する事後予測標本
[パス。やはりコード例がついている]
-----------
　フガフガとめくっただけだけど、R-INLAの設計思想みたいなのがだんだんわかってきた。formulaでは観察方程式の各項の指定を並べ、各項の指定の際にその項に対応する状態方程式の指定を埋め込む、という感じなのか。だから、3.5節の場合、状態方程式に定数が入っているというより、観察方程式に時点の項が入っていると考えないといけない。要するに、dlmパッケージやKFASパッケージを使うときと同様、ユーザは自分の組みたいモデルをどうにかしてきれいな状態空間表現(状態方程式には定数項がない)に落としこまなきゃいけない、ってことかな？
　しかも、状態の指定はそれほど柔軟にはできないので、複数の状態変数の間に依存性が生じるような場合には工夫が必要になる。3.6節の場合、dlmやKFASなら素直に書けるモデルだと思うのだけれど、R-INLAの場合は、ダミーデータをくっつけて状態変数間の関係を無理やり観察方程式の側で表現してやるというようなトリックが必要になるわけだ。

　うむむむ。R-INLA、意外にハードル高いなあ... このハードルの高さは、R-INLAの設計に由来するのだろうか、それとも、「状態はGMRFに従う」というINLAの制約に由来するのだろうか。

2023/09/30追記: 3.8節以降を追加。前半と後半の2本にわけるつもりだったんだけど、この1本にまとめた。