随机采样之接受拒绝采样

创作时间:

作者:

@小白创作中心

随机采样之接受拒绝采样

引用

CSDN

https://blog.csdn.net/u011426016/article/details/143592888

逆变换采样（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 xxx，有f ( x ) ≤ M ⋅ g ( x ) f(x) \leq M \cdot g(x)f(x)≤M⋅g(x)，其中f ( x ) f(x)f(x)是目标分布，M MMM是一个正常数。
抽样：从提议分布g ( x ) g(x)g(x)中抽取样本x xxx和从均匀分布U ( 0 , 1 ) U(0, 1)U(0,1)中抽取样本u uuu。
接受-拒绝条件：如果u ≤ f ( x ) M ⋅ g ( x ) u \leq \frac{f(x)}{M \cdot g(x)}u≤M⋅g(x)f(x) ，则接受x xxx作为目标分布f ( x ) f(x)f(x)的一个样本；否则拒绝x xxx。

接受拒绝采样可以使用下图进行表示（图片来源：【数之道】马尔可夫链蒙特卡洛方法是什么？十五分钟理解这个数据科学难点）。

二、接受拒绝采样证明

要证明接收-拒绝采样确实产生服从目标分布f ( x ) f(x)f(x)的样本，我们需要证明对于所有的x xxx，有：

P ( X = x ) = f ( x ) (1) P(X=x) = f(x)\tag1P(X=x)=f(x)(1)

其中P ( X = x ) P(X=x)P(X=x)是样本x xxx被接受的概率。

证明：

接受概率：样本x xxx被接受的概率是f ( x ) M ⋅ g ( x ) \frac{f(x)}{M \cdot g(x)}M⋅g(x)f(x) ，因为u uuu是从U ( 0 , 1 ) U(0, 1)U(0,1)中抽取的。
联合概率：样本x xxx从提议分布g ( x ) g(x)g(x)中抽取的概率是g ( x ) g(x)g(x)，并且u uuu在[ 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}\tag2P(X=x,U≤M⋅g(x)f(x) )=g(x)⋅M⋅g(x)f(x) =Mf(x) (2)

边缘概率：现在我们需要计算X XXX的边缘概率P ( X = x ) P(X=x)P(X=x)，即样本x xxx被接受的总概率。由于u uuu是均匀分布的，我们可以将联合概率在u uuu的所有可能值上积分：

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}\tag3P(X=x)=∫01 P(X=x,U=u)du=∫01 Mf(x) du=Mf(x) ⋅∫01 du=Mf(x) (3)

归一化常数：由于M MMM是使得f ( x ) ≤ M ⋅ g ( x ) f(x) \leq M \cdot g(x)f(x)≤M⋅g(x)对所有x xxx成立的最小常数，我们可以将上式中的M MMM移到f ( x ) f(x)f(x)的定义中，从而得到：

P ( X = x ) = f ( x ) (4) P(X=x) = f(x)\tag4P(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)\tag5f(x)=0.3exp(−(x−0.3)2)+0.7exp(−(x−2)2/0.3)(5)

其归一化常数为Z = 1.2113 Z = 1.2113Z=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.5M=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的图像，蓝色区域是采样结果，可见采样结果跟真实分布几乎一致。