蒙特卡罗方法

19.1 蒙特卡罗法

19.1.1 随机抽样

统计学和机器学习的目的是基于数据对概率分布的特征进行推断，蒙特卡罗法要解决的问题是：假设概率分布的定义已知，通过抽样获得概率分布的随机样本，并通过得到的随机样本对概率分布的特征进行分析。比如，从样本得到经验分布，从而估计总体分布；或者从样本计算出样本均值，从而估计总体期望。所以蒙特卡罗法的核心是随机抽样（random sampling）。

一般的蒙特卡罗法有直接抽样法、接受-拒绝抽样法、重要性抽样法等。接受-拒绝抽样法、重要性抽样法适用于概率密度函数复杂（如密度函数含有多个变量，各变量相互不独立，密度函数形式复杂）的情况。

这里介绍接受-拒绝抽样法（accept-reject sampling method）。假设有机变量 $x$，取值 $x\in \mathcal{X}$，其概率密度函数为 $p(x)$。目标是得到该概率分布的随机样本，以对这个概率分布进行分析。

接受-拒绝法的基本想法如下。假设 $p(x)$ 不能直接抽样。找一个可以抽样的分布，称为建议分布（proposal distribution）。假设 $q(x)$ 是建议分布的概率密度函数，并且有 $q(x)$ 的 $c$ 倍一定大于等于 $p(x)$，其中 $c>0$，如图 19.1 所示。按照 $q(x)$ 进行抽样。假设抽样得到的结果是 $x^{\ast }$，再按照 $\frac{p(x^{\ast })}{cq(x^{\ast })}$ 的比例随机决定是否接受 $x^{\ast }$。直观上，落到 $p(x^{\ast })$ 范围内的就接受（绿色），落到 $p(x^{\ast })$ 范围外的就拒绝（红色）。接受-拒绝法实际是按照 $p(x)$ 的涵盖面积（或涵盖体积）占 $cq(x)$ 的涵盖面积（或涵盖体积）的比例进行抽样。

接受-拒绝抽样的具体算法如下。

算法 19.1（接受-拒绝法）

输入：抽样的目标概率分布的概率密度函数 $p(x)$。
输出：概率分布的随机样本 $x_1,x_2,\dots ,x_n$。
参数：样本数 $n$。

1. 选择概率密度函数为 $q(x)$ 的概率分布，作为建议分布，使其对任一 $x$ 满足 $cq(x)\ge p(x)$，其中 $c>0$。
2. 按照建议分布 $q(x)$ 随机抽样得到样本 $x^{\ast }$，再按照均匀分布在 (0, 1) 范围内抽样得到 $u$。
3. 如果 $u\le \frac{p(x^{\ast })}{cq(x^{\ast })}$，则将 $x^{\ast }$ 作为抽样结果；否则，回到步骤 (2)。
4. 直至得到 $n$ 个随机样本，结束。

接受-拒绝法的优点是容易实现，缺点是效率可能不高。如果 $p(x)$ 的涵盖体积占 $cq(x)$ 的涵盖体积的比例很低，就会导致拒绝的比例很高，抽样效率很低。注意，一般是在高维空间进行抽样，即使 $p(x)$ 与 $cq(x)$ 很接近，两者涵盖体积的差异也可能很大（与我们在三维空间的直观不同）。

19.1.2 数学期望估计

一般的蒙特卡罗法如直接抽样法、接受-拒绝抽样法、重要性抽样法，也可以用于数学期望估计（estimation of mathematical expectation）。假设有随机变量 $x$，取值 $x\in \mathcal{X}$，其概率密度函数为 $p(x)$，$f(x)$ 为定义在 $\mathcal{X}$ 上的函数，目标是求函数 $f(x)$ 关于概率密度函数 $p(x)$ 的数学期望 $E_{p(x)}[f(x)]$。

针对这个问题，蒙特卡罗法按照概率分布 $p(x)$ 独立地抽取 $n$ 个样本 $x_1,x_2,\dots ,x_n$，比如用以上的抽样方法，之后计算函数 $f(x)$ 的样本均值 ${\hat{f}}_n$：
$${\hat{f}}_n=\frac{1}{n}\sum _{i=1}^{n}f(x_i)\text{\hspace{1em}(19.1)}$$

