蒙特卡洛采样之拒绝采样（Reject Sampling）

引子

蒙特卡洛（Monte Carlo）方法是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为基础的数值计算方法。它的核心思想就是使用随机数（或更常见的伪随机数）来解决一些复杂的计算问题。

当所求解问题可以转化为某种随机分布的特征数（比如随机事件出现的概率，或者随机变量的期望值等）时，往往就可以考虑使用蒙特卡洛方法。通过随机抽样的方法，以随机事件出现的频率估计其概率，或者以抽样的数字特征估算随机变量的数字特征，并将其作为问题的解。这种方法多用于求解复杂的高维积分问题。

实际应用中，我们所要面对的第一个问题就是如何抽样？在统计学中， 抽样（或称采样）是指从目标总体中抽取一部分个体作为样本的过程。

例如，我们想知道一所大学里所有男生的平均身高。但是因为学校里的男生可能有上万人之多，所以为每个人都测量一下身高可能存在困难，于是我们从每个学院随机挑选出100名男生来作为样本，这个过程就是抽样。

但是在计算机模拟时，我们所说的抽样，其实是指从一个概率分布中生成观察值（observations）的方法。而这个分布通常是由其概率密度函数（PDF）来表示的。而且，即使在已知PDF的情况下，让计算机自动生成观测值也不是一件容易的事情。从本质上来说，计算机只能实现对均匀分布（Uniform distribution）的采样。

具体来说，我们可能要面对的问题包括：

• 计算机只能实现对均匀分布的采样，但我们仍然可以在此基础上对更为复杂的分布进行采样，那具体该如何操作呢？
• 随机分布的某些数字特征可能需要通过积分的形式来求解，但是某些积分可能没有（或者很难求得）解析解，彼时我们该如何处理呢？
• 在贝叶斯推断中，后验概率的分布是正?于先验和似然函数之积的，但是先验和似然函数的乘积形式可能相对复杂，我们又该如何对这种形式复杂的分布进行采样呢？

欢迎关注白马负金羁的博客 http://blog.csdn.net/baimafujinji，为保证公式、图表得以正确显示，强烈建议你从该地址上查看原版博文。本博客主要关注方向包括：数字图像处理、算法设计与分析、数据结构、机器学习、数据挖掘、统计分析方法、自然语言处理。

Inverse CDF 方法

比较简单的一种情况是，我们可以通过PDF与CDF之间的关系，求出相应的CDF。或者我们根本就不知道PDF，但是知道CDF。此时就可以使用Inverse CDF的方法来进行采样。这种方法又称为逆变换采样（Inverse transform sampling）。

如果你对PDF和CDF的概念有点模糊，我们不妨先来一起回顾一下它们的定义。对于随机变量 X，如下定义的函数 F：

F(x)=P{X≤x},?∞<x<∞

称为 X 的累积分布函数（CDF，Cumulative Distribution Function）。对于连续型随机变量 X 的累积分布函数 F(x)，如果存在一个定义在实数轴上的非负函数 f(x)，使得对于任意实数 x，有下式成立：

F(x)=∫∞?∞f(t)dt

则称 f(x) 为 X 的概率密度函数（PDF，Probability Density Function）。显然，当概率密度函数存在的时候，累积分布函数是概率密度函数的积分。

所以，通常我们可以通过对PDF（如下图中的左图所示为正态分布的PDF）进行积分来得到概率分布的CDF（如下图中的右图所示为正态分布的CDF）。然后我们再得到CDF的反函数 F?1(u)，如果你想得到 m 个观察值，则重复下面的步骤 m 次：

• 从 Uniform(0,1) 中随机生成一个值（前面已经说过，计算机可以实现从均匀分布中采样），用 u 表示。
• 计算 F?1(u) 的值 x，则 x 就是从 f(x) 中得出的一个采样点。

以下图为例，如果从 Uniform(0,1) 中随机生成的值 u=0.8413，则可以算得 F?1(u)=1，则此次从正态分布中生成的随机数就是 1。

你可能会好奇，面对一个具有复杂表达式的函数， Inverse CDF 方法真的有效吗？来看下面这个例子。假设现在我们希望从下面这个PDF中抽样：

f(x)=?????8x83?83x0,if 0≤x<0.25,if 0.25≤x≤1,otherwise

