1.予測理論この記事では、状態空間モデルをStanで推定するときの収束を良くするコツを説明します。コードはGitHubから参照できます。状態空間モデルは説明能力が高く、データに合わせて柔軟に構造を変えることができます。あまりに複雑な構造を指定すると、結果が収束しないこともしばしばあります。収束が悪い時には、弱情報事前分布を指定したり、MCMC実行時の設定を変える(iterやwarmupを増やす等)で対応することが多いと思います。こちらの方法で多くの場合は解決しますが、複雑なモデルですと、そもそものStanコードの実装の方法から変えた方が良いかもしれません。この記事では、状態空間モデルの収束を良くするために、Stanコードの書き方を工夫するやり方を紹介します。modeling』という記事を見つけました。上記の記事で提案されている方法を参考にしています(コードは少し変えました)。この記事の方法を使うと、収束が良くなることがあるようですが、正直コードの可読性は下がってしまいます。ミスが発生する可能性が増えるので注意してください。万が一、この記事の内容に誤りがあった場合は、コメントやメールなどでご連絡いただければ幸いです。この記事では、読者がすでにR言語やStan、そして状態空間モデルの基本について知っていることを前提としています。ベイズ統計モデリングによるデータ分析入門』などを参照してください。一部のコードは、この本からの引用となります。2019の14日目の記事となります。2019」を参照してください。このリンクをたどると、Stanに関する様々な記事を読むことができます。基本方針ベクトル化したローカルレベルモデル再パラメータ化したローカルレベルモデル再パラメータ化したローカル線形トレンドモデル1.シミュレーションデータの作成まずは、分析対象となるデータを作成します。ローカル線形トレンドモデルに従う乱数を生成します。コードはGitHubから参照できます。sd_v)}シミュレーションデータの折れ線グラフを描きます。"シミュレーションデータ")最初は全体的に減少傾向にあるデータですが、そのあとで上昇トレンドに変化していることがわかります。最後に、後で使いやすくするためにlistにまとめておきます。先ほど作成したシミュレーションデータを対象として、まずは素朴な実装を行ったStanコードでモデルを推定します。ローカルレベルモデルまずはローカルレベルモデルを推定します。シミュレーションデータはローカル線形トレンドモデルに従っているので、ローカルレベルモデルではちょっと力不足なのですが、コードが短くて簡単なため載せておきます。ベイズ統計モデリングによるデータ分析入門』からの引用です。『local-level.stan』というファイル名称で保存されているとします。}}こちらでMCMCを実行すると、収束に関するワーニングが出ました。MCMCの結果は以下のようになっています。Rhatは1.1未満ですが、ワーニングがあるので、改善の余地はありそうです。1.00ローカル線形トレンドモデル続いてローカル線形トレンドモデルです。ローカル線形トレンドモデルは、説明力の高い優れたモデルですが、経験上、収束しにくいように思います。ベイズ統計モデリングによるデータ分析入門』のサンプルコードを一部改変したものです。『local-linear-trend.stan』というファイル名称で保存されているとします。}}こちらのコードを実行すると、ワーニングが多く出ます。Rhatも1.1を超えており、収束したとは言えない状況です。1.01トレースプロットを見ても、問題があるのがわかります。"s_v"))誤解の無いように記しておくと、先のStanコードに大きな問題があるわけではないです。素朴な実装ですし、初学者の方が最初に状態空間モデルの実装をする場合は、こちらの素朴な実装から入るべきだと思います。データを作成したシミュレーションのコードを比較してもらうとわかるように、こちらの実装は「データの生成過程をそのまま」コードで表現できています。ベイズ統計モデリングの勉強という意味では、むしろこちらの素朴な実装を使う方が好ましいと思います。収束を良くするためのコードは、それなりにテクニカルな実装になります。今回はあえてワーニングがたくさん出るように、何の工夫もせずに実行させてあります。上記の素朴なStanコードであっても、iterやwarmup、thinを増やしたり、弱情報事前分布を指定したりすることで、収束を改善させることができます(拙著ではその方法を使っていました)。いえ、この記事のテーマはStanコードの書き方を工夫することで収束を良くするということです。iterやwarmupはそのままにして、収束を改善することを試みます。基本方針Stanコードを修正する基本的な方針は2つあります。ベクトル化をすることです。ベクトル化というのは少し難しいように聞こえる用語ですが、forループをやめて実装するイメージです。後ほど具体例を挙げます。ベクトル化しても、収束がすぐに良くなるわけではありません。そのうえで再パラメータ化をします。再パラメータ化がキモで、これで収束を良くします。再パラメータ化をするとサンプリングが遅くなることもあるので、その点は注意が必要です。再パラメータ化も、言葉では説明が難しいので、実際のコードを後ほど紹介します。ベクトル化したローカルレベルモデルまずはローカルレベルモデルのStanコードをベクトル化します。『local-level-vec.stan』というファイル名で保存しておきます。s_v);}modelブロックのみ変更があります。muの2時点目は、平均値が『muの1時点目』である正規分布から得られます。muのT時点目は、平均値が『muのT-1時点目』である正規分布から得られます。s_w);としています。観測方程式はさらに短いコードになっています。長さが同じなので、角カッコを使って添え字を指定する必要もありません(添え字を入れても別に良い)。ベクトル化をすることで、コードが短くなりました。実効速度が速くなることもあります。収束に関するワーニングは引き続き出てきているようです。再パラメータ化したローカルレベルモデル続いて再パラメータ化をしたローカルレベルモデルのコードを載せます。『local-level-remodeling.stan』という名称で保存しておきます。s_v);}ベクトル化したコードと比較すると、dataブロック以外はすべて変わっています。parametersブロックから、状態推定値muがなくなりました。その代わりにmu_errという「水準成分の増減量」があります。parametersブロックはいったん飛ばして、modelブロックを見ます。mu_err[1]は、無情報事前分布を指定しました。modeling』という、最初に挙げた参考資料と異なったやり方です。いえ、初期状態が0近辺にあるという保証はないので、ここは無情報事前分布を使うほうが自然かと思いました。ここはいろいろの実装があり得そうです。mu_err[2]以降は、標準偏差1の正規乱数とみなしています。観測方程式は、ベクトル化したバージョンのコードと同じです。parametersブロックです。1時点目のmuはmu_errと同じとしています。2時点目以降が重要です。forループの中です。標準偏差1の正規乱数すなわちmu_errに対して、過程誤差の標準偏差s_wを掛け合わせることで、実際の過程誤差としています。このように「増減量を外だし」してあげることで、モデルの収束が良くなるようです。このような事例はGLMMなどの階層モデルでしばしばみられますね。実行すると、ワーニングが減りました。adapt_deltaを変えるよう促されています。problemsワーニングメッセージの通りに、adapt_deltaやmax_treedepthを増やして再実行させます。すると、きれいにワーニングがなくなりました。)再パラメータ化したローカル線形トレンドモデル最後に、ローカル線形トレンドモデルのStanコードを修正します。『local-trend-remodeling.stan』という名称で保存しておきます。s_v);}基本的な方針は、ローカルレベルモデルと同じです。増減量をmu_errやdelta_errとして外だしします。これらを標準偏差1の正規分布だとみなし、これにs_wやs_zという標準偏差を表すパラメータを掛け合わせることで、実際の状態の推定値を得ています。上記のコードを使って実行した結果はこちらです。大きく改善されました。いえワーニングはまだ少し残っているようです。problemsワーニングに従い、adapt_deltaやmax_treedepthを増やすと、ワーニングは一切出ません。)推定結果を確認すると、Rhatは1になっているようです。1状態推定値も含めて、すべての推定値のRhatを確認しても、問題ないようです。収束の確認mcmc_rhat(rhat(mod_llt_remodeling_2))トレースプロットなどを得ると、とてもきれいな形になりました。ベイズ統計モデリングによるデータ分析入門このブログの管理人が書いたベイズ統計モデリングの入門書です。RとStanを使ってベイズ統計モデリングをするための基本をまとめています。ベイズ統計モデリングをやってみようという方にお勧めします。本記事は、この書籍を補足する意味も持たせています。今後も、書籍をサポートする内容を充実させていきたいと思います。サポートページはこちらです。StanとRでベイズ統計モデリングこちらはより応用的な内容が載っているベイズ統計モデリングの教科書です。Stanのコードも豊富にあります。コメントをキャンセルメールアドレスが公開されることはありません。コメントデータの処理方法の詳細はこちらをご覧ください。technology.。