0x01 蒙特卡罗与马尔科夫链

参考的资料:

MCMC(一)蒙特卡罗方法

MCMC(二)马尔科夫链

MCMC(三)MCMC采样和M-H采样

MCMC(四)Gibbs采样

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

简单易学的机器学习算法—马尔可夫链蒙特卡罗方法MCMC

为什么要使用MCMC方法?

面对的场景

对于一个随机变量$X$, 我们想求关于这个随机变量的一些性质, 相关值, 统计学特性等等, 如:

• 关于$x$的函数$f(x)$的积分

• 相关的面积

• 后验概率, 期望等

这里的随机变量$X$可以是连续的, 即有无限个状态, 也可以是离散的, 即有限个状态.

对于这些问题, 一个关键的问题就是必须知道随机变量$X$的概率分布$P(X)$. 然而当这个$P(X)$十分复杂, 或者难以用显式的方式表示时, 又要怎么办呢?

蒙特卡罗方法的引入

这时, 我们可以引入蒙特卡罗思想, 通过随机采样的方法, 即接受-拒绝采样的方法, 来对随机变量$X$的分布概率进行估计.

首先, 我们看一个随机采样的例子:

假设, 我们想估计下图中圆圈的面积:

我们在正方形内随机地采样(可以对x, y坐标分别进行均匀分布的采样)得到N个点, 计算落在园内点的比例, 再乘以正方型的面积, 我们就能估计出园的面积.

对于更复杂的图形, 我们都可以通过这种随机采样的方法进行估计.

那么, 对于随机变量的分布概率, 要怎么通过采样进行呢?

对于随机变量$X$本身的分布概率$P(X)$, 本身太复杂导致无法采样. 那么我们设定一个可采样的分布$q(X)$, 例如高斯分布:

image

使用这个方便采样的分布$q(x)$, 以及使用一个常量$k$, 使得真正的概率分布$p(x)$总在$kq(x)$的下方.

• 首先根据$q(x)$进行采样, 得到一个$x_0$值
• 从均匀分布$[0, kq(x)]$中采样得到一个值$u$
• 若$u>p(x_0)$, 即在灰色区域, 拒绝这次抽样, 否则接受这次抽样
• 重复以上过程得到$n$个接受的样本$x_0,x_1,\cdots,x_{n-1}$
• 通过这$n$个点, 我们就得到了随机变量$X$整体的分布情况, 可以继续进行别的运算, 如积分, 求期望等

以上采样是对连续情况的采样, 但这种采样只能在我们能够显式地得到概率分布的是否才能使用, 但实际上, 我们会遇到很多情况, 无法得到显式的概率分布, 比如:

• 概率分布无法被显示的计算

• 高维情况下, $p(x_1,x_2,...,x_n)$, 去寻找合适的$q(x)$和$k$是非常困难的

  • 即使可行, 每个维度都要进行采样, 当把这些维度放在一起, 需要采样的数量是指数级的, 这不是现有计算能力能解决的.

因此, 蒙特卡洛方法需要解决在复杂概率分布(现实中多数是高维情况)时对应的采样困难问题, 这时, 就需要马尔科夫链来解决这个问题了.

马尔科夫链的引入

马尔科夫链是如何解决这个问题的呢? 马尔科夫链每个时刻都会有一个分布, 对于离散的就是每个状态的概率, 连续的就是一个概率分布函数.

我们将这个分布和我们的目标采样分布$p(x)$结合起来, 利用马尔科夫链的平稳分布进行采样来代替真实的分布.

我们采样的样本, 就是马尔科夫链经过转移达到平稳分布, 之后每个时刻采样一个样本, 从这个时刻的分布从抽样一个状态, 即一个样本.

关于马尔科夫链为何能解决高维空间样本组合爆炸的问题, 移步马尔可夫链蒙特卡洛(MCMC)采样中的第三部分.

现在我们要做的构造一个具有平稳分布$p(x)$的马尔科夫链了.

马尔科夫链的构造

1. 马尔科夫链

   马尔科夫链有一个很重要的性质, 状态转移的概率只依赖于当前时刻, 具体来说对于序列状态$...X_{t-2}, X_{t-1}, X_{t}, X_{t+1},...$, 我们在时刻$X_{t+1}$的状态转移概率只依赖于$X_{t}$, 即:

   $P(X_{t+1}=x|X_t, X_{t-1}, \cdots) =P(X_{t+1}=x|X_t)$