可以算得相应的CDF为

F(x)=???????04x283x?43x2?131,if x<0,if 0≤x<0.25,if 0.25≤x≤1,if x>1

对于 u∈[0,1]，它的反函数为

F?1(u)=?????u√21?3(1?u)￣￣￣￣￣￣￣￣√2,if 0≤u<0.25,if0.25≤u≤1

下面我们在R中利用上面的方法采样10000个点，并以此来演示抽样的效果。

m <- 10000
u <- runif(m, 0, 1)
invcdf.func <- function(u) {
+ if (u >= 0 && u < 0.25)
+         sqrt(u)/2
+     else if (u >= 0.25 && u <= 1)
+         1 - sqrt(3 * (1 - u))/2
+ }
x <- unlist(lapply(u, invcdf.func))

curve(8*x, from = 0, to = 0.25, xlim=c(0,1), ylim=c(0,2), col = "red",xlab = "", ylab="")
par(new=TRUE)
curve((8/3)-(8/3)*x, from = 0.25, to = 1, xlim=c(0,1), ylim=c(0,2),
+ col = "red",xlab = "", ylab="")
par(new=TRUE)
plot(density(x), xlim=c(0,1), ylim=c(0,2), col = "blue", xlab = "x", ylab="density")

从下图中你可以发现我的采样点与原始分布非常吻合。

下面再举一个稍微复杂一点的例子，已知分布的PDF如下：

h(x)=2m2(1?m2)x3,x∈[m,1]

可以算得相应的CDF为

H(x)=∫x?∞$h(t)dt=?????011?m2?m2(1?m2)x21,if x<m,if x∈[m,1],if x>1

对于 u∈[0,1]，它的反函数为

H?1(u)=m21?(1?m2)u￣￣￣￣￣￣￣￣￣￣￣￣￣￣√

同样，我们给出R中的示例代码如下：

invcdf <- function(u, m) {
    return(sqrt(m^2/(1 - (1 - m^2) * u)))
}

sample1 <- sapply(runif(1000), invcdf, m = .5)

下面这段代码利用R中提供的一些内置函数实现了已知PDF时基于Inverse transform sampling 方法的采样，我们将新定义的函数命名为 samplepdf() 。当然，对于那些过于复杂的PDF函数（例如很难积分的），samplepdf() 确实有力所不及的情况。但是对于标准的常规PDF，该函数的效果还是不错的。

endsign <- function(f, sign = 1) {
    b <- sign
    while (sign * f(b) < 0) b <- 10 * b
    return(b)
}

samplepdf <- function(n, pdf, ..., spdf.lower = -Inf, spdf.upper = Inf) {
    vpdf <- function(v) sapply(v, pdf, ...)  # vectorize
    cdf <- function(x) integrate(vpdf, spdf.lower, x)$value
    invcdf <- function(u) {
        subcdf <- function(t) cdf(t) - u
        if (spdf.lower == -Inf)
            spdf.lower <- endsign(subcdf, -1)
        if (spdf.upper == Inf)
            spdf.upper <- endsign(subcdf)
        return(uniroot(subcdf, c(spdf.lower, spdf.upper))$root)
    }
    sapply(runif(n), invcdf)
}

下面我们就用 samplepdf() 函数来对上面给定的 h(x) 来进行采样，然后再跟之前所得之结果进行对比。

h <- function(t, m) {
    if (t >= m & t <= 1)
        return(2 * m^2/(1 - m^2)/t^3)
    return(0)
}

sample2 <- samplepdf(1000, h, m = .5)

plot(density(sample1), xlim=c(0.4, 1.1), ylim=c(0, 4), col = "red", xlab = "", ylab="", main="")
par(new=TRUE)
plot(density(sample2), xlim=c(0.4, 1.1), ylim=c(0, 4), col ="blue",
+ xlab = "x, N=1000", ylab="density", main="")
text.legend = c("my_invcdf","samplepdf")
legend("topright", legend = text.legend, lty = c(1,1), col = c( "red", "blue"))

代码执行结果如下图所示。

拒绝采样（Reject Sampling）

我们已经看到 Inverse CDF 方法确实有效。但其实它的缺点也是很明显的，那就是有些分布的 CDF 可能很难通过对 PDF 的积分得到，再或者 CDF 的反函数也很不容易求。这时我们可能需要用到另外一种采样方法，这就是我们即将要介绍的拒绝采样。

