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

基礎からのベイズ統計学 Skype読書会(13): 参加メモ

勉強会

この勉強会に参加させていただきました: https://github.com/7shi/study/wiki/bayes_hmc
読んでいる本(出典): 基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門 | 豊田 秀樹 | 本 | Amazon.co.jp
前回: 基礎からのベイズ統計学 Skype読書会(10): 参加メモ - クッキーの日記

読んだ範囲: 115~125ページ
11、12回を欠席してしまったので流れがよくわかっていないけど後から読み返してメモ。

  • ハミルトニアンモンテカルロ法112ページ
    • 事後分布にしたがう乱数がほしい。
    • 事後分布と、それとは独立な標準正規分布との同時分布にしたがう乱数でも別に構わない。
    • 事後分布と標準正規分布の同時分布を書き下すと以下のようになる(章末問題の 8 問目を見据えて、尤度が試行回数  n の二項分布、事前分布が母数  \alpha,\,\beta のベータ分布の場合にしてみる)。
         \displaystyle f(\theta, p|x)=\exp \biggl[\log \bigl( f(\theta|x) \bigr) +\log \bigl( f(p) \bigr) \biggr]
             \displaystyle \propto \exp\biggl[\log \bigl( \theta^x (1-\theta)^{n-x} \theta^{\alpha-1} (1-\theta)^{\beta-1} \bigr)+\log \bigl( \exp(-\frac{1}{2}p^2) \bigr) \biggr]
             \displaystyle = \exp\biggl[ (\alpha + x -1) \log \theta + (\beta+n-x-1)\log (1-\theta)-\frac{1}{2}p^2 \biggr]
             \displaystyle = \exp\bigl[ -\mathcal{H}(\theta, p) \bigr]
      最後の解釈は、
      位置-運動量空間母数-標準正規分布にしたがう確率変数空間
      位置 \theta母数 \theta
      運動量 p標準正規分布にしたがう確率変数 p
      ハミルトニアン同時分布の対数の-1倍
      ポテンシャルエネルギー事後分布の対数の-1倍
      位置エネルギー標準正規分布対数の-1倍
    • \alpha + x=10.2, \; \beta + n - x =5.8 のとき  f(\theta, p|x) を図示してみると以下。
      • 赤いところ: ハミルトニアンが大きくて、同時確率が低い。
      • 青いところ: ハミルトニアンが小さくて、同時確率が高い。52ページによると事後確率最大値は  \hat{\theta}_{\rm MAP} = 0.657 で、下の絵の青っぽいところもだいたいその辺にある気はする。
        f:id:cookie-box:20160626181743p:plain:w450
    •  (p_0, \, \theta_0)=(0, \, 0.1) からリープフロッグ法で時間を何ステップか進めてみる(章末問題 8, 9)。
      •  h(\theta) = - (\alpha + x -1) \log \theta - (\beta+n-x-1)\log (1-\theta) なので、 h'(\theta) = - (\alpha + x -1) / \theta + (\beta+n-x-1) / (1-\theta) を111ページの式にあてはめればよい。
      • 章末問題 9 で指示されている方の  \varepsilon=0.01 をつかって、 L=15 ステップ動いた図が以下。上図では横のプロット範囲が狭かったので少し広げた。 \varepsilon=0.05 では軌跡がもっとカクカクしていた。
        f:id:cookie-box:20160626185828p:plain:w450
    •  L ステップ動いたらその点の受容/棄却を判断する(ほぼ受容されるはず)。新しい点  (p_1, \, \theta_1) を受容したら、正規乱数  \tilde{p} を発生させて  (\tilde{p}, \, \theta_1) にワープする。そこからまた  L ステップ動く。繰り返し。

 \thetaEAP推定値まで出してないけどお腹すいたので終わる。スクリプトは以下。

library(fields)

p <- seq(-5.5, 5.5, length.out=100)
theta <- seq(0.01, 0.99, length.out=100) # log(0) を避ける
H <- c()

for (i in 1:length(p)) {
  for (j in 1:length(theta)) {
    H <- c(H, 9.2 * log(theta[[j]]) + 4.8 * log(1 - theta[[j]]) - 0.5 * p[[i]] * p[[i]])
  }
}
H <- exp(H)

mycolors <- c(rgb(1,(84:256)/256,(84:256)/256), rgb((255:84)/256,(255:84)/256, 1)) # 赤~青グラデーション
image.plot(p, theta, matrix(H, ncol=100, byrow=TRUE), col=mycolors)
abline(v=c(-5.55), h=c(0.995)) # image.plot すると左辺と上辺の黒線が見えなくなるので描き足す苦肉の策

### リープフロッグ法

L <- 15
epsilon <- 0.01
p <- 0
theta <- 0.1
p.trace <- c(p)
theta.trace <- c(theta)

for (i in 2:L) {
  p <- p - 0.5 * epsilon * (-9.2/theta + 4.8/(1-theta))
  theta <- theta + epsilon * p
  p <- p - 0.5 * epsilon * (-9.2/theta + 4.8/(1-theta))
  p.trace <- c(p.trace, p)
  theta.trace <- c(theta.trace, theta)
}

#図示
points(p.trace[1], theta.trace[1], cex=1.5, pch=19)
lines(p.trace, theta.trace, lwd=2)
points(p.trace[L], theta.trace[L], cex=1.5, pch=17)