本文に掲載していないコードとその解説

以下，本文に掲載していないコードとその解説です．

コードは（本文掲載のものとテストデータを含めて）zipでまとめてあります． ここからダウンロード (2017/06/21版）

これらのコードと解説は，本書の一部ではなく，正式の記事と比較して，必ずしも十分吟味されていない内容も含まれます．利用は自己責任でお願いします．

誤りや問題点についてはtwitter の岩波DSアカウント@iwanamiDSまでご連絡頂ければ幸いです．

PF_smooth.R

固定ラグ平滑化のシンプルな実装例です．線形ガウスモデルについてはdlm及びKFASによる平滑化と結果を比較します．

平滑化の場合，データの時系列がジャンプする点より手前でも誤差が大きくなることがわかります．ラグLおよび粒子数を変えて実行してみると興味深いですが，乱数によるばらつきがかなりあることに注意する必要があるでしょう．

PF_res.R

DS6の解説の末尾注で述べた２つの事項（重みの計算の注意(logsumexp），整数部分を乱数を使わずに確保する手法）の実装例です．

PF_like.R, PF_like_n.R

尤度計算の実装例です．対数尤度を計算してKFASの結果と比較します．

PF_like_n.Rでは末尾注の「重みの計算の注意（logsumexp）」を実装しています．この数値例では結果は同一でした．

ここでは，観測ノイズを固定して，いろいろなシステムノイズの値について尤度を計算しています．

システムノイズが大きいときには，粒子フィルタで求めた値はKFASとよく一致します．僅かな系統誤差がありますが，この原因ははっきりしません（初期値の関係？見落とした定数項がある？）

システムノイズが小さいときには，粒子フィルタで求めた値はKFASと一致しませんが，これは本文で述べた「データが飛躍する部分で粒子が事後密度の大きい領域にほとんど到達できない」という現象のためではないかと思われます．この場合も，粒子の数をうんと増やせば正しい答えが得られるはずですが，なかなか収束しない場合もあります．いまの場合，もともと尤度が低い領域で起きているので，そこで計算が合わなくても問題はないという考え方もできるかもしれません．