※現在、ブログ記事を移行中のため一部表示が崩れる場合がございます。
順次修正対応にあたっておりますので何卒ご了承いただけますよう、お願い致します。
2017年 06月 23日
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` は一直線に最尤推定値に向かうのではなく、時々戻ったりするのが分かります。