直接采样、拒绝采样与重要性采样

【摘要】蒙特卡洛(Monte Carlo method)是一种以概率统计理论为指导的重要数值计算方法。它使用随机数来解决随机变量(或随机函数)的期望值积分求解、仿真模拟等非常棘手的计算问题,特别适用于无解析形式的复杂概率分布。根据对蒙特卡洛方法的理解,会发现其中最为核心的部分是如何在给定一个复杂分布时,按照概率随机、高效地获得样本,即采样方法问题。本文将介绍其中最为基础和直觉的几种早期方法,分别是基于 CDF 的直接采样、拒绝采样和重要性采样。

1 直接采样

直接采样的思想是:计算机适合于随机的均匀采样,如果能够把任意概率分布的采样转化成对均匀分布的采样,就可以解决采样问题。

假设 yy 服从某项分布 p(y)p(y),其累积分布函数( CDF )为 h(y)h(y),现有均匀分布的样本 zUniform(0,1)z \sim \operatorname{Uniform}(0,1),令 z=h(y)z = h(y),即 y=h1(z)y = h^{-1}(z),结果 yy 即为对分布 p(y)p(y) 的采样。

图 1: 直接采样算法流程

举个例子:

图 2:直接采样示例

直接采样要求分布的累积分布函数 h(x)h(x) 及其拟函数存在解析解,因此只适用于形式已知、具有解析形式的概率分布,对于形态复杂的概率密度函数(或概率质量函数)无能为力。

2 接受-拒绝采样

接受-拒绝采样的思想是:对于形态比较复杂的分布 p(z)p(z) ,很难进行直接采样时,可以考虑构建一个形式更为简单的、能够进行直接采样的新分布 q(z)q(z),然后通过对 q(z)q(z) 的采样结果进行筛选,来获得 p(z)p(z) 的样本。那么我们怎么可以通过 q(z)q(z) 采样得到 p(z)p(z) 的样本呢?

先看下图:

图 3:接收-拒绝采样原理示意图

红色的复杂分布是 p(z)p(z) , 蓝色的简单分布是 q(z)q(z),我们对 q(z)q(z) 乘一个参数 kk,使得 kq(z)kq(z) 正好能包住p(z)p(z),那么对于每一个采用直接方法从 q(z)q(z) 得到的样本 z0z_0,我们采用一定的概率去接受它,概率大小为 p(z0)/kq(z0)p(z_0) / kq(z_0)。很容易就能看出来,在 p(z)p(z)kq(z)kq(z) 相切的地方采样,接受概率为 11。可能有人会问,接受概率是统计结果,虽然有理论设计,但对于某个具体的样本 z0z_0,该如何判断接受还是拒绝呢?我们有 uUniform[0,1]u \sim Uniform\left[0,1\right],对于样本 z0z_0,我们抽取一个样本 u0u_0,如果 u0p(z0)/kq(z0)u_0 \leq p(z_0) / kq(z_0) 就接受,否则拒绝。重复此过程,得到的样本就服从分布 p(z)p(z)

图 4:接受-拒绝采样算法流程

当然,q(z)q(z) 的选取要有一定规则:一是 q(z)q(z)p(z)p(z) 的外形要相近,二是 q(z)q(z) 采样比较容易。

现在有一个问题,kk 怎么求?

其实也很简单,看图有 kq(z)p(z)kq(z) \geq p(z),那么 kp(z)/q(z)k \geq p(z) / q(z),我们只需要求 p(z)/q(z)p(z) / q(z) 的最大值,即可以确定 kk

我们以截断正态分布的接受-拒绝采样为例来加以说明。 截断正态分布的意思就是对于正态分布 N(a,b)N(a, b)xx 属于 [0,4]\left[0, 4\right],其在 [0,4]\left[0, 4\right] 上的积分为 11,而不是在负无穷到正无穷上的积分为 11。截断正态分布不是正态分布。维基上对截断正态分布给出了如下定义:

图 5:截断正态分布的维基百科定义

截断正态分布定义式中的分子中,ϕ\phi 是标准正态分布的概率密度函数,分母中的 Φ\Phi 是标准正态分布的累积分布函数。我们有 p(z)p(z) 服从 N(1,1)N(1, 1)(0<=x<=4)(0 <= x <= 4),令 q(z) U[0,4]q(z)~U\left[0, 4\right],根据上图公式,p(z)p(z) 就有了,q(z)=1/4q(z) = 1/4,所以 k=max(p(z)/q(z))k = max(p(z) / q(z)),在 z=1z = 1(均值)的时候 p(z)/q(z)p(z) / q(z) 取最大,所以得到kk

这个例子比较巧,分母没有 zz,我们可以直接判断在 z=1z=1 时,p(z)/q(z)p(z)/q(z) 取最大,如果 p(z)p(z), q(z)q(z) 都有zz,那么就需要通过求导方式求 kk 了。

3 重要性采样

重要性采样与 “接受-拒绝采样” 有异曲同工之妙。接受-拒绝采样时,通过接受/拒绝的方式对通过 q(z)q(z) 得到的样本进行筛选,使得最后得到的样本服从分布 p(z)p(z),所有接受的样本没有高低贵贱之分,一视同仁。而重要性采样的思想是:对于通过 q(z)q(z) 得到的样本,我全部接受!全部接受的话会有一个问题:最后样本点的分布不服从 p(z)p(z) 分布。为了矫正此偏差,我们给每个样本赋予一个重要性权重,例如,对于 p(z0)/q(z0)=1p(z_0)/q(z_0)=1 的样本,权重大一点,p(z0)/q(z0)=0.1p(z_0)/q(z_0)=0.1 的样本,权重小一点。但是这个权重怎么算呢?其实,p(z0)/q(z0)p(z0)/q(z0) 就是一个很好的权重标识。

重要性采样取得到的是带有重要性权重的、服从 q(z)q(z) 分布的样本,这个权重乘以样本之后的结果其实就是服从 p(z)p(z) 分布的。对于这种特殊的样本,我们怎么用呢?通常我们基于分布来求均值等积分量。

图 6:重要性采样的使用

应当注意的是,图中 p(z)p(z)f(z)f(z) 的关系,p(z)p(z) 是一种分布,是相对于 zz 轴的采样点而言的,比如在红色的两个峰值处,zz 的取点比较多,在其他地方 zz 的取点就比较少,这叫样本分布服从 p(z)p(z)。而 f(z)f(z) 通常是我们希望得到期望值的某种映射,将 zz 值映射到其他维度。比如我们熟悉的 y=f(x)y = f(x),将 xx 映射到 yy, 我们期望的到 yy 的期望(即 f(x)f(x) 的均值)。

对比接受-拒绝采样和重要性采样两种方法,前者按照 p(z)p(z) 的概率密度来控制不同区域样本的数量,进而直接为蒙特卡洛积分提供有效样本;而后者则在不同区域按照 q(z)q(z) 的概率密度控制不同区域的样本数量,在蒙特卡洛积分时,通过重要性系数对冲样本数量带来的误差。