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

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

今日読んだ範囲: 94~104ページ
8、9回を欠席してしまったので流れがよくわかっていないけど雑多なメモ。

  • p値は、初心者からしばしば誤解されます(97ページ)。
    • p値は「p値が有意水準 α=0.05 より小さいとき、帰無仮説を棄却する」などのようにつかう。
    • 誤: p値は、帰無仮説が正しい確率。それが小さいときは、帰無仮説を棄却する。
    • 正: p値は、帰無仮説が正しいときにそのデータが観測される確率。p値が小さいときは「帰無仮説が本当に正しいとしたらそのようなデータが観測される確率は低い」ということなので帰無仮説を棄却する。
      • 例えば、コインを7回投げたら7回とも裏だったとする。
      • このとき「帰無仮説: このコインは裏が出る確率も表が出る確率も等しい」を検定する。
      • 本当に裏が出る確率も表が出る確率も等しいコインだったら、7回とも裏がでる確率は 0.0078125 で、1% にも満たないので、「偶然7回とも裏が出た」と考えるにはちょっと低い。
      • ちょっと低い、だと基準が曖昧なので、「帰無仮説が正しいとしたら起きる確率が5%以下になるくらい珍しいこと」を「偶然それが起きたと考えるには低い」と定義する。この帰無仮説に対しては、「表が出すぎ」側も「表が出なさすぎ」側もどちらも起こりにくいことなので、珍しいエリアは両端から 2.5% ずつになる。
      • 0.0078125 < 0.025 から「7回とも裏が出ることは5%も起きない珍しいこと」といえるので、帰無仮説を棄却する。
      • もし、α=0.01 だったら、0.0078125 > 0.005 なので、「7回とも裏が出ることは1%も起きない珍しいこと」とはいえず、帰無仮説を棄却できない。というより7回投げたくらいではそんな精度は出ない。α=0.01 の有意水準で検定したかったらもっと投げるべき。
  • 98ページ表4.4「夕食」は、船さんが夕食を用意しなければならない確率。表だけみるとびっくりする。
  • 5章はハミルトニアンの意味からはじめて、ハミルトニアンモンテカルロ法の説明らしい。ちゃんと読んでいないけど、MH法では次のサンプル点候補をどちらの方角に取るかは適当だったのが、HMC法ではサンプリング対象の空間にポテンシャルをしきつめてなんかそれをつかって方角も上手いこと決めるということなのか。というか採択率がよくなるのがいいところなのか。



4章が何だったのかのついでに章末問題をみておく。

1~2)

  • 解くのは省略。

3~5)

  • 解くのは省略。
  • この「前日がカレー/ラーメン/焼きそばだったときに今日がカレー/ラーメン/焼きそばになる遷移確率」は、詳細釣り合い条件を満たさないが定常分布に収束する例らしい(カレーのリピート率高い…)。
  • 詳細釣り合い条件はなんだったかというと(82ページ)、これをみたせば絶対定常分布に収束する。だからMHアルゴリズムの遷移確率はこれを満たすようにする。詳細釣り合い条件は目的に照らし合わせると強すぎる条件だけど、取り扱いが容易だからこれを満たすようにしている?という認識。
  • 「カレー/ラーメン/焼きそば」が詳細釣り合い条件を満たしていないのに定常分布に収束するのは、「カレーとラーメン」「ラーメンと焼きそば」のように2つずつの組で流出入が釣り合っているのではなくて、「カレーとラーメンと焼きそば」の3つ合わせて流出入が釣り合っているから、のはず。

5~8)
ガンマ分布から独立MH法でサンプリングする問題だけ。以下のスクリプトはおそらく何か駄目。

set.seed(1234)

# 提案分布
mn.prop <- 1.0
sd.prop <- sqrt(0.5)
GetPointFromDistProp <- function(){return(rnorm(mean=mn.prop, sd=sd.prop, 1))}
GetValueOfDistProp <- function(x){return(dnorm(mean=mn.prop, sd=sd.prop, x))}
# 目標分布
GetValueOfDistTarget <- function(x){return(dgamma(shape=11.0, scale=1.0/13.0, x))}

GetNextPoint <- function(theta) {
  # 提案分布から候補点を発生します
  a <- GetPointFromDistProp()
  # 候補点を受容するかどうか判定します
  denom <- GetValueOfDistProp(a) * GetValueOfDistTarget(theta)
  numer <- GetValueOfDistProp(theta) * GetValueOfDistTarget(a)
  if (denom > numer) {
    r <- numer / denom
    if (runif(1) < r) {
      return(a) # 受容
    } else {
      return(theta) # 留まる
    }
  } else {
    return(a) # 受容
  }
}

vec.x <- c()
x <- 0
vec.x <- c(vec.x, x)
for (trial in 2:100000) {
  x <- GetNextPoint(x)
  vec.x <- c(vec.x, x)
}
vec.x <- vec.x[1001:100000] # 最初の方は burn-in 期間として捨てる
print(mean(vec.x))
[1] 0.8477589

点数のわりに理論値 0.8461538 とずれる。解答例だとちゃんと 0.846 になっている。何か乱数の扱いが駄目なのかもしれない。