Rで試すMCMC(ランダムウォーク連鎖の場合)

前回は詳細つり合い条件について確認しました。
bananarian.hatenablog.com

今回はとうとうMCMCに移っていきます。
提案分布は色々考えられますが、具体例として最も簡単なランダムウォーク連鎖について説明します。

ランダムウォーク連鎖

$x^{(t)}=x$ であるとき、候補として次のようなものを考える。
$y=x+\epsilon,$ $\epsilon ～N(0,\sigma^2)$ ←別に一様分布でも良い

ここで $\sigma$ はステップサイズと呼ばれる。

$u$ を一様分布 $U(0,1)$ で発生させ、次のように選ぶ。

$(u≦\alpha(x^{(t)},y)であるとき)→x^{(t+1)}=y,$
$(その他)→x^{(t+1)}=x^{(t)}$

ただしこの時採択確率 $\alpha$ は目標分布 $\pi$ を用いて次のように表される。

$\alpha(x,y)=min(1,\frac{\pi(y)}{\pi(x)})$

Rで試す

それでは、 $\pi(x)=\frac{1}{\sqrt(2\pi)}exp(-\frac{(x-5)^2)}{2}$ という目標分布からサンプリングする状況を考えます。これは普通の正規分布であり、分散1,期待値5であることが式の形からすぐわかります。

#目標分布を((sqrt(2*pi))^(-1))*exp(-((x-5)^2)/2)とする。

target=function(x){
	((sqrt(2*pi))^(-1))*exp(-((x-5)^2)/2)
}

# 初期値(別に何でもよいです)
s0 = 0.5
# 初期値を目標分布へつっこむ
pi0 =target(s0)
s = c(s0)
for(n in 1:10000){
　# 正規分布でランダムウォーク連鎖を作成
　 s1 = s0 + rnorm(1)
    # 一様分布
    u = runif(1)
    # 候補の値の代入
    pi1 = target(s1)
    # 採択確率を求める
    alpha = min(1, pi1/pi0)
    #採択について
    if(alpha>u){
      # s1をs0に上書き
      s0 = s1
      # pi0に候補値の尤度を代入
      pi0 = pi1
    }
    # s0をサンプル集合sに追加
    s <- c(s, s0)
}