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

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

読んだ範囲： 126～132ページ

126ページのカタログ刷新問題をやる。
正規分布の場合の事後分布、ポテンシャルエネルギー、その母数での偏微分は以下のようになるはず。

　 $\displaystyle f(\mu, \sigma^2, p_1, p_2|x)=\exp \biggl[\log \bigl( f(\mu, \sigma^2|x) \bigr) +\log \bigl( f(p_1, p_2) \bigr) \biggr]$
　　　　　　　　　　 $\displaystyle \propto \exp \biggl[\log \prod_{i=1}^n \biggl( \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \Bigl( - \frac{(x_i - \mu)^2}{2 \sigma^2} \Bigr) \biggr) + \log \biggl( \exp \bigl( -\frac{1}{2}(p_1^2 + p_2^2) \bigr) \biggr) \biggr]$
　　　　　　　　　　 $\displaystyle \propto \exp \biggl[ -\frac{n}{2} \log \sigma^2 - \sum_{i=1}^n \frac{(x_i - \mu)^2}{2 \sigma^2} - \frac{1}{2}(p_1^2 + p_2^2) \biggr]$
　 $\displaystyle h(\mu, \sigma^2) = \frac{n}{2} \log \sigma^2 + \sum_{i=1}^n \frac{(x_i - \mu)^2}{2 \sigma^2}$
　 $\displaystyle \frac{\partial h(\mu, \sigma^2)}{\partial \mu} = - \sum_{i=1}^n \frac{(x_i - \mu)}{\sigma^2}$
　 $\displaystyle \frac{\partial h(\mu, \sigma^2)}{\partial \sigma^2} = \frac{n}{2\sigma^2} - \sum_{i=1}^n \frac{(x_i - \mu)^2}{2 \sigma^4}$

これで $\varepsilon=0.1, \; L=100, \; T=100000$ としてHMC法を適用すると $\mu$ と $\sigma^2$ のトレースは以下のようになる。
特に $\sigma^2$ の方収束しているのかが大いにあやしいけど目をつぶる。本当は付録2を読んで適切な $\varepsilon, \; L, \; T$ を決めないといけない。128ページには $T=11000$ として $T=1000$ までを破棄してと書いてあるけど $\varepsilon$ と $L$ の値がない。 $T=11000$ でも多少 $\varepsilon$ と $L$ を試したが以下の図よりもっと悲惨だったのでとりあえずサンプリング数を増やした。 $\varepsilon$ と $L$ の具体的な数値を教えてくれていないのは理論を理解して自分で決定するようにという教育的配慮なんだろうか。
f:id:cookie-box:20160703130039p:plain:w450

f:id:cookie-box:20160703130039p:plain:w450

それで求まった10万点のサンプルからEAP推定値を出すと128ページの表と近い。ただ事後標準偏差が小さい。
事後標準偏差が小さいなりに95%信頼区間もちょっと狭くなっている。
「カタログBの平均購入金額が従来のものより高い確率」も事後標準偏差が小さいなりに教科書よりちょっと高い。
「カタログBの平均購入金額が3000円より高い確率」も事後標準偏差が小さいなりに教科書よりちょっと低い。
「効果量が0.8より大きい確率」も事後標準偏差が小さいなりに教科書よりちょっと低い。
でもまあ雰囲気的には近いからいいや。

> mean(mu.trace)
[1] 2925.961
> sd(mu.trace)
[1] 158.6847
> quantile(mu.trace, probs=c(0.025, 0.975))
    2.5%    97.5% 
2611.585 3238.228 
> length(mu.trace[mu.trace > 2500])/length(mu.trace)
[1] 0.99646
> length(mu.trace[mu.trace > 3000])/length(mu.trace)
[1] 0.32493
> length(mu.trace[(mu.trace - 2500)/sqrt(sigma2.trace) > 0.8])/length(mu.trace)
[1] 0.13042

スクリプトは以下。収束がいまいちなのは $\varepsilon$ と $L$ の調整の問題だと思いたい。
ちなみに受容率は 100% だった。

# カタログBを使用した顧客20人の購入金額データ
x <- c(3060, 2840, 1780, 3280, 3550, 2450, 2200, 3070, 2100, 4100,
       3630, 3060, 3280, 1870, 2980, 3120, 2150, 3830, 4300, 1880)

p_1 <- rnorm(1)
p_2 <- rnorm(1)
mu  <- 2926.5      # mu の初期値： 標本平均にしておく。
sigma2 <- 572855.5 # sigma2 の初期値： 不偏分散にしておく。

p_1.trace <- c(p_1)
p_2.trace <- c(p_2)
mu.trace  <- c(mu)
sigma2.trace <- c(sigma2)

epsilon <- 0.1
L <- 100

# ポテンシャルエネルギーの mu 偏微分
p_h_p_mu <- function(mu, sigma2) {
  return(-sum(x - mu) / sigma2)
}
# ポテンシャルエネルギーの sigma2 偏微分
p_h_p_sigma2 <- function(mu, sigma2) {
  return(10.0 / sigma2 - 0.5 * sum((x - mu) * (x - mu))/(sigma2 * sigma2))
}
# ハミルトニアン
H <- function(p_1, p_2, mu, sigma2) {
  return(10.0 * log(sigma2) + 0.5 * sum((x - mu) * (x - mu))/sigma2 + 0.5 * (p_1 * p_1 + p_2 * p_2))
}

# 10万点サンプリング
for (i in 2:200000) { # 一応棄却されることを考えて多めにループ
  p_1_save <- p_1
  p_2_save <- p_2
  mu_save <- mu
  sigma2_save <- sigma2

  for (j in 1:L) { # リープフロッグ法で次の候補点を得る
    p_1 <- p_1 - 0.5 * epsilon * p_h_p_mu(mu, sigma2)
    p_2 <- p_2 - 0.5 * epsilon * p_h_p_sigma2(mu, sigma2)
    mu     <- mu + epsilon * p_1
    sigma2 <- sigma2 + epsilon * p_2
    p_1 <- p_1 - 0.5 * epsilon * p_h_p_mu(mu, sigma2)
    p_2 <- p_2 - 0.5 * epsilon * p_h_p_sigma2(mu, sigma2)
  }
  
  # 需要・棄却の判定
  r <- exp( H(p_1, p_2, mu, sigma2) - H(p_1_save, p_2_save, mu_save, sigma2_save) )

  if (r < runif(1)) { # 棄却
    p_1 <- p_1_save
    p_2 <- p_2_save
    mu <- mu_save
    sigma2 <- sigma2_save
    next
  }
  
  # 受容
  p_1.trace <- c(p_1.trace, p_1)
  p_2.trace <- c(p_2.trace, p_2)
  mu.trace  <- c(mu.trace, mu)
  sigma2.trace <- c(sigma2.trace, sigma2)
  
  if (length(p_1.trace) >= 100000) { # 10万点たまったら抜ける
    print(i) # ここでループのインデックスも 100000 なら1回も棄却なし
    break
  }
  
  # 次の探索の準備 -- 正規乱数をセットしておく
  p_1 <- rnorm(1)
  p_2 <- rnorm(1)
}