下面这张图很好地阐释了拒绝采样的基本思想。假设我们想对 PDF 为 p(x) 的函数进行采样，但是由于种种原因（例如这个函数很复杂），对其进行采样是相对困难的。但是另外一个 PDF 为 q(x) 的函数则相对容易采样，例如采用 Inverse CDF 方法可以很容易对对它进行采样，甚至 q(x) 就是一个均匀分布（别忘了计算机可以直接进行采样的分布就只有均匀分布）。那么，当我们将 q(x) 与一个常数 M 相乘之后，可以实现下图所示之关系，即 M?q(x) 将 p(x) 完全“罩住”。

然后重复如下步骤，直到获得 m 个被接受的采样点：

• 从 q(x) 中获得一个随机采样点 x(i)
• 对于 xi 计算接受概率（acceptance probability）
  
  α=p(xi)Mq(xi)
• 从 Uniform(0,1) 中随机生成一个值，用 u 表示
• 如果 α≥u，则接受 xi 作为一个来自 p(x) 的采样值，否则就拒绝 xi 并回到第一步

你当然可以采用严密的数学推导来证明Reject Sampling的可行性。但它的原理从直观上来解释也是相当容易理解的。你可以想象一下在上图的例子中，从哪些位置抽出的点会比较容易被接受。显然，红色曲线和绿色曲线所示之函数更加接近的地方接受概率较高，也即是更容易被接受，所以在这样的地方采到的点就会比较多，而在接受概率较低（即两个函数差距较大）的地方采到的点就会比较少，这也就保证了这个方法的有效性。

我们还是以本文最开始给出的那个分段函数 f(x) 为例来演示 Reject Sampling 方法。如下面图所示，参考分布我们选择的是均匀分布（你当然可以选择其他的分布，但采用均匀分布显然是此处最简单的处理方式）。而且令常数 M=3。

下面给出R中的示例代码，易见我们此处采样点数目为10000。

f.x <- function(x) {
+ if (x >= 0 && x < 0.25)
+ 8 * x
+ else if (x >= 0.25 && x <= 1)
+ 8/3 - 8 * x/3
+ else 0
+ }

g.x <- function(x) {
+ if (x >= 0 && x <= 1)
+ 1
+ else 0
+ }

M <- 3
m <- 10000
n.draws <- 0
draws <- c()
x.grid <- seq(0, 1, by = 0.01)

while (n.draws < m) {
+ x.c <- runif(1, 0, 1)
+ accept.prob <- f.x(x.c)/(M * g.x(x.c))
+ u <- runif(1, 0, 1)
+ if (accept.prob >= u) {
+     draws <- c(draws, x.c)
+     n.draws <- n.draws + 1
+     }
+ }

同样，我们还是用图形来展示一下Reject Sampling的效果如何。绘图的代码如下：

curve(8*x, from = 0, to = 0.25, xlim=c(0,1), ylim=c(0,2),
+ col = "red", xlab = "", ylab="", main="")
par(new=TRUE)
curve((8/3)-(8/3)*x, from = 0.25, to = 1, xlim=c(0,1),
+ ylim=c(0,2), col = "red",xlab = "", ylab="", main="")
par(new=TRUE)
plot(density(draws), xlim=c(0,1), ylim=c(0,2), col = "blue",
+ xlab = "x, N=10000", ylab="density")

text.legend = c("Target Density","Rejection Sampling")
legend("topright", legend = text.legend, lty = c(1,1), col = c( "red", "blue"))

上述代码的执行结果如下图所示，可见采样结果是非常理想的。

最后给出的这个额外的例子来自文献[3]，这里同样采用uniform分布作为参考分布。而需要被采样的分布则是我们在之前的文章《先验概率、后验概率以及共轭先验》中介绍过的Beta分布。如果你读过前面的博文，或者对机器学习比较了解，应该知道Beta分布的PDF表达式非常复杂。而且，你应该能够看出，跟前面不太一样的地方是，这的 Mq(x)，所取之值就是Beta分布的极值。它的图形应该是与Beta(3,6)的极值点相切的一条水平直线。

sampled <- data.frame(proposal = runif(50000,0,1))
sampled$targetDensity <- dbeta(sampled$proposal, 3,6)

