缘起

如果一个函数 f(x)f(x)f(x) 在整个坐标轴上的积分有界：
∫−∞∞f(x)dx=M<∞\int^{\infty}_{-\infty} f(x) dx = M < \infty ∫−∞∞f(x)dx=M<∞
那么就可以用它来定义一个分布 f^\hat{f}f^：
f^(x)≡f(x)M⇒∫−∞∞f^(x)dx=∫−∞∞f(x)dxM=1\begin{aligned} & \hat{f}(x) \equiv \frac{f(x)}{M} \\ \Rightarrow & \int^{\infty}_{-\infty} \hat{f}(x) dx =\frac{\int^{\infty}_{-\infty} f(x) dx}{M} = 1 \end{aligned} ⇒f^(x)≡Mf(x)∫−∞∞f^(x)dx=M∫−∞∞f(x)dx=1
如何从这个以 f^\hat{f}f^ 为概率密度函数的分布中采样呢？
很多时候由于f(x)f(x)f(x)的复杂性，我们甚至不能计算出上述的归一化因子MMM，那么如何采样，可以使得样本的分布正比于 f(x)f(x)f(x) 呢？

举个例子，
f(x)=e−∣x−4∣+1(1+x2)2+log⁡(1+∣x∣)(1+∣x+2∣)2f(x) = e^{-|x-4|} + \frac{1}{(1+x^2)^2} + \frac{\log(1+|x|)}{(1+|x+2|)^2} f(x)=e−∣x−4∣+(1+x2)21+(1+∣x+2∣)2log(1+∣x∣)
形状有3个尖峰：

import numpy as np
import matplotlib.pyplot as pltdef f(x):return np.exp(-np.abs(x-4)) + 1/np.square(1 + np.square(x))+ np.log(1+np.abs(x))/(1+np.abs(x+2))**2x = np.linspace(-100,100,1000)
y = f(x)
plt.plot(x,y)
plt.xlim([-20,20])

大家感兴趣可以试试把积分即归一化常数 MMM 算出来
为了方便之后方案的对比验证，我这里直接用 scipy 的数值积分函数，计算出该分布的期望为：0.1575

def expectedValue(f):from scipy import integrateM, error = integrate.quad(f, -np.inf, np.inf)assert(error < 1e-6)expected, error = integrate.quad(lambda x: x*f(x), -np.inf, np.inf)assert(error < 1e-6)expected /= Mreturn expected, MEV, M0 = expectedValue(f)
# 期望，归一化常数
# (0.1574930705005623, 6.317327048449105)

拒绝采样

接收-拒绝算法是利用符合某个分布的函数 ggg 去模拟目标概率密度函数 fff，且 f/gf/gf/g 是有界的，即：
f(x)/g(x)<M<+∞,∀xf(x)/g(x)< M < +\infty, \quad \forall x f(x)/g(x)<M<+∞,∀x

标准的接收-拒绝算法为：¹

从ggg 中生成样本XiX_{i}Xi
从 [0,1][0,1][0,1] 区间的均匀分布生成 UiU_{i}Ui
如果 Ui≤f(Xi)/Mg(Xi)U_{i}\leq f(X_{i})/Mg(X_{i})Ui≤f(Xi)/Mg(Xi)，则接收 XiX_{i}Xi 为函数的采样点
如果不满足上述条件，则 i+1i+1i+1

例如，我们选用正态分布作为 g(x)∼N(0,52)g(x) \sim \mathbb{N}(0, 5^2)g(x)∼N(0,52)，且选择 M=20M = 20M=20

from scipy.stats import norm
colors = [ "#348ABD", "#A60628", "#7A68A6", ]g = norm(0, 5)
M = 20plt.plot(x, M *g.pdf(x), label='Mg(x)', color=colors[1])
plt.fill_between(x, M *g.pdf(x), color=colors[1], alpha=.33)
plt.plot(x, y, label='f(x)', color=colors[0])
plt.fill_between(x, y, color=colors[0], alpha=.33)
plt.xlim([-20,20])
plt.legend()
plt.show()

虽然对于主要区域，Mg(x)Mg(x)Mg(x) 可以将 f(x)f(x)f(x) 覆盖住，但事实上，并不存在有限的MMM使得 f<Mgf < Mgf<Mg恒成立，因为正态分布曲线在两侧的下降速度要快于f(x)f(x)f(x) 的，可以看看上图的局部放大情况：

下面为拒绝采样的代码实现

code1

按照基本逻辑实现

def rejectSampling_slow(f, g, M, N_try=1):samples = []for i in range(N_try):x = g.rvs()u = np.random.random()if u < f(x) / M*g.pdf(x):samples.append(x)return samples

code2

使用 numpy array向量化加速，效率很高哦！

def rejectSampling(f, g, M, N_try=1):samples = g.rvs(N_try)U = np.random.random(N_try)return samples[U < f(samples) / M*g.pdf(samples)]

采样并计算样本均值、样本接受率

N_try = 100000000
samples = rejectSampling(f, g, M, N_try)
# 计算样本均值、样本接受率
sample_mean, accept_rate = np.mean(samples), len(samples)/N_try
# (0.26507388249562314, 0.00126865)

