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

ill-identified diary

所属組織の業務や見解とは一切無関係なアフィリエイト付きメモ帳。所属とは関係ないけどここを見て所属先にも興味を持っていただけると喜びます。

[STAN] [R] STAN の出力加工方法2 DIC の計算

概要

Spiegelhalter et al. (2002) で提案された DIC, デビアンス情報量規準はベイズ統計でモデルの選択に用いられる指標である. 詳しい説明は元論文に任せて, この記事では簡単な説明だけに留める. Gelman et al. (2013) の Ch. 7, takehiko-i-hayashi.hatenablog.com, あるいは 小西 (2008) にも少しだけ言及がある*1.

2016/8/18 不勉強だったので知らなかったのだが, 最近は DIC は理論的背景が不適切で好ましくないという解釈になっているらしく, ベイズ統計など複雑になりがちなモデルに対しては WAIC を使うのがよいらしい. もう少し自分なりに理解してから最近の情報量規準について続編を書く予定.

DIC の考え方

DIC とは, によって提案された情報量規準の一種である. AIC など情報量規準には罰則項としてパラメータの数を計算に含むが, 階層ベイズモデルのようなベイジアンの複雑なモデル, たとえば

\begin{align}
y_{i}\sim & \mathcal{N}(\theta_{i},\,\tau_{i}^{-1}),\\
\theta_{i}\sim & \mathcal{N}(\psi,\,\lambda^{-1})\end{align}

のようなモデルでは, パラメータの数をどう数えれば良いのかが不明瞭である. 1つめの式から  \theta_{i},\,\tau_{i} の数がパラメータ数なのか,

\begin{align}
y_{i}|\theta_{i}\sim & \mathcal{N}(\psi,\,\tau_{i}^{-1}+\lambda^{-1})\end{align}

というモデルと解釈するかでパラメータの数が変わってしまう. これが DIC という概念を生んだモチベーションの1つで, DIC はこれを解決するために「有効パラメータ数」という指標を含んでいる. 有効パラメータ数 (effective number of parameters)  p_{D}

\begin{align}
p_{D}= & \overline{D(\theta)}-D(\bar{\theta})\end{align}

で表される.  D(\cdot) デビアンスを表す.  \overline{D(\theta)} は (事後分布によって評価された) デビアンスの期待値 (平均デビアンス) であり,  D(\bar{\theta}) はパラメータの事後平均のデビアンスである. パラメータ数と言いながらあきらかに整数にはならないので, この指標は一見すると突飛に見えるが, 元論文では「伝統的な」統計学の理論を拡張して導いていることが主張されている. この有効パラメータ数を用いて DICは

\begin{align}
\mathrm{DIC}= & p_{D}+\overline{D(\theta)},\\
= & D(\bar{\theta})+2p_{D}\end{align}

と表される.

DIC の計算方法

MCMC の結果から計算する場合, MCサンプルとは事後分布から生成されたものであるから, 平均デビアン \overline{D(\theta)} は MCサンプル単位で計算したデビアンスの平均値となり, パラメータの事後平均のデビアン D(\bar{\theta}) はMCサンプルから先にパラメータの事後平均を求め, そのデビアンスを計算することになる. デビアンスはふつう, マイナス2掛ける対数尤度と定義されるから, MCサンプリングの試行回数を  S とすると,

\begin{align}
\overline{D(\theta)} & =-2\times\frac{1}{S}\sum_{i=1}^{S}\ln L(\theta_{i}),\\
D(\bar{\theta}) & =-2\times\ln L\left(\frac{1}{S}\sum_{i}\theta_{i}\right)\end{align}

と表せる. そして  \overline{D(\theta)},\,D(\bar{\theta}) が得られれば DIC もただちに計算できることになる. 前提としては, stan 側で generated quantities ブロックで対数尤度を計算する. デビアンスを直接計算するより対数尤度で出したほうが他の診断基準の計算にも使いまわせるからである. WAIC ではさらにデータの個体ごとに対数尤度を計算していたが, DIC の計算ではそれは不要で, その試行におけるサンプルをモデルに当てはめた場合の対数尤度あるいはデビアンスだけが必要になる. サンプリングと違って generated — での計算はずっと早いため, ここでの計算が増えてもあまり気にする必要はないため, 両方とも計算するようなプログラムでも問題はない.

前回書いたように, Stan の結果はサンプリングの試行単位で得られるので, Stan だけで DIC を計算することは出来ない. そこで R などで DIC の計算プログラムを書くことになる. R なら  \overline{D(\theta)} の計算は簡単だが, モデルをどう設定するかはケースバイケースなので,  D(\bar{\theta}) の計算はその都度自分で関数を用意する必要がある. ここでは stan 側で loglike という変数名で対数尤度 を計算したという前提なので下のような DIC() 関数を作成した. また, 事後平均のデビアンスは, 事後平均パラメータのベクトルと入力のデータフレームを引数に取る関数 dev を与えることで計算するという前提にした*2. R ならば dnorm() など確率密度を計算する関数が豊富なのでさほど複雑にはならないだろう. だいたいの場合はデビアンスを計算する関数は

function(parms, df.input){
    return(-2* sum(with(df.input, dnorm(y, param[1], param[2], log=T))))
}

のような感じになるだろう*3. これら d~~() 関数は log=T という引数を与えれば対数確率を返してくれるというのもポイントである. 計算結果は  \overline{D(\theta)} ,  D(\bar{\theta}) ,  p_{D} ,  \mathrm{DIC} の4要素を返す numeric 型のベクトルになる. プログラムは以下のようになる. 前回の stanfit をデータフレームにする処理を使いまわしているが、別々の関数に分けたほうが使いやすいかもしれない.


gistad42d075729bc80374d9d6dc0ba761e3



参考文献

Gelman, Andrew, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, and Donald B. Rubin. 2013. Bayesian Data Analysis. 3rd ed. CRC Press.

Bayesian Data Analysis, Third Edition (Chapman & Hall/CRC Texts in Statistical Science)

Bayesian Data Analysis, Third Edition (Chapman & Hall/CRC Texts in Statistical Science)

  • 作者: Andrew Gelman,John B. Carlin,Hal S. Stern,David B. Dunson,Aki Vehtari,Donald B. Rubin
  • 出版社/メーカー: Chapman and Hall/CRC
  • 発売日: 2013/11/05
  • メディア: ハードカバー
  • この商品を含むブログを見る

Spiegelhalter, David J, Nicola G Best, Bradley P Carlin, and Angelika van der Linde. 2002. “Bayesian Measures of Model Complexity and Fit.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 64 (4): 583–639. doi:http://dx.doi.org/10.1111/1467-9868.00353.

小⻄貞則・越智 義道・大森祐浩. 2008. 計算統計学の方法 -ブートストラップ・EMアルゴリズムMCMC-. 朝倉書店.

計算統計学の方法―ブートストラップ・EMアルゴリズム・MCMC (シリーズ予測と発見の科学 5)

計算統計学の方法―ブートストラップ・EMアルゴリズム・MCMC (シリーズ予測と発見の科学 5)

*1:このシリーズに情報量規準をあつかったものがあるが, 私は持っていない

*2:stan に与える場合は list 型にする場合が多いが, それ以前の加工の段階ではたいていの人は data.frame を使うから問題ないだろう

*3:マイナス2を掛けるのを忘れずに