maxDens = max(sampled$targetDensity, na.rm = T)
sampled$accepted = ifelse(runif(50000,0,1) < sampled$targetDensity / maxDens, TRUE, FALSE)

hist(sampled$proposal[sampled$accepted], freq = F, col = "grey", breaks = 100, main="")
curve(dbeta(x, 3,6),0,1, add =T, col = "red", main="")

上述代码的执行结果如下图所示，可见我们的采样结果与目标分布相当吻合。

Reject Sampling方法确实可以解决我们的问题。但是它的一个不足涉及到其采样效率的问题。就本文所给出的两个例子而言，第二个的处理方法是更好的选择。尽管两个例子都采用uniform来作为参考分布，但是第一个例子中的 Mq(x) 其实离目标分布还有一定距离，这样在采用过程中被拒绝的点就会更多。而第二个例子中，我们选择了离目标函数最近的参考函数，就uniform而言，已经不能有更进一步的方法了。但即使这种，在这个类似钟型的图形两侧其实我们仍然会拒绝掉很多很多采样点，这种开销其实相当浪费。最理想的情况下，参考分布应该跟目标分布越接近越好，从图形上来看就是包裹的越紧实越好。但是这种情况的参考分布往往又不那么容易得到。在满足某些条件的时候也确实可以采用所谓的改进方法，即Adaptive Rejection Sampling。

自适应的拒绝采样（Adaptive Rejection Sampling）

前面我们已经分析了，拒绝采样的弱点在于当被拒绝的点很多时，采样的效率会非常不理想。同时我们也支持，如果能够找到一个跟目标分布函数非常接近的参考函数，那么就可以保证被接受的点占大多数（被拒绝的点很少）。这样一来便克服了拒绝采样效率不高的弱点。如果函数是 log-concave 的话，那么我们就可以采样自适应的拒绝采样方法。什么是 log-concave 呢？还是回到我们之前介绍过的 Beta 分布的PDF，我们用下面的代码来绘制 Beta(2, 3) 的函数图像，以及将 Beta(2, 3) 的函数取对数之后的图形。

> integrand <- function(x) {(x^1)*((1-x)^2)}
> integrate(integrand, lower = 0, upper = 1)
0.08333333 with absolute error < 9.3e-16
> f<-function(x,a,b){log(1/0.08333)+(a-1)*log(x)+(b-1)*log(1-x)}
> curve(f(x, 2, 3))
> curve(dbeta(x, 2, 3))

上述代码的执行结果如下所示。其中左图是Beta(2, 3) 的函数图像，右图是将 Beta(2, 3) 的函数取对数之后的图形，你可以发现结果是一个凹函数（concave）。那么Beta(2, 3) 就满足log-concave的要求。

然后我们在对数图像上找一些点做图像的切线，如下图所示。因为对数图像是凹函数，所以每个切线都相当于一个超平面，而且对数图像只会位于超平面的一侧。

我同时给出用以绘制上图的代码，要知道R语言的一个强项就是绘图！

log_f <- function(x,a,b){log(1/0.08333)+(a-1)*log(x)+(b-1)*log(1-x)}
g <- function(x,a,b){(a-1)/x-(b-1)/(1-x)}

log_f_y1 <- log_f(0.18, 2, 3)
log_f_y2 <- log_f(0.40, 2, 3)
log_f_y3 <- log_f(0.65, 2, 3)
log_f_y4 <- log_f(0.95, 2, 3)

g1 <- g(0.18, 2, 3)
b1 <- log_f_y1 - g1*0.18
y1 <- function(x) {g1*x+b1}

g2 <- g(0.40, 2, 3)
b2 <- log_f_y2 - g2*0.40
y2 <- function(x) {g2*x+b2}

g3 <- g(0.65, 2, 3)
b3 <- log_f_y3 - g3*0.65
y3 <- function(x) {g3*x+b3}

g4 <- g(0.95, 2, 3)
b4 <- log_f_y4 - g4*0.95
y4 <- function(x) {g4*x+b4}

curve(log_f(x, 2, 3), col = "blue", xlim = c(0,1), ylim = c(-7, 1))
curve(y1, add = T, lty = 2, col = "red", to = 0.38)
curve(y2, add = T, lty = 2, col = "red", from = 0.15, to=0.78)
curve(y3, add = T, lty = 2, col = "red", from = 0.42)
curve(y4, add = T, lty = 2, col = "red", from = 0.86)

