新規作成日:2015年12月6日
最終更新日:2016年9月22日

Stanを使って、ローカルレベルモデルを推定しましょう。
今回はナイル川の流量データを対象として、ローカルレベルモデルを推定します。

なお、Stanとベイズ推定の基礎に関しては、以下の記事をご覧ください
ベイズ統計学基礎
ベイズと統計モデルの関係
ベイズとMCMCと統計モデルの関係
Stanによるベイズ推定の基礎

状態空間モデルの基礎に関しては、以下の記事をご覧ください
状態空間モデルとは
dlmの使い方(Stanを使わないで状態空間モデルを推定する方法について書いてあります)
ローカルレベルモデル
季節とトレンド 

コードをまとめたものはこちらにおいておきます。コピペする際はこちらをお使いください。



スポンサードリンク

 

目次

1.ローカルレベルモデルとは
2.状態空間モデルの「左端」の問題
3.Stanによるローカルレベルモデルの推定
4.推定結果の取得と利用方法
おまけ1.ggplot2を使ったきれいなグラフを作る
おまけ2.乱数の一覧を取得する方法
おまけ3.結果の保存

 

1.ローカルレベルモデルとは

ローカルレベルモデルについては、それだけで一つの記事を書いていますので、こちらをご覧ください。
ここでは簡単な紹介をするにとどめます。

状態空間モデルはその名の通り「状態」を推定します。
私たちが手に入れたデータは「観測」結果として「状態」とは明確に区別します。

状態の変化の仕方を表す部分を「状態方程式」と呼びます。
状態から観測データが得られるプロセスを表す部分を「観測方程式」と呼びます。

ローカルレベルモデルは、「状態」の変化の仕方をランダムウォークだと仮定します。
「観測」の結果は、「状態」に正規分布に従うノイズが加わったものだと仮定します。
なので、ローカルレベルモデルは、別名「ランダムウォーク・プラス・ノイズ モデル」と呼ばれます。
ランダムウォークに従って遷移する「状態」にノイズが加わって「観測」データが得られると仮定しているからです。

ランダムウォークとは、「一期前の値からスタートしてノイズを付け加えていった結果」だと思っていただければ結構です。
一回のノイズはあまり大きな値をとらないかもしれません。でも、ノイズが累積されていくと、大きな変化を生み出す時があります。
例えばノイズが「1,2,1,-1,3」だとしましょう。
「一期前の値からスタートしてノイズを付け加えていった結果」は以下のようになります


1+2 = 3
     ↓
     3+1 = 4
          ↓
          4-1 = 3
               ↓
               3+3 = 6
というわけで、ノイズの累積和は「1,3,4,3,6」ですね。
一つ一つのノイズを見ると最大値は3なのですが、累積和をとるとどんどん大きくなっていくこともあります。
もちろん、負のノイズが累積されると、どんどん小さな値になっていくこともあり得ます。

私たちの目に見えない「状態」にランダムウォークを仮定し、その「状態」にノイズが付け加わったものが「観測」データです。
「状態」に付け加わるノイズのことを「過程誤差」と呼びます。
「観測」が得られる際に入るノイズのことを「観測誤差」と呼びます。
今回は、「過程誤差」も「観測誤差」も平均値が0の正規分布に従うホワイトノイズだと仮定します。
よって、「過程誤差」と「観測誤差」の分散の値が推定できれば、ローカルレベルモデルが推定できることになります。

 

2.状態空間モデルの「左端」の問題

状態空間モデルを推定するにあたって、悩ましいのが「左端」の存在です。「左端」は私の造語ですが、ほかの言い方が難しいですのでそのまま使っています。よい用語があれば教えてください。

ランダムウォークとは、「一期前の値からスタートしてノイズを付け加えていった結果」です。
今回はナイル川の流量データをモデリングします。
ナイル川の流量データは1871年から1970年まで、100年間分のデータがあります。
ランダムウォークしていると仮定する以上、1970年の「状態」は、1969年の「状態」に過程誤差が加わって生成されますね。
1969年の「状態」は1968年の「状態」に過程誤差が加わって生成されます。
それをずーっと戻していきます。
1872年の「状態」は1871年の「状態」に過程誤差が加わって生成されます。
では、初年度である1871年の「状態」は、いったいどこから推定すればよいのでしょう。1870年のデータは存在しないのに。

