TIM Labs

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

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

この章では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)
  if (ql < pl) {
    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 に更新する、という方法です。 develop 関数の中で q0.01 足すか引くかして q を変化させています。 コードを実行すると次のような出力が得られます。

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
      1       3       4       4       5       6
Analyticaly, q = 0.5

start q = 0.35
q = 0.36 likelihood -41.69397
q = 0.37 likelihood -40.76192
q = 0.38 likelihood -39.90848
q = 0.39 likelihood -39.13129
q = 0.39 likelihood -39.13129
q = 0.39 likelihood -39.13129
q = 0.39 likelihood -39.13129
q = 0.4 likelihood -38.42821
q = 0.4 likelihood -38.42821
q = 0.41 likelihood -37.79737
q = 0.42 likelihood -37.23712
q = 0.42 likelihood -37.23712
q = 0.43 likelihood -36.74602
q = 0.44 likelihood -36.32282
q = 0.45 likelihood -35.96647
q = 0.45 likelihood -35.96647
q = 0.45 likelihood -35.96647
q = 0.46 likelihood -35.67609
q = 0.47 likelihood -35.45097
q = 0.48 likelihood -35.29055
q = 0.49 likelihood -35.19445
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
last q = 0.5

この例の場合は、最尤推定値を解析的に求めることができるので、まず最初にその値を出力しています。 q をランダムに変化させていくと、最終的にはその値に近づいていくのが分かります。

トラックバック(0)

トラックバックURL: http://labs.timedia.co.jp/mt/mt-tb.cgi/647

コメントする

このブログ記事について

このページは、sato.yuichiroが2017年6月20日 12:00に書いたブログ記事です。

ひとつ前のブログ記事は「ChainerがMicrosoft Azure, Windows に対応」です。

次のブログ記事は「並列に Deep Learning して、高速に賢くなれるか」です。

最近のコンテンツはインデックスページで見られます。過去に書かれたものはアーカイブのページで見られます。