进行了 10810^8108 次采样，样本接收率仅为 0.1%，而且样本均值 0.2651 与理论期望值 0.1575 相差挺大

来看看样本的密度曲线与期望分布的差距：

# _ = plt.hist(samples, bins=100, density=True, label='sample density', color=colors[0])
density = gaussian_kde(samples)
plt.plot(x, density(x), label='sample kde', color=colors[1])
plt.fill_between(x, density(x), color=colors[1], alpha=.33)
plt.plot(x, y/M0, label=r'f(x)/$M_0$', color=colors[0])
plt.fill_between(x, y/M0, color=colors[0], alpha=.33)
plt.xlim([-20,20])
plt.legend()
plt.show()

可以看出，实际采集的样本比预期分布得更集中，主要原因为：在偏离中心（原点）的两侧区域，实际采集的样本比例要低于理论值，导致实际的样本密度曲线向中心集中。

拒绝采样的缺点

看完上述案例，想必就能理解拒绝采样的缺点：要找到一个能完全覆盖函数 f(x)f(x)f(x) 且样本拒绝率低的函数 Mg(x)Mg(x)Mg(x) 是很难的

要想在整个坐标轴上覆盖住f(x)f(x)f(x)，一般会导致常数 MMM 选取得很大，而一旦 MMM 很大，会导致整体上样本拒绝率高，有效采样率很低。

MCMC

拒绝采样先从分布 g(x)g(x)g(x) 中采样，每次采样都是独立进行的。试想，如果某一次采样 XiX_iXi 没有被拒绝，意味着 fff 在XiX_iXi 附近很可能是相对较大的，那么如果继续在这附近采样，那被拒绝的概率应当也不会太高，这样每次都从上一次样本的附近采样是否就能真题提高样本的接收率呢？

这就是 MCMC 的出发点，每次采样都是基于上一步采样结果进行的，且待下回分解。

参考资料

https://www.jiqizhixin.com/graph/technologies/c84b6d7e-447f-4bcb-b4a5-c807d7b8a5f7 ↩︎

从拒绝采样（Reject Sampling）到马尔科夫链蒙特卡洛（MCMC）相关推荐

蒙特卡洛分析_随机模拟：马尔科夫链蒙特卡洛采样MCMC与EM算法「2.3」
最近学习了机器学习中的马尔科夫链蒙特卡洛(Markov Chain Monte Carlo, 简称MCMC) 相关的知识. 主要内容包括: [1]蒙特卡洛原则,及其应用于采样的必要性(已经发布在头条) ...
MCMC原理解析(马尔科夫链蒙特卡洛方法)
马尔科夫链蒙特卡洛方法(Markov Chain Monte Carlo),简称MCMC,MCMC算法的核心思想是我们已知一个概率密度函数,需要从这个概率分布中采样,来分析这个分布的一些统计特性,然而 ...
马尔科夫链蒙特卡洛_蒙特卡洛·马可夫链
马尔科夫链蒙特卡洛 A Monte Carlo Markov Chain (MCMC) is a model describing a sequence of possible events wher ...
第十五课.马尔科夫链蒙特卡洛方法
目录 M-H采样 Metropolis-Hastings采样原理 M-H采样步骤 Gibbs方法 Gibbs核心流程 Gibbs采样的合理性证明 Gibbs采样实验在第十四课中讲述了马尔科夫链与其 ...
马尔科夫链和马尔科夫链蒙特卡洛方法
前言译自:<Training Restricted Boltzmann Machines: An Introduction > 马尔科夫链在RBM的训练中占据重要地位,因为它提供了从复杂 ...
matlab mcmc工具箱,马尔科夫链蒙特卡洛模拟（MCMC）matlab工具箱
马尔科夫链蒙特卡洛模拟(MCMC)matlab工具箱 matlab 2021-2-10 下载地址 https://www.codedown123.com/64660.html 马尔科夫链蒙特卡洛模拟( ...
马尔科夫链蒙特卡洛（MCMC）
在以贝叶斯方法为基础的机器学习技术中,通常需要计算后验概率,然后通过最大后验概率(MAP)等方法进行参数推断和决策.然而,在很多时候,后验分布的形式可能非常复杂,这个时候寻找其中的最大后验估计或者对后 ...
马尔科夫链与MCMC方法
马尔科夫链概述基本思想: 过去所有的信息都已经被保存到了现在的状态,基于现在就可以预测未来. Example: 假如每天的天气是一个状态的话,那今天是不是晴天只依赖于昨天的天气,而和前天的天气没有任 ...
Python用MCMC马尔科夫链蒙特卡洛、拒绝抽样和Metropolis-Hastings采样算法
最近我们被客户要求撰写关于MCMC的研究报告,包括一些图形和统计输出. 我们将研究两种对分布进行抽样的方法:拒绝抽样和使用 Metropolis Hastings 算法的马尔可夫链蒙特卡洛方法 (MC ...

从拒绝采样（Reject Sampling）到马尔科夫链蒙特卡洛（MCMC）

文章目录

缘起