[Stan] 状態空間モデルでシステムノイズに非ガウス分布
2014-12-31
北川本「時系列解析入門」14.4節の例題をやってみましたので途中経過の作業ログとして残します。まずは失敗例から。
16行目でシステムノイズにCauchy分布を使っていますが、このモデルは全く収束しません。
これを解決するヒントはStanのマニュアルの「19. Optimizing Stan Code」（の「Reparameterizing the Cauchy」節）にあります。Cauchy分布では幹の部分をうろちょろするのと比べて裾の重い部分をうろちょろするのは相対的にステップサイズが小さくなるため、rejectが増えて時間がとてもかかるようになります。こういう状況を解消する方法がこの章にはいくつか載っています。これをうけての改良版が以下。
・9,16,17行目: 累積分布関数を使う方法でCauchy分布に従う乱数を発生させています。この一様分布に従う乱数（mu_unif）を代わりに推定します。
キックするRコードは以下です。
・4-8行目: データを作成しています。
・2,11行目: hoxo_m氏作成のRの並列化ライブラリpforeachを使っています。超使いやすいです。
計算時間は1chainあたり数分でした。Rhatはすべて<1.1です。推定結果は以下のとおりです。
細い黒線で振動が激しいのが実データ、帯は外側から順にMCMCサンプルの99.9%信用区間、98%信用区間、85%信用区間です。真ん中の黒線はMCMCサンプルの中央値です。
さらに書籍に載っているコーシー分布より裾の重いピアソン分布族で実行したいです。これもモデルの組み方に工夫が必要でした。参考になったのは以下の書籍の9章です。web版を公開してくださっており大変感謝です。
Luc Devroye, "Non-Uniform Random Variate Generation", (SpringerVerlag, New York, 1986)
この本によりますと欲しいピアソン分布族（VII）に従う乱数は以下の手順で生成できます（書籍の b-1/2 を b と置き換えています）。
これをそのまま実装したのが以下になります。
・4行目: ピアソン分布族のパラメータです。B=0.5の時がCauchy分布に対応します。
書籍の方ではB=0.1の結果があったのでそれを実行しました。計算時間は1chainあたり数分でRhatはすべて<1.1です。推定結果は以下のとおりです。
ちなみにBをパラメータとして推定させることもできました。その時はBの事後分布の中央値は0.81・平均値は0.79になりました。
- 関連記事