作为数学期望 $E_{p(x)}[f(x)]$ 的近似值。

根据大数定律，当样本容量增大时，样本均值以概率 1 收敛于数学期望：
$${\hat{f}}_n\to E_{p(x)}[f(x)],n\to \infty \text{\hspace{1em}(19.2)}$$

这样就得到了数学期望的近似计算方法：
$$E_{p(x)}[f(x)]\approx \frac{1}{n}\sum _{i=1}^{n}f(x_i)\text{\hspace{1em}(19.3)}$$

19.1.3 积分计算

一般的蒙特卡罗法也可以用于定积分的近似计算，称为蒙特卡罗积分（Monte Carlo integration）。假设有一个函数 $h(x)$，目标是计算该函数的积分：
$${\int }_{\mathcal{X}}h(x)dx$$

如果能够将函数 $h(x)$ 分解成一个函数 $f(x)$ 和一个概率密度函数 $p(x)$ 的乘积的形式，那么就有
$${\int }_{\mathcal{X}}h(x)dx={\int }_{\mathcal{X}}f(x)p(x)dx=E_{p(x)}[f(x)]\text{\hspace{1em}(19.4)}$$

于是函数 $h(x)$ 的积分可以表示为函数 $f(x)$ 关于概率密度函数 $p(x)$ 的数学期望。实际上，给定一个概率密度函数 $p(x)$，只要取 $f(x)=\frac{h(x)}{p(x)}$，就可得式 (19.4)。也就是说，任何一个函数的积分都可以表示为某一个函数的数学期望的形式，而函数的数学期望又可以通过函数的样本均值估计。于是，就可以利用样本均值来近似计算积分，这就是蒙特卡罗积分的基本想法。
$${\int }_{\mathcal{X}}h(x)dx=E_{p(x)}[f(x)]\approx \frac{1}{n}\sum _{i=1}^{n}f(x_i)\text{\hspace{1em}(19.5)}$$

例 19.1

用蒙特卡罗积分法求 ${\int }_0^1{e}^{-x^2/2}dx$。

解：令 $f(x)=e^{-x^2/2}$，

$$p(x)=1(0<x<1)$$

也就是说，假设随机变量 $x$ 在 (0, 1) 区间遵循均匀分布。

使用蒙特卡罗积分法，如图 19.2 所示，在 (0, 1) 区间按照均匀分布抽取 10 个随机样本 $x_1,x_2,\dots ,x_{10}$。计算样本的函数均值 ${\hat{f}}_{10}$：

$${\hat{f}}_{10}=\frac{1}{10}\sum _{i=1}^{10}{e}^{-x_i^2/2}=0.832$$

也就是积分的近似。随机样本数越大，计算就越精确。

例 19.2

用蒙特卡罗积分法求 ${\int }_{-\infty }^{\infty }x\frac{1}{\sqrt{2\pi }}\exp \left(-\frac{x^2}{2}\right)dx$。

解：令 $f(x)=x$，

$$p(x)=\frac{1}{\sqrt{2\pi }}\exp \left(-\frac{x^2}{2}\right)$$

$p(x)$ 是标准正态分布的密度函数。

使用蒙特卡罗积分法，按照标准正态分布在区间 $(-\infty ,\infty )$ 抽样 $x_1,x_2,\dots ,x_n$，取其平均值，就得到要求的积分值。当样本增大时，积分值趋于 0。

本章介绍的马尔可夫链蒙特卡罗法也适用于概率密度函数复杂、不能直接抽样的情况，旨在解决一般的蒙特卡罗法，如接受-拒绝抽样法、重要性抽样法，抽样效率不高的问题。一般的蒙特卡罗法中的抽样样本是独立的，而马尔可夫链蒙特卡罗法中的抽样样本不是独立的，样本序列形成马尔可夫链。

---

例19.1 p(x)=1(0＜x＜1) 公式解析
例19.1 计算蒙特卡洛积分的具体过程
例 19.2：用蒙特卡洛积分法求解具体过程