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

六本木で働くデータサイエンティストのブログ

元祖「銀座で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木

Rでベイジアン動的線形モデルを学ぶ(2.5):最尤法でパラメータ推定してみる

前回サクッとローカルレベル・モデルを推定してみたわけですが、そう言えばパラメータ推定は何もしなかったのでした。既に線形モデルも一般化線形モデルもこのブログで見てきている以上最小二乗法や最尤法やMCMCでパラメータ推定するというのは常識なわけですが、今回テキストにしているPetris本とCommandeur本の2冊とも


Rによるベイジアン動的線形モデル (統計ライブラリー)

Rによるベイジアン動的線形モデル (統計ライブラリー)

状態空間時系列分析入門

状態空間時系列分析入門

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


何故か最尤法によるパラメータ推定の話題をかなり後の方に持ってきていてこのままだとそこまで到達するのに時間がかかるので*1、先にさっさと最尤法でのパラメータ推定を{dlm}でどうやってしまうのかを取り上げておくことにしておきます。


なお、今回の話題はPetris本のpp.147-150にも載っていますが、前回も沖本本の時も参考にさせていただいたLogics of Blueさんの記事も分かりやすくて良いと思うので是非是非どうぞ。


dlmMLE関数を使う


基本的にはdlmMLE関数に適当な初期値を与えてやることで最尤推定ができます。ということで、前回の記事までやったところの続きとして以下のようにやってみましょう。

> build1<-function(param){dlmModPoly(order=1,dV=exp(param[1]),dW=exp(param[2]))}
# 最尤推定したいモデルを先に関数定義の形で作っておく必要がある
# dlmModPolyで普通にランダムウォーク・プラス・ノイズモデルを指定する
# なお、基本的に対数尤度に直してから計算している都合があるので、
# パラメータは先にexpしておく必要がある。。。らしい

> fit<-dlmMLE(d2,c(1,1),build1)
# 初期値にdV=1, dW=1を入れたつもり
# build1関数を引数として与える

> d2.model2<-build1(fit$par)
# ここでは多段階最適化を行っていることになっているので、
# (※注:そうではないらしいです下記参照)
# 推定したfitの中のパラメータ部$parをbuild1関数に与えることで、
# 最尤推定値に基づくモデルd2.model2を改めて作る

> d2.filt2<-dlmFilter(d2,d2.model2)
# これをカルマンフィルタで推定する

> matplot(cbind(d2,d2.filt1$m[-1],d2.filt2$m[-1]),type='l',lwd=c(2,3,3),col=c('red','blue','purple'),lty=c(1,2,3))
# 後は前回の記事で推定したd2.filt1$mと併せてプロットするだけ

> par(mfrow=c(2,1))
> hist(d2-d2.filt1$m[-1],col='blue',main='Not MLE')
> hist(d2-d2.filt2$m[-1],col='red',main='MLE')
# 残差の分布も見てみる

f:id:TJO:20140927155406p:plain

f:id:TJO:20140927155909p:plain


ということで、最尤法でパラメータ推定した方が残差の分布もより正規性が強くなっていることが何となく見て取れます。


注意点


Logics of Blueさんもご指摘の通りで、一般の最尤法と同じく初期値にかなり依存するケースがあるので、推定値がおかしいと思ったら毎回初期値を変えてやり直す必要があります。そのためにも、プロットなどで可視化できる場合は可視化し、できない場合でも生の推定値を目視しておかしくないかどうか確認した方が良いでしょう。


追記1


@先生から以下の如くツッコミが。



ということでもうちょっと勉強を。。。


追記2


愚かにもCommandeur本の内容をRでやってみましたというPDF本があるじゃないですかー。

あーもうやる気なくなったわ。。。じゃなくて最低でもこれはなぞってやってみます。。。

*1:Petris本だと4章まで待つ必要がある上にCommandeur本では最尤法のサの字も出てこない