1. Home
  2. 技術者ブログ
  3. データ解析のための統計モデリング入門 マルコフ連鎖モンテカルロ(MCMC)法とベイズ統計モデル 読書メモ2

データ解析のための統計モデリング入門 マルコフ連鎖モンテカルロ(MCMC)法とベイズ統計モデル 読書メモ2

2017年 06月 23日

このブログ記事は『データ解析のための統計モデリング入門』（久保拓弥　著、岩波書店）という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。
今回は第８章、マルコフ連鎖モンテカルロ(MCMC)法とベイズ統計モデルについてのまとめの二回目です。

この章ではMCMCやメトロポリス法について説明がされています。
なので、実際にメトロポリス法を動かしてみるコードを用意しました。
コードはRで書きました。

dataSize <- 20
N <- 8

data <- rbinom(dataSize, N, 0.45)
summary(data)
cat("Analytically, q =", sum(data) / (dataSize * N), "\n\n")

likelihood <- function (q) {
  sum(log(dbinom(data, N, q)))
}

develop <- function(q) {
  if (runif(1, 0, 1) < 0.5) {
    q + 0.01
  } else {
    q - 0.01
  }
}

q <- round(0.3 + runif(1, 0, 0.3), 2)
ql <- likelihood(q)
i <- 0
cat("start q =", q, "\n")
while(i < 30) {
  p <- develop(q)   
  pl <- likelihood(p)
  r <- log(runif(1, 0, 1))
  if (ql < pl || r < pl - ql) {
    q <- p
    ql <- pl
  }
  cat("q =", q, "likelihood", ql, "\n")
  i <- i + 1
}
cat("last q =", q, "\n")

サンプルデータとしては、本の中でやっているのと同じように、事象の発生確率が 0.45 で、事象の試行回数は 8 の二項分布から、20 個のサンプルを取って用意しました。
このデータに対して、事象の発生確率 q が分かっていないものとして推定していきます。

メトロポリス法ではランダムに q を変化させた後、尤度が低くなる場合でも、変化の前後で尤度の比を計算して、その確率で q を更新します。
今の場合、計算の途中で使っているのは対数尤度なので、比ではなくて差になります。
コードで言うと以下の部分です。

  r <- log(runif(1, 0, 1))
  if (ql < pl || r < pl - ql) {
    q <- p
    ql <- pl
  }

コードを実行すると次のような出力が得られます。

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    2.0     3.0     4.0     3.7     4.0     6.0
Analytically, q = 0.4625

start q = 0.31
q = 0.31 likelihood -40.44572
q = 0.32 likelihood -39.35181
q = 0.32 likelihood -39.35181
q = 0.32 likelihood -39.35181
q = 0.33 likelihood -38.34881
q = 0.33 likelihood -38.34881
q = 0.34 likelihood -37.43294
q = 0.35 likelihood -36.60087
q = 0.36 likelihood -35.84958
q = 0.37 likelihood -35.17642
q = 0.38 likelihood -34.579
q = 0.39 likelihood -34.05522
q = 0.4 likelihood -33.60322
q = 0.39 likelihood -34.05522
q = 0.4 likelihood -33.60322
q = 0.41 likelihood -33.22138
q = 0.42 likelihood -32.90828
q = 0.42 likelihood -32.90828
q = 0.43 likelihood -32.66271
q = 0.43 likelihood -32.66271
q = 0.44 likelihood -32.48365
q = 0.43 likelihood -32.66271
q = 0.44 likelihood -32.48365
q = 0.43 likelihood -32.66271
q = 0.44 likelihood -32.48365
q = 0.45 likelihood -32.37025
q = 0.46 likelihood -32.32184
q = 0.45 likelihood -32.37025
q = 0.44 likelihood -32.48365
q = 0.43 likelihood -32.66271
last q = 0.43

この例の場合は、最尤推定値を解析的に求めることができるので、まず最初にその値を出力しています。
メトロポリス法は確率的に尤度が低くなる変化もするので、q は一直線に最尤推定値に向かうのではなく、時々戻ったりするのが分かります。

< 前の記事へ 次の記事へ >

＃人気の記事

＃タグ

＃アーカイブ