MCMC：伪随机数采样与近似积分

伪随机数采样

假设计算机能够生成服从于$[0,1]$ 均匀分布的随机变量/随机数（random numbers/variables, r.v.'s），那么所谓生成随机样本（generating random sample）的任务就是根据上述假设寻找一个快速的生成其他非均匀分布的随机变量（non-uniform r.v.'s）的方法。该任务也被称之为采样或伪随机数采样（pseudo-random number sampling）。

逆变换采样方法

逆变换采样（inverse transform sampling）是一种给定想要生成的随机变量$X$ 的累计分布函数（Cumulative Distribution Function, CDF）$F_X$ 后，通过对其进行逆变换（反函数）实现随机样本生成的方法。其中$F_X$ 是一个连续且严格单调递增的函数。

具体来说，我们需要找到一个变换/映射$T:[0,1]\mapsto\mathbb R$ 使得$U\sim\text{Unif}[0,1]$ 变换为服从$F_X$ 的新随机变量$X$，即$T(U)=X$。不难发现，有：

$$\begin{aligned} F_X(x)&=\Pr(X\leq x)=\Pr(T(U)\leq x)\\ &=\Pr(U\leq T^{-1}(x))\\ &=T^{-1}(x) \end{aligned}$$

也就是说，根据任意一个随机数$u\in[0,1]$ ，我们可以生成另一个随机数$x=T(u)=F_X^{-1}(u)$，其中$F_X^{-1}(u):=\inf\{x|F_X(x)\geq u\}$，$\inf$ 表示 infimum，即下确界。这么定义该反函数的原因是 CDF 是单调且右连续的。

Example #1

指数分布

标准拉普拉斯分布

逆变换采样网格法

对于$X\not\in\mathbb R$ 的随机变量$X\sim F_X(x)$，设其 support 为$\mathcal S_X$ ，即所有可能的取值的集合。为了生成$X$，我们可以选择适当的网格$\{x_i\}^d_{i=1}$ 将$\mathcal S_X$ 覆盖。
则$X$ 的概率密度函数（Probability Density Function, PDF）可以使用离散随机变量$\{x_i\}^d_{i=1}$ 的离散概率分布进行近似：

$$p_i=\frac{f_X(x_i)}{\sum_{j=1}^df_X(x_j)}\qquad i=1,2,...,d$$

进而采用离散分布的 CDF 使用 inverse sampling 生成。

显然，The Grid Method 同样可以适用于一些非标准的分布函数，乃至部分很难求“反”的函数。

Example #2

截断正态分布

【文献阅读 02】截断(truncated)正态分布 - 知乎 (zhihu.com)

接受-拒绝采样方法

Inverse CDF 采样方法固然有效，但对于那些 CDF 很难通过对 PDF 的积分得到，或者 CDF 的反函数不容易求得的分布，我们可能需要用到另外一种采样方法。而接受-拒绝采样方法（Acceptance-Rejection sampling），通常也直接称为拒绝采样（Rejection sampling）就是一种可以解决上述问题的方法之一。

参考上图，假设$p(x)$ 是目标采样分布$X\sim p(\cdot)$ 的 PDF，它在自身的支撑集$\mathcal S_X$ 内可以被另一个 PDF 为$q(x)$ 的分布$Y\sim q(\cdot)$ 包络起来（$Y$ 的分布称为 提议分布 proposal distribution）。即：

$$p(x)\leq k·q(x)\qquad \forall x\in\mathcal S_X$$

其中$k\geq1$ 为包络常数（envelope constant），$q(x)$ 是提议分布的 PDF，也称为包络密度函数（envelope density）。

Rejection sampling 的算法步骤如下：

1. 独立地生成两个随机变量$U\sim\text{U}[0,1]$ 和$Y\sim q(\cdot)$；
2. 如果$U\leq p(Y)/kq(Y)$ 则$Y$ 即为所要的服从于$p(\cdot)$ 的随机变量，接受$Y$，返回$X=Y$；否则拒绝$Y$ 并返回到第一步。

感性认知

从图像上来理解，$p(z_0)/kq(z_0)$ 越趋近于1，说明此时$p(z_0)$ 和$kq(z_0)$ 的距离越小，图中竖直线段上处于灰色阴影部分的线段长度就越短，此时在线段上进行均匀采样，就有极大的概率使采样点$u_0$ 落在非阴影区域，以该概率判断是接受还是拒绝$z_0$ 。

显然 rejection sampling 是符合直觉的认知策略，下面还将在理论上进行推导验证该方法的可行性。

理论证明

令$h(Y)=p(Y)/kq(Y)$ ，从而：

