拒绝接受采样(R语言实现)

与其他采样最不同的点，拒绝接受采样无需利用积分计算的思想。比如有一个复杂分布f(x)，我们可以从好采样的分布中采样（比如均匀分布），然后依据一定的原则去拒绝或接受好采样的分布产生的一些样本值。使得接受的样本值可以合理的来自于分布f(x),其思想称之为拒绝接受采样。

假定有一分布$h(x)$使得$h(x)$恒大于$f(x)$，如下图所示：

随机取得点x1，再给定第二个维度，使得对于点x1的第二个维度为0到h（x1）的均匀分布如下：

若点x1在第二个维度上低于f(x)(绿色)则接受该样本点，否则拒绝该样本点(红色)

因此接受x1的机率与f(x1)的高度成正比，即随机变量X在点x1处的密度：

接下来可以生成一系列样本点，即可求得f(x)的密度

接下来对h(x)进行修正：

因为如上的拒绝接受形式效率不高。

因此可以找到更一般的均匀分布g(x)

定义m使得

$$m=max\left(\frac{f(x)}{g(x)}\right)$$

令

$$h(x)=mg(x)$$

能确保

$$h(x)\ge f(x)\text{对于任意的点x}$$

（注：在实际应用中并不一定要求m取到精确的$m=max\left(\frac{f(x)}{g(x)}\right)$，比最大值要大也可，只是效率没有精确的最大值来的高）

总结下拒绝接受采样生成来自于f(x)的随机变量的步骤：

1.从g(x)中采样生成随机样本x1。

2.从均匀分布U(0,mg(x1))生成第二个维度的随机变量u。

3.若mg(x1)<f(x1)，则落在f(x)内（绿色的点），此时接受点x1，认为此点来自于f(x)，否则拒绝，之后重复第一个步骤，生成一系列来自于f(x)的点。

若令U为0~1上的均匀分布，则等价于

$$U<\frac{f(X_1)}{mg(X_1)}\text{时接受}$$

最终可算得：

$$P(acc|X)=P(U<\frac{f(X)}{mg(X)}|X)=\frac{f(X)}{mg(X)}$$

那么因此总的接受概率为：

$$\sum _{x}P(acc|x)P(X=x)=\sum _x\frac{f(x)}{mg(x)}g(x)=\frac{1}{m}$$

m越小，总的接受概率越大则由g(x)产生的点更容易被视作为来自于f(x)，因此采样效率越高。

R语言实践

从beta(2,2)中采样$$f(x)=6x(1-x),0<x<1$$

令

$$g(x)\sim U(0,1)$$

此时m=3/2

则当：

$$4x(1-x)>u\text{时接受样本x}$$

#像产生来自beta(2,2)样本数
n <- 1000
#接受次数
k <- 0
#迭代次数
j <- 0
y <- numeric(n)

while(k<n){
  u <- runif(1)
  j <- j+1
  #从g(x)中采样
  x <- runif(1)
  if(4*x*(1-x)>u){
    #接受样本，认为从建议分布中生成的样本等价于从目标分布中生成的样本
    k <- k+1
    y[k] <- x
  }
}
cat("当m=3/2时，从建议分布中产生beta(2,2)的1000个样本共需迭代次数:",j)

比较拒绝接受采样与实际采样之间的差异

p <- seq(.1,.9,.1)
#拒绝接受的实际采样点
Qhat <- quantile(y,p)
#理论
Q <- qbeta(p,2,2)
round(rbind(Qhat, Q), 3)

可以看到拒绝接受采样获得数据与实际beta(2,2)的数据相差不大，并且当我们调整样本量n更大，会产生更为精确的拒绝接受采样的样本。

当然若取得的m比精确最大值m还大，也是可以的，只是会影响采样效率，比如我们取得m=6

此时
$$x(1-x)>u\text{时接受样本x}$$

#样本数
n <- 1000
#接受次数
k <- 0
#迭代次数
j <- 0
y <- numeric(n)

while(k<n){
  u <- runif(1)
  j <- j+1
  #从g(x)中采样
  x <- runif(1)
  if(x*(1-x)>u){
    #接受样本
    k <- k+1
    y[k] <- x
  }
}
j

p <- seq(.1,.9,.1)
#实际
Qhat <- quantile(y,p)
#理论
Q <- qbeta(p,2,2)
round(rbind(Qhat, Q), 3)

随m的增加，采样的效率变低