蒙特卡洛采样之拒绝采样（RejectSampling ）
引⼦
蒙特卡洛（Monte Carlo ）⽅法是⼆⼗世纪四⼗年代中期由于科学技术的发展和电⼦计算机的发明，⽽被提出的⼀种以概率统计理论为基础的数值计算⽅法。它的核⼼思想就是使⽤随机数（或更常见的伪随机数）来解决⼀些复杂的计算问题。
当所求解问题可以转化为某种随机分布的特征数（⽐如随机事件出现的概率，或者随机变量的期望值等）时，往往就可以考虑使⽤蒙特卡洛⽅法。通过随机抽样的⽅法，以随机事件出现的频率估计其概率，或者以抽样的数字特征估算随机变量的数字特征，并将其作为问题的解。这种⽅法多⽤于求解复杂的⾼维积分问题。
实际应⽤中，我们所要⾯对的第⼀个问题就是如何抽样？在统计学中， 抽样（或称采样）是指从⽬标总体中抽取⼀部分个体作为样本的过程。
例如，我们想知道⼀所⼤学⾥所有男⽣的平均⾝⾼。但是因为学校⾥的男⽣可能有上万⼈之多，所以为每个⼈都测量⼀下⾝⾼可能存在困难，于是我们从每个学院随机挑选出100名男⽣来作为样本，这个过程就是抽样。
普通话测试成绩查询但是在计算机模拟时，我们所说的抽样，其实是指从⼀个概率分布中⽣成观察值（obrvations ）的⽅法。⽽这个分布通常是由其概率密度函数（PDF ）来表⽰的。⽽且，即使在已知PDF 的情况下，让计算机⾃动⽣成观测值也不是⼀件容易的事情。从本质上来说，计算机只能实现对均匀分布（Uniform distribution ）的采样。具体来说，我们可能要⾯对的问题包括：
计算机只能实现对均匀分布的采样，但我们仍然可以在此基础上对更为复杂的分布进⾏采样，那具体该如何操作呢？
随机分布的某些数字特征可能需要通过积分的形式来求解，但是某些积分可能没有（或者很难求得）解析解，彼时我们该如何处理呢？在贝叶斯推断中，后验概率的分布是正⽐于先验和似然函数之积的，但是先验和似然函数的乘积形式可能相对复杂，我们⼜该如何对这种形式复杂的分布进⾏采样呢？
Inver CDF ⽅法
⽐较简单的⼀种情况是，我们可以通过PDF 与CDF 之间的关系，求出相应的CDF 。或者我们根本就不知道PDF ，但是知道CDF 。此时就可以使⽤Inver CDF 的⽅法来进⾏采样。这种⽅法⼜称为逆变换采样（）。
如果你对PDF 和CDF 的概念有点模糊，我们不妨先来⼀起回顾⼀下它们的定义。对于随机变量 ，如下定义的函数 ：
称为 的累积分布函数（CDF ，Cumulative Distribution Function ）。对于连续型随机变量 的累积分布函数 ，如果存在⼀个定义在实数轴上的⾮负函数 ，使得对于任意实数 ，有下式成⽴：
则称 为 的概率密度函数（PDF ，Probability Density Function ）。显然，当概率密度函数存在的时候，累积分布函数是概率密度函数的积分。
所以，通常我们可以通过对PDF （如下图中的左图所⽰为正态分布的PDF ）进⾏积分来得到概率分布的CDF （如下图中的右图所⽰为正态分布的CDF ）。然后我们再得到CDF 的反函数 ，如果你想得到 个观察值，则重复下⾯的步骤 次：
X F F(x)=P{X ≤x},−∞<x <∞
X X F(x)f(x)x F(x)=f(t)dt
∫∞−∞
f(x)X (u)F −1m m u
从 Uniform(0,1) 中随机⽣成⼀个值（前⾯已经说过，计算机可以实现从均匀分布中采样），⽤ 表⽰。计算 的值 ，则 就是从 中得出的⼀个采样点。以下图为例，如果从 Uniform(0,1) 中随机⽣成的值 ，则可以算得 ，则此次从正态分布中⽣成的随机数就是
1。
你可能会好奇，⾯对⼀个具有复杂表达式的函数， Inver CDF ⽅法真的有效吗？来看下⾯这个例⼦。假设现在我们希望从下⾯这个PDF 中抽样：
可以算得相应的CDF 为
对于 ，它的反函数为 下⾯我们在R 中利⽤上⾯的⽅法采样10000个点，并以此来演⽰抽样的效果。
从下图中你可以发现我的采样点与原始分布⾮常吻合。
u (u)F −1x x f(x)u =0.8413(u)=1F −1f(x)=⎧⎩⎨⎪⎪
⎪⎪8x −x
83830
,
if 0≤x <0.25,if 0.25≤x ≤1,otherwi
F(x)=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪04x 2
x −−8343x 21
31
,if x <0
,if 0≤x <0.25,if 0.25≤x ≤1,if x >1
u
∈[0,1]
(u)=F −1⎧⎩⎨⎪⎪
⎪⎪
u √21−3(1−u)−−
−−−−−√2
,if 0≤u <0.25,if0.25≤u ≤1
m <- 10000u <- runif(m, 0, 1)
invcdf.func <- function (u) {+ if (u >= 0 && u < 0.25)+ sqrt (u)/2
+ el 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")
12345678910111213141516
下⾯再举⼀个稍微复杂⼀点的例⼦，已知分布的PDF 如下：
可以算得相应的CDF 为
对于 ，它的反函数为 同样，我们给出R 中的⽰例代码如下：
下⾯这段代码利⽤R 中提供的⼀些内置函数实现了已知PDF 时基于Inver transform sampling ⽅法的采样，我们将新定义的函数命名为samplepdf() 。当然，对于那些过于复杂的PDF 函数（例如很难积分的），samplepdf() 确实有⼒所不及的情况。但是对于标准的常规PDF ，该函数的效果还是不错的。
下⾯我们就⽤ samplepdf() 函数来对上⾯给定的 来进⾏采样，然后再跟之前所得之结果进⾏对⽐。
h(x)=,
x ∈[m,1]
2m 2
(1−)m 2x 3
数学悖论
H(x)=$h(t)dt =∫
x −∞
⎧⎩⎨⎪⎪0
−
11−m 2m 2(1−)m 2x 2
1
,if x <m ,if x ∈[m,1],if x >1
u ∈[0,1]
(u)=H
−1
m 2
1−(1−)u
m 2−−
−−−−−−−−−−√
invcdf <- function(u, m) {
return(sqrt(m^2/(1 - (1 - m^2) * u)))}
sample1 <- sapply(runif(1000), invcdf, m = .5)
七年级上册地理期中试卷12345
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)}
12345678910111213141516171819
h(x)
代码执⾏结果如下图所⽰。
拒绝采样（Reject Sampling ）
我们已经看到 Inver CDF ⽅法确实有效。但其实它的缺点也是很明显的，那就是有些分布的 CDF 可能很难通过对 PDF 的积分得到，再或者 CDF 的反函数也很不容易求。这时我们可能需要⽤到另外⼀种采样⽅法，这就是我们即将要介绍的拒绝采样。
下⾯这张图很好地阐释了拒绝采样的基本思想。假设我们想对 PDF 为 的函数进⾏采样，但是由于
种种原因（例如这个函数很复杂），对其进⾏采样是相对困难的。但是另外⼀个 PDF 为 的函数则相对容易采样，例如采⽤ Inver CDF ⽅法可以很容易对对它进⾏采样，甚⾄ 就是⼀个均匀分布（别忘了计算机可以直接进⾏采样的分布就只有均匀分布）。那么，当我们将 与⼀个常数
相乘之后，可以实现下图所⽰之关系，即 将 完全“罩住”。
然后重复如下步骤，直到获得 个被接受的采样点：
从 中获得⼀个随机采样点 对于 计算接受概率（acceptance probability ）
从 Uniform(0,1) 中随机⽣成⼀个值，⽤ 表⽰如果 ，则接受 作为⼀个来⾃ 的采样值，否则就拒绝 并回到第⼀步
你当然可以采⽤严密的数学推导来证明Reject Sampling 的可⾏性。但它的原理从直观上来解释也是相当容易理解的。你可以想象⼀下在上图的例⼦中，从哪些位置抽出的点会⽐较容易被接受。显然，红⾊曲线和绿⾊曲线所⽰之函数更加接近的地⽅接受概率较⾼，也即是更容易被接受，所以在这样的地⽅采到的点就会⽐较多，⽽在接受概率较低（即两个函数差距较⼤）的地⽅采到的点就会⽐较少，这也就保证了这个⽅法的有效性。
我们还是以本⽂最开始给出的那个分段函数 为例来演⽰ Reject Sampling ⽅法。如下⾯图所⽰，参
考分布我们选择的是均匀分布（你当然可以选择其他的分布，但采⽤均匀分布显然是此处最简单的处理⽅式）。⽽且令常数 。
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"))
1234567891011121314
p(x)q(x)q(x)q(x)M M ⋅q(x)p(x)m q(x)x(i)
x i α=
p()x i M q()时间过的快的成语
x i u α
≥u x i p(x)x i f(x)M中国传统麻将
=3
下⾯给出R 中的⽰例代码，易见我们此处采样点数⽬为10000。
同样，我们还是⽤图形来展⽰⼀下Reject Sampling 的效果如何。绘图的代码如下：
上述代码的执⾏结果如下图所⽰，可见采样结果是⾮常理想的。
最后给出的这个额外的例⼦来⾃⽂献[3]，这⾥同样采⽤uniform 分布作为参考分布。⽽需要被采样的分布则是我们在之前的⽂章《 》中介绍过的Beta 分布。如果你读过前⾯的博⽂，或者对机器学习⽐较了解，应该知道Beta 分布的PDF 表达式⾮常复杂。⽽且，你应该能够看出，跟前⾯不太⼀样的地
⽅是，这的 ，所取之值就是Beta 分布的极值。它的图形应该是与Beta(3,6)的极值点相切的⼀条⽔平直线。
f .x <- function(x) {
+ if (x >= 0 && x < 0.25)+ 8 * x
+ el if (x >= 0.25 && x <= 1)+ 8/3 - 8 * x/3+ el 0+ }
g .x <- function(x) {
+ if (x >= 0 && x <= 1)+ 1+ el 0+ }
M <- 3m <- 10000n .draws <- 0draws <- c()
x .grid <- q(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+ }+ }
1234567891011121314151617181920212223242526272829
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"))
1234567891011
M q(x)
