Rで粒子フィルタ：RcppSMCのテスト

経緯

Rで粒子フィルタ（逐次モンテカルロ）ってできるんだっけ、とふと思ったので少し調べてみました。すると、Rのパッケージでは、例えば

• https://cran.r-project.org/web/packages/SMC/SMC.pdf : {SMC}
• https://cran.r-project.org/web/packages/RcppSMC/RcppSMC.pdf : {RcppSMC}

というパッケージがあるようです。{RcppSMC}の方は、Rcppで開発されているパッケージです。まずは、簡単にテストできそうな{RcppSMC}を使ってみることにしました。

参考記事

ネットで探しているときに目についた方のブログを参考にしました。
というかほぼ丸パクリです。

d.hatena.ne.jp
ito-hi.blog.so-net.ne.jp

また、粒子フィルタに関する書籍も載せておきます。自分が読んだことあるものは、例えば以下です。ざっくり概要を知るのは樋口さんの本のほうがよさげです。
www.kspub.co.jp
時系列解析入門https://www.iwanami.co.jp/.BOOKS/00/4/0054550.htmlwww.iwanami.co.jp

テスト

パッケージのインストール

if(!require(RcppSMC)){
  install.packages("RcppSMC", dependencies = TRUE)
}
require(RcppSMC)

テストデータの生成

線形ガウス状態空間モデルのシミュレーションデータを生成します。データ点は特になんの思想もなく100点にしました。

n <- 100
sim <- simGaussian(len = n)

simはlistで、状態変数(state)と、観測データ(data)が格納されています。

プロット

無駄に{ggplot2}で描きます。

# パッケージのインストール
if(!require(ggplot2)){
  install.packages("ggplot2", dependencies = TRUE)
}
require(ggplot2)

# data.frameに変換 
df <- data.frame(iteration = 1:length(sim$state), 
                   state = sim$state, 
                   data = sim$data)

# プロット
ggplot(data = df) + 
  geom_line(mapping = aes(x = iteration, y = state), 
            colour = "magenta") + 
  geom_point(mapping = aes(x = iteration, y = data), 
             shape = 21, size = 3, colour = "blue")

観測データから状態を推定

ここでは、逆に与えられた観測データから状態変数を推定します。{RcppSMC}では、blockpfGaussianOpt関数を使用して行います。plot = TRUEとすると、自動的に図を出力してくれます。ここで、粒子数は1000としています。

blockpfGaussianOpt(sim$data, 
                   particle = 1000, 
                   lag = 5, 
                   plot = TRUE)

推定値を取り出してもとのシミュレーションデータと比較

# 推定値などを格納
res <- blockpfGaussianOpt(sim$data, 
                          particle = 1000, 
                          lag = 5, 
                          plot = FALSE)

# 推定値を取り出す
estimate <- apply(res$values, 2, mean)

# 前に作っていたdata.frameに追加
df <- cbind(df, estimate)

# プロット
ggplot(data = df) + 
  geom_line(mapping = aes(x = iteration, y = state), 
            colour = "magenta") + 
  geom_line(mapping = aes(x = iteration, y = estimate), 
            colour = "green") + 
  geom_point(mapping = aes(x = iteration, y = data), 
             shape = 21, size = 3, colour = "blue")

緑の線が推定値で、そこそこ推定できているとも言えなくもない結果がでてきました。

おわりに

他にもpfLineartBS、pfNonlinBSといった関数が用意されていますが、マニュアルに載っている論文の結果を再現するだけです。他のパッケージとしては、{SMC}もありますが、こちらはパッケージが2011年から更新されていないようです。こちらも後日時間があったらテストしてみたいと思っています。