ML——MCMC采样

本文介绍MCMC采样法和他的两个常用方法:MH(Metropolis-Hastings)采样法和Gibbs采样法
马尔科夫蒙特卡洛(Markov Chain Monte Carlo, MCMC)采样法

  • 对于高维空间中的随机向量,拒绝采样和重要性采样经常难以找到合适的参考分布,容易导致采样效率低下(样本的接受概率太小或者重要性权重太低),此时可以考虑马尔科夫蒙特卡洛采样法(MCMC)
  • MCMC中常见的有两种,MH(Metropolis-Hastings)采样法和Gibbs采样法

MCMC概述

  • 马尔科夫蒙特卡洛(Markov Chain Monte Carlo, MCMC)采样法可分为两个部分(两个MC)描述

蒙特卡洛法

  • 蒙特卡洛法(Monte Carlo)是指: 基于采样的数值型近似求解方法

马尔科夫链

又称离散时间马尔可夫链(discrete-time Markov chain，缩写为DTMC),或者马氏链

  • 马尔科夫链(Markov Chain)是指: 状态空间中经过从一个状态到另一个状态的转换的随机过程, 该随机过程满足马尔科夫性
  • 马尔科夫性(Markov property)是指: 当一个随机过程在给定现在状态及所有过去状态情况下，其未来状态的条件概率分布仅依赖于当前状态；换句话说，在给定现在状态时，它与过去状态（即该过程的历史路径）是条件独立的，那么此随机过程即具有马尔可夫性质。具有马尔可夫性质的过程通常称之为马尔可夫过程
    • 马尔科夫性的简单理解: “无记忆性”: 下一状态的概率分布只能由当前状态决定，在时间序列中它前面的事件均与之无关

MCMC基本思想

  • 针对采样的目标分布,构造一个马尔科夫链,使得该马尔科夫链的平稳分布就是目标分布
  • 从任何一个初始状态出发,沿着马尔科夫链进行状态转移,直到马尔科夫链收敛(到达平稳状态)
  • 继续在马尔科夫链上进行状态转移,收敛后继续采样得到的样本就是原始目标分布的样本
    • burn-in处理: 现实应用中,我们需要丢弃收敛前的样本,只保留收敛后的样本
    • burn-in 原意为”老化,定型”之意,在这里表示我们只取后面马氏链定型(收敛)后采样得到的样本
    • 假设采样到收敛用了n次采样,那么服从原始分布的k个样本为\((x^{n+1}, x^{n+2},,,, x^{n+k})\)有时候为了得到近似独立的样本,可以间隔每r次再取出其中一个样本\((x^{n+r}, x^{n+2r},,,, x^{n+kr})\)
    • 真正独立同分布的k个样本可用多条k条不同的收敛后的马氏链得到,不同马氏链采样得到的样本是独立同分布的
  • 核心: 马尔科夫链的构造,也就是确定马尔科夫链的状态转移概率

常见的MCMC方法

MH(Metropolis-Hastings)采样法

  • 对于原始目标分布\(p(x)\)
  • 选择一个参考条件分布\(q(x^{\star}|x)\), 定义接受概率\(A(x,x^{\star})\):
    • (注意这里是参考条件分布,因为是马尔科夫链,所以每个状态都由上一个状态转移而来,需要定义的参考分布应该是条件分布,不是一般拒绝采样中的普通参考分布)
      $$
      A(x,x^{\star}) = min\left ( 1, \frac{p(x^{\star})q(x|x^{\star})}{p(x)q(x^{\star}|x)} \right )
      $$
  • MH采样法构建满足平稳分布就是目标分布\(p(x)\)的秘诀就是让每次采样时,当前状态以一定概率停留在上一状态
    • 与拒绝采样对应: 接受意味着跳转到下一状态,拒绝意味着停留在当前状态
