RcppSMCパッケージを使って逐次モンテカルロ（Sequential Monte Carlo）・粒子フィルタ（Particle Filter)を実行する
前フリ
逐次モンテカルロ・粒子フィルタというのは状態推定*1をシミュレーションに基づいて実行する技でして、最近は「データ同化」なんてワードと絡めて有名になってきている気がしています。アイディアとしてはカルマンフィルタの非線形拡張＆解析ベースではなくシミュレーションベースでの状態推定という感じ。
入門書としては
の本がお勧めで、バッチリ6章に粒子フィルタの解説が載っている。
RcppSCMパッケージ
このライブラリの計算エンジン自体はC++で実装されてて「 Johansen, 2009, J Statistical Software, 30:6」で提案されたC++ライブラリを使っている模様。そしてそのC++のコードをRcpp経由でRから叩けるようにしているのがこのパッケージだぞとそういうことです。
パッケージのインストール
このパッケージを使うためにはRcppSMCパッケージはもちろん、Rcppパッケージも必要なのですが、そのインストールがちょいと面倒。昔書いた記事Rcpp・inlineパッケージを使ってC++とRを連携させる - My Life as a Mock Quantがあるのでこれを参考にするか適当にググって調べましょう。
RcppSMCパッケージ本体は他のパッケージ同様に
install.packages("RcppSMC")
を実行すればOK。
サンプルデータの作成
まずは状態推定をおこなうためのサンプルデータの作成を行う。そのための関数である「simGaussian関数」はその名の通り「線形ガウシアン状態空間モデル」に従うシミュレーション値を返却してくれるもので、式で書くと
状態方程式：
観測方程式：
に従って生成されるx, yを返却してくれる。e,fはそれぞれ独立な標準正規分布に従う乱数。
実行するとこんな感じで状態変数(x)と観測データ(y)が取得できる。
> sim <- simGaussian(len=10) > sim $state [1] -0.09427418 2.82961623 4.70898832 5.06479708 5.84263024 6.49104229 [7] 6.42444557 5.72142367 6.69790596 6.69520952 $data [1] -1.253615 4.063993 5.038656 5.710719 5.553496 5.908212 7.038596 5.616786 [9] 7.454828 8.131660
stateがxでdataがyに相当してて、ここで取得したyのみの値を使ってxを推定したい！というそんな状況設定です。
次にここで取得できたdata(y)の値を使ってstateの値を推定するためにblockpfGaussianOpt関数を使う。この関数で実装されてる手法はまだ読めてないけど以下の論文によるもの（らしい）
実際に以下のように実行してみると(データは少なすぎるので数を増やして再作成）
sim <- simGaussian(len=250) res <- blockpfGaussianOpt(sim$data,lag=5,plot=TRUE) lines(sim$state, lty=1, lwd = 3, col="red")
こんな感じの図が得られる。青が推定値で赤が実際の状態の値でなかなかよく推定出来ているように見える。予測値の誤差は多分こんな感じで計算できて・・・
error <- as.numeric(t(res$weight) %*% res$values / sum(res$weight)) - sim$state plot(error, ylab ="Error")
となる。最大で３程度のズレはあるもののそれなりな当てはまりか。
パッケージ自体がかなり開発途上っぽくて、残りのpfLinear・BSpfNonlinBS関数についても論文にあるサンプルを再現できるだけで、自分でモデルを突っ込んで云々はできないので一旦ここでやめ。また続報があり次第書く。
【参考】
RcppSMC
