三、MCMC 采样

1. 概率图模型中最常用的采样技术是马尔可夫链蒙特卡罗方法Markov Chain Monte Carlo:MCMC。

   给定连续随机变量 的概率密度函数 ，则 在区间 中的概率可以计算为：

   

   如果函数 ， 则可以计算 的期望： 。

   • 如果 不是一个单变量，而是一个高维的多元变量 ，且服从一个非常复杂的分布，则对于上式的求积分非常困难。为此，MCMC先构造出服从 分布的独立同分布随机变量 ， 再得到 的无偏估计：

     

   • 如果概率密度函数 也很复杂，则构造服从 分布的独立同分布随机变量也很困难。MCMC 方法就是通过构造平稳分布为 的马尔可夫链来产生样本。

3.1 MCMC 算法

1. MCMC 算法的基本思想是：先设法构造一条马尔可夫链，使其收敛到平稳分布恰好为 。然后通过这条马尔可夫链来产生符合 分布的样本。最后通过这些样本来进行估计。

   这里马尔可夫链的构造至关重要，不同的构造方法将产生不同的MCMC算法。Metropolis-Hastings:MH算法是MCMC的重要代表。

2. 假设已经提供了一条马尔可夫链，其转移矩阵为 。目标是另一个马尔科夫链，使转移矩阵为 、平稳分布是 。

   通常 ，即 并不满足细致平稳条件不成立。但是可以改造已有的马尔可夫链，使得细致平稳条件成立。

   引入一个函数 ，使其满足： 。如：取 ，则有：

   

   令： ，则有 。其中 构成了转移矩阵 。而 恰好满足细致平稳条件，因此它对应的马尔可夫链的平稳分布就是 。

3. 在改造 的过程中引入的 称作接受率。其物理意义为：在原来的马尔可夫链上，从状态 以 的概率跳转到状态 的时候，以 的概率接受这个转移。

   • 如果接受率 太小，则改造马尔可夫链过程中非常容易原地踏步，拒绝大量的跳转。这样使得马尔可夫链遍历所有的状态空间需要花费太长的时间，收敛到平稳分布 的速度太慢。

   • 根据推导 ，如果将系数从 1 提高到 ，则有：

     

     于是： 。因此，即使提高了接受率，细致平稳条件仍然成立。

4. 将 同比例放大，取： 。

   • 当 时： ，此时满足细致平稳条件。
   • 当 时： ，此时满足细致平稳条件。
   • 当 时： ，此时满足细致平稳条件。

5. MH 算法：

   • 输入：

     • 先验转移概率矩阵
     • 目标分布

   • 输出： 采样出的一个状态序列

   • 算法：

     • 初始化

     • 对 执行迭代。迭代步骤如下：

       • 根据 采样出候选样本 ，其中 为转移概率函数。

       • 计算 ：

         

       • 根据均匀分布从 内采样出阈值 ，如果 ，则接受 ， 即 。否则拒绝 ， 即 。

     • 返回采样得到的序列

   注意：返回的序列中，只有充分大的 之后的序列 才是服从 分布的采样序列。

3.2 Gibbs 算法

1. MH算法不仅可以应用于一维空间的采样，也适合高维空间的采样。

   对于高维的情况，由于接受率 的存在（通常 ），MH算法的效率通常不够高，此时一般使用 Gibbs 采样算法。

2. 考虑二维的情形：假设有概率分布 ，考察状态空间上 坐标相同的两个点 ，可以证明有：

   

   于是 。则在 这条平行于 轴的直线上，如果使用条件分布 作为直线上任意两点之间的转移概率，则这两点之间的转移满足细致平稳条件。

   同理：考察 坐标相同的两个点 ， 有 。在 这条平行于 轴的直线上，如果使用条件分布 作为直线上任意两点之间的转移概率，则这两点之间的转移满足细致平稳条件。

   可以构造状态空间上任意两点之间的转移概率矩阵 ： 对于任意两点 ， 令从 转移到 的概率为 ：

   • 如果 ，则 。
   • 如果 ，则 。
   • 否则 。

   采用该转移矩阵 ，可以证明：对于状态空间中任意两点 ，都满足细致平稳条件：

   

   于是这个二维状态空间上的马尔可夫链将收敛到平稳分布 ，这就是吉布斯采样的原理。

3. Gibbs 算法：

   • 输入：目标分布 ，其中

   • 输出： 采样出的一个状态序列

   • 算法步骤：

     • 初始化：随机初始化 。

     • 执行迭代，迭代步骤如下：

       • 随机或者以一定次序遍历索引 。遍历过程为（设当前索引为 ）：

         • 将 保持不变，计算条件概率： 。

           该条件概率就是状态转移概率

         • 根据条件概率 对分量 进行采样，假设采样结果为 ，则获得新样本 。

         • 令 ，继续遍历。

     • 最终返回一个状态序列 。

4. 吉布斯采样Gibbs sampling 有时被视作MH算法的特例，它也使用马尔可夫链获取样本。