1870年の「観測」データは存在しないんですけど、1870年の「状態」は計算できなきゃ困ります。
というわけで、1870年の「状態」をこの記事では「左端」と呼んで、別途推定することにします。
この「左端」もStanを使えばちゃんとベイズ推定できますので心配はいりません。
とはいえ、コードを理解するためには、「左端」の存在を理解しなければいけないですので、ここで補足しておきました。

 

3.Stanによるローカルレベルモデルの推定

Stanを使って、ローカルレベルモデルを推定してみます。
コードをまとめたものはこちらにおいておきます。コピペする際はこちらをお使いください。

Stanの基礎やインストールの方法についてはこちらを参照してください。

以下の5つのステップを踏んで、モデルを推定します。
ステップ1:Stanを使用する準備
ステップ2:統計モデルの記述
ステップ3:データの指定
ステップ4:ベイズ推定の実行
ステップ5:うまく推定されているかをチェック~修正

順にみていきます。

ステップ1:Stanを使用する準備

Stanがインストールできていれば、ライブラリを読み込むだけでOKです。
計算の並列化も試してみます。

ステップ2:統計モデルの記述

ローカルレベルモデルは、以下のように記述します。

上記のコードをもう少し詳しく見ていきます。

dataブロック

データのブロックについては、特に複雑な指定はありません。
サンプルサイズと、ナイル川の流量データを指定しているだけです。
「intってなんやねん」という方はこちらをご覧ください。

parameters

ローカルレベルモデルの場合、推定するのは大きく2つあります。
1つは、「状態」の値そのもの。
1つは、過程誤差や観測誤差の分散の大きさです。

2,3行目が「状態」を表すパラメタです。
muZeroが「左端」すなわち1871年からしかないナイル川流量データの1870年時点の「状態」を表すパラメタです。
次のmuが100年間分の「状態」の推定値です。

4,5行目が過程誤差や観測誤差の分散です。
sigmaVが観測誤差、sigmaWが過程誤差の分散です。
VやWといった命名はdlmパッケージに合わせました。

model

modelブロックは少々複雑なので、上から順を追って解説します。

まずはナイル川の流量の初年度のデータ、1871年の話から始めましょう。
最初の状態、すなわち1番目のmuは、1870年の「状態」であるmuZero(「左端」のことですね)に過程誤差が加わって生成されます。
よって、4行目では「期待値muZero、分散sigmaWに従う正規分布」に従って1871年の状態(mu[1])が得られると指定しています。

「分散sigmaWの過程誤差が加わる」動作が「期待値muZero、分散sigmaWの正規分布に従う乱数を発生させる」ことと同じだということはイメージがわくでしょうか。
過程誤差は、期待値0のホワイトノイズです。
加えられる過程誤差の期待値が0なので、「過程誤差が加わった結果の期待値」は、過程誤差が加わる前の値(muZero)と変わりがありません。
ということで、期待値に1年前の状態であるmuZeroが入るのでした。

8行目に書かれた次の状態の遷移は、初年度(1871年)の状態推定においては無視します。
だって、muZeroから状態がすでに遷移しているから。
forループの中が(i in 2:n)と「2」からスタートしていることに注意してください。

そして13行目に移り、初年度の「状態」であるmu[1]に観測誤差が加わって、初年度の「観測」であるNile[1]が得られます。
これも、観測誤差の期待値が0であるため、「期待値mu[1]、分散sigmaVに従う正規分布」に従って「観測」データが得られると記述してあることに注意してください。

次は2年目です。
2年目の「状態」の推定値mu[2]は、1年前の状態の推定値mu[1]に過程誤差が加わって生成されます。
これは、8行目に記述されています。
「mu[i] ~ normal(mu[i-1], sqrt(sigmaW));」により、「 i 年目の状態は、期待値が『 i-1 年目の状態』で分散が sigmaW の正規分布に従う」ことを指定しています。1年前の状態から当年の状態へ遷移しています。これが2年目以降ずっと続きます。

「観測」については、初年度と変わらず、「 i 年目の観測データは、期待値が『同じ年(i 年)の状態』で分散が sigmaV の正規分布」に従います。

これを100年間ループ(サンプルサイズnには100を指定する予定です)すれば、すべての年で「状態」が推定できます。

ステップ3:データの指定

データには、Rにもともと入っているナイル川流量データを使います。

これをStanに渡すために整形します。

ステップ4:ベイズ推定の実行

stan関数を使って、ベイズ推定を実行します。

