「マルコフ連鎖モンテカルロ法」1章のメトロポリス・ヘイスティング・アルゴリズムによるガンベル分布の乱数生成

「マルコフ連鎖モンテカルロ法 (統計ライブラリー)」PP.11-14のメトロポリス・ヘイスティング・アルゴリズムによる二重指数分布からの乱数生成（サンプリング）をやってみた．

この問題は本書にも実験結果が掲載されているので，自分が行った実験結果と見比べることができる．結論から言うと，一応大まかな結果としては同じ傾向を示したが，乱数を用いた実験ということもあって多少の値のズレがでている．また，パラメータによっては結果の値/分布が安定しないものもあるので，何度か実験を行なって結果を吟味している部分もあるが，結果を曲解するためのものではないということを付け加えておく．

なお，この文章ではメトロポリス・ヘイスティング・アルゴリズムを以降MHと表記する．

問題

ガンベル分布（二重指数分布）からの乱数生成を考える．今回はガンベル分布からの直接的な乱数生成ができないとしたうえで，MHによる乱数生成を行う．ガンベル分布は $\frac{1}{2b}\exp(\frac{-|x-\mu|}{b})$ で表され，今回は $\mu=0$ ， $b= 3$ の分布を考える．なお，ガンベル分布の平均と分散は解析的に得られ，それぞれ $\mu$ と $2b^2$ である．

実験条件

今回は，MHの中でも酔歩過程/ランダムウォーク（Random walk）と独立連鎖（Independent chain）の2種類の方法を試す．それぞれのMHにおける推移カーネル $K$ は

酔歩過程： $K(x^{(t)},x^{(t+1)}) = \mathcal{N}(x^{(t)},c_{ran}^2)$
独立連鎖： $K(x^{(t)},x^{(t+1)}) = \mathcal{N}(0,c_{ind}^2)$

となる．ここで， $c_{ran}^2$ および $c_{ind}^2$ は調整母数/チューニングパラメータであり，今回はそれぞれ{0.25, 1, 5, 15, 150}の5パターンで実験を行った．

また，生成した乱数の系列において，系列初期の $t=1\dots 1000$ までの結果は初期値に依存した値となっているため，今回は結果として用いていない．

結果：以下の図の見方

以下に示す図は，各行にチューニングパラメータの値，各列に系列のプロット，系列の分布，系列の平均の推移，系列の自己相関プロットを示したものである．図のタイトルに，チューニングパラメータの値，生成した乱数の平均と分散，そして採択率を示している（便宜的にタイトルの部分に情報を書き加えているだけで，タイトルと図は直接的には関連していない）．

結果：酔歩過程/ランダムウォーク

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")
}

$c_{ran}^2=0.25$ の場合

$c_{ran}^2=1$ の場合

$c_{ran}^2=5$ の場合

$c_{ran}^2=15$ の場合

$c_{ran}^2=150$ の場合

結果：独立連鎖の場合

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)
}

$c_{ind}^2=0.25$ の場合

$c_{ind}^2=1$ の場合

$c_{ind}^2=5$ の場合

$c_{ind}^2=15$ の場合

$c_{ind}^2=150$ の場合

解釈

さて，このように色々とプロットしてみたわけだが，酔歩過程の場合も独立連鎖の場合もチューニングパラメータの値が0.25や150といった非常に小さい/大きい場合においては，乱数の分布や収束速度などに問題があることがわかる．

酔歩過程においては，マルコフ連鎖の候補点の生成は少しずつ値を移動させていくという点で，分布の表面をなぞる形となる．そのため，チューニングパラメータの値である分散が小さい場合では候補点の移動範囲が狭いため，分布全体を探索するのに非常に時間がかかる． $c_{ran}^2=0.25$ の図を見ると，系列のプロットではあまり大幅な移動（縦のぶれ）が少なく，生成した乱数の分布もかなりいびつである．また，値があまり変動しないという点で，自己相関係数も非常に大きな値を取っていることがわかる．一方チューニングパラメータの値が大きい場合においては，採択率が非常に小さな値を取り，系列のプロットからも値の変動しない部分が他と比べて多い．

一方，独立連鎖においては，マルコフ連鎖の候補点の生成は現在の状態に依存しない．そのため，酔歩過程のように値を移動させていくということは見られないが，一方で採択確率の部分でチューニングパラメータの値が影響していることがわかる．