采样过程
  • 随机选取初始样本x^{0}
  • for t = 1,2,3,…:
    • 参考条件分布采样\(x^{\star} \sim q(x^{star}|x^{t-1})\)
    • 均匀分布采样\(u \sim U(0,1)\)
    • 判断是否接受: 如果\(u < A(x^{t-1}, x^{\star})\),则接受,令: \(x^{t} = x^{\star}\), 否则拒绝,令: \(x^{t}=x^{t-1}\)
  • burn-in处理: 丢弃采样到平稳分布前的样本, 只保留平稳分布后的样本即为服从原始分布\(p(x)\)的样本
    • 假设采样到收敛用了n次采样,那么服从原始分布的k个样本为\((x^{n+1}, x^{n+2},,,, x^{n+k})\)有时候为了得到近似独立的样本,可以间隔每r次再取出其中一个样本\((x^{n+r}, x^{n+2r},,,, x^{n+kr})\)
    • 真正独立同分布的k个样本可用多条k条不同的收敛后的马氏链得到,不同马氏链采样得到的样本是独立同分布的
  • 采样次数一般来说是凭经验选择一个足够大的值,现实是现实可以使用一些参数变化量这类的指标来判断采样是否收敛,参考NLP——LLDA的Gibbs采样实现
与拒绝采样的区别
  • MH采样基于拒绝采样来逼近平稳分布
  • 拒绝采样中: 如果样本某一步被拒绝,那么该步不会产生新的样本,需要重新对当前步进行采样
  • MH中: 每一步都会产生一个样本,被拒绝后,就令当前样本和上一个样本相同即可
    • 因为这里是为了使得每个状态的转出概率等于转入概率,所以拒绝就意味着当前采样步骤状态不跳转
    • MH采样法最核心的思想就是一定概率停留在上一个状态来实现对马尔科夫链的构建的
MH采样法正确性证明

  • MH采样法构造的马尔科夫链(状态转移概率矩阵)是正确的吗?

  • 细致平稳条件, 如果非周期的马氏链的状态转移矩阵P和分布\(\pi(x)\)满足下面的式子对任意的\(i,j\)都成立:
    $$\pi(x^{i})P_{ij} = \pi(x^{j})P_{ji}$$

    • 上式为细致平稳分布条件(detailed balance condition)
    • 其中\(\pi(x)\)为马氏链的平稳分布,在这里等于我们的原始分布\(p(x)\)

  • 证明\(\pi(x)\)为马氏链的平稳分布:
    $$
    \begin{align}
    \sum_{i=1}^{\infty}\pi(x^{i})P_{ij} = \sum_{i=1}^{\infty}\pi(x^{j})P_{ji} = \pi(x^{j})\sum_{i=1}^{\infty}P_{ji} = \pi(x^{j}) \\
    => \pi(x) P = \pi(x)
    \end{align}
    $$

    • 由于\(\pi(x)\)为方程\(\pi(x) P = \pi(x)\)的解,所以\(\pi(x)\)是状态转移矩阵P对应的马氏链的平稳分布

  • 在MH采样法中

    • 参考条件分布函数本对应状态转移矩阵的一个元素,\(q(x^{i}|x^{j}) = P_{ij}\)(注意: 实际上一般不相等)
    • 但是很难构造这样方便采样的函数,于是我们使用一个接受率来修正\(q(x^{i}|x^{j})\alpha(x^{j}, x^{i}) = P_{ij}\)
      • \(\alpha(x^{j}, x^{i})\)表示从\(x^{j}\)跳转到\(x^{i}\)的接受率, 其值可如下求得:
        $$
        \begin{align}
        \pi(x^{i}) P_{ij} &= \pi(x^{j})P_{ji} \\
        \pi(x^{i}) q(x^{j}|x^{i})\alpha(x^{i}, x^{j}) &= \pi(x^{j})q(x^{i}|x^{j})\alpha(x^{j}, x^{i})
        \end{align}
        $$
    • 显然,直接取:
      $$
      \begin{align}
      \alpha(x^{i}, x^{j}) &= \pi(x^{j})q(x^{i}|x^{j})\\
      \alpha(x^{j}, x^{i}) &= \pi(x^{i})q(x^{j}|x^{i})\\
      \end{align}
      $$
    • 即可

  • 但是由于\(\alpha(x^{j}, x^{i})\)一般来说太小,所以我们考虑将\(\alpha(x^{j}, x^{i})\)和\(\alpha(x^{i}, x^{j})\)同时扩大M倍,使得其中大的那个为1,即可得到最大的接受率

    • 使用原始接受率的方法称为一般MCMC方法
    • 使用扩大M被接受率的方法称为MCMC的改进方法: MH方法

  • 改进后的接受率为从\(x^{i}\)跳转到\(x^{j}\)的接受率:
    $$
    \begin{align}
    A(x^{i}, x^{j}) = min\left ( 1, \frac{p(x^{j})q(x^{i}|x^{j})}{p(x^{i})q(x^{j}|x^{i})} \right )
    \end{align}
    $$

    • 理解:
      $$
      \begin{align}
      p(x^{j})q(x^{i}|x^{j}) &> p(x^{i})q(x^{j}|x^{i})时: A(x^{i}, x^{j}) = 1 \\
      p(x^{j})q(x^{i}|x^{j}) &< p(x^{i})q(x^{j}|x^{i})时: A(x^{i}, x^{j}) = \frac{p(x^{j})q(x^{i}|x^{j})}{p(x^{i})q(x^{j}|x^{i})}
      \end{align}
      $$

  • 在MH中表现为从\(x\)跳转到\(x^{\star}\)的接受率:
    $$
    A(x,x^{\star}) = min\left ( 1, \frac{p(x^{\star})q(x|x^{\star})}{p(x)q(x^{\star}|x)} \right )
    $$

