Metropolis-Hastings(MH)乱数の自己相関

  MHアルゴリズムによるN(0,1)からのサンプリングを行う。提案分布は一様分布とする。乱数の自己相関の状況について、特に採択率やthinningを変えた場合の影響を例示する。

１．Rコードの一例
thin <- 0 #サンプリング間隔（0のとき間引きなし）
thin1 <- thin+1
Ne <- 1000 # 有効サンプル数
N <- thin1*Ne # 総サンプル数
samp <- numeric(Ne)
count <- 0 # 採択数のカウント
x0 <- 4.5 # 初期値
a <- 4 # 採択率の調整
samp[1] <- x0
for (i in 1:N) {
cand <- runif(1, -a, a) # 提案分布（一様分布）から乱数生成
x1 <- x0 + cand # random walk
alpha <- min(1, dnorm(x1)/dnorm(x0)) # 採択確率の計算、提案分布は対称
u <- runif(1)
if (u < alpha){
x0 <- x1
count <- count + 1
}
if(i %% thin1 == 0){
j <- (i+thin1)/thin1
samp[j] <- x0
}
}
count/N #採択率
mean(samp)
sd(samp)
par(mfrow=c(2,1))
plot(samp,type="l")
hist(samp) # 密度関数を描くときは plot(density(samp))
par(mfrow=c(1,1))
par(ask=TRUE)
NN <- N/2
for(i in 1:NN){
j <- 2*i-1
plot(samp[j],samp[j+1],xlim=c(-5,5),ylim=c(-5,5),col="blue")
par(new=T)
}
par(new=F)
acf(samp) # 自己相関係数
# http://aoki2.si.gunma-u.ac.jp/R/acf.html には自己相関係数を計算する
# 関数acf2(x,k)のRコードが掲載されている。

２．採択率と調整パラメータの関係など（有効サンプル数=10000）

a	thin	採択率	サンプル平均	サンプル標準偏差
4	0	0.399	0.0192	0.9855
4	3	0.388	-0.0119	1.0055
4	6	0.389	-0.0037	1.0025
1	0	0.800	0.0806	0.9881
1	3	0.805	-0.0092	1.0058
1	6	0.802	0.0059	0.9966

３．結果の表示（有効サンプル数=1000)
・thin=0(a=4)

 

・thin=3(a=4)

・thin=10(a=4)

・thin=0(a=1)

・thin=3(a=1)

・thin=10(a=1)

採択率が0.3～0.4であれば、thin=3でも自己相関の低下にはかなり効果的と思われる。