$$h(x)\cdot q(x)=\frac1k\cdot p(x)$$

根据贝叶斯定理，有：

$$\begin{aligned} f_Y(x|U\leq h(Y))&=\frac{\Pr\{U\leq h(Y)|Y=x\}\cdot q(x)}{\Pr\{U\leq h(Y)\}}\\ &=\frac{\Pr\{U\leq h(x)\}\cdot q(x)}{\Pr\{U\leq h(Y)\}}\\ &=\frac{h(x)\cdot q(x)}{\Pr\{U\leq h(Y)\}} \end{aligned}$$

而分母：

$$\begin{aligned} \Pr\{U\leq h(Y)\}&=\int_{\mathcal S_X}\Pr\{U\leq h(Y)|Y=x\}\cdot q(x)\;\mathrm dx\\ &=\int_{\mathcal S_X} h(x)\cdot q(x)\;\mathrm dx\\ &=\frac1k\int_{\mathcal S_X}p(x)\;\mathrm dx\\ &=\frac1k \end{aligned}$$

综上可得，$f_Y(x|U\leq h(Y))=p(x)$。

也就是说，在$U\leq h(Y)$ 的条件下，如果一个随机变量$x$ 服从于$f_Y(\cdot)=q(\cdot)$ ，那么它也服从于$p(\cdot)$。

算法有效性

从图像上来看，过大的$k$ 值会让$kq(x)$ 的图像远在$p(x)$ 上方，阴影部分最多，拒绝的概率也就增大了。此外根据前面的推导我们知道，$\Pr\{U\leq h(Y)\}=1/k$。也就是说我们随机生成的能满足接受条件的随机变量的有效性与包络常数$k$ 成反比。

更数学化地，令算法迭代次数为$N$ ，则获得一次被接受的随机数对$(U,X)$ 的迭代数$n$ 同样是一个随机变量，它服从于几何分布：

$$\Pr\{N=n\}=p(1-p)^{n-1},\; n\geq1$$

其中，$p=\Pr\{U\leq h(Y)\}=1/k$.
该几何分布的期望$\mathbb E(N)=1/p=k$，表明算法第一次获得被接受的随机数时的迭代次数的期望。

总而言之，$k$ 尽可能地小，算法的迭代次数就能尽可能地少，每次迭代生成的随机变量就尽可能多地符合条件并被接受。由$p(x)\leq k·q(x)$ 可推出一个最优的$k$ 应该满足：

$$k=\max_{x\in\mathcal S_X}\frac{p(x)}{q(x)}$$

实际上，在实践中往往需要同时选取一个合适的包络函数，即从一个函数族$\{q_{\theta}(x):\theta\in\Theta\}$ 中选择 candidate envelopes。所以我们一般是对所有的$\theta$ 都找到一个两个分布距离最远的坐标$\hat{x}_\theta$：

$$\hat{x}_\theta=\arg\max_{x\in\mathcal S_X}\frac{p(x)}{q_\theta(x)}=\arg\max_{x\in\mathcal S_X}[\log p(x)-\log q_\theta(x)]$$

然后从中寻得最小的$k$：

$$k=\min_{\theta\in\Theta}k_\theta=\min_{\theta\in\Theta}\frac{p(\hat x)}{q_\theta(\hat x)}$$

Example #3

对数凹密度函数

截断正态分布

数值积分

拉普拉斯近似

虽然前面提到的采样方法已经能够应对非标准的分布，但是寻找合适的建议分布（包络函数）是较为困难的。此外，在我们对后验概率分布（posterior distribution）等这种概率密度难以给出解析式的分布感兴趣时，如何对其进行合理的近似就是一个关键。

拉普拉斯近似（Laplace Approximation）就是这样一种通过使用高斯分布来近似复杂分布的技术。显然它只能用于近似连续型分布。

核心思想

对于一个非负可积函数$f$ ，记其积分为$I(f)$：

$$I(f)=\int_{\mathcal X}f(x)\;\mathrm dx$$

拉普拉斯近似方法定义函数$h(x)$ 如下：

$$h(x):=\frac1n\log f(x)$$

其中，$h(x)$ 存在一个最值点$\tilde{x}$，且有$\nabla h(\tilde{x})=0$。 将其在$x=\tilde{x}$ 处泰勒展开到第三项：

$$h(x)=h(\tilde{x})+\frac{(x-\tilde{x})^2}{2!}\nabla^2h(\tilde x)+\frac{(x-\tilde x)^3}{3!}\nabla^3h(\tilde x)+R_3(x)$$

又由于$e^y=1+y+y^2/2!+R_3'$。从而有：

