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

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

例 6.1 メトロポリス・ヘイスティング・アルゴリズムでベータ分布の乱数を生成する

例2.7では受理・棄却法を用いてベータ分布の乱数を生成したが，本問題ではメトロポリス・ヘイスティング・アルゴリズムを使ってベータ分布からサンプリングを行うことを目的とする．

a <- 2.7
b <- 6.3
c <- 2.669
Nsim <- 5000
X <- rep(runif(1), Nsim)

for(i in 2:Nsim){
  Y <- runif(1)
  rho <- dbeta(Y, a, b)/dbeta(X[i-1], a, b)
  X[i] <- X[i-1] + (Y - X[i-1])*(runif(1)<rho)
}

plot(1:Nsim, X, type="l", xlab="Iterations", main="Sequence X^t")
plot(4500:4800, X[4500:4800], type="l", xlab="Iterations", main="Sequence X^t for t = 4500 - 4800")

まず，以下の図は系列 ${X^{(t)}}$ をプロットしたものである．

そして，その中からシミュレーション回数が4500〜4800回までの系列を抜き出したのが以下の図となっている．このように一部分を拡大して見てみると，系列を表す線が変化せずにx軸に並行に遷移している部分が見られる．これは $\rho(x^{(t)},Y_t)$ の確率で棄却されて値が変化していないことを表している．これは対応する $Y_t$ に依存し，マルコフ連鎖の性質から次の状態にも影響するため，値が変化しない状態が一定期間続くことがあるからである．

また，得られた系列がベータ分布の乱数となっているかを確かめるため，左側に今回のメトロポリス・ヘイスティング・アルゴリズムで作成した乱数，右側にRの組込み関数rbeta()で生成した乱数をヒストグラムで図示し，今回実験で設定したベータ分布 $\mathcal{Be}(2.7,6.3)$ を赤の曲線で表した．

par(mfrow=c(1,2))
hist(X, nclass=100, freq=F, xlim=c(0,1.0), main="Metropolis-Hastings")
curve(dbeta(x,a,b), add=T, col="red", lwd=2)
hist(rbeta(Nsim,a,b), freq=F, nclass=100, xlim=c(0,1.0), main="Direct Generation")
curve(dbeta(x,a,b), add=T, col="red", lwd=2)

この図を見ると，メトロポリス・ヘイスティング・アルゴリズムを使った場合でも，受理・棄却法の時と同じく乱数の生成に成功していることがわかる．

例 6.2 メトロポリス・ヘイスティング・アルゴリズムでコーシー分布の乱数を生成する

例6.1のベータ分布の乱数生成における提案分布は $\mathcal{U}_{0,1}$ だったが，コーシー分布の乱数生成における提案分布には，

標準正規分布： $\mathcal{N}(0,1)$
自由度0.5のスチューデントt分布： $\mathcal{T}_{1/2}$

を使うことができる．このような提案分布は独立提案分布であり，マルコフ連鎖ではない．

なお，本問題ではこれらの分布を表す単語として候補分布と提案分布という2つの言葉が出てくるが，今回は表記として提案分布で統一する．

# Standard Normal Distributionp
Nsim <- 10^4
X <- c(rt(1,1))
for(t in 2:Nsim){
  Y <- rnorm(1)
  rho <- min(1, dt(Y, 1)*dnorm(X[t-1]) / (dt(X[t-1],1)*dnorm(Y)))
  X[t] <- X[t-1] + (Y - X[t-1])*(runif(1) < rho)
}

# t-distribution (0.5 degree of freedom)
X2 <- c(rt(1,1))
for(t in 2:Nsim){
  Y <- rt(1, 0.5)
  rho <- min(1, dt(Y, 1)*dt(X2[t-1], 0.5) / (dt(X2[t-1],1)*dt(Y, 0.5)))
  X2[t] <- X2[t-1] + (Y - X2[t-1])*(runif(1) < rho)
}

まずはそれぞれの系列をプロットして特徴を見ていく．左側に提案分布を標準正規分布とした場合，右側に提案分布をスチューデントt分布とした場合で，上から順に系列のプロット，系列のヒストグラム，系列の自己相関プロットを図示したのが以下の図になっている．

まずは標準正規分布を使った場合からみていく．系列のプロットは例6.1と同じような形を示しているが，イテレーションが8500を過ぎたあたりで値が一定期間変化しない部分がある．これは乱数から生成する正規候補 $Y_t$ が大きい（または小さい）値を取ってしまい，連鎖で長い間引きずることになったからである．この現象は，系列のヒストグラムの右端において値が少し大きくなっていることからも確認することができる．

一方，スチューデントt分布を使った場合においては，系列のプロットは0付近で集中しており，たまに正規候補 $Y_t$ が非常に大きい値を取ることはあるものの，それが継続して続くわけではないことがわかる．系列のヒストグラムを見ても，標準正規分布の場合と比べて分布がなだらかでコーシー分布に当てはまっていることがわかる．ただし，実際には両端の外れ値を除外して系列のヒストグラムをプロットしているので，このヒストグラムの描写範囲から外れている値があることには注意が必要である．

次に，2つの提案分布におけるメトロポリス・ヘイスティング・アルゴリズムによって生成したコーシー乱数の累積収束プロット，すなわち乱数の平均値がどのように収束していくかを図示する．

黒線が提案分布に標準正規分布を使った場合，黄色の線が提案分布にスチューデントt分布を使った場合となっている．今回シミュレートしているコーシー分布から求めた正確な値は0.896なので，スチューデントt分布はイテレーション回数を重ねるごとに真の値に収束していくことがわかる．一方で標準正規分布においては，今回のイテレーション回数(N=5000)では真の値に収束していないようにみえる．しかし，平均値の傾向を見ると値が徐々に上昇しては一気に下がり…といった傾向が見られ，理論的にも標準正規分布を提案分布として用いた場合でも真の値に収束することがわかっている．なので，今回の場合はイテレーション回数が少なすぎるために，収束の様子が見られない．