King, A.A., Ngyen, D., Ionides, E.L. (2016) Statistical Inference for Partially Observed Markov Processes via the R Package pomp. Journal of Statistical Software, 69(12), 1–43.
仕事の都合でマルコフ・スイッチング・モデルについて調べているんだけど、CRAN Task ViewsのTime Seriesで挙げられているマルコフ・スイッチング・モデルに関係しそうな雰囲気のパッケージとしてNTS, MSwM, depmixS4, pompがある。このうちpompってやつが一体なにやってんのか全然わかんなかったので、試しにvignetteをめくってみた。おそらく同内容の論文がJ. Statistical Softwareに載っている(上記)。
1. イントロダクション
部分的に観察されたマルコフ過程(POMP)モデルというのは、観察されていない潜在的なマルコフ過程の不完全でノイジーな観察からなる。
pompパッケージが実装しているアルゴリズムは任意のPOMPモデルに適用可能なはずである。もっとも、追加的なモデル構造について検討するための方法は実装していない。特に、線形でガウシアンな近似が適切な時や、潜在過程が少数の離散集合の値をとるときは、アンサンブルカルマンフィルタとか正確隠れマルコフモデルのような追加的な仮定を置く方法があるのだが、まだ実装していない。pompパッケージが焦点を当てるのは、非線形で非ガウシアンで状態空間が大きいPOMPモデルである。
[ここからよくわからないので細かくメモ:] POMPモデルはマルコフ過程の遷移密度と測定密度で特徴づけられるのだが、遷移密度に基づくシミュレーションだけが必要になるような手法もあれば、密度の評価が必要な手法もある。また、手法の中には線形化のような近似をしないといけないのもある。シミュレータだけを使って動的モデルを特定するようなアルゴリズムのことをplug-and-play法という。多くのPOMPモデルがこれでいける。さて、ある候補モデルがtractableな遷移確率を持っていたとしても、遷移確率がtractableでないような別のモデルについて考えたくなるかもしれない。plug-and-play法の場合、モデルにおける変動分析をのためには、シミュレータのコードを数行書き換えるだけでよいことが多い。plug-and-play法はこのように柔軟だが、その代償は計算コストである。
[…]
2. POMPモデルとpompにおけるその表現
パラメータを\(\theta \in \mathbb{R}^D\)とする。\(\{X(t; \theta), t \in T\}\)がマルコフ過程で、\(X(t; \theta)\)が\(\mathbb{R}^{dim(X)}\)の値をとるとする。\(T \subset \mathbb{R}\)は間隔でも離散でもよい。
観察できるのは\(\{t \in T, i = 1, \ldots, N\}\)における\(X(t;\theta)\)とする。初期時間を\(t_0 \in T\)とし、\(t_0, t_1, \ldots, t_N\)は順に並んでいるとする(\(t_0と t_1\)の間のみ\(t_0 \leq t_1\)で、あとは等号がつかない)。以下、\(X_i = X(t_i; \theta), X_{i:j} = (X_i, X_{i-1}, \ldots, X_j)\)と略記する。
観察できるのは\(Y_{1:N}\)だけで、それは\(X_{0:N}\)の下で条件付き独立な確率変数だとする[ええええ。さっき\(X\)が観察できるって言ったじゃん～]。データ\(y^*_{1:N}\)は固定された値とし、この観察過程の実現値としてモデル化する。
POMP構造とはこういうのである。\(X_{0:N}, Y_{1:N}\)は同時密度\(f_{X_{0:N}, Y_{1:N}} (x_{0:N}, y_{1:N}; \theta)\)を持つとする。この同時密度が初期密度\(f_{X_0}(x_0; \theta)\)と条件付き遷移確率密度\(f_{X_n|X_{n-1}}(x_n | x_{n-1}; \theta)\)と観察密度\(f_{Y_n|X_n}(y_n|x_n; \theta)\)で決まると考える。すなわち、$$ f_{X_{0:N}, Y_{1:N}} (x_{0:N}, y_{1:N}; \theta) = f_{X_0}(x_0;\theta) \prod_{n=1}^N f_{X_n|X_{n-1}}(x_n|x_{n-1}; \theta) f_{Y_n|X_n}(y_n|x_n; \theta)$$
[ぶひー。難しいぞ。ここまでを平たく言い換えよう。まずパラメータ\(\theta\)があります。そのもとでの、連続状態マルコフ過程\(X_i\)があります。\(X_0\)の密度を\(f_{X_0}\)とし、時点\(n\)の条件付き密度を\(f_{X_n|x_{n-1}}\)とします。で、このマルコフ過程にiidなノイズが乗ったのが\(Y_i\)で、その条件付き密度を\(f_{Y_i|x_i}\)とします。
\(X_{0:N}\)と\(Y_{1:N}\)の同時密度関数はこうなります。初期密度\(f_{X_0}\)、かける、時点1の条件付き密度\(f_{X_1|X_0}\)かける観察密度\(f_{Y_1|X_1}\)、かける、時点2の条件付き密度\(f_{X_2|X_1}\)かける観察密度\(f_{Y_2|X_2}\)、かける…。
要は\(Y\)が観察変数、\(X\)が状態変数で、状態変数はAR(1)ないしVAR(1)みたいなもんなわけね。あれ？ 観察方程式は時変すんの？ つまり\(f_{Y_n|X_n}(y_n|x_n; \theta)\)は\(n\)に依存すんの？ どうやらそうみたいね。ぶひー、一般的すぎる～～]
2.1 POMPモデルの実装
pompパッケージでは、rprocessメソッドで\(f_{X_n|X_{n-1}}\)からのシミュレーション、dprocessメソッドで\(f_{X_n|X_{n-1}}\)の評価をする。またrmeasureメソッドで\(f_{Y_n|X_n}\)からのシミュレーション、dmeasureメソッドで評価をする。plug-and-play法ではdprocessを使わない。[あ、plug-and-playって遷移確率密度を決め打ちするってことなの？]
[以下、力尽きたのでパス。小節の見出しのみメモする]
2.2 初期状態
2.3 共変量 [あ、共変量時系列をいれられるんだ… \(f_{X_n|x_{n-1}})\)も\(f_{Y_n|X_n}\)も共変量に依存させられるらしい]
3. POMPモデルの方法論
[もう無理です、許してください、というわけでまるごとパス。小節の見出しのみメモする]
3.1 尤度関数と系列的モンテカルロ [粒子フィルタのことらしい。知らんけど]
3.2 反復フィルタリング
3.3 粒子MCMC [なにそれ…]
3.4 要約統計量のsynthetic尤度
3.5 近似的ベイジアン計算(ABC)
3.6 非線形予測
4. モデル構築とデータ分析: 単純な例
[気を取り直してちょっとだけ読んでみよう]
4.1 最初の例: Gompertzモデル
人口成長のGompertzモデルについて考えよう。このモデルによれば、時点\(t + \Delta t\)の人口はこうなる。人口容量を\(K\)、正のパラメータを\(r\)として$$ X_{t+\Delta t} = K^{1-\exp(-r \Delta t)} X_t^{\exp(-r \Delta t)} \epsilon_t$$ $$ \log \epsilon \sim_{iid} N(0, \sigma^2)$$ さらに観察誤差を考えて $$ \log Y_t \sim N(\log X_t, \tau^2) $$ とする。なお、このモデルは線形ガウシアンなのでカルマンフィルタでいける。plug-and-play法はいらない。
pompパッケージの使い方。(1) \(X_t\)の更新を表す関数を定義する。(2)\(Y_t\)を表す関数を定義する。(3)観察密度を表す関数を定義する(ここではdlnorm()を使えばよい)。で、(4)pomp()をコール。引数dataにデータを、rprocessに(1)を、dmeasureに(2)を渡す。[ちゃんと読んでないけど、要はモデルを関数の形で自力で書いて渡せってことね… ]
[めんどくさくなってきたのであとは見出しのみメモ]
4.2 SMCによる尤度計算
4.3 反復フィルタリングによる最尤推定
4.4 PMCMCによる完全情報ベイズ推論
4.5 第二の例: Rickerモデル
4.6 feature-based synthetic尤度最大化
4.7 ABCによるベイジアンfeature matching [なんだそれは]
4.8 simulated擬似尤度によるパラメータ推定
5. より複雑な例: 連続時間の疫学
[もちろんパス。小節の見出しのみメモする]
5.1 確率的な季節SIRモデル
5.2 pompでの実装
5.3 さらなる複雑性の統合
6. 結論
[略]
——————-
ひいいい。難しくてよくわかんなかった。
でも、まあ、どんな感じのパッケージなのかはなんとなくわかったのでよしとしよう。要するに、時系列について、それが連続マルコフ過程にノイズが乗ったものだと考えてモデル化したりする、んですよね。それって、要は状態空間表現できるようなわりかし普通のモデルであって(そうですよね？)、線形ガウシアンなら別にカルマンフィルタで推定できる。でもそうじゃない場合のために、なんかその… こう、いろいろと… 工夫？とか？ があるわけだ。知らんけど。
えーと、これ、マルコフ・スイッチングとはちょっと違う話ということで合ってますでしょうか… (読了にしておくけれど全然わかってない)