「マルコフ連鎖モンテカルロ法 (統計ライブラリー)
この問題は本書にも実験結果が掲載されているので，自分が行った実験結果と見比べることができる．結論から言うと，一応大まかな結果としては同じ傾向を示したが，乱数を用いた実験ということもあって多少の値のズレがでている．また，パラメータによっては結果の値/分布が安定しないものもあるので，何度か実験を行なって結果を吟味している部分もあるが，結果を曲解するためのものではないということを付け加えておく．
なお，この文章ではメトロポリス・ヘイスティング・アルゴリズムを以降MHと表記する．
問題
ガンベル分布（二重指数分布）からの乱数生成を考える．今回はガンベル分布からの直接的な乱数生成ができないとしたうえで，MHによる乱数生成を行う．ガンベル分布は
実験条件
今回は，MHの中でも酔歩過程/ランダムウォーク（Random walk）と独立連鎖（Independent chain）の2種類の方法を試す．それぞれのMHにおける推移カーネル
酔歩過程：
独立連鎖：
となる．ここで，{0.25, 1, 5, 15, 150} の5パターンで実験を行った．
また，生成した乱数の系列において，系列初期の
結果：以下の図の見方
以下に示す図は，各行にチューニングパラメータの値，各列に系列のプロット，系列の分布，系列の平均の推移，系列の自己相関プロットを示したものである．図のタイトルに，チューニングパラメータの値，生成した乱数の平均と分散，そして採択率を示している（便宜的にタイトルの部分に情報を書き加えているだけで，タイトルと図は直接的には関連していない）．
結果：酔歩過程/ランダムウォーク
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
par( mfrow= c( 5 , 3 ))
for ( i in c( 0.25 , 1 , 5 , 15 , 150 )){
c.ran <- i^ 2
alpha <- numeric( 0 )
Niter <- 10000
x <- rnorm( 1 , mean= 0 , sd= c.ran)
for ( t in 2 :Niter){
z <- rnorm( 1 , mean= 0 , sd= c.ran)
alpha[ t- 1 ] <- min( 1 , exp( - ( 1 / 3 ) * ( abs( x[ t- 1 ] + z) - abs( x[ t- 1 ]))))
x[ t] <- x[ t- 1 ] + z* ( runif( 1 ) < alpha[ t- 1 ])
}
X <- x[ 1001 :Niter]
stat <- table( alpha== 1 )
plot( 1001 :Niter, X, type= "l" , main= paste( "c_ran=" , i), xlab= "Iterations" , ylab= "" )
plot( density( X), main= paste( "mean=" , round( mean( X), digit= 2 ), " var=" , round( var( X), digit= 2 )))
plot( cumsum( X) / 1 :length( X), type= "l" , main= paste( "accept ratio=" , round( 100 * stat[ 2 ] / stat[ 1 ]), "%" ))
abline( h= 0 , lwd= 2 , col= "red" )
}
結果：独立連鎖の場合
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
par( mfrow= c( 5 , 4 ))
for ( i in c( 0.25 , 1 , 5 , 15 , 150 )){
c.ind <- i
alpha <- numeric( 0 )
Niter <- 10000
x <- rnorm( 1 , mean= 0 , sd= c.ind)
for ( t in 2 :Niter){
z <- rnorm( 1 , mean= 0 , sd= c.ind)
alpha[ t- 1 ] <- min( 1 , exp( - ( 1 / 3 ) * ( abs( z) - abs( x[ t- 1 ]))
- ( 1 / 2 ) * (( x[ t- 1 ] / c.ind) ^ 2 - ( z/ c.ind) ^ 2 )
)
)
if ( runif( 1 ) < alpha[ t- 1 ]){
x[ t] <- z
} else{
x[ t] <- x[ t- 1 ]
}
}
X <- x[ 1001 :Niter]
stat <- table( alpha== 1 )
plot( 1001 :Niter, X, type= "l" , main= paste( "c_ind=" , i), xlab= "Iterations" , ylab= "" )
plot( density( X), main= paste( "mean=" , round( mean( X), digit= 2 ), " var=" , round( var( X), digit= 2 )))
plot( cumsum( X) / 1 :length( X), type= "l" , main= paste( "accept ratio=" , round( 100 * stat[ 2 ] / stat[ 1 ]), "%" ), ylab= "" , xlab= "Iterations" )
abline( h= 0 , lwd= 2 , col= "red" )
acf( X)
}
解釈
さて，このように色々とプロットしてみたわけだが，酔歩過程の場合も独立連鎖の場合もチューニングパラメータの値が0.25や150といった非常に小さい/大きい場合においては，乱数の分布や収束速度などに問題があることがわかる．
酔歩過程においては，マルコフ連鎖の候補点の生成は少しずつ値を移動させていくという点で，分布の表面をなぞる形となる．そのため，チューニングパラメータの値である分散が小さい場合では候補点の移動範囲が狭いため，分布全体を探索するのに非常に時間がかかる．
一方，独立連鎖においては，マルコフ連鎖の候補点の生成は現在の状態に依存しない．そのため，酔歩過程のように値を移動させていくということは見られないが，一方で採択確率の部分でチューニングパラメータの値が影響していることがわかる．
参考
」PP.11-14のメトロポリス・ヘイスティング・アルゴリズムによる二重指数分布からの乱数生成（サンプリング）をやってみた．
