Visser, I., & Speekenbrink, M. (2010). depmixS4: An R Package for Hidden Markov Models. Journal of Statistical Software, 36(7), 1–21.
RパッケージdepmixS4のvignette。上記はJSSの論文だが、実際にはパッケージ添付のvignetteを読んだ。
実戦投入検討前のセレモニーとして大急ぎでめくった奴なので、ほとんどの部分をちゃんと読んでいない。
1. イントロダクション
[時間がないのでパス]
2. 従属混合モデル
長さ\(T\)の\(m\)変量時系列 $$ \mathbf{O}_{1:T} = (O^1_1, \ldots, O^m_1, O^1_2, \ldots, O^m_2, \ldots, O^1_T, \ldots, O^m_T) $$ について考える。このうち時点\(t\)の部分を\(\mathbf{O}_t\)と略記する。
多変量時系列について潜在マルコフモデルを考えることは多いけど、推定方法はたいてい長い時系列に向かない。いっぽう隠れマルコフモデルは長い時系列にしか使われない。これらを関連付けるために従属混合モデルという新しい名前を使う。
従属混合モデルの基本的想定はこうだ。いかなる時点でも、観察は\(n\)個の要素(状態)の混合として分布する。観察の時間依存性は、混合要素の間の時間依存性、つまり要素間の遷移確率によって生じる[え、じゃあ状態と無関係なAR(1)攪乱項とかは考えないってこと？]。この依存性は一次マルコフ過程に従う。
というわけで、潜在状態を\(\mathbf{S}_{1:T} = (S_1, \ldots, S_T)\)、共変量を\(\mathbf{z}_{1:T} = (\mathbf{z}_1, \ldots, \mathbf{z}_T)\)として、同時尤度は… $$ P(\mathbf{O}_{1:T}, \mathbf{S}_{1:T} | \mathbf{\theta}, \mathbf{z}_{1:T}) = \pi_i (\mathbf{z}_1) \mathbf{b}_{S_1} (\mathbf{O}_1 | \mathbf{z}_1) \prod_{t=1}^{T-1} a_{ij}(\mathbf{z}_t) \mathbf{b}_{S_t}(O_{t+1} | z_{t+1}) $$ えーと、まず\(S_t\)ってのは状態ね。\(n\)個あることにする。
\(\pi_i(\mathbf{z}_1)\)ってのはつまり\(P(S_1 = i | \mathbf{z}_i)\)のこと。時点1における状態の確率である。
\(a_{ij}(\mathbf{z}_t) = P(S_{t+1} = j | S_t = i, \mathbf{z}_t)\)。つまり状態\(i\)から\(j\)への遷移確率。
最後に、\(\mathbf{b}_{S_t}\)というのは、その要素が観察密度 \(b^k_j (\mathbf{z}_t) = P(O^k_t | S_t = j, \mathbf{z}_t) \)であるようなベクトル [ハナから多変量で、しかも共変量を入れて定式化しているせいですごくややこしくなっているが、要するに、潜在状態の下での観察の条件付き確率のことだ]。\(b^k_j\)の分布は、正規でもベルヌーイでもなんでもよい。
2.1 尤度
[時間ないのでパス]
2.2 パラメータ推定
[時間ないのでパスするけど、デフォルトではEMアルゴリズムで最尤推定。ほかにRsolnpパッケージ、Rdonlp2パッケージというのをコールして直接的に推定することもできるんだそうな。とにかくベイズ推定ではないわけね]
2.3 欠損データ
[パス]
3. depmixS4パッケージの使い方
[以下、すごく粗いメモ…]
まずdepmix()でモデルを指定してfit()でフィッティングする、というのが基本。前者は、潜在クラスモデル・有限混合モデルの場合はmix()というのを使ってもよい。
3.1 データ例: speed
[反応時間の速度-正確性トレードオフの話が説明されている。いずれちゃんと読もうっと]
3.2 単純なモデル
depmix()のresponse引数にフォーミュラを私、nstates引数に状態数、trstart引数とかに遷移パラメータの初期値を渡す。で、できたオブジェクトをfit()に渡す。
3.3 遷移パラメータに対する共変量
[読まずにパスしたけど、遷移行列の説明変数をいれられるみたいね]
3.4 多変量データ
[えっ、多変量もいけるの？ へーそうなんだ。じゃあマルコフ・スイッチングVAR(1)とかもできるわけ？ いずれ読み直そう]
3.5 パラメータの固定・制約
[ここはいま関心があるので細かくメモする]
パラメータに一般的な線形等値制約や非等値制約をかけることができる。その場合はRsolnp, Rdonlp2パッケージがコールされる。
制約はfit()のconpat引数で指定する。それぞれのパラメータについて、固定したかったら0, 自由推定したかったら1以上の値を指定する。等値制約はパラメータに同じ値を与えれば指定できる。固定だけならconpat引数じゃなくてfixed引数でもよい。
いずれにせよ、パラメータの番号を知る必要がある。
パラメータ数はnpar(mod)でわかる。パラメータの順序は setpars(mod, valud = 1:npar(mod))でわかる。事前モデルパラメータ、遷移モデルパラメータ、状態ごとの反応モデルパラメータ、の順に並んでいる。
[…] というわけで、ユーザは好きなように制約を掛けられる。自分の責任でやるように。
3.6 事前確率に対する共変量の追加
[ふーん… パス]
4. depmixS4の拡張
[反応の分布をユーザが細かく指定できる機能があるそうな。パス]
5. 結論と将来の課題
ほかにflexmixパッケージというのもあるよ。云々。
————-
ふーん… 要するに、glm()の全パラメータをレジームスイッチできるパッケージなわけね。MSwMパッケージやNTSパッケージなんかとちがって、パラメータごとにこまかく制約を掛けられたり、多変量でも扱えたり、遷移行列の説明変数をいれたりできるらしい。
あれ？ ちょっと待てよ、このパッケージ、測定モデルの攪乱項にダイナミクスをいれることはできるの？ 攪乱項が自己回帰するとか？ たぶん駄目なんだろうな… うーん…