model_codeはStep2で作成されたモデル式を、dataにはStep3で整形されたナイル川流量データを指定します。
1回のシミュレーションの繰り返し数(iter)は2500とし、最初の1000個のシミュレーション結果は使用せずに捨てます(warmup=1000)。
間引きはせず(thin=1)、2500回のシミュレーションを3回繰り返します(chains=3)。

ステップ5:うまく推定されているかをチェック~修正

推定されたパラメタを表示するには、以下のコードを実行します。

今回は状態を100年間も推定したので、100以上パラメタがあります。少々長いので、ここでは省略します。

過程誤差と観測誤差の分散が正しく推定されているかを確認します。

引数「pars」を指定すると、特定のパラメタだけを表示させることができます。
推定結果が多い時はここを指定しないと、ほしい結果が出てきません。

Stanによるローカルレベルモデル_あまりよくない結果

グラフを見ると、難しいところですがsigmaWがあまりきれいではありませんね。
きれいかどうかの判別って難しいんですが、各々のChainがまんべんなく混ざり合っている状態が理想です。
Rhatは1.1以下であり、そこまで問題はないかと思いますが、念のためシミュレーションをやり直します。

繰り返し数(iter)を増やしたうえで、シミュレーション結果の間引きを行いました。
繰り返し数を増やすと計算に時間がかかるのが難ですが、実務の際は解析にかけられる時間などスケジュールと相談して実行してください。
以前にも書きましたが、ベイズ推定を行う際は短いスケジュールでやるべきではないと思います。
(ローカルレベルモデルくらいなら数分で終わりますが)

過程誤差と観測誤差の分散をもう一度確認します。

Stanによるローカルレベルモデル_修正結果

こっちのほうが、各chainごとにきれいに混ざり合っていてよさそうです。
ちなみにRhatはぴったり1になっています(無理にぴったり1にする必要はないのですが)。

長くて恐縮なのですが、こちらは結果をすべて表示させておきます。

 

4.推定結果の取得と利用方法

せっかく状態が推定できたのですから、それを取り出したり図示したりしたくなります。
しかし、stan関数で得られた結果から直接中身を取り出すのは若干面倒です。
少々変則的ですが、いったんsummary関数を介してから結果を取り出すと楽でした。

まず、推定結果にsummary関数を直接適用するなら、以下のコードを実行します。

summary(Nile_LocalLevelModel_2)

この結果は、各chainごとの結果と3つのchainをまとめた結果が各々出力されます。
3つのchainをまとめた結果だけを抽出するには、以下のコードを実行します。

summary(Nile_LocalLevelModel_2)$summary

最初の6行だけ取り出してみます。
(右端が切れているときはコードの中で右スクロールしてください。隠れているだけで載っています。)

ここまでくれば、例えば1列目のEAP推定量一覧を取り出すことは難しくありません。
やってみます。

鍵かっこを使う際は[行番号, 列番号]と指定すれば、ほしい結果だけ抽出できます。
「状態」であるmuだけを取り出しましょう。muは2行目から101行目までを占めていることに注意します。
(下のコードの2行目で指定しています)

EAP推定量と95%信用区間を取り出して保管します。

これで、推定結果の抽出ができました。

次は、図示してみましょう。
パラメタを見るだけだと結果がよく理解できません。推定結果はぜひ、実データと合わせて図示されることをお勧めします。これを見て「明らかにヘン」な結果になっているとすれば、モデルを修正することを検討します。

図示にはいろいろな方法がありますが、まずは、もともとのナイル川流量データと合わせて、推定値を時系列データに変換させるやり方を使ってみます。

開始年や終了年はナイル川の流量データに合わせました。

後はplotするだけです。

Stanによるローカルレベルモデル_ローカルレベルモデル図示
(画像はクリックすると拡大します)

パッと見問題なさそうです。
これで、ひとまず解析完了でしょうか。
お疲れ様でした。

 

おまけ1.ggplot2を使ったきれいなグラフを作る

おまけその1です。
さっきのやり方で図示はできるのですが、論文に出したり本に書いたりするにはちょっとカッコ悪い図ですね。
もっときれいなグラフを描くには「ggplot2」という外部パッケージを使うのが便利です。

まずはRを管理者として実行したうえで「install.packages("ggplot2")」を実行します。
これでパッケージをインストールした後で、以下を実行してください。

これで準備完了。
きれいなグラフを作りましょう。

