事後確率分布に従う乱数

森林学会のポスター発表で話していて、MCMC推定した事後確率分布に従うような乱数を得られないか、という話になった。で、MCMCサンプルから値を得てくればどうか、というような結論になったのだが、うまくいくのかちょっと試してみた。

ワーキングディレクトリにCODAchain1.txtとCODAindex.txtとがあるとする。

post <- read.openbugs("")
plot(post[,"sigma"], trace = FALSE)

こういう関数を定義して、

rpost <- function(n, samples) {
  ind <- floor(runif(n, min = 1, 
                     max = length(samples) + 0.999))
  samples[ind]
}

使ってみる。

samples <- post[[1]][,"sigma"]
hist(rpost(500, samples))

これを「乱数」というのは間違っているような気もするが、MCMCのサンプルサイズが十分に大きければそれなりに役に立つような気もする。

タグ：MCMC R