Slice Smapling算法

最新推荐文章于 2021-07-19 19:47:21 发布

原创最新推荐文章于 2021-07-19 19:47:21 发布 · 1w 阅读

11 ·

CC 4.0 BY-SA版权

文章标签：

#算法 #function #plot #list #c

本文介绍了一种名为Slice Sampling的新采样技术，并详细解释了其算法流程。通过与Adaptive Rejection Sampling和Gibbs抽样的结合使用，可以有效解决从复杂分布中抽样的问题。文章还提供了使用R语言进行模拟实验的代码及结果。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

今天看到一种新的Sampling技术，颇为兴奋。配上adaptive rejection sampling, Gibbs抽样就可以“无敌”了。

大致说说Slice sampling的算法：

我们想从一个密度分布为f(x)的函数中抽样

初始化：任意选择x0，步长w

for i= 1:n

固定xi,在[0,f(xi)]区间上均匀抽样，得到yi

固定yi，在函数值大于yi的x区间（{x: f(x)>g(yi)}）均匀抽样。

重复步骤，最后得到{xi,yi}的联合分布，从里面直接抽取到{xi}即得到x的边际模拟抽样结果。

上述算法实现的难点在第二步，如何从{x: f(x)>g(yi)}上均匀抽样。一种简便的算法如下：

1.对于第i步的xi，任意选取包含xi的长度为w的区间。

2.该区间的左右两个端点分别向左和向右延拓，知道端点的值都小于g（xi）

3.在上述区间中均匀抽样得到y，但是我们以一定概率接受它；如果g(y)<g(yi)则拒绝，以g(y)为新的端点重新抽样；如果g(y)>g(yi)则接受yi=y.

以下是我simulation的结果：

真实分布是一个高斯的混合分布，用slice sampling模拟的密度估计如图所示：

附代码(R):

p=function(x,u1,sig1,u2,sig2){
(1/3)*(1/(sqrt(2*pi)*15)*exp(-0.5*(x-70)^2/15^2)+1/(sqrt(2*pi)*sig1)*exp(-0.5*(x-u1)^2/sig1^2)+1/(sqrt(2*pi)*sig2)*exp(-0.5*(x-u2)^2/sig2^2))
}

sample=function(init,time,w,u1,sig1,u2,sig2){
X=NULL
Y=NULL
X[1]=init

for (i in 1:time){
Y[i]=runif(1,0,p(X[i],u1,sig1,u2,sig2))
xl=X[i]-runif(1,0,w)
xr=xl+w
while(!((p(xl,u1,sig1,u2,sig2)<Y[i])&(p(xr,u1,sig1,u2,sig2)<Y[i]))){
xl=xl-w
xr=xr+w
}
x=runif(1,xl,xr)
while(p(x,u1,sig1,u2,sig2)<Y[i]){
if(x>X[i]){xr=x}else{xl=x}
x=runif(1,xl,xr)
}
if(p(x,u1,sig1,u2,sig2)>Y[i]){X[i+1]=x} else{print('error')}
}
return(list(X,Y))
}