Slice sampling（スライスサンプリング）の簡単な例をRでやってみた。

テラモナギさんの記事と，PRMLの11.4節を参考にしてみた。
https://teramonagi.hatenablog.com/entry/20130203/1359888495#fn-b0191759
本当はもとの論文を見るべきだが，大体あっていそうなのでOKなような気がする。

（1変数の確率密度からの）スライスサンプリングのアルゴリズムまとめ

1. サンプリングのための変数（x）の初期化（x_oldとする。）
2. 目標とする確率密度（規格化していないものでOK）をfとしたとき，一様分布(min=0, max=f(x_old) )からスライスする位置uをサンプリングする。
3. {x: u < f(x)}となる，領域から一様にx_newをサンプリングして，x_newをサンプルとして保存し，x_old = x_newとする。
4. 2~3を繰り返す。

3を実行するときに，lower_xとupper_xを適当に決める。このときに，f(lower_x), f(upper_x) < uになるように，lower_x, upper_xを引き伸ばす必要がある（steppin out）。
さらに，その後，一様分布(min=lower_x, max=upper_x )からx_newをサンプリングし，f(x_new) < uとなってしまっていたら，lower_x, upper_xをx_newに置き換えて（shrinkage），再度x_newを一様分布からサンプリングする。
f(x_new) >= uであれば，x_newを新しいMCMCサンプルとして保存する。
x_newはスライスの範囲，つまり，{x: u < f(x)}に入っていなければならないので，この条件に入らないx_newはサンプルとして採用できない。
この辺のややめんどい手続きが必要になるが，これは，{x: u = f(x)}となるxの点を数値的に見つけないといけないから。
しかし，メトロポリスヘイスティングアルゴリズムのような，採択・棄却サンプリングのように，提案分布の調整がいらないので，とても楽。
しかも，確率密度関数のカーネルだけわかっていればよいので，汎用性が高い。
多変量の場合も，サンプリングしたいパラメタ以外のパラメタなどをすべて条件づけたものを利用すればよい。
（確率密度の部分については，対数をとったものを考えるというのも手。数値がアンダーフローになるのを防ぐため。）

#
# Slice サンプリング
#

#
# サンプリングを目指したい混合分布
# 
f <- function(x, 
              mixing=c(0.7, 0.3),
              mu=c(-1, 4),
              sigma=c(0.5,2)){
  y <- 0
  for(i in seq_along(mixing)){
    y <- y + mixing[i]*exp(-(x-mu[i])^2/(2*(sigma[i]^2))) 
  }
  y 
}

#
# 目標の分布
# 
curve(f, from=-3, to = 10)

# 
set.seed(111)

is_f_x_in_s <- function(u, f_x) u <= f_x

slice_one_step <- function(x_old, width = 1.0){
  
  
  #x_old <- rnorm(1)
  
  f_x_old <- f(x_old)
  u <- runif(1, min = 0, max=f_x_old)
  #f_u <- f(u)
  
  
  lower_x <- x_old - width*runif(1)
  upper_x <- x_old + width*runif(1)
  
  f_lower_x <- f(lower_x)
  f_upper_x <- f(upper_x)
  
  
  #
  # Stepping out
  # 
  while(is_f_x_in_s(u, f_lower_x)){
    lower_x <- lower_x - width 
    f_lower_x <- f(lower_x)
  }
  while(is_f_x_in_s(u, f_upper_x)){
    upper_x <- upper_x + width 
    f_upper_x <- f(upper_x)
  }
  
  #
  # Shrinkage
  # 
  repeat{
    x_new <- runif(1, min=lower_x, max=upper_x)
    f_x_new <- f(x_new)
    if(is_f_x_in_s(u, f_x_new)){break}
    
    if(x_new < x_old){
      lower_x <- x_new
    } else {
      upper_x <- x_new
    }
  }
  
  x_new
}


n_iter = 1000
slice_sampling <- function(n_iter = 1000){
  x_samples <- rep(0, n_iter)
  x_samples[1] <- x_old <- rnorm(1)
  
  for(i in 2:n_iter){
    cat("i = ",i,"\n")
    x_new <- slice_one_step(x_old, width = 0.5)
    x_samples[i] <- x_new
    x_old <- x_new
  }
  x_samples
}



res_x <- slice_sampling(n_iter = 10000)

hist(res_x,freq = F,ylim=c(0,0.4),
     breaks = seq(min(res_x)-0.5, max(res_x)+0.5,by=0.1),
     col=rgb(0,0,1,alpha=0.7),
     ylab="")
f_const <- integrate(f, lower=-10, upper=10)

#
# 数値的に規格化した密度
# 
f2 <- function(x)f(x)/f_const$value
curve(f2, from = min(res_x)-0.5 ,
      to= max(res_x)+0.5,
      col=rgb(1,0,0,alpha=0.8),lwd=2, 
      xlim=c(min(res_x)-0.5, max(res_x)+0.5),
      ylim=c(0,0.4),
      ylab="",
      add=T)

サンプリングの結果とサンプリングをしたい密度