$$\begin{aligned} I(f)&=\int_{\mathcal X}f(x)\;\mathrm dx=\int_{\cal X} e^{nh(x)}\mathrm dx\\ &=\int_{\cal X}e^{nh(\tilde x)}e^{0.5n(x-\tilde x)^2\nabla^2h(\tilde x)}e^{n(x-\tilde x)^3\nabla^3h(\tilde x)/6}\mathrm dx\\ &\approx e^{nh(\tilde x)}\int_{\cal X}e^{0.5n(x-\tilde x)^2\nabla^2h(\tilde x)}\left[1+\frac{n(x-\tilde x)^3}{6}\nabla^3h(\tilde x)+\frac{n^2(x-\tilde x)^3}{72}(\nabla^3h(\tilde x))^2\right]\mathrm dx \end{aligned}$$

考虑一阶近似，我们可以得到：

$$I_1(f)=e^{nh(\tilde x)}\int_{\cal X}e^{0.5n(x-\tilde x)^2\nabla^2h(\tilde x)}\mathrm dx$$

其中的被积函数正好是一个$\mu=\tilde x,\sigma^2=-1/n\nabla^2h(\tilde x)$ 的正态分布的密度函数。

由于使用到了二阶梯度，故得名为拉普拉斯近似。

多维版本

拉普拉斯近似方法 (Laplace Approximation) - 知乎 (zhihu.com)

局限性

根据推导不难得出 Laplace approximation 只能对单值函数$f$ 进行近似。并且虽然可以通过增大近似阶数来增加近似程度，但是不可避免地需要计算复杂的 Hesse 矩阵，在计算机中的代价较大。

后验分布推导

贝叶斯推断之拉普拉斯近似 - 大熊猫同学 - 博客园 (cnblogs.com)

经典蒙特卡洛积分

经典蒙特卡洛积分 （Classical Monte Carlo integration）将积分式作如下变换：

$$I(f)=\int_{\cal X}\frac{f(x)}{q(x)}q(x)\;\mathrm dx=\mathbb E_x\left[\frac{f(x)}{q(x)}\right]:=\mathbb E_x[h(x)]$$

其中$q(x)$ 是某个随机变量的 PDF.

根据辛钦大数定律，我们可以用一系列服从$q(x)$ 的独立同分布 (independent and identically distributed, i.i.d.) 的样本$\{X_i\}_{i=1}^m$ 的均值来近似此期望，即：
X

$$I(f)\approx\frac1m\sum_{i=1}^m\frac{f(X_i)}{q(X_i)}:=\bar\mu_m$$

原始版本

相关性质

数值分析：数值积分（蒙特卡洛积分） - 知乎 (zhihu.com)

黎曼和估计法

对于一元函数来说，假设取$\mathcal X=[a,b]$ ，并取$m+1$ 个点$a_i\quad(i=1,2,...,m+1)$ 使得$a=a_1\lt a_2\lt\cdots\lt a_{m+1}=b$ 则当$m\to\infty$ 时，黎曼和趋近于积分值：

$$\lim_{m\to\infty}\sum_{i=1}^mf(a_i)[a_{i+1}-a_i]=I(f)=\int_a^bf(x)\;\mathrm dx=\int_a^bh(x)q(x)\;\mathrm dx$$

其中$h(x),q(x)$ 的定义与之前一致。

现在我们将确定点$\{a_i\}_{i=1}^{m+1}$ 更换为一组从$q(x)$ 中抽取的 i.i.d. 的随机样本点，并排好序为$\{X_i\}_{i=1}^{m+1}$ ，则积分的估计值为：

$$I(f)\approx\sum_{i=1}^mf(X_i)[X_{i+1}-X_i]:=\mu^R$$

这个值也被称为黎曼和估计（Riemannian sum estimator）.

可以证明该方法在一元函数上的的收敛速度为$O(m^{-2})$ 比经典蒙特卡洛积分方法的$O(m^{-1})$ 更优。

重要性采样方法

加权估计

参考

1. Inverse transform sampling - Wikipedia
2. 用逆变换采样方法构建随机变量生成器 - MyEncyclopedia
3. 【文献阅读 02】截断(truncated)正态分布 - 知乎
4. 接受拒绝采样（Acceptance-Rejection Sampling） - 知乎
5. Laplace Approximation-拉普拉斯近似 - 知乎
6. 贝叶斯推断之拉普拉斯近似 - 大熊猫同学 - 博客园
7. 拉普拉斯近似方法 (Laplace Approximation) - 知乎
8. 数值分析：数值积分（蒙特卡洛积分） - 知乎
9. 蒙特卡洛积分与重要性采样详解 - 烈日行者 - 博客园