par(new=TRUE)

xs = c(0.18, 0.40, 0.65, 0.95)
ys = c(log_f_y1, log_f_y2, log_f_y3, log_f_y4)
plot(xs, ys, col="green", xlim=c(0,1), ylim=c(-7, 1), xlab="", ylab="")

再把这些切线转换回原始的Beta(2, 3)图像中，显然原来的线性函数会变成指数函数，它们将对应用下图中的一些曲线，这些曲线会被原函数的图形紧紧包裹住。特别是当这些的指数函数变得很多很稠密时，以彼此的交点作为分界线，我们其实相当于得到了一个分段函数。这个分段函数是原函数的一个逼近。用这个分段函数来作为参考函数再执行Reject Sampling，自然就完美的解决了我们之前的问题。

下面是我用来画出上图的R语言代码。

e_y1 <- function(x) {exp(g1*x+b1)}
e_y2 <- function(x) {exp(g2*x+b2)}
e_y3 <- function(x) {exp(g3*x+b3)}
e_y4 <- function(x) {exp(g4*x+b4)}

curve(dbeta(x, 2, 3), col="blue", ylim=c(0, 2.0))
curve(e_y1(x), add=T, lty=3, col = "red")
curve(e_y2(x), add=T, lty=3, col = "red", from = 0.2, to = 0.75)
curve(e_y3(x), add=T, lty=3, col = "red", from=0.48)
curve(e_y4(x), add=T, lty=3, col = "red", from=0.86)

这无疑是一种绝妙的想法。而且这种想法，我们在前面其实已经暗示过。在上一部分最后一个例子中，我们其实就是选择了一个与原函数相切的uniform函数来作为参考函数。我们当然回想去选择更多与原函数相切的函数，然后用这个函数的集合来作为新的参考函数。只是由于原函数的凹凸性无法保证，所以直线并不是一种好的选择。而ARS（Adaptive Rejection Sampling）所采用的策略则非常巧妙地解决了我们的问题。当然函数是log-concave的条件必须满足，否则就不能使用ARS。

下面给出了一个在R中进行Adaptive Rejection Sampling的例子。显然，这个例子要比之前的代码简单许多。因为R中ars包为我们提供了一个现成的用于执行Adaptive Rejection Sampling的函数，即ars()。关于这个函数在用法上的一些细节，读者还可以进一步参阅R的帮助文档，我们这里不再赘言。此次我们需要指出：ars()函数中两个重要参数，一个是对原分布的PDF取对数，另外一个则是对PDF的对数形式再进行求导（在求导时我们忽略了前面的系数项），其实也就是为了确定切线。

f<-function(x,a,b){(a-1)*log(x)+(b-1)*log(1-x)}
fprima<-function(x,a,b){(a-1)/x-(b-1)/(1-x)}
mysample<-ars(20000, f, fprima, x=c(0.3,0.6), m=2, lb=TRUE, xlb=0,
+ ub=TRUE, xub=1, a=1.3, b=2.7)
hist(mysample, freq=F)
curve(dbeta(x,1.3,2.7), add=T, col="red")

上述代码的执行结果如下图所示。

至此，我们已经对蒙特卡洛采样法中的拒绝采样进行了完整的介绍。在后续的文章中，我们将会对包括重要性采样（Importance Sampling）和MCMC在内的其他基于蒙特卡洛思想设计的采样方法再做介绍。但最后我们仍然需要指出，这些采样方法在现代机器学习技术中占据着非常重要的地方。因为这些采样方法可以作为很多模型参数估计的基础，这一点我们后续的文章中也会有进一步的说明。

参考文献与推荐阅读

[1] http://blog.quantitations.com/tutorial/2012/11/20/sampling-from-an-arbitrary-density/

[2] http://www.people.fas.harvard.edu/~plam/teaching/methods/sampling/sampling.pdf

[3] http://www.r-bloggers.com/a-simple-explanation-of-rejection-sampling-in-r/

[4] Gilks, W.R., P. Wild. (1992) Adaptive Rejection Sampling for Gibbs Sampling, Applied Statistics, 41:337–348.

[5] 同时推荐悉尼科大徐亦达博士的机器学习公开课