vec 演算子とトレース
行列を以下の様に定義する。
ここで は縦ベクトルとする。
このとき
また基本的な操作として
がある。
{rpart}{partykit}のプロットで日本語を使用する方法
自分用メモ。実行例は後日載せるかも。
fit.dt <- rpart(y ~., data) plot(as.party(fit.dt), gp = gpar(fontfamily = "Osaka"))
pythonによる粒子フィルタの実装
2階差分トレンド+季節(週)トレンドを考慮した以下の状態空間モデル(線形ガウス状態空間モデル)
について、粒子フィルターをpythonで実装しました。
ここではそれぞれ時刻における観測値、トレンド(平均)、季節トレンドです。
状態空間モデルの詳しい説明については
予測にいかす統計モデリングの基本―ベイズ統計入門から応用まで (KS理工学専門書)
- 作者: 樋口知之
- 出版社/メーカー: 講談社
- 発売日: 2011/04/07
- メディア: 単行本(ソフトカバー)
- 購入: 9人 クリック: 180回
- この商品を含むブログ (12件) を見る
が分かりやすいです。
行列とベクトルの表記を使って綺麗に状態空間モデルのモデリングや種々のアルゴリズムの導出を行っていますし、P.89には粒子フィルタの疑似コードが載っています。
ただしこの本に従って実際に実装するとなると、モデルによっては膨大な量の行列とベクトルを実装しなければなりません。
特に拡大状態ベクトルを用いて固定ラグ平滑化を導出していますが、いざ固定ラグ平滑化を実装しようとなると拡大状態ベクトルの実装が面倒だと感じます。また固定ラグ平滑化については実装例が載っていません。
そこで、固定ラグ平滑化まで含めて粒子フィルタをpythonで実装しましたので参考にしていただけたらと思います。
使用したデータは上記書籍のサポートページ
http://daweb.ism.ac.jp/yosoku/
からダウンロードできます。
固定ラグ平滑化の概要は以下であると解釈しています。予測とフィルタリングは普通の粒子フィルタです。縦に5つ並んでいる丸が粒子で、この図では5つの粒子が描かれています。は時点でのデータで条件付けられた時点での番目の粒子ということを表しています。粒子フィルタはシステムモデルによって時点のデータから時点の値を"予測"し、実際に時点でのデータが手に入った時点で観測モデルの尤度でもって状態をフィルタリング(更新)します。
さて、時点での最新のデータが手に入った時点で例えばやの過去のデータを更新するというのが平滑化という操作で、平滑化の方法はMCMCといった方法で実行することもできますが、粒子フィルターの枠組みでは固定ラグ平滑化という方法が一般に用いられます。その固定ラグ平滑化の概要は上図のように、最新のデータでもって現在の状態をフィルタリングした後に、そのフィルタリングされた粒子は前の時点ではどの粒子から派生したものかというのを逆向きに辿ることで平滑化ができる(と私は解釈しています)。
青い四角で囲まれた数字は、前の時点のどの(上から何番目の)粒子から派生したかのインデクスを表しています。上図ではラグ幅が1の固定ラグ平滑化を実行していることになっていますが、プログラムを実装する上では上図のように各時点の粒子がどの粒子から派生したものかのインデックスを記憶しておいて、それをデータが手に入ってフィルタリングをするたびに記憶しておいたインデックスを頼りに過去に辿っていけば平滑化ができるということです。
今回は冒頭で紹介した書籍を参考に、pythonで実際に粒子フィルタを実装をしてみたわけですが、これをStanで実装した例として@berobero11さんのスライドやブログ
は大変参考になります。
間違いや質問等ございましたら@kefismまで遠慮なくご連絡ください。