2. 状态转移概率

   状态转移概率是指随机变量在一个时刻由状态$s_i$转移到下一个时刻状态$s_j$的概率, 即:

   $P\left ( i\rightarrow j \right ):=P_{i,j}=P\left ( X_{t+1}=s_j\mid X_t=s_i \right )$

   记$\pi _k^{\left ( t \right )}$为随机变量$X$在$t$时刻取值为$s_k$的概率, 则随机变量$X$在$t+1$时刻为$s_i$的概率为的概率为:

   $$\begin{aligned} \pi i^{\left ( t+1 \right )} &=P\left ( X{t+1}=si \right ) \ &= \sum{k}P\left ( X{t+1}=s_i\mid X{t}=sk \right )\cdot P\left ( X{t}=sk \right )\ &= \sum{k}P_{k,i}\cdot \pi _k^{\left ( t \right )} \end{aligned}

   假设状态的数目有$$n$$种, 这时则有: $$\left ( \pi _1^{\left ( t+1 \right )},\cdots ,\pi _n^{\left ( t+1 \right )} \right )=\left ( \pi _1^{\left ( t \right )},\cdots ,\pi _n^{\left ( t \right )} \right )\begin{bmatrix} P_{1,1} & P_{1,2} & \cdots & P_{1,n}\\ P_{2,1} & P_{2,2} & \cdots & P_{2,n}\\ \vdots & \vdots & & \vdots \\ P_{n,1} & P_{n,2} & \cdots & P_{n,n} \end{bmatrix}

   我们把$\left ( \pi _1^{\left ( t \right )},\cdots ,\pi _n^{\left ( t \right )} \right )$称为$t$时刻的状态向量(随机变量共有$n$个离散的状态), 也是这个时刻下随机变量的分布情况; 而右侧的矩阵被称为转移概率矩阵$P$, 决定了一个马尔科夫链的表现.

   ==注意==

   在高维离散的情况下, 这里的状态向量将是一个多维联合分布, 概率转移矩阵也将是一个多维的条件分布.

3. 连续情况的状态转移

   以上是对于离散的情况, 即状态为有限多个. 那么对于状态为无限多个, 也就是连续的情况又该如何定义呢?

   首先, 状态向量$\pi$是一个连续函数, 即连续随机变量$X$的概率函数. 若进行抽样就是根据此时的概率函数进行抽样.

   而转移概率矩阵此时也将不是一个矩阵, 而是一个转移核, 也可以理解为一个条件概率分布函数$P(y|x)$, 其中$x, y$分别代表随机过程在两个不同时态(此时的时态仍然是离散的)上的取值. 转移核的构建有很多方法, 比如使用核函数表示两个值之间的距离. 在PyMC包中的贝叶斯推断方法, 会将训练样本的真实值考虑进来, 来构建这个转移核, 等等.

   将以上两者相乘(即将概率函数和转移核相乘), 就得到了随机过程在下一个时间点上的连续分布函数. 我们在这个函数上抽样, 就得到下一个时间点上的某一状态了.

4. 马尔科夫链的平稳分布

   我们还是用状态是离散的情况来继续考虑.

   如果马尔科夫链满足以下两条性质:

   • 非周期性: 周期性指的是马尔科夫链的状态转换是循环的, 即经过有限次的状态转移, 又回到了自身. 这里要求马尔科夫链没有周期性

   • 不可约性: 即两两状态之间都能够相互转移

     如果一个马尔科夫链满足以上两条性质, 那么对于任何的初始状态向量$\pi ^{\left ( 0 \right )}$, 经过若干次的转移, 随机变量$X$的分布都能够收敛到唯一的平稳分布$\pi ^{\ast }$, 即:

     $\underset{t\rightarrow \infty }{lim}\pi ^{\left ( 0 \right )}\mathbf{P}^t=\pi ^{\ast }$

     且这个平稳分布满足:

     $\pi ^{\ast }\mathbf{P}=\pi ^{\ast }$

     也就是说在达到平稳分布以后的时刻, 再进行状态转移, 其分布也不会改变, 仍未$\pi ^{\ast }$

     因此, 这个平稳分布$\pi ^{\ast }$就是我们想要想要进行采样的分布$p(x)$, 这个平稳分布就可以近似代替随机变量的真实分布了, 我们在马尔科夫链收敛到平稳分布之后, 每个时刻抽样一个样本即可, 就相当与在真实分布上进行抽样.
