バナナでもわかる話

計量経済学・統計学が専門の大学院生です。統計学・経済学・投資理論・マーケティング等々に関する勉強・解説ブログ。ときどき趣味も。極力数式は使わずイメージで説明出来るよう心掛けていますが、時々暴走します。

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

sをプロットしてみるとこんな感じ
f:id:bananarian:20180902032535p:plain

見にくいですが、よくよく見ると初期値の最初の方は初期値に依存して小さい値をとっています。この初期値に依存している部分を取り除きましょう。この初期値に依存している時間をバーンイン期間と呼びます。

s = tail(s, length(s)-300)

再びプロットするとこんな感じ
f:id:bananarian:20180902032909p:plain


それでは取り出せたサンプルの算術平均と不偏分散を求めてみます。
f:id:bananarian:20180902033104p:plain

真の値にかなり近い値を得られていることが確認出来ました。



正規分布は別にこんな面倒な事をしなくてもrnormで簡単にサンプリングできるわけですが、もう少し複雑な分布だと重宝します。