Question

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

   给定连续随机变量 $ MathJax-Element-184 $ 的概率密度函数 $ MathJax-Element-185 $ ，则 $ MathJax-Element-102 $ 在区间 $ MathJax-Element-186 $ 中的概率可以计算为：

   $ P(\mathbb A)=\int_\mathbb A\tilde p(x)dx $

   如果函数 $ MathJax-Element-191 $ ， 则可以计算 $ MathJax-Element-192 $ 的期望： $ MathJax-Element-208 $ 。

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

     $ \mathbb E_{\vec X\sim \tilde p(\mathbf {\vec x})}[f(\vec X)]=\frac 1N \sum_{i=1}^{N}f(\mathbf {\vec x}_i) $

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

3.1 MCMC 算法

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

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

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

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

   引入一个函数 $ MathJax-Element-260 $ ，使其满足： $ MathJax-Element-265 $ 。如：取 $ MathJax-Element-266 $ ，则有：

   $ \tilde p(i)Q_{i,j}\alpha(i,j)=\tilde p(i)Q_{i,j}\tilde p(j)Q_{j,i}=\tilde p(j)Q_{j,i}\tilde p(i)Q_{i,j}=\tilde p(j)Q_{j,i}\alpha(j,i) $

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

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

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

   • 根据推导 $ MathJax-Element-266 $ ，如果将系数从 1 提高到 $ MathJax-Element-280 $ ，则有：

     $ \alpha^{*}(i,j)=K\tilde p(j)Q_{j,i}=K\alpha (i,j)\\ \alpha^{*}(j,i)=K\tilde p(i)Q_{i,j}=K\alpha (j,i) $

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

4. 将 $ MathJax-Element-289 $ 同比例放大，取： $ MathJax-Element-290 $ 。

   • 当 $ MathJax-Element-292 $ 时： $ MathJax-Element-293 $ ，此时满足细致平稳条件。
   • 当 $ MathJax-Element-294 $ 时： $ MathJax-Element-295 $ ，此时满足细致平稳条件。
   • 当 $ MathJax-Element-296 $ 时： $ MathJax-Element-297 $ ，此时满足细致平稳条件。

5. MH 算法：

   • 输入：

     • 先验转移概率矩阵 $ MathJax-Element-257 $
     • 目标分布 $ MathJax-Element-210 $

   • 输出： 采样出的一个状态序列 $ MathJax-Element-314 $

   • 算法：

     • 初始化 $ MathJax-Element-171 $

     • 对 $ MathJax-Element-300 $ 执行迭代。迭代步骤如下：

       • 根据 $ MathJax-Element-301 $ 采样出候选样本 $ MathJax-Element-302 $ ，其中 $ MathJax-Element-304 $ 为转移概率函数。

       • 计算 $ MathJax-Element-305 $ ：

         $ \alpha( x^{*} \mid x_{t-1})=\min \left(1,\frac{\tilde p( x^{*})Q( x _{t-1} \mid x^{*})}{\tilde p( x _{t-1})Q( x^{*} \mid x _{t-1})} \right) $

       • 根据均匀分布从 $ MathJax-Element-306 $ 内采样出阈值 $ MathJax-Element-307 $ ，如果 $ MathJax-Element-308 $ ，则接受 $ MathJax-Element-311 $ ， 即 $ MathJax-Element-310 $ 。否则拒绝 $ MathJax-Element-311 $ ， 即 $ MathJax-Element-312 $ 。

     • 返回采样得到的序列 $ MathJax-Element-314 $

   注意：返回的序列中，只有充分大的 $ MathJax-Element-173 $ 之后的序列 $ MathJax-Element-315 $ 才是服从 $ MathJax-Element-210 $ 分布的采样序列。

3.2 Gibbs 算法

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

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

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

   $ \tilde p(x_1,y_1)\tilde p(y_2\mid x_1)=\tilde p(x_1)\tilde p(y_1\mid x_1)\tilde p(y_2\mid x_1)\\ \tilde p(x_1,y_2)\tilde p(y_1\mid x_1)=\tilde p(x_1)\tilde p(y_2\mid x_1)\tilde p(y_1\mid x_1) $

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

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

   可以构造状态空间上任意两点之间的转移概率矩阵 $ MathJax-Element-257 $ ： 对于任意两点 $ MathJax-Element-414 $ ， 令从 $ MathJax-Element-415 $ 转移到 $ MathJax-Element-416 $ 的概率为 $ MathJax-Element-417 $ ：

   • 如果 $ MathJax-Element-420 $ ，则 $ MathJax-Element-421 $ 。
   • 如果 $ MathJax-Element-423 $ ，则 $ MathJax-Element-424 $ 。
   • 否则 $ MathJax-Element-425 $ 。

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

   $ \tilde p(A)Q(A\rightarrow B)=\tilde p(B)Q(B\rightarrow A) $

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

3. Gibbs 算法：

   • 输入：目标分布 $ MathJax-Element-462 $ ，其中 $ MathJax-Element-635 $

   • 输出： 采样出的一个状态序列 $ MathJax-Element-641 $

   • 算法步骤：

     • 初始化：随机初始化 $ MathJax-Element-511 $ 。

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

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

         • 将 $ MathJax-Element-670 $ 保持不变，计算条件概率： $ MathJax-Element-715 $ 。

           该条件概率就是状态转移概率 $ MathJax-Element-417 $

         • 根据条件概率 $ MathJax-Element-715 $ 对分量 $ MathJax-Element-680 $ 进行采样，假设采样结果为 $ MathJax-Element-581 $ ，则获得新样本 $ MathJax-Element-730 $ 。

         • 令 $ MathJax-Element-655 $ ，继续遍历。

     • 最终返回一个状态序列 $ MathJax-Element-641 $ 。

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