読者です 読者をやめる 読者になる 読者になる

StatsFragments

Python, R, Rust, 統計, 機械学習とか

R で 状態空間モデル: 状態空間時系列分析入門を {rstan} で再現したい

前の記事でもリンクさせていただいているが、サイト 「状態空間時系列分析入門」をRで再現する では以下のテキストを {dlm}, {KFAS} で再現されており非常にありがたい。これらのパッケージの使い方については リンク先を読めば困らない感じだ。

自分も勉強のために似たことやりたい、、でも同じことやるのもなあ、、と考えた結果 同テキストの内容 {rstan} を使ってやってみた。

補足 Stan には状態空間表現用の関数 gaussian_dlm_obs ( 利用例 ) があるのだが、自分は使ったことがない。7章までのモデルは全て漸化式で表現されているため、それらを Stan のモデルとして記述した。

状態空間時系列分析入門

状態空間時系列分析入門

  • 作者: J.J.F.コマンダー,S.J.クープマン,Jacques J.F. Commandeur,Sime Jan Koopman,和合肇
  • 出版社/メーカー: シーエーピー出版
  • 発売日: 2008/09
  • メディア: 単行本
  • 購入: 2人 クリック: 4回
  • この商品を含むブログを見る

github リポジトリ

スクリプト / モデルは GitHub にあげた。現在 7 章までのモデルと kick する R スクリプト, markdown ファイルを置いている。

github.com

モデルの一覧

各章のモデルについて、その実行結果を RPubs にあげている。各モデルの概要とリンクを以下に記載する。自分のモデリング能力不足のため いくつか課題があると考えており、特に怪しいものには * をつけた。

かんたんな説明

必要パッケージ

まず、実行には以下のCRAN非公開パッケージが必要である。後者はどうでもよいが {pforeach} はいろいろ捗るのでインストールしていない方はただちに入れるべきだ。

データの読み込み

データは著者のサポートサイト、An Introduction to State Space Time Series Analysis から入手できる。各データは適宜 ts インスタンス化しておく。

Stan の実行

Stan の呼び出しは全て以下のように並列化して行う。

model_file <- 'モデルを記述したファイルのパス'

# モデルのコンパイルを実行
stan_fit <- stan(file = model_file, chains = 0)

# Stan を並列実行する。実行結果に sflist2stanfit を適用し、stanfit インスタンスを得る
fit <- pforeach(i = 1:4, .final = sflist2stanfit)({
  stan(fit = stan_fit, data = standata, chains = 1, seed = i)
})

また、モデルが収束したかどうかは 全ての Rhat が 1.1 未満となっているかどうかで確認する。そのままでは収束しない / しにくいモデルには、初期値 / 事前分布 / 制約を入れて収束させた。

is.converged <- function(stanfit) {
  summarized <- summary(stanfit)  
  all(summarized$summary[, 'Rhat'] < 1.1)
}

stopifnot(is.converged(fit))

また、テキストに最尤推定の結果が記載されている場合は、以下のようにして近い値になっているかを確認した。値がずれる場合は stopifnot を外し、差異を出力する。

# is.almost.fitted の定義は各リンク先を参照
stopifnot(is.almost.fitted(mu[[1]], 7.4157))
結果のプロット

テキストのグラフを再現する形で行う。使っているパッケージの説明はこちら。

課題

以下の2つについてはもう少し Stan に習熟できたら直したい。

複数の確率的状態をもつモデルをうまく収束させたい

このようなモデルは 収束しにくく、また 最尤推定した状態推定値 ( テキスト、{dlm}{KFAS}の結果 ) とずれやすい。特に 現状 制約を入れて収束させているモデルについて、事前分布 / 初期値 / 尤度関数の範囲でなんとかしたい。

モデルの一覧で * をつけたものがこれにあたる。

テキストの尤度関数を再現したい

テキストに記載の式 (概要は以下記事参照) に揃えようとすると、自分の記述では収束しにくくなってしまう。

sinhrks.hatenablog.com