「Rによるモンテカルロ法入門」読書ノート：6章ランダム・ウォーク・メトロポリス・ヘイスティング・アルゴリズム

読書ノートアーカイブ：http://yagays.github.com/blog/2012/10/20/archive-introducing-monte-carlo-methods-with-r/

例 6.4 ランダム・ウォークのパラメータ設定によるサンプリングへの影響を，正規分布の乱数生成を例に確かめる

摂動が $\varepsilon_t \sim \mathcal{U}_{[-\delta,\delta]}$ となるようなランダム・ウォークにおいて，正規分布 $\mathcal{N}(0,1)$ をサンプリングするときのメトロポリス・ヘイスティング・アルゴリズムの受理確率は

$\rho(x^{(t)},y_t) = \min \big[1, \exp{\frac{x^{(t)2}-y_t^2}{2} } \big]$

で表される．

今回は，一様候補のパラメータを $\delta = 0.1, 1, 10$ と変更して，メトロポリス・ヘイスティング・アルゴリズムによるローカルな探索がどういう挙動を示すのかを調べる．

par(mfcol=c(3,3))
for(d in c(0.1,1,10)){
  Niter <- 2500
  x <- 0
  for(i in 2:Niter){
    y <- x[i-1] + runif(1, min=-d, max=d)
    if (runif(1) < min(1, exp(-(1/2)*(y^2 - x[i-1]^2)) )){
      x[i] <- y
    }else{
      x[i] <- x[i-1]
    }
  }
  plot(1:Niter, x, type="l", xlab="Iterations", ylab="",main=paste("d = ",d))
  hist(x, nclass=150, freq=F, xlim=c(-4,4), main="")
  curve(dnorm(x), lwd=2, add=T, col="red")
  acf(x)
}

それぞれのパラメータにおける系列の変動，分布，自己相関プロットを示したのが以下の図である．それぞれのサンプリング結果を見比べてみると，

一様候補のパラメータが小さい場合（）
- ランダム・ウォークで移動する幅が狭い
- 分布が偏っている
- 自己相関が高い（前の値からあまり変わらない）
一様候補のパラメータが大きい場合（）
- ランダム・ウォークで移動する幅が大きい
- これも同様に分布が安定しない

という結果となり，ランダム・ウォーク・メトロポリス・ヘイスティング・アルゴリズムにおいては摂動（ $\varepsilon_t$ ）のパラメータ設定が重要であることがわかる．

例 6.10 ランダム・ウォーク・メトロポリス・ヘイスティング・アルゴリズムの受理率を調べる

メトロポリス・ヘイスティング・アルゴリズムにおけるアルゴリズム効率化の調整は，本書の2章で扱った受理・棄却アルゴリズムの時の最大受理率に基づいた調整の場合とは異なる．まず，上の例 6.4で示した正規分布からのサンプリングにおける受理率を見てみる．パラメータごとの受理率は以下のようになった（Rのコードは省略）．

$\delta = 0.1$ ： $r = 0.9792$
$\delta = 1$ ： $r = 0.8132$
$\delta = 10$ ： $r = 0.1716$

この結果から，受理率が高ければサンプリングの結果も良くなるわけではないことがわかる．このことから，ランダム・ウォーク・メトロポリス・ヘイスティング・アルゴリズムにおいては，受理率を高めるようにアルゴリズムを調節してはいけないと判断できる．このアルゴリズムにおいて，どのような受理率を取ればよいかといった疑問に関しては，以下の練習問題6.14のc.で詳しく見ていく．

練習問題 6.14 ランダム・ウォーク候補がサンプリング結果に与える影響について，複数の候補分布を用いてより詳細に調べる

a. P.205の図6.7を再現する

図6.7の再現と $\delta$ の解釈は，上の例6.4の解答を参考．

b. コーシー候補とラプラス候補を検討する

今度は，一様分布の代わりにコーシー分布やラプラス分布を候補分布として使った場合の，ランダム・ウォーク・メトロポリス・ヘイスティング・アルゴリズムの挙動を調べる．

コーシー候補

ランダム・ウォーク候補は $\varepsilon_t \sim \mathcal{C}(x_0,\gamma)$ とし，コーシー分布のパラメータは $x_0 = 0$ ， $\gamma = 0.01,\ 1,\ 10$ を考える．

# Cauchy candidate
par(mfcol=c(3,3))
for(g in c(0.01,1,10)){
  Niter <- 2500
  x <- 0
  for(i in 2:Niter){
    y <- x[i-1] + rcauchy(1, location=0, scale=g)
    if (runif(1) < min(1, exp(-(1/2)*(y^2 - x[i-1]^2)) )){
      x[i] <- y
    }else{
      x[i] <- x[i-1]
    }
  }
  plot(1:Niter, x, type="l", xlab="Iterations", ylab="",main=paste("Cauchy candidate,","gamma=",g))
  hist(x, nclass=150, freq=F, xlim=c(-4,4), main="")
  curve(dnorm(x), lwd=2, add=T, col="red")
  acf(x)
}

ラプラス候補

ランダム・ウォーク候補は $\varepsilon_t \sim \mathcal{L}(\mu,b)$ とし，ラプラス分布のパラメータは $\mu=0$ ， $b = 0.1,\ 1,\ 10$ を考える．

# Laplace candidate
par(mfcol=c(3,3))
for(a in c(0.1,1,10)){
  Niter <- 2500
  x <- 0
  for(i in 2:Niter){
    y <- x[i-1] + ifelse(runif(1)>0.5, 1, -1)*rexp(1)/a
    if (runif(1) < min(1, exp(-(1/2)*(y^2 - x[i-1]^2)) )){
      x[i] <- y
    }else{
      x[i] <- x[i-1]
    }
  }
  plot(1:Niter, x, type="l", xlab="Iterations", ylab="",main=paste("Laplace candidate,","b=",a))
  hist(x, nclass=150, freq=F, xlim=c(-4,4), main="")
  curve(dnorm(x), lwd=2, add=T, col="red")
  acf(x)
}

以上の結果から，コーシー候補やラプラス候補においても，ランダム・ウォークの値の取り方によって適切なサンプリングができるかどうかが決まってくることがわかる．また，どちらも自己相関が低い値の時にサンプリング結果が一番良いことが見て取れる．これは例6.4の一様候補の場合とは違う結果となったが，恐らくコーシー分布やラプラス分布はともに一峰型の正規分布のような形であることが効いていて，ランダム・ウォークの移動が適度に調節されているからだろう．

c. パラメータの調節により受理確率を0.25に近づける

さて，例6.10で述べたように，ランダム・ウォークを使ったメトロポリス・ヘイスティング・アルゴリズムにおいては，受理率の最大化がアルゴリズムの効率化に繋がらないことがわかった．受理率はどれくらいを目安にすればよいかという疑問に関しては，本書では「高次元では0.25，低次元では0.5程度」というRoberts et al.(2007)の提案を引用している．そこで，先ほどの3つのランダム・ウォーク候補について受理率を0.25に近づけるために，パラメータの値を0から10まで0.1刻みで設定したときの，一様候補・コーシー候補・ラプラス候補を用いたアルゴリズムにおける受理率の変化を以下の図でプロットした．