之前提到的逆变换采样(Inverse Transform Sampling)是一种生成随机样本的方法,它利用累积分布函数(CDF)的逆函数来生成具有特定分布的随机变量。以下是逆变换采样的缺点:
- 计算复杂性:对于某些分布,找到累积分布函数(CDF)的逆函数可能是困难的,甚至是不可能的。
- 效率问题:对于具有重尾分布的随机变量,逆变换采样可能非常低效,因为CDF的逆可能需要大量的计算。
- 数值稳定性:在数值计算中,由于浮点数的精度限制,逆变换采样可能会引入误差,尤其是在CDF的值接近1时。
一、接受拒绝采样
接受-拒绝采样(Accept-Reject Sampling)方法是一种更为通用的采样方法,它可以用来生成具有任意分布的随机样本。这种方法不要求我们知道CDF的逆,而是利用一个简单的概率分布(称为提议分布)来生成样本,然后以一定的概率接受或拒绝这些样本。
接收-拒绝采样的基本步骤:
-
选择提议分布 g ( x ) g(x) g(x):选择一个容易从中抽样的分布 g ( x ) g(x) g(x),并且确保对于所有的 x x x,有 f ( x ) ≤ M ⋅ g ( x ) f(x) \leq M \cdot g(x) f(x)≤M⋅g(x),其中 f ( x ) f(x) f(x)是目标分布, M M M是一个正常数。
-
抽样:从提议分布 g ( x ) g(x) g(x)中抽取样本 x x x和从均匀分布 U ( 0 , 1 ) U(0, 1) U(0,1)中抽取样本 u u u。
-
接受-拒绝条件:如果 u ≤ f ( x ) M ⋅ g ( x ) u \leq \frac{f(x)}{M \cdot g(x)} u≤M⋅g(x)f(x),则接受 x x x作为目标分布 f ( x ) f(x) f(x)的一个样本;否则拒绝 x x x。
接受拒绝采样可以使用下图进行表示(图片来源:【数之道】马尔可夫链蒙特卡洛方法是什么?十五分钟理解这个数据科学难点)。
二、接受拒绝采样证明
要证明接收-拒绝采样确实产生服从目标分布
f
(
x
)
f(x)
f(x)的样本,我们需要证明对于所有的
x
x
x,有:
P
(
X
=
x
)
=
f
(
x
)
(1)
P(X=x) = f(x)\tag1
P(X=x)=f(x)(1)
其中 P ( X = x ) P(X=x) P(X=x)是样本 x x x被接受的概率。
证明:
-
接受概率:样本 x x x被接受的概率是 f ( x ) M ⋅ g ( x ) \frac{f(x)}{M \cdot g(x)} M⋅g(x)f(x),因为 u u u是从 U ( 0 , 1 ) U(0, 1) U(0,1)中抽取的。
-
联合概率:样本 x x x从提议分布 g ( x ) g(x) g(x)中抽取的概率是 g ( x ) g(x) g(x),并且 u u u在 [ 0 , f ( x ) M ⋅ g ( x ) ) [0, \frac{f(x)}{M \cdot g(x)}) [0,M⋅g(x)f(x))区间的概率是 f ( x ) M ⋅ g ( x ) \frac{f(x)}{M \cdot g(x)} M⋅g(x)f(x)。因此,联合概率是:
P ( X = x , U ≤ f ( x ) M ⋅ g ( x ) ) = g ( x ) ⋅ f ( x ) M ⋅ g ( x ) = f ( x ) M (2) P(X=x, U \leq \frac{f(x)}{M \cdot g(x)}) = g(x) \cdot \frac{f(x)}{M \cdot g(x)} = \frac{f(x)}{M}\tag2 P(X=x,U≤M⋅g(x)f(x))=g(x)⋅M⋅g(x)f(x)=Mf(x)(2)
-
边缘概率:现在我们需要计算 X X X的边缘概率 P ( X = x ) P(X=x) P(X=x),即样本 x x x被接受的总概率。由于 u u u是均匀分布的,我们可以将联合概率在 u u u的所有可能值上积分:
P ( X = x ) = ∫ 0 1 P ( X = x , U = u ) d u = ∫ 0 1 f ( x ) M d u = f ( x ) M ⋅ ∫ 0 1 d u = f ( x ) M (3) P(X=x) = \int_0^1 P(X=x, U=u) \, du = \int_0^1 \frac{f(x)}{M} \, du = \frac{f(x)}{M} \cdot \int_0^1 du = \frac{f(x)}{M}\tag3 P(X=x)=∫01P(X=x,U=u)du=∫01Mf(x)du=Mf(x)⋅∫01du=Mf(x)(3)
-
归一化常数:由于 M M M是使得 f ( x ) ≤ M ⋅ g ( x ) f(x) \leq M \cdot g(x) f(x)≤M⋅g(x)对所有 x x x成立的最小常数,我们可以将上式中的 M M M移到 f ( x ) f(x) f(x)的定义中,从而得到:
P ( X = x ) = f ( x ) (4) P(X=x) = f(x)\tag4 P(X=x)=f(x)(4)
这就证明了接收-拒绝采样确实产生了服从目标分布 f ( x ) f(x) f(x)的样本。
三、接受拒绝采样模拟
借用作者anshuai_aw1的例子,设我们需要采样的pdf为:
f
(
x
)
=
0.3
exp
(
−
(
x
−
0.3
)
2
)
+
0.7
exp
(
−
(
x
−
2
)
2
/
0.3
)
(5)
f(x)=0.3 \exp \left(-(x-0.3)^{2}\right)+0.7 \exp \left(-(x-2)^{2} / 0.3\right)\tag5
f(x)=0.3exp(−(x−0.3)2)+0.7exp(−(x−2)2/0.3)(5)
其归一化常数为
Z
=
1.2113
Z = 1.2113
Z=1.2113, 参考分布为
g
(
x
)
=
N
(
μ
=
1.4
,
σ
2
=
(
1.
2
2
)
)
g(x) =N(\mu=1.4,\sigma^2=(1.2^2))
g(x)=N(μ=1.4,σ2=(1.22)),
M
=
2.5
M=2.5
M=2.5, 以确保
M
⋅
g
(
x
)
≥
f
(
x
)
M \cdot g(x) \geq f(x)
M⋅g(x)≥f(x)。采样的代码如下:
import numpy as np
import matplotlib.pyplot as plt
def f(x):
return (0.3*np.exp(-(x-0.3)**2) + 0.7* np.exp(-(x-2.)**2/0.3))/1.2113
x = np.arange(-4.,6.,0.01)
plt.plot(x,f(x),color = "red")
size = int(1e+07)
mu = 1.4
sigma = 1.2
M = 2.5
x = np.random.normal(loc = mu,scale = sigma, size = size)
g_x = 1/(np.sqrt(2*np.pi)*sigma)*np.exp(-0.5*(x-mu)**2/sigma**2)
u = np.random.uniform(low = 0, high = M*g_x, size = size) #在[0,M*g_x]中均匀采样
fx = 0.3*np.exp(-(x-0.3)**2) + 0.7* np.exp(-(x-2.)**2/0.3)
sample = x[u <= fx] # u < fx(x)
plt.hist(sample,bins=150, density=True, edgecolor='black')
plt.show()
结果如下,其中红色曲线的是公式(5)所示pdf的图像,蓝色区域是采样结果,可见采样结果跟真实分布几乎一致。
参考资料:
[1]【数之道】马尔可夫链蒙特卡洛方法是什么?十五分钟理解这个数据科学难点
[2] 逆采样(Inverse Sampling)和拒绝采样(Reject Sampling)原理详解