Gibbs采样法

Gibbs采样是MH采样法的一个特例,每次只更新样本的一个维度

  • 针对维度很高的多维向量,同时采样多个维度难度较高,且接受率很小
  • 使用Gibbs采样每次采样一个维度可以解决这个问题
准备
  • 求得已知其他维度下,每一维度的条件概率\(p(x_{i}|x_{1},,,x_{i-1}, x_{i+1},,,x_{d})\)
采样过程

  • 随机选择初始状态\(x^{0} = (x_{1}^{0}, x_{2}^{0}, x_{3}^{0},,,, x_{d}^{0})\)

  • for t = 1,2,3,…:

    • 基于前一次产生的样本:\(x^{t-1} = (x_{1}^{t-1}, x_{2}^{t-1}, x_{3}^{t-1},,,, x_{d}^{t-1})\), 依次采样和更新每个维度的值,即依次按照:
      • \(x_{1}^{t} \sim p(x_{1}|x_{2}^{t-1},,,x_{d}^{t-1})\)
      • \(x_{2}^{t} \sim p(x_{2}|x_{1}^{t}, x_{3}^{t-1},,,x_{d}^{t-1})\)
      • \(x_{i}^{t} \sim p(x_{i}|x_{1}^{t},,,x_{i-1}^{t}, x_{i+1}^{t-1},,,x_{d}^{t-1})\)
      • \(x_{d}^{t} \sim p(x_{d}|x_{1}^{t}, x_{2}^{t},,,x_{d-1}^{t})\)
    • 得到新的一个样本:\(x^{t} = (x_{1}^{t}, x_{2}^{t}, x_{3}^{t},,,, x_{d}^{t})\)

  • burn-in处理: 丢弃采样到平稳分布前的样本, 只保留平稳分布后的样本即为服从原始分布\(p(x)\)的样本

    • 假设采样到收敛用了n次采样,那么服从原始分布的k个样本为\((x^{n+1}, x^{n+2},,,, x^{n+k})\), 有时候为了得到近似独立的样本,可以间隔每r次再取出其中一个样本\((x^{n+r}, x^{n+2r},,,, x^{n+kr})\)
    • 真正独立同分布的k个样本可用多条k条不同的收敛后的马氏链得到,不同马氏链采样得到的样本是独立同分布的
    • 注意: 采样完成部分维度生成的中间样本如\((x_{1}^{t},,,x_{i-1}^{t}, x_{i}^{t}, x_{i+1}^{t-1},,,x_{d}^{t-1})\)是不能作为最终样本的,因为他们之间(同一轮次的多个中间样本)相互依赖性太强,不具有独立同分布的(虽然完全采样完成的也不能视为具有独立同分布,但是可近似的认为是独立同分布的)

  • 采样次数一般来说是凭经验选择一个足够大的值,现实是现实可以使用一些参数变化量这类的指标来判断采样是否收敛,参考NLP——LLDA的Gibbs采样实现