5. 基于马尔科夫链的采样过程

   采样过程具体为:

   1. 初始化: 输入马尔科夫链的状态转移矩阵$P$, 设定状态转移次数阈值$n_1$(达到平稳过程需要的状态转移次数), 需要的样本个数$n_2$.
   2. 初始化: 从先验分布$\pi^0$采样得到初始的状态值$x_0$
   3. 进行状态转移, 从条件概率分布$P(x|x_t)$中采样得到下一时刻的状态$x_{t+1}$, 这个过程重复$n_1+n_2$次
   4. 样本集合$(x_{n_1}, x_{n_1+1},..., x_{n_1+n_2-1})$就是我们从平稳分布中采集到的逼近于真实分布的样本.

      上面可以看到, 我们每个时刻只会考虑一个状态, 而不会对所有状态进行考虑. 因为某一个马尔科夫链的性质就是完全由它的状态转移矩阵决定的.

      ==注意==

      由此来讲, 如果随机变量$X$共有$n$种离散的状态, 我们不必考虑每个时刻所有状态的转移(这样我们需要考虑$n^2$种情况). 在实现代码的时候, 我们只需要考虑实现当前状态$x_i$到下一时刻的状态$x_j$的条件概率$P(x_j|x_i)$的接口, 并不需要一个完整的概率矩阵.

      ==实际实现==

   5. 实际过程中, 我们只需要为当前状态和目标状态定义一个统一的转移概率计算方式就可以了

   6. 此外, 在高维(离散)空间中, 我们还可以提前限定好一个状态下一时刻所能到达的可能状态数, 即每次转移只允许某一个维度上变化.

      通过以上两种方式, 就可以完全避免多维产生的维度爆炸问题.

      ==本质==

      通过MCMC来做抽样本质上就是将样本纵向在空间上的分布转化成为了样本横向在时间上的等价分布, 所以MCMC具备避免空间爆炸的能力是必然的.

   但是, 如何构造马尔科夫链中的转移概率矩阵$P$仍未解决. 那么具体如何采样仍不得知.

   接下来讨论MCMC如何采样的问题.

6. 细致平稳条件

   要解决MCMC采样的问题, 我们从细致平稳条件来切入.

   马尔科夫链还有一个重要的细致平稳条件, 即马尔可夫链在某一时刻的分布$\pi$为平稳态的一个充分不必要条件为:

   $\pi_{i}\times p_{ij} = \pi_j\times p_{ji}$

   这是在离散情况下的公式, $\pi_{i}$表示随机变量处于状态$x_i$(即上文的$s_i$)的概率, $p_{ij}$为状态$x_i$转移到$x_j$的概率, 即$P(x_{j}|x_{i})$.

   证明:

   $\sum_{i}\pi_{i}p_{ij} = \sum_{i}\pi_{j}p_{ji} = \pi_j$

   (第一个等式是将一直平稳条件带入得到的.)

   将上式由单个状态$x_i$的概率$\pi_{i}$扩充成所有状态, 将上式写成矩阵形式, 就会得到$\pi \mathbf{P}=\pi$的形式, 因此此时的状态分布是平稳的.

   因此, 我们要找到一个马尔科夫链(实际上就是寻找状态转移矩阵$P$), 使得平稳分布$\pi$满足上述的性质, 这是我们继续进行的目标.
