StatModeling Memorandum

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

「The BUGS Book : A Practical Introduction to Bayesian Analysis」 David Lunn et al.

★★★★★の良書です。WinBUGS, OpenBUGSの作者らによる初の教科書です。登場遅すぎですよ。 非常によくまとまっており、久保先生の緑本の次に読むべき本と言えそうです。買いの一択です。

The BUGS Book (Chapman & Hall/CRC Texts in Statistical Science)

The BUGS Book (Chapman & Hall/CRC Texts in Statistical Science)

  • 作者:Lunn, David
  • 発売日: 2012/11/01
  • メディア: ペーパーバック

目次はここから。サポートページやサンプルコードはここから。

内容は理論と実践(BUGSコード)がほどよいバランスです。 何がいいってBUGSコードがよい。シンプルで習得しやすい。さすがは北斗神拳を作った人たちです。

まず2章のRepair問題で読者をうならせます。 問題「1回の修理費は平均=100円, 標準偏差=50円のガンマ分布に従うとする。予算1000円で何回修理できそうか?」に対するBUGSコードが以下(サンプルコードより引用)。

model {
   for (i in 1:20) {
      Y[i] ~ dgamma(4, 0.04)
   }
   cum[1] <- Y[1]
   for (i in 2:20) {
      cum[i] <- cum[i - 1] + Y[i]
   }
   for (i in 1:20) {
      cum.step[i] <- i*step(1000 - cum[i])
   }
   number <- ranked(cum.step[], 20)  # maximum number in cum.step
   check  <- equals(cum.step[20], 0) # always 1 if I=20 big enough
}

美しいです。解説の詳細は書籍参考。

11章のZIP modelもさすがです(サンプルコードより引用)。ZIP とは Zero-inflated Poisson の略でゼロがインフレ気味(過剰気味)のポアソン分布です。

model {
   for (i in 1:N) {
      y[i] ~ dpois(m[i])
      m[i] <- group[i] * mu   # mean is 0 if group is 0
      group[i] ~ dbern(p)
   }
   # proportion of claims that could be positive
   p ~ dunif(0,1)
   mu ~ dgamma(0.5, 0.0001) # approximate Jeffreys prior
}
list(N=35,y=c(0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,4,4,5))

その他、個人的に参考となった部分は8章、9章、10章です。

8章はモデルチェックの部分がよかったです。実データがモデルからどれぐらい乖離しているか見積もるために mid p-value というものを算出して判断しています。これは乖離しちゃ困るような統計量を(工夫して)作って、モデルから生成した擬似データと実データをstep関数で比較してp-valueを求めています。

また、WinBUGSの結果に出てくるpD(実効的なパラメータ数)やDICをどれほど重要視したらよいかの感覚も文面から読み取れます。

9章は各トピックは短いですが要点を抑えていて有用です。 例えば変数の範囲の制限のかけ方の部分です。theta1とtheta2がtheta1-theta2平面でドーナツ状に分布するような制限がある時は以下のコードで実現できます(サンプルコードより引用)。

model {
   theta[1]    ~ dnorm(0, 1)
   theta[2]    ~ dnorm(0, 1)
   z          <- 1
   z           ~ dbern(constraint)
   constraint <- step(theta[1]*theta[1] + theta[2]*theta[2] - 1)
}

10章はモデルチェック部分、収束を早める方法、cut関数の使い方など勉強になりました。

11章の時系列データの解析は紙面も少なく無難で微妙です。ARとARMAが軽く紹介されている程度です。 PK(化合物の血中濃度)の微分方程式をWinBUGSの拡張パッケージのode関数で解いている例はびっくりしました。内部的にルンゲクッタで数値積分しているようです。いつかフォローしたいと思います。