Metropolis Hastings サンプリングする

# 値の台
xx <- seq(from=-5,to=20,length=1000)

# 適当な(密度関数)
dd <- dnorm(xx) + dnorm(xx,1,2) + dnorm(xx,10,2)

# 確率密度関数っぽい形
plot(xx,dd)

# 
library(mcmc)

# 確率密度に比例したあたいの対数を返す関数を書く
h <- function(x){
	log(dnorm(x) + dnorm(x,1,2) + dnorm(x,10,2))
}
# サンプリングする
# ２番目の引数は、始めに試す乱数の値
# 
out <- metrop(h, 0, 10000)
out$accept

# out[[2]]には試した値が格納される
# out$accept.batchには、0/1が格納され、採択されたものには1が立つのでそれを選ぶ

hist(out[[2]][which(out$accept.batch==1)])