Slice sampling

  • ????????????????????

ベイズ統計などではよくマルコフ連鎖モンテカルロ(MCMC)を使います。MCMCというと Metropolis Hastings とか Gibbs Sampler を思い浮かべる人が多いと思います。しかし、MCMCというのはこの2つだけではないようです。ここでは最近、MCMCを実装する中で知った slice sampling という方法を実装してテストしてみた結果について書いてみます。参照したのは

Slice sampling, Radford M. Neal, Source: Ann. Statist. Volume 31, Number 3 (2003), 705-767.

です。

次のような状況を考えてみます。

  • 一変数の確率分布を考える
  • f(x) はサンプリングしたい分布 p(x) に比例するような関数(f はわかっているけど p はわからない)
  • f(x) (および p(x) は単峰:モードがひとつしか無い)

このような仮定のもと、p(x) からサンプルを得るための slice sampler は R を使って以下のように書くことができます。

slice.sample <- function(n, x0, fun, w=1.0, m=5, log.dens=F){

  if(!w>0)
    error("w must be positive")

  if(!m>0)
    error("m must be positive")

  if(!log.dens){
    log.fun <- function(x) log(fun(x))
  } else {
    log.fun <- fun
  }

  cnt=1
  samples=numeric(n)

  repeat { 
    z <- log.fun(x0)-rexp(1)
    U <- runif(1)
    L <- x0-w*U
    R <- L+w
    J <- floor(m*runif(1))
    K <- (m-1)-J

    # Finding an interval around x0
    while( J>0 && z<log.fun(L) ){
      L <- L-w
      J <- J-1
    }

    while( K>0 && z<log.fun(R) ){
      R <- R+w
      K <- K-1
    }

    # Sampling the next point
    repeat {
      x1 <- L+runif(1)*(R-L)

      if( z<log.fun(x1) ) break
      if( x1<x0 )
        L <- x1
      else
        R <- x1
    }

    samples[cnt] <- x1

    if(cnt==n)  break

    # Modifying w
    w <- (w+abs(x1-x0))/2.0
    cnt <- cnt+1
    x0 <- x1
  }
  samples
}  

使い方は、以下のように。引数は、1番目;サンプル数、2番目;初期値、3番目;f(x)、4番目;ステップ幅(内部で自動調整される)、5番目;スライスの幅の最大値を定める。大きいほどベター(ただし計算は遅くなる)??

数値実験

たとえば「標準正規分布の0.5乗に比例する確率分布」を考えてみます。こんな嫌な感じの分布でも手軽にサンプリングできます。

s=slice.sample(50000,0,function(x) 4*dnorm(x,3)^0.5,w=1,m=5)

これをプロットしてみると、以下のようになります。

ヒストグラム

サンプルパス

なんだかそれらしい分布を得ることができました。

結果の確認

えられた分布が本当にもとの分布から得られたかチェックしてみましょう。標準正規分布よりも裾の厚い分布が得られているはずですね。

plot(sort(s),(1:50000)/50000,type="l",xlab="X",ylab="CDF",xlim=c(-6,6))
x=seq(-10,10,0.001)
y=cumsum(dnorm(seq(-10,10,0.001))^0.5)
y=y/y[20001]  # normalize
lines(x,y,col=2)
lines(x,pnorm(x),col=3)   # 標準正規分布との比較
legend("topleft",legend=c("MCMC result","True (normal^0.5)","Normal distribution"),col=1:3,lw=1)


ということでOKですね。

はてなブックマーク - Slice sampling
Pocket