「Rによるモンテカルロ法入門」読書ノート：6章マルコフ連鎖とランダム・ウォーク

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

練習問題 6.1 自己回帰モデルAR(1)が定常分布を持つことを実験的に示す

この問題はいわゆるAR(1)というもので，本書の1章 P.34に同様のコードおよび作図が載っている．なお，本問題の詳細はMCMC-UseR.pdfのP.130を参照．

$X^{(t+1)} = \rho X^{(t)}+\varepsilon_t$ ， $\varepsilon_t \sim \mathcal{N}(0,1)$ で定義されるマルコフ連鎖において，系列 ${X^{(t)}}$ の定常分布は $\mathcal{N}(0, \frac{\sigma^2}{1-\rho^2})$ で表される．今回は $t=10^4$ までのシミュレーション結果をヒストグラムとしてプロットし，定常分布を赤い曲線で示した．

rho <- 0.9
x <- rnorm(1)
for(t in 2:10^4){
  x[t] <- x[t-1]*rho + rnorm(1)
}

hist(x, nclass=150, freq=F)
curve(dnorm(x, mean=0, sd=(1/sqrt(1-rho^2))), add=T, col="red")

練習問題 6.2 ランダム・ウォークが定常分布を持たないことを実験的に示す

$X^{(t+1)}=X^{(t)}+\varepsilon_t$ で表されるランダム・ウォークが，シミュレーション回数を増やしても定常分布に収束しないことを実験的に確かめる．今回は初期値 $X^{(0)}=0$ として， $t=10^4$ と $t=10^6$ までマルコフ連鎖をシミュレートする．この問題は練習問題6.1における $\rho = 1$ の場合と同じになるので，上で書いたコードのパラメータを変更するだけで再現できる．ただ，今回は別の書き方としてcumsum()関数を用いてランダム・ウォークの系列を計算している．

init <- 0
x1 <- cumsum(rnorm(10^4-1))
x2 <- cumsum(rnorm(10^6-1))

par(mfrow=c(2,2))
hist(x1, nclass=150, main="t=10^4")
plot(1:10^4, c(init,x1), type="l", xlab="t", main="t=10^4")
hist(x2, nclass=150, main="t=10^6")
plot(1:10^6, c(init,x2), type="l", xlab="t", main="t=10^6")

シミュレーション回数を $t=10^4$ と $t=10^6$ としてランダム・ウォークを実行した結果が上図である．左図に系列のヒストグラム，右図に試行回数を横軸に取った系列をプロットしている．上のように図示すると，たとえ10⁶ 回試行したとしても何らかの分布に収束しないことがわかる．

雑感

練習問題 6.1と練習問題6.2を比較すると，ランダム・ウォークにおいて先行する状態に $0 \leq \rho < 1$ のような値を掛け合わせるだけで定常分布を取るというのは，なんとも面白い結果となっている．

補足：MCMC-UseR.pdfのP.132

練習問題の $\rho$ と以下図の $\theta$ は同じ
pdfにも書いてあるがスケール（x軸とy軸の描写範囲）に注意．いかに $\theta$ が効いているかがわかる
左上はマルコフ連鎖ではない

par(mfrow=c(2,2))
for(rho in c(0, 0.4, 0.8, 0.95)){
  x <- rnorm(1)
  y <- rnorm(1)
  for(t in 2:1000){
    x[t] <- x[t-1]*rho + rnorm(1)
    y[t] <- y[t-1]*rho + rnorm(1)
  }
  plot(x,y, type="l", main=paste("theta=",rho))
}