ggplot2を使う最初のステップは、「描画したいデータをまとめる」ことです。
描画にしか使わないので、雑な名前を付けておきました。
d1にはもともとのナイル川の流量データが、d2には、推定された「状態」が入っています。これは時系列データ形式ではなくデータフレーム形式(data.frame)にしておきます。

次は描画する作業に移ります。
描画は、「グラフを重ねていく」イメージで作ります

NileGraphという名前でグラフ描画に必要な情報を格納します。

まずは、グラフの外側を作ります。
次に、さっき作ったNileGraphにデータ点を追加します。

11行目の「NileGraph <- NileGraph + ~~~」の部分が大事です。 2~8行目作ったグラフの情報NileGraphにさらに情報を追加していくのです。 今回はデータ点を付け加えるためにgeom_point()を足してやりました。次に、Stanで推定された「状態」を追加します。

「stat="identity"」を指定しないと、勝手に別の方法で平滑化された結果がプロットされるので気を付けてください。 ほかはデータの指定や線の太さの指定(lwd=1.2)です。グラフの情報を作ったので、これを描画します。


Stanによるローカルレベルモデル_ggplot2で図示

きれいなグラフができました。

 

おまけ2.乱数の一覧を取得する方法

今まではsummaryの結果を使って、パラメタのEAP推定量や95%信用区間を求めていました。
しかし、シミュレーション結果はそのまま残っていますので、手作業で求めることも可能です。

まずは、シミュレーション結果の一覧を取得します。extract関数を使います。

このシミュレーション結果は何個あるでしょうか。
例えば「左端」すなわち1870年の「状態」のシミュレーション結果の個数を出します。

この8700が出る計算式はこちらです。
(1回のシミュレーション回数 - 捨てられたシミュレーションの個数)÷ 間引き数 × chain数
=((iter - warmup) /thin)*chain

ということで確認してみます。

8700個あるシミュレーション結果を使って平均値や95%区間を求めても、summary関数と同じ結果が出てきます。

 

おまけ3.結果の保存

最後に、解析手法とは関係ないのですが、実務に使う際に覚えておくとよいだろうことを書いておきます。

MCMCは計算に時間がかかります。
とても簡単なローカルレベルモデルでさえ数分かかります。
これがポアソン分布になりランダムエフェクトが入り、外生要因が入り、季節変動が入り……と複雑になっていくと、計算にものすごく時間がかかることになります。
その計算を何度も繰り返したくないですね。

というわけで、MCMCした結果は、確実に保存するようにしましょう。
そして、その推定結果を二次解析するというのが効率の良いやり方かと思います。

結果の保存は、Rのコンソール画面から行います。
エディタではなく、確実に「コンソール」すなわち、計算結果がずらずら表示される画面をクリックしてから行ってください。
コンソールをクリックしてから、上のほうにあるメニューの「ファイル」~「作業スペースの保存」を選択します。
例えば「localLevelModel.RData」という名前で保存します。
今回の計算例では大体20Mバイトありました。

次回以降は、いちいちstan関数を回さなくてもすぐに結果が手に入ります。
先ほど保存した「localLevelModel.RData」をダブルクリックして開きます。
そして、ライブラリの読み込みだけやっておきます。

ちなみに、library関数とrequire関数は平文で使う場合は対して違いがないので、気にしないでください。

あとは、推定された結果をそのまま使うことができます。

stanを回すやり方と、二次解析のやり方を覚えておき、作業を分けられるようになればよいかと思います。
今回の記事では「4.推定結果の取得と利用方法」の章からが二次解析になります。
二次解析に関しては、実行に時間がかからないので、保存する必要はないでしょう。コードだけ残しておけば十分です。

 

参考文献


岩波データサイエンス Vol.1
特集「ベイズ推論とMCMCのフリーソフト」


 
何度も紹介している本ですね。
Stanを使って状態空間モデルを推定する方法が書かれています。
基礎を押さえた後は、こちらの本を読まれて、どんどん応用されていくのがよいかと思います。
 


状態空間時系列分析入門

 
状態空間モデルの基礎を押さえたい方にお勧めの本です。
一番簡単な状態空間モデルの入門書です。
 

リンク

この記事はベイズ推定を応用して状態空間モデルを推定する一連の記事の一つです。
記事の一覧とそのリンクは以下の通りです。
ベイズ統計学基礎
ベイズと統計モデルの関係
ベイズとMCMCと統計モデルの関係
Stanによるベイズ推定の基礎
Stanで推定するローカルレベルモデル



スポンサードリンク