概要

前回の(R) Stan の出力加工方法 - ill-identified diaryのおまけ的な形で書いた.
BUGS には DIC を計算する機能があるらしいが, rstan にはないので書いてみた.

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

2017/6/9 追記 不勉強だったのでよく理解していなかったのだが, DIC は計算に事後分布の平均を利用しているため、正則モデルに対してのみ有効である. ベイズ統計に特有の複雑な構造のモデルはしばしば正則モデルにならず, そのような場合は 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}$

と表される.

有効パラメータ数という名前は, AIC の罰則項に含まれるパラメータ数 $k$ の拡張であることに由来する. 線形モデルでは有効パラメータ数はAICの $k$ に漸近的に一致する. しかし, このような式で $k$ を近似できるというのもピンとこない. そこで, Gelman, et al. (2013) にかかれている簡単な例を参考に解説する. 今, $y\sim\mathcal{N}(\theta, 1)$ というモデルと $\theta\sim\mathcal{U}(0,\infty)$ という事前分布のベイズモデルを例に考える. もし観測値 $y$ の値がゼロに近いなら, 尤度と事前分布がそれぞれ半分づつ事後分布の説明に寄与しているといえる.
よって有効パラメータ数は1/2 になる. 一方で $y$ が非常に大きい正の値なら, 事前分布は事後分布に対して意味をなさないため, 尤度だけが事後分布を説明することになり, 有効パラメータ数は1になり, 1/2 の場合よりもモデルが好ましいという判定になる.

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)

作者: 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-. 朝倉書店.