Monte-Carlo Sampling

基础采样方法

在机器学习中，我们常常遇到这样的问题：设 $X$ 为一随机变量，其 PDF 为 $p(x)$；又设 $f$ 为关于 $X$ 的函数，考虑如下期望： $$\mathbb{E}[f(X)]=\int f(x)p(x)\mathrm{d}x$$ 当 $p(x)$ 或 $f(x)$ 比较复杂时，上述积分往往是无法计算的。蒙特卡洛采样的思想是用随机样本近似期望： $$\hat{f}(X^{(1)},\dots ,X^{(n)})=\frac{1}{N}\sum _{n=1}^{N}f(X^{(n)}),\text{where }{X}^{(n)}\sim p(x)$$ 容易知道，该统计量是真实结果的无偏估计，即： $$\mathbb{E}[\hat{f}]=\mathbb{E}[\frac{1}{N}\sum _{n=1}^{N}f(X^{(n)})]=\mathbb{E}[f(X)]$$ 而估计方差为： $$\text{var}(\hat{f})=\text{var}(\frac{1}{N}\sum _{n=1}^{N}f(X^{(n)}))=\frac{1}{N}\text{var}(f(X))$$ 可见样本量 $N$ 越大，方差越小，估计就越准确。直接使用蒙特卡洛方法要求我们能够从 $p(x)$ 中采样，下面介绍三种基础的采样方法——变量替换、拒绝采样和重要性采样。

变量替换 (change of variables)

定理：设 $p_X(x)$ 为随机变量 $X$ 的 PDF，单调函数 $f$ 将 $X$ 映射到 $Y$，即 $Y=f(X)$，则： $$p_Y(y)=p_X(x)\left|\frac{\mathrm{d}x}{\mathrm{d}y}\right|$$ 对于多维变量有类似结论： $$p_Y(y_1,\dots ,y_M)=p_X(x_1,\dots ,x_M)\left|\frac{\partial (x_1,\dots ,x_M)}{\partial (y_1,\dots ,y_M)}\right|$$ 其中 ${\displaystyle \frac{\partial (x_1,\dots ,x_M)}{\partial (y_1,\dots ,y_M)}}$​ 表示 Jacobian 行列式。

例子：Box-Muller 方法是从二维正态分布中采样的方法。首先生成 $x_1,x_2\sim U(0,1)$，然后设： $$y_1=\sqrt{-2\ln x_1}\cos (2\pi x_2),y_2=\sqrt{-2\ln x_1}\sin (2\pi x_2)$$ 那么可以证明 $(y_1,y_2)$ 服从各分量独立的标准二维正态分布： $$\begin{array}{rl}{p}_Y(y_1,y_2)& =p_X(x_1,x_2)\left|\frac{\partial (x_1,x_2)}{\partial (y_1,y_2)}\right|={\left|\frac{\partial (y_1,y_2)}{\partial (x_1,x_2)}\right|}^{-1}\\ & ={|det\left(\begin{array}{cc}-\frac{\cos (2\pi x_2)}{x_1\sqrt{-2\ln x_1}}& -2\pi \sqrt{-2\ln x_1}\sin (2\pi x_2)\\ -\frac{\sin (2\pi x_1)}{x_1\sqrt{-2\ln x_1}}& 2\pi \sqrt{-2\ln x_1}\cos (2\pi x_2)\end{array}\right)|}^{-1}\\ & ={|-\frac{2\pi }{x_1}{\cos }^2(2\pi x_2)-\frac{2\pi }{x_1}{\sin }^2(2\pi x_2)|}^{-1}\\ & =\frac{x_1}{2\pi }=\frac{1}{2\pi }\exp (-\frac{y_1^2+y_2^2}{2})\\ & =[\frac{1}{\sqrt{2\pi }}\exp (-y_1^2/2)][\frac{1}{\sqrt{2\pi }}\exp (-y_2^2/2)]\end{array}$$ 直观上看，可以认为 $\sqrt{-2\ln x_1}$ 在采样向量的模长，而 $2\pi x_2$ 在采样向量的幅角。

注：上述过程生成的是标准正态随机变量 $y\sim \mathcal{N}(0,I)$，若要获取 $\mathcal{N}(\mu ,\mathrm{\Sigma })$，只需对 $\mathrm{\Sigma }$ 做 Cholesky 分解 $\mathrm{\Sigma }=L{L}^T$，然后做变换 $y^{\prime }=\mu +Ly$ 即可。

例子：对于可写出 CDF 及其反函数的简单分布 $p_Y(y)$，有一种从均匀分布变换到该分布的方法，一些书籍称之为概率积分变换。设 $x\sim U(0,1)$，我们希望找到一个变换函数 $f$ 使得 $y=f(x)$ 服从 $p_Y(y)$. 假设 $f$ 单调递增并且取值在 $(0,1)$ 之间，则： $$F_Y(y)=P(Y\le y)=P(f(X)\le y)=P(X\le f^{-1}(y))=f^{-1}(y)⟹f(y)=F_Y^{-1}(y)$$ 这说明我们要找的 $f$ 就是随机变量 $Y$ 的 CDF 的逆 $F_Y^{-1}$. 例如，对于指数分布： $$p_Y(y)=\lambda \exp (-\lambda y)$$ 其 CDF 为： $$F_Y(y)={\int }_0^y\lambda \exp (-\lambda t)\mathrm{d}t=1-\exp (-\lambda y)$$ 其逆为： $$F_Y^{-1}(y)=-\frac{1}{\lambda }\ln (1-y)$$ 因此对于 $x\sim U(0,1)$，只需要做变换 $y=-\frac{1}{\lambda }\ln (1-x)$，就有 $y\sim p_Y(y)$.

该方法要求目标分布的 PDF 或 CDF 具有解析形式，因此只对比较简单的分布有效。

拒绝采样 (rejection sampling)

假设要采样的分布 $p(x)$ 在任意一点 $x$ 处的值（在不考虑归一化常数下）是可计算的，即设： $$p(x)=\tilde{p}(x)/Z_p$$ 其中 $\tilde{p}(x)$ 是已知的，$Z_p$ 为未知的归一化常数，那么我们可以使用拒绝采样方法进行采样。

首先引入一个提议分布 (proposal distribution) $$，满足该提议分布是容易采样的，并且存在常数 $$ 使得 $$，如图所示：

拒绝采样方法先从 $$ 中采样一个 $$，然后以 $$ 的概率接受该采样结果。直观来看，就是拒绝掉上图中的灰色区域。如此，$$ 最终被采出来的概率就是： $$ 这验证了拒绝采样的正确性。容易知道，拒绝采样的总接受率为： $$ 因此 $$ 越小，总接受率越大，算法效率越高。然而，$$ 小也意味着 $$ 本身就要与 $$ 比较相似，对于复杂的 $$ 而言寻找到一个合适的 $$ 非常困难的。

重要性采样 (importance sampling)

重要性采样不从 $$ 中采样出样本，而是直接近似 $$. 假设 $$ 在任意一点 $$ 处的值都是可计算的，仍然引入一个提议分布 $$，满足该提议分布是容易采样的。重要性采样的思路是将对 $$ 的采样转化为对 $$ 的采样。具体而言，有：

$$ 于是只需要从 $$ 中采样即可近似上式： $$ 其中系数 $$​​ 称作重要性权重 (importance weights)。

进一步地，假设 $$ 且我们只能计算 $$ 而归一化系数未知，类似地假设 $$ 并且只能计算 $$ 而归一化系数未知，那么： $$ 其中系数 $$ 可以通过同一套样本估计： $$ 代入得： $$ 与拒绝采样一样，重要性采样的效果与提议分布 $$ 同 $$ 的接近程度紧密相关。当 $$ 比较复杂时，选择合适的 $$ 是非常困难的。

马尔可夫链蒙特卡洛 (MCMC)

在高维场景下，通常难以为拒绝采样和重要性采样找到合适的提议分布，导致采样效率低下。此时我们可以考虑马尔可夫链蒙特卡洛 (MCMC) 方法。顾名思义，MCMC 指基于马尔可夫链的随机采样方法，主要包括 Metropolis-Hastings 采样和 Gibbs 采样。

马尔可夫链

首先回顾一下有关马尔可夫链的基本知识。

马尔可夫链：若离散随机变量序列 $$ 满足如下条件（马尔可夫性质）： $$ 其中 $$ 为随机变量的状态空间，称 $$ 为马尔可夫链。

一步转移概率：$$ 时刻从状态 $$ 转移到状态 $$ 的概率： $$ 齐次马尔可夫链：若一步转移概率与时刻 $$ 无关，则称该马尔可夫链为齐次马尔可夫链： $$ 一步转移矩阵：齐次马尔可夫链的一步转移概率构成矩阵： $$ 转移矩阵的每一行之和为 1.

平稳分布：若状态空间上的分布 $$ 满足方程： $$ 则称 $$ 为马尔可夫链的平稳分布（不变分布）。值得注意的是，平稳分布不一定存在，存在也不一定唯一。下面给出平稳分布存在的一个充分条件。

细致平衡 (detailed balance) 条件：若分布 $$ 满足如下条件，则 $$ 是该马尔可夫链的平稳分布： $$ 证明很简单： $$ 称满足细致平衡条件的马尔可夫链为可逆的 (reversible) 马尔可夫链。

极限分布：若 $$ 对任意状态 $$ 都存在，则称 $$ 为马尔可夫链的极限分布。极限分布如果存在，则是唯一的，且与初始分布无关。

平稳分布与极限分布的关系：一般情形下二者没有必然的联系；但在一些约束条件（两两联通、非周期、正常返）下，齐次马尔可夫链存在极限分布，且该极限分布就是平稳分布。

MCMC 基本思想

MCMC 方法的基本思想是：构建一个马尔可夫链，使其极限分布就是要采样的分布。这样，只需要从任一初始分布出发，采样 $$，然后采样 $$……当迭代次数足够多时，$$ 开始服从极限分布，那么 $$ 就是从极限分布中采样出来的样本了。不过需要注意的是，相邻两次采样得到的 $$ 显然并不独立，如果有独立性要求，可以在两次采样之间间隔一定的步数。

现在，问题的关键是如何构建这样的马尔可夫链，也就是设计转移矩阵 $$ 使其极限分布为给定的分布。在极限分布就是平稳分布的假设下，我们只需要让 $$ 满足细致平衡条件即可，但这依旧不是一件容易的事。受到拒绝采样的启发，Metropolis-Hastings 采样通过引入一个提议转移矩阵 $$ 和接受率 $$ 解决该问题。

Metropolis-Hastings 采样

如上一节所说，现在我们的目标是找到一个满足细致平稳条件的转移矩阵 $$. 为此，首先随机找一个转移矩阵 $$，那么对不相同的两个状态 $$，它大概率并不满足细致平稳条件： $$ 注：$$ 时等式恒成立，因此这里只考虑 $$.

为了让细致平稳条件成立，引入接受率 $$，容易验证： $$ 因此只需要令： $$ 就找到了要求的 $$. 等价地，这样的一步转移相当于先按照 $$ 从 $$ 转移到 $$，如果 $$，则以 $$ 的概率接受转移，否则不转移。

不过上述方法有一个问题：接受率 $$ 可能会很小，影响算法效率。注意到等比放大 $$ 并不影响细致平稳条件的成立，因此做出如下改进： $$ 简单来说，就是等比放大 $$​ 直到其中一个变成 1 为止，这样任意两个状态之间至少有一个方向是始终能转移成功的。

综上，M-H 采样的算法伪代码如下：

def MHSampling(p, Q):
    // p is the target distribution
    // Q is the proposal transfer matrix
    sample s from any distribution
    while True:
        sample t from Q[s,t]
        sample u from U(0,1)
        alpha = min((p[t]*Q[t,s]) / (p[s]*Q[s,t]), 1)
        if u <= alpha:
            s = t
        if sufficiently large time:
            yield s
PSEUDOCODE

Gibbs 采样

Gibbs 采样可以视为 M-H 采样的一种特殊形式。在 M-H 采样中，如果每个状态都是高维向量，那么状态数量十分巨大，选择合适的状态转移矩阵或状态转移概率 $$ 比较困难。为了解决这个问题，Gibbs 采样限制每次转移时只改动一个维度。假设当前状态为： $$ 转移时只改动第 $$​ 维，转移到新的状态： $$ 其中 $$ 表示 $$，即除了第 $$ 维以外的其他维度。Gibbs 采样用目标分布的条件概率来定义这个转移概率： $$ 如此接受率为： $$ 我们发现 Gibbs 采样的接受率为 1，因此在实际应用中有着较高的效率，是最广为使用的 MCMC 方法之一。

Gibbs 采样每次转移只改动一个维度，所以使用时需要指定维度顺序，具体顺序与具体问题相关，一般依次选择即可。

综上，Gibbs 采样的算法流程如下：

参考资料

1. Bishop, Christopher. Pattern recognition and machine learning. ↩︎
2. 刘建平Pinard. MCMC(一)蒙特卡罗方法. http://www.cnblogs.com/pinard/p/6625739.html ↩︎
3. 刘建平Pinard. MCMC(二)马尔科夫链. http://www.cnblogs.com/pinard/p/6632399.html ↩︎
4. 刘建平Pinard. MCMC(三)MCMC采样和M-H采样. https://www.cnblogs.com/pinard/p/6638955.html ↩︎
5. 刘建平Pinard. MCMC(四)Gibbs采样. http://www.cnblogs.com/pinard/p/6645766.html ↩︎
6. 机器学习-白板推导系列(十三)-MCMC(Markov Chain Monte Carlo). https://www.bilibili.com/video/BV1dW411z7qs ↩︎