7. 基础的MCMC采样

   我们的目标现在是满足细致平稳条件, 即寻找一个状态转移矩阵$P$, 使得:

   $\pi(i)P(i,j) = \pi(j)P(j,i)$

   但是对于任意一个容易采样的分布$Q$, 很难满足:

   $\pi(i)Q(i,j) \neq \pi(j)Q(j,i)$

   注意这里的$P(i,j)$和$Q(i,j)$等都是条件概率.

   我们加以改造, 对上式引入$\alpha(i,j)$, 使得上式可以恒等:

   $\pi(i)Q(i,j)\alpha(i,j) = \pi(j)Q(j,i)\alpha(j,i)$ ​ 这要求$\alpha(i,j)$满足:

   $\alpha(i,j) = \pi(j)Q(j,i)$

   $\alpha(j,i) = \pi(i)Q(i,j)$

   这样, 我们就得到了马尔科夫链的转移概率矩阵, 满足:

   $P(i,j) = Q(i,j)\alpha(i,j)$

   ==理解==

   马尔科夫链的状态转移矩阵$P$通过一个容易采样的条件概率矩阵$Q$和一个接受率$\alpha(i,j)$组合获得. 接受率$\alpha(i,j)$的取值在$[0, 1]$之间, 表示着是否接受这个采样的概率值.

   这个形式与前面蒙特卡洛方法中的接受-拒绝采样很相似.

   • 那里是以一个常用分布通过一定的接受-拒绝概率得到一个非常见分布

   • 这里是以一个常见的马尔科夫链状态转移矩阵$Q$通过一定的接受-拒绝概率得到目标转移矩阵$P$, 两者的解决问题思路是类似的

   现在MCMC采样的过程为:

   • 选定一个常用的马尔科夫链状态转移矩阵$Q$, 设定状态转移次数阈值$n_1$, 需要的样本个数$n_2$.
   • 从先验分布$\pi^0$采样得到初始的状态值$x_0$
   • for t=0 to $n_1$ + $n_2$ - 1:

     • 从条件概率分布$Q(x|x_t)$中采样得到样本$x_*$
     • 从均匀分布中采样得到$u \sim uniform[0,1]$
     • 若$u < \alpha(x_t,x_{*}) = \pi(x_{*})Q(x_{*},x_t)$, 这里的$\alpha(x_t,x_{*})$就是从$$x

       t$转移到$x$的接受率, 则接受转移:$x{t+1}= x{}$$
     • 否则不接受转移

   • 样本集合$(x_{n_1}, x_{n_1+1},..., x_{n_1+n_2-1})$就是我们从平稳分布中采集到的逼近于真实分布的样本.
8. Metropolis-Hasting采样(M-H采样)

   上面的MCMC采样方法在实际中是难以使用的, 原因为:

   接受率$\alpha(x_t,x_{*}) = \pi(x_{*})Q(x_{*},x_t)$往往是非常小的, 例如0.1, 这会导致大部分的采样都会被拒绝转移, 这样收敛过程以及采样的效率都会非常的低, 在高维的情况下情况更加严重.

   为了解决接受率过小的问题, 我们可以尝试着将等式两边的接受率同时方法相同的比例, 即:

   $\pi(i)Q(i,j)[n*\alpha(i,j)] = \pi(j)Q(j,i)[n*\alpha(j,i)]$

   这样等式依然是成立的, 但接受率被放大之后, 被拒绝的情况就少了.

   放大到何种程度呢, 我们将$\alpha(i,j)$和$\alpha(j,i)$两者中的最大值放大到1, 即一定接受的程度, 又由于:

   $\alpha(i,j) = \pi(j)Q(j,i)$

   $\alpha(j,i) = \pi(i)Q(i,j)$

   当$\alpha(i,j)$为最大值或$\alpha(j,i)$为最大值时就会有着不同的放大倍数, 而综合起来, $\alpha(i,j)$最终会被放大为:

   $\alpha(i,j) = min\{ \frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)},1\}$

   证明当$P(i,j) = Q(i,j)\alpha(i,j)$且$\alpha(i,j) = min\{ \frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)},1\}$时, 满足细致平稳条件:

   $\begin{aligned}\pi(i)P(i,j) &= \pi(i)Q(i,j)\alpha(i,j) \\ &= \pi(i)Q(i,j)\min\{ \frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)},1\} \\ &= \min\{\pi(j)Q(j,i), \pi(i)Q(i,j)\} \\ &= \pi(j)Q(j,i)\min\{\frac{\pi(i)Q(i,j)}{\pi(j)Q(j,i)}, 1\} \\ &= \pi(j)Q(j,i)\alpha(j,i) \\ &= \pi(j)P(j,i) \end{aligned}$

   证明完毕.

   经过改造后的M-H采样过程如下:

   • 选定一个常用的马尔科夫链状态转移矩阵$Q$, 设定状态转移次数阈值$n_1$, 需要的样本个数$n_2$.
   • 从先验分布$\pi^0$采样得到初始的状态值$x_0$
   • for t=0 to $n_1$ + $n_2$ - 1:

     • 从条件概率分布$Q(x|x_t)$中采样得到样本$x_*$
     • 从均匀分布中采样得到$u \sim uniform[0,1]$
     • 若$u < \alpha(x_t,x_{*}) = min\{ \frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)},1\}$, 这里的$\alpha(x_t,x_{*})$就是从$$x

       t$转移到$x$的接受率, 则接受转移:$x{t+1}= x{}$$
     • 否则不接受转移

   • 样本集合$(x_{n_1}, x_{n_1+1},..., x_{n_1+n_2-1})$就是我们从平稳分布中采集到的逼近于真实分布的样本.

   M-H采样存在的问题:

   1. 高维情况下, 计算接受率计算式$\alpha(x_t,x_{*}) = min\{ \frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)},1\}$需要的时间会很多, 导致算法的效率很低. 而且由于$\alpha\lt{1}$, 很多转移提议被否决, 这就浪费了计算力, 能否做到不拒绝转移呢?
   2. 高维(特征很多)情况下, 甚至很难求出目标的各个特征的联合分布, 但可以方便求出各个特征的之间的条件概率分布, 能不能使用条件概率分布进行采样呢

9. Gibbs采样

   为了解决上面的问题, 我们使用Gibbs采样

   Gibbs的思路是:

   • 限定下一个时刻所能达到的可能的状态数, 即从当前时刻的状态, 只能转移到有限的状态上, 而不是全部的状态

   • 具体来说, 对于高维空间, 每次转移只允许在某一个维度上进行变化

   二维空间下的细致平稳条件:

   $\pi(x_1,x_2)$是一个二维联合分布, 观测两个点$A(x_1^{(1)},x_2^{(1)})$和$B(x_1^{(1)},x_2^{(2)})$, 两个点在第一个维度上是相等的, 则以下两式成立:

   • 从点$A$进行状态转移, 且第一维度保持不变, 第二维度上变为$x_2^{(2)}$的概率:

     $\pi(x_1^{(1)},x_2^{(1)}) \pi(x_2^{(2)} | x_1^{(1)}) = \pi(x_1^{(1)})\pi(x_2^{(1)}|x_1^{(1)}) \pi(x_2^{(2)} | x_1^{(1)})$
   • 从点$B$进行状态转移, 且第一维度保持不变, 第二维度上变为$x_2^{(1)}$的概率

     $\pi(x_1^{(1)},x_2^{(2)}) \pi(x_2^{(1)} | x_1^{(1)}) = \pi(x_1^{(1)}) \pi(x_2^{(2)} | x_1^{(1)})\pi(x_2^{(1)}|x_1^{(1)})$

     两式的右端是相等的, 因此两式左端相等, 即:

     $\pi(x_1^{(1)},x_2^{(1)}) \pi(x_2^{(2)} | x_1^{(1)}) = \pi(x_1^{(1)},x_2^{(2)}) \pi(x_2^{(1)} | x_1^{(1)})$

     也就是:

     $\pi(A) \pi(x_2^{(2)} | x_1^{(1)}) = \pi(B) \pi(x_2^{(1)} | x_1^{(1)})$

     这就是Gibbs采样中的二维情况下的细致平稳条件:

   • 在$x_1 = x_1^{(1)}$这条直线上, 使用条件概率分布$\pi(x_2| x_1^{(1)})$作为状态转移概率, 则直线上的任意两点满足细致平稳条件.
   • 因此, 在任意一个维度的具体值确定的直线, 直线上的任意两点之间的条件概率分布就能够满足细致平稳条件的状态转移概率

     因此, 我们构造如下的状态转移概率矩阵$P$:

     $P(A \to B) = \pi(x_2^{(B)}|x_1^{(1)})\;\; if\; x_1^{(A)} = x_1^{(B)} =x_1^{(1)}$

     $P(A \to C) = \pi(x_1^{(C)}|x_2^{(1)})\;\; if\; x_2^{(A)} = x_2^{(C)} =x_2^{(1)}$

     $P(A \to D) = 0\;\; else$

     以上的转移概率矩阵$P$对于平面上任意两点$E$, $F$, 都能满足细致平稳条件:

     $\pi(E)P(E \to F) = \pi(F)P(F \to E)$

   二维Gibbs采样方法:

   • 输入平稳分布$\pi(x_1,x_2)$, 设定状态转移次数阈值$n_1$, 需要的样本个数$n_2$
   • 随机初始化初始状态值$x_1^{(0)}$, $x_2^{(0)}$
   • for $t=0$ to $n_1 +n_2-1$:

     • 从条件概率分布$P(x_2|x_1^{(t)})$中采样得到样本$x_2^{t+1}$
     • 从条件概率分布$P(x_1|x_2^{(t+1)})$中采样得到样本$x_1^{t+1}$
   • $\{(x_1^{(n_1)}, x_2^{(n_1)}), (x_1^{(n_1+1)}, x_2^{(n_1+1)}), ..., (x_1^{(n_1+n_2-1)}, x_2^{(n_1+n_2-1)})\}$样本集就是我们需要的平稳分布对应的样本集.

   多维Gibbs采样方法:

   将二维的情况推广到多维, 我们每次转移的时候只需要对一个维度进行转移, 其他维度不变, 通过这种思路来构建状态转移概率矩阵.

   比如一个$n$维的概率分布$\pi(x_1,x_2,...x_n)$, 通过对$n$个坐标轴上轮换采样, 来得到新的样本. 对于轮换到的任意一个坐标轴$x_i$上的转移, 状态转移概率为$P(x_i|x_1,x_2,...,x_{i-1},x_{i+1},...,x_n)$, 即固定$n-1$个坐标轴, 在某一个坐标轴上移动.

   具体过程如下:

   • 输入平稳分布$\pi(x_1,x_2,...,x_n)$, 或对应的所有特征的条件概率分布. 设定状态转移次数阈值$n_1$, 需要的样本个数$n_2$
   • 随机初始化初始状态值$(x_1^{(0)},x_2^{(0)},...,x_n^{(0)})$
   • for $t=0$ to $n_1 +n_2-1$:

     • 从条件概率分布$P(x_1|x_2^{(t)}, x_3^{(t)},...,x_n^{(t)})$中采样得到$x_1^{t+1}$
     • 从条件概率分布$P(x_2|x_1^{(t+1)}, x_3^{(t)}, x_4^{(t)},...,x_n^{(t)})$中采样得到$x_2^{t+1}$
     • 依次从每个维度中采样

   • $\{(x_1^{(n_1)}, x_2^{(n_1)},..., x_n^{(n_1)}), ..., (x_1^{(n_1+n_2-1)}, x_2^{(n_1+n_2-1)},...,x_n^{(n_1+n_2-1)})\}$样本集就是我们需要的平稳分布对应的样本集.
10. MCMC抽样方法的细节

    不管是用哪种抽样方法(M-H抽样, Gibbs抽样, 以及其他等等抽样), 我们在抽样过程中都不是直接使用$\pi(x)$本身, 因为随机变量$X$可能有很多状态, 甚至无限多的状态(连续情况下), 而且在高维情况下状态更加爆炸.

    在采样过程中, 每个时间状态只采样出一个样本, 然后下一个时刻的状态只于这一个具体的状态有关. 因此, 这里的每轮的采用都将使用条件概率分布$\pi(y_{i}|x_{-i})$代替原来每个时刻的$\pi(x)$, 这样就解决了分布本身难以计算的问题.

    对于M-H采样, 每一步提议新的样本是从提议矩阵$Q$代表的条件概率分布$Q(x|x^t)$采样得到的, 即对于每个时刻, 原本应当从目标分布$\pi$中采样的下一时刻的状态, 现在改为对提议矩阵代表的条件概率分布$Q(x|x^t)$进行采样.

    因此在计算接受率$\alpha$中使用到的$\pi$中的值可以用对应的条件概率分布$Q(x|x^t)$的值来代替.