机器学习-白板推导系列(十三)- 马尔可夫链&蒙特卡洛方法(MCMC, Markov Chain&Monte Carlo)之概念及MH算法

13. 马尔可夫链&蒙特卡洛方法(MCMC)

13.1 采样方法介绍

• 上一章讲贝叶斯角度主要是计算$P(Z|X)$，大致采用两种思路来解决这个问题，也就是精确推断和近似推断。精确推断无法达到我们想要的结果时，就会采用近似推断的方法。而近似推断中我们又可以分成两大类，即为确定性近似(VI)和随机近似(MCMC)。

• Monte Carlo Method是一种基于采样的随机近似算法。目标是求解后验概率$P(Z|X)$，其中$Z$为Latent data，$X$为Observed data。知道分布后，通常的目标是求解：
  $${\mathbb{E}}_{Z|X}[f(Z)]={\int }_{Z}P(Z|X)f(Z)dZ\approx \frac{1}{N}\sum _{i=1}^{N}f(z_i)\text{\hspace{1em}(13.1.1)}$$

  那么知道了后验分布$P(Z|X)$，怎么去采样呢？即：如何通过采样得到$z^{(1)},z^{(2)},\cdots ,z^{(N)}\sim P(Z|X)$。这一节将要主要介绍三种采样方法：概率分布采样，拒绝采样和重要性采样。

13.1.1 概率分布采样

1. 缩写强调

   • $c.d.f(概率分布函数)$；
   • $p.d.f(概率密度函数)$；
   • $i.i.d(独立同分布)$。

2. 介绍

   • 已知：$p.d.f(概率密度函数)$ ；求：采样 { x ( 1 ) , x ( 2 ) , ⋯   , x ( N ) } \{ x^{(1)},x^{(2)},\cdots,x^{(N)} \} {x(1),x(2),⋯,x(N)} $N$个样本点。
   • 思路：
     • 由$p.d.f(概率密度函数)$，求$c.d.f(概率分布函数)$；
     • 由于$c.d.f(概率分布函数)$的值一定是$[0,1]$之间，所以可以在均匀概率密度分布$U(0,1)$上采样，，得到$u^{(i)}\sim U(0,1)$；
     • 求$x^{(i)}\sim cd{f}^{-1}(u^{(i)})$就可以计算得到采样的 { x ( 1 ) , x ( 2 ) , ⋯   , x ( N ) } \{ x^{(1)},x^{(2)},\cdots,x^{(N)} \} {x(1),x(2),⋯,x(N)} $N$个样本点。
       

3. 缺陷

   • 很多情况根本不知道$p.d.f(概率密度函数)$的具体表现形式。
   • 很多时候$c.d.f(概率分布函数)$也并不好求。

   所以很多情况下，概率分布采样并没有那么的好求。

13.1.2 拒绝采样(Rejection Sampling)

1.引入概念

• 提议分布（proposal distribution）$q(z)$：较简单的分布；
• 接受率：$\alpha =\frac{P(z^{(i)})}{M\cdot q(z^{(i)})}$，$0\le \alpha \le 1$，实际是上图中绿色的部分。
  

  绿色的部分称为$拒绝区域$，蓝色部分称为$接受区域$。采样点的概率要到$蓝色区域内$，所以这个采样方法就是拒绝采样就是这样来的。

2. 介绍
   对于较复杂的概率分布$p(z)$，引入简单的$提议分布(proposaldistribution)$$q(z)$，具体的采样方法步骤为：
   • ①选择概率密度函数为$q(z)$作为提议分布，使其 $\forall z$满足$Mq(z_i)\ge p(z_i)$，其中$M>0$；
   • ②按照$提议分布q(z)$随机抽样得到$样本z_i$，再按照$均匀分布在(0,1)$范围内抽样得到$u_i$；
   • ③如果$u_i\le \frac{p(z_i)}{Mq(z_i)}$，则将$z_i$作为抽样结果；否则，返回步骤②；
   • ④直到获得$N$个样本，结束。
3. 缺陷
   $q(z)$很难选取；采样效率可能不高。
   • 如果$M\cdot q(z)$比$p(z)$大很多的话，那么采样老是是失败的，这样采样效率低下。
   • 当$M\cdot q(z)=p(z)$的时候，$\alpha =1$，我们每次采样的结果都是接受的。
   • 当$p(z)$的分布形式非常的复杂，根本没有办法得到准确的结果(特别是采样cost非常高)，经常性的采样失败带来的损失是很大的。

13.1.3 重要性采样（Importance Sampling）

1. 重要性采样（Importance Sampling）
   • 重要性采样在我们的强化学习(PPO)中的应用非常的多。重要性采样并不是直接对概率进行采样，而是对概率分布的期望进行采样。引入另一个分布$q(z)$：
     $$\begin{array}{ll}{\mathbb{E}}_{p(z)}[f(z)]& =\int p(z)f(z)dz\\ & =\int \frac{p(z)}{q(z)}q(z)f(z)dz\\ & =\int f(z)\frac{p(z)}{q(z)}q(z)dz\\ & \approx \frac{1}{N}{\sum }_{i=1}^{N}f(z_i)\cdot \underset{weight}{\underbrace{\frac{p(z_i)}{q(z_i)}}}\\ & z_i\sim q(z),i=1,2,\cdots ,N\end{array}\text{\hspace{1em}(13.1.2)}$$
     这里的$\frac{p(z_i)}{q(z_i)}$也就是$Weight$，用来平衡不同的概率密度值之间的差距。采样在$q(z)$中采样，并通过权重计算和。
   • 当两个分布之间的差距太大了话，总是采样采不到重要的样本，采的可能都是实际分布概率值小的部分。也就是采样效率不均匀的问题。在这个基础上，我们进一步提出了Sampling Importance Resampling。
2. 重要性重采样(Sampling Importance Resampling)
   经过重要性采样后，得到了$N$个样本点以及对应的权重。重要性重采样(Sampling Importance Resampling) 用权重来作为采样的概率，重新测采样出$N$个样本。也就是如下图所示：
   
   通过二次采样可以降低采样不平衡的问题。至于为什么呢？
   • $\frac{p(z_i)}{q(z_i)}$是Weight，如果Weight比较大的话，说明$p(z_i)$比较大而$q(z_i)$比较的小，也就是通过$q(z_i)$采出来的数量比较少。
   • 那么按权重再来采一次，就可以增加采到重要性样本的概率，成功的弥补了重要性采样带来的缺陷，有效的弥补采样不均衡的问题。

---

13.2 马尔可夫链(Markov Chain)

13.2.1 基本概念

1. 随机过程
   Random Process研究的对象是一个随机变量的序列 X = { X 0 , X 1 , ⋯   , X t , ⋯   } X=\left \{X_{0},X_{1},\cdots ,X_{t},\cdots \right \} X={X0​,X1​,⋯,Xt​,⋯}。随机过程就是一个序列，这个序列中的每一个元素都是一个随机变量。
2. Markov Chain
   • Markov Chain就是一个特殊的随机过程，它的$时间$和$状态$都是$离散的$。
   • Markov Chain满足Markov性质，即：$未来$和$过去$是$无关的$。用数学的语言表达是：
     $$P(X_{t+1}=x|X_1,X_2,\cdots ,X_t)=P(X_{t+1}|X_1,X_2,\cdots ,X_{t-m})\text{\hspace{1em}(13.2.1)}$$
     公式(13.2.1)是$m阶马尔可夫性质$。
   • 齐次马尔可夫链
     当公式(13.2.1)的$m=0$时，就得到了齐次(一阶)马尔可夫链，即满足：
     $$P(X_{t+1}=x|X_1,X_2,\cdots ,X_t)=P(X_{t+1}|X_t)\text{\hspace{1em}(13.2.2)}$$
     条件概率分布$P(X_{t+1}|X_t)$称为马尔可夫链的转移概率分布。
   • 时间齐次的马尔可夫链（Time Homogenous Markov Chain）
     当转移概率分布$P(X_t|X_{t-1})$与$t$无关，也就是说不同时刻的转移概率是相同的，则称该马尔可夫链为时间齐次的马尔可夫链（Time Homogenous Markov Chain），形式化的表达是：
     $$P(X_{t+s}|X_{t-1+s})=P(X_t|X_{t-1}),t=1,2,\cdots ;s=1,2,\cdots \text{\hspace{1em}(13.2.3)}$$
3. 转移概率矩阵和状态分布
   • 转移概率矩阵
     如果马尔可夫链的随机变量$X_t(t=0,1,2,\cdots )$定义在$离散空间$，则转移概率分布可以由矩阵表示。若马尔可夫链在时刻$t-1$处于状态$j$，在时刻$t$移动到状态$i$，将转移概率记作：
     $$p_{ij}=(X_t=i|X_{t-1}=j),i=1,2,\cdots ;j=1,2,\cdots \text{\hspace{1em}(13.2.4)}$$
     满足：
     $$p_{ij}\ge 0,\sum _i{p}_{ij}=1\text{\hspace{1em}(13.2.5)}$$
     马尔可夫链的转移概率可以由矩阵表示：
     $$P=\left[\begin{array}{cccc}{p}_{11}& p_{12}& p_{13}& \cdots \\ p_{21}& p_{22}& p_{23}& \cdots \\ p_{31}& p_{32}& p_{33}& \cdots \\ \cdots & \cdots & \cdots & \cdots \end{array}\right](p_{ij}\ge 0,\sum _i{p}_{ij}=1)\text{\hspace{1em}(13.2.6)}$$
   • 状态分布
     • 考虑马尔可夫链 X = { X 0 , X 1 , ⋯   , X t , ⋯   } X=\left \{X_{0},X_{1},\cdots ,X_{t},\cdots \right \} X={X0​,X1​,⋯,Xt​,⋯}在时刻$X_t(t=0,1,2,\cdots )$的概率分布，称为时刻$t$的状态分布，记作：
       $$\pi (t)=\left[\begin{array}{c}{\pi }_1(t)\\ {\pi }_2(t)\\ ⋮\end{array}\right]\text{\hspace{1em}(13.2.7)}$$
       其中${\pi }_i(t)$表示时刻$t$状态为$i$的概率$P(X_i=i)$，即：
       $${\pi }_i(t)=P(X_i=i),i=1,2,\cdots \text{\hspace{1em}(13.2.8)}$$
     • 对于马尔可夫链的初始状态分布：
       $$\pi (0)=\left[\begin{array}{c}{\pi }_1(0)\\ {\pi }_2(0)\\ ⋮\end{array}\right]\text{\hspace{1em}(13.2.9)}$$
       其中${\pi }_i(0)$表示时刻$0$状态为i的概率$P(X_0=i)$。通常初始分布${\pi }_i(0)$的向量只有一个分量是$1$，其余分量为$0$，表示马尔可夫链是从一个具体状态开始的。
     • 马尔可夫链在时刻$t$的状态分布，可以由在时刻$t-1$的状态分布以及转移概率分布决定：
       $$\pi (t)=P\pi (t-1)\text{\hspace{1em}(13.2.10)}$$
       其中：
       $$\begin{array}{ll}{\pi }_i(t)& =P(X_t=i)\\ & ={\sum }_{m}P(X_t=i|X_{t-1}=m)P(X_{t-1}=m)\\ & ={\sum }_m{P}_{im}{\pi }_m(t-1)\end{array}\text{\hspace{1em}(13.2.11)}$$
       马尔可夫链在时刻t的状态分布，可以通过递推得到：
       $$\pi (t)=P\pi (t-1)=P(P\pi (t-2))=P^2\pi (t-2)\text{\hspace{1em}(13.2.12)}$$
       递推得到：
       $$\pi (t)=P^t\pi (0)\text{\hspace{1em}(13.2.13)}$$

       1. 这里的递推式表明马尔可夫链的状态分布由$初始分布$和$转移概率分布$决定。
       2. 这里的$P^t$称为$t步转移概率矩阵$：
          $$P_{ij}^t=P(X_t=i|X_0=j)\text{\hspace{1em}(13.2.14)}$$

     • 给定一个马尔可夫链$X$，状态分布为$\pi =({\pi }_1,{\pi }_2,\cdots {)}^T$为$X$的$充要条件$为$\pi$是下列方程组的解：
       $$\{\begin{array}{l}①x_i={\sum }_j{p}_{ij}{x}_j,i=1,2,\cdots \\ ②x_i\ge 0,i=1,2,\cdots \\ ③{\sum }_i{x}_i=1\end{array}\text{\hspace{1em}(13.2.15)}$$

     1. 必要性证明
        假设$\pi$是平稳分布，显然满足①和②，又因为${\pi }_i={\sum }_j{p}_{ij}{\pi }_j$，所以满足③。
     2. 充分性证明
        由②和③知$\pi$是一概率分布。假设$\pi$是$X_t$的分布，则有：
        $$P(X_t=i)={\pi }_i=\sum _j{p}_{ij}P(X_{t-1}=j)\text{\hspace{1em}(13.2.15)}$$
        又因为$\pi$满足①，所以有：
        $$P(X_t=i)={\pi }_i=\sum _j{p}_{ij}{\pi }_j\text{\hspace{1em}(13.2.16)}$$
        综合两式则可得出$P(X_{t-1}=j)={\pi }_j,j=1,2,\cdots$，即$\pi$也是$X_{t-1}$的概率分布。事实上这对任意$t$都成立，所以$\pi$是马尔可夫链的平稳分布。

     5. 平稳分布
        对于一个马尔可夫链$X$，如果其状态空间上存在一个分布：
        $$\pi =\left[\begin{array}{c}{\pi }_1\\ {\pi }_2\\ ⋮\end{array}\right]\text{\hspace{1em}(13.2.17)}$$
        使得：$\pi =P\pi$
        则称$\pi$为马尔可夫链$X$的平稳分布。直观地来看，如果以该平稳分布作为初始分布，面向未来进行随机状态转移，之后任何一个时刻的状态分布都是该平稳分布。具体内容在（13.2.2）中会描述。
4. 连续状态马尔可夫链
   1. 概率转移核
      • 连续状态马尔可夫链的转移概率分布由概率转移核或转移核（transition kernel）表示。
      • 在连续状态空间$S$上，对任意的$x\in S,A\subset S$（$A$可以理解为一个区间），转移核$P(x,A)$定义为：
        $$P(x,A)={\int }_{A}p(x,y)\mathrm{d}y\text{\hspace{1em}(13.2.18)}$$
        其中$p(x,\cdot )$为概率密度函数，满足$p(x,\cdot )\ge 0,P(x,S)={\int }_{S}p(x,y)\mathrm{d}y=1$。转移核$P(x,A)$表示从$x\sim A$的转移概率：
        $$P(X_t=A|X_{t-1}=x)=P(x,A)\text{\hspace{1em}(13.2.19)}$$

        有时也将概率密度函数$p(x,\cdot )$称为转移核。

      2. 若马尔可夫链的状态空间$S$上的概率分布$\pi (x)$满足条件：
         $$\pi (y)=\int p(x,y)\pi (x)\mathrm{d}x,\forall y\in S\text{\hspace{1em}(13.2.20)}$$
         则$\pi (x)$为该马尔可夫链的平稳分布。等价地：
         $$\pi (A)=\int p(x,A)\pi (x)\mathrm{d}x,\forall A\in S\text{\hspace{1em}(13.2.21)}$$
         或简写为：
         $$\pi =P\pi \text{\hspace{1em}(13.2.22)}$$
5. 马尔可夫链的性质
   • 不可约
     在状态空间$S$中对于任意状态$i,j\in S$，如果存在一个时刻$t(t>0)$满足：
     $$P(X_t=i|X_0=j)>0\text{\hspace{1em}(13.2.23)}$$
     也就是说，时刻$0$从状态$j$出发，时刻$t$到达状态$i$的概率$大于0$，则称此马尔可夫链是不可约的（irreducible），否则称马尔可夫链是可约的（reducible）。

     直观上：$一个不可约的马尔可夫链，从任意状态出发，当经过充分长时间后，可以到达任意状态。$

     • 不可约：
       
     • 可约：
       

   2. 非周期
      在状态空间$S$中对于任意状态$i\in S$，如果时刻$0$从状态$i$出发，$t$时刻返回状态的所有时间长 { t : P ( X t = i ∣ X 0 = i ) > 0 } \color{blue}\left \{t:P(X_{t}=i|X_{0}=i)> 0\right \} {t:P(Xt​=i∣X0​=i)>0}的最大公约数是$1$，则称此马尔可夫链是非周期的（aperiodic），否则称马尔可夫链是周期的（periodic）。

      直观上，一个非周期性的马尔可夫链，不存在一个状态，从这一个状态出发，再返回到这个状态时所经历的时间长呈一定的周期性，也就是说$非周期性的马尔可夫链的任何状态都不具有周期性$。

      • 非周期：
        
      • 周期：
        

   3. 正常返
      对于任意状态$i,j\in S$，定义概率$p_{ij}^t$为时刻$0$从状态$j$出发，时刻$t首次$转移到状态$i$的概率，即$p_{ij}^t=P(X_t=i,X_s\ne i,s=1,2,\cdots ,t-1|X_0=j),t=1,2,\cdots$。若对所有状态$i,j$都满足${\lim }_{t\to \infty }{p}_{ij}^t>0$，则称马尔可夫链是正常返的（positive recurrent）。

      • 直观上，一个正常返的马尔可夫链，其中任意一个状态，从其他任意一个状态出发，当时间趋于无穷时，首次转移到这个状态的概率不为$0$。
      • 定理：$不可约、非周期且正常返的马尔可夫链，有唯一平稳分布存在。$

   4. 遍历定理

      • 定理：
        若马尔可夫链是不可约、非周期且正常返的，则该马尔可夫链有唯一平稳分布$\pi =({\pi }_1,{\pi }_2,\cdots {)}^T$，并且转移概率的极限分布是马尔可夫链的平稳分布：
        $$\underset{t\to \infty }{\lim }P(X_t=i|X_0=j)={\pi }_i,i=1,2,\cdots ,j=1,2,\cdots \text{\hspace{1em}(13.2.24)}$$

        也就是：$$\underset{t\to \infty }{\lim }P(X_t=i)={\pi }_i,i=1,2,\cdots$$

        若$f(X)$是定义在状态空间上的函数，$E_{\pi }[f(X)]<\infty$，则：
        $$P\left\{{\hat{f}}_t\to E_{\pi }[f(X)]\right\}=1\text{\hspace{1em}(13.2.25)}$$
        此处：

        • ${\hat{f}}_t=\frac{1}{t}{\sum }_{s=1}^{t}f(x_s)$；
        • $E_{\pi }[f(X)]={\sum }_{i}f(i){\pi }_i$是$f(X)$关于平稳分布$\pi =({\pi }_1,{\pi }_2,\cdots {)}^T$的数学期望；
        • $P\left\{{\hat{f}}_t\to E_{\pi }[f(X)]\right\}=1$表示${\hat{f}}_t\to E_{\pi }[f(X)],t\to \infty$几乎处处成立或以概率$1$成立。
      2. 意义
         • 遍历定理的直观解释：满足相应条件的马尔可夫链，当时间趋于无穷时，马尔可夫链的状态分布趋近于平稳分布，随机变量的函数的样本均值以概率收敛于该函数的数学期望。
         • 样本均值可以认为是时间均值，数学期望是空间均值。遍历定理表述了遍历性的含义：当时间趋于无穷时，时间均值等于空间均值。
         • 遍历定理的三个条件：不可约、非周期、正常返，保证了当时间趋于无穷时达到任意一个状态的概率不为。
      3. 实际应用
         理论上并不知道经过多少次迭代，马尔可夫链的状态分布才能接近平稳分布，在实际应用遍历定理时，取一个足够大的整数$m$，经过$m$次迭代之后认为状态分布就是平稳分布，这时计算从第$m+1$次迭代到第$n$次迭代的均值，即：
         $$E_{\pi }[f(X)]=\frac{1}{n-m}\sum _{i=m+1}^{n}f(x_i)\text{\hspace{1em}(13.2.26)}$$
         称为遍历均值。
6. 可逆马尔可夫链

   • 定义
     对于任意状态i,j\in S，对任意一个时刻t满足：
     $$P(X_t=i|X_{t-1}=j){\pi }_j=P(X_{t-1}=j|X_t=i){\pi }_i,i,j=1,2,\cdots \text{\hspace{1em}(13.2.27)}$$
     或简写为：
     $$p_{ij}{\pi }_j=p_{ji}{\pi }_i,i,j=1,2,\cdots \text{\hspace{1em}(13.2.28)}$$
     则称此马尔可夫链为可逆马尔可夫链（reversible Markov chain），上式又被称作细致平衡方程（detailed balance equation）。13.2.2中会讲解到。

     直观上，如果有可逆马尔可夫链，那么以该马尔可夫链的平稳分布作为初始分布，进行随机状态转移，无论是面向过去还是面向未来，任何一个时刻的状态分布都是该平稳分布。

   • 定理
     满足细致平衡方程的状态分布$\pi$就是该马尔可夫链的平稳分布，即满足$P\pi =\pi$。

     证明：
     $(P\pi {)}_i={\sum }_j{p}_{ij}{\pi }_j={\sum }_j{p}_{ji}{\pi }_i={\pi }_i\underset{=1}{\underbrace{\sum _j{p}_{ji}}}={\pi }_i,i=1,2,\cdots$

     该定理说明，可逆马尔可夫链一定有唯一平稳分布，给出了一个马尔可夫链有平稳分布的充分条件（不是必要条件）。也就是说，可逆马尔可夫链满足遍历定理的条件。

---

13.2.2 平稳分布

1. Markov Chain 理解
   一个Markov Chain可以用下图来表示：
   
   此图就是一个时间序列，$x_i$就表示在第$i$时刻的状态，而每一个状态都是一个随机变量。而${\pi }_i$描述的就是第$i$个随机变量的分布。对于一个马氏链来讲，它在第$t+1$时刻的概率分布，可以被我们表达为：
   $${\pi }_{t+1}(x^{\ast })=\int {\pi }_t(x)\cdot P(x\mapsto x^{\ast })dx\text{\hspace{1em}(13.2.29)}$$
   此处$x^{\ast }$表示$t+1$时刻随机变量的某一个值，$x$表示$t$时刻随机变量。即：所有可能转移到状态$x_{t+1}$的概率分布的和。
2. 平稳分布定义
   假如存在一个$\pi$，这里的$\pi$和前面的${\pi }_t$和${\pi }_{t+1}$都$没有关系$。假如$\pi$是一个概率分布，那么它写成一个无限维向量的形式：
   $$\pi =[\pi (1),\pi (2),\cdots ,\pi (t),\cdots ],\sum _{i=1}^{\infty }\pi (i)=1\text{\hspace{1em}(13.2.30)}$$
   如果式(13.2.16)成立：
   $$\pi (x^{\ast })=\int \pi (x)\cdot P(x\mapsto x^{\ast })dx\text{\hspace{1em}(13.2.31)}$$
   我们就称 { π ( k ) } k = 1 ∞ \{ \pi(k) \}_{k=1}^\infty {π(k)}k=1∞​是马氏链 { x k } \{x_{k}\} {xk​}的$平稳分布$。

   用通俗的话讲就是：对于一个马氏链，每个时刻都符合一个概率分布，如果每一个时刻的概率的分布都是一样的，都是$\pi (k)$，那么就可以称这个马氏链满足一个平稳分布。

3. 引入平稳分布的意义
   • 本章开头，我们引入MCMC的目标是求$P(Z)$，此处$P(Z)$可以看成是一个平稳分布$\pi (k)$，那么就可以通过构建一系列的 { x 1 , x 2 , ⋯   , x t , ⋯   } \{ x_1,x_2,\cdots,x_t,\cdots \} {x1​,x2​,⋯,xt​,⋯}的马氏链，让其来逼近这个平稳分布。
   • 构建的马氏链，包括随机变量和转移矩阵，如果它满足平稳分布的条件，确实是可以收敛到平稳分布的。那么，我就可以让构建出来的这个马氏链收敛到平稳分布来求得$P(Z)$。
   • 如何构建平稳分布
     引入一个条件，也就是Detailed Balance：
     $$\pi (x)\cdot P(x\mapsto x^{\ast })=\pi (x^{\ast })\cdot P(x^{\ast }\mapsto x)\text{\hspace{1em}(13.2.32)}$$
     即：对于任意两个状态之间，使用概率分布称为转移概率得到的结果都是可逆的，那么这两个状态之间的分布一定是一样的。
     $$DetailedBalance\Rightarrow 平稳分布(充分不必要条件)\text{\hspace{1em}(13.2.33)}$$

   如果一个分布满足Detailed Balance那么它一定可以是一个平稳分布，但是反过来并不能成立。而证明过程如下：
   $$\begin{array}{ll}\int \pi (x)\cdot P(x\mapsto x^{\ast })dx& =\int \pi (x^{\ast })\cdot P(x^{\ast }\mapsto x)dx\\ & =\pi (x^{\ast })\underset{{\sum }_{j=1}^{\infty }{P}_{ij}=1}{\underbrace{\int P(x^{\ast }\mapsto x)}}dx\\ & =\pi (x^{\ast })\end{array}$$

---

13.3 MH Algorithm

13.3.1 引子

• Monte Carlo Method
  Monte Carlo Method是一种基于采样的随机近似算法。目标是求解后验概率$P(Z|X)$，其中$Z$为Latent data，$X$为Observed data。知道分布后，通常的目标是求解：
  $${\mathbb{E}}_{Z|X}[f(Z)]={\int }_{Z}P(Z|X)f(Z)dZ\approx \frac{1}{N}\sum _{i=1}^{N}f(z_i)\text{\hspace{1em}(13.3.1)}$$
  常用的方法有：概率分布采样，拒绝采样和重要性采样。
• 马氏链（平稳分布）
  1. 对于遍历的马尔可夫链，当其转移次数足够多时，将会进入稳态分布$\pi (x)$，即各状态的出现概率趋于稳定。
  2. 进一步，假设变量$X$所有可能的取值，构成某遍历马氏链的状态空间。
     则无论从什么状态出发，只要转移步数足够大，该马氏链将趋于稳定，即逐渐开始依次输出符合$p(x)$分布的状态$X_i$。
  3. 这样，通过收集马氏链收敛后产生的$X_i$，我们就得到了符合分布$p(x)$的样本。
     此时稳态分布即：$$\pi (x)=p(x)\text{\hspace{1em}(13.3.2)}$$
• MH 算法
  1. MH 算法是利用Monte Carlo Method中的拒绝采样和马氏链中的平稳分布的$采样方法$。
  2. MH算法是Metropolis 算法的一个常用的改进变种。
  3. 由于马氏链的收敛性质主要由转移矩阵$P$决定, 所以基于马氏链做采样的关键问题是：如何构造转移矩阵$P$，使平稳分布恰好是我们要的分布。

13.3.2 概念

1. Detailed Balance
   上一节中我们讲解了Detailed Balance，这是平稳分布的充分必要条件。Detailed Balance为：
   $$\pi (x)P(x\mapsto x^{\ast })=\pi (x^{\ast })P(x^{\ast }\mapsto x)\text{\hspace{1em}(13.3.3)}$$
   这里的$P(x\mapsto x^{\ast })$实际上就是条件概率$P(x^{\ast }|x)$，这样写只是便于理解。换一种写法可以是：
   • 定理——细致平稳条件
     若非周期马氏链的转移矩阵$P$和分布$\pi (x)$满足：
     $$\pi (i)P_{ij}=\pi (j)P_{ji},\forall i,j\text{\hspace{1em}(13.3.4)}$$
     $i,j$表示相应时间的状态，则分布$\pi (x)$是马氏链的平稳分布。
   • 这说明：
     $$\pi (x)P=\pi (x)\text{\hspace{1em}(13.3.5)}$$
     因此$\pi (x)$是平稳分布。
   • 物理意义：
     对于任意两个状态$i,j$，从状态$i$流失到状态$j$的概率质量，等于反向流动的概率质量，因此是状态概率是稳定的。
2. 目标
   推断问题的核心目标：求的是后验概率分布$P(Z)$。求$P(Z)$主要是为了求在$P(Z)$概率分布下的一个相关函数的期望，也就是：
   $${\mathbb{E}}_{P(Z)}[f(Z)]\approx \frac{1}{N}\sum _{i=1}^{N}f(z^{(i)})\text{\hspace{1em}(13.3.6)}$$
   • 我们的目标是 采 样 来 得 到 P ( Z ) ∼ { z ( 1 ) , z ( 2 ) , ⋯   , z ( N ) } 样 本 点 \color{red}采样来得到P(Z) \sim \{ z^{(1)},z^{(2)},\cdots, z^{(N)} \}样本点 采样来得到P(Z)∼{z(1),z(2),⋯,z(N)}样本点。
   • 如果$\pi (x)$是最终的平稳分布($满足DetailedBalance条件$)，求$P(Z)$则可以转换为：求出概率转移矩阵$P_{ij}$。
   • 如果是平稳分布，对不同时刻进行采样，利用马氏链转换到需要的平稳时刻。最终就可以转换到$P(Z)$，MH算法则是如何构则是如何采样，得到想要的关于$P(Z)$的$N$个样本。

13.3.3 M-H算法实现

1. M-H算法实现

   • 假设我们已经有了一个转移矩阵$Q$。一般情况下，该转移矩阵不满足细致平稳条件：
     $$P(i)Q_{ij}\ne P(j)Q_{ji}\text{\hspace{1em}(13.3.7)}$$
     我们引入接受率$\alpha (i,j)$，满足：
     $$P(i)Q_{ij}\alpha (i,j)=P(j)Q_{ji}\alpha (j,i)\text{\hspace{1em}(13.3.8)}$$
     其中：

     • $P(i)$是状态$i$出现的概率，是$假设马氏链满足复杂分布时输出状态的概率$。
     • $Q$是初始转移概率矩阵，$满足任意一个给定的简单分布$，便于抽样即可。

   • 显然，接受率最简单的构造方法是：
     $$\begin{array}{l}\alpha (i,j)=P(j)Q_{ij}\\ \alpha (j,i)=P(i)Q_{ji}\end{array}\text{\hspace{1em}(13.3.9)}$$
     此时，细致平稳条件成立，$该马氏链输出的状态满足稳态分布p(x)！$

   • 注意
     综上，有几个要点必须要理解，有助于我们编程实现：

     • 我们引入接受率，使得该马氏链$以转移概率Q_{ij}\alpha (i,j)从状态i转移到状态j$。
     • 该转移概率的实现方式是：我们$先按Q_{ij}生成新状态，再按\alpha (i,j)的概率接受转移结果$，乘积即为$Q_{ij}\alpha (i,j)$转移概率。
     • 一定要理解这个接受率！
       • 在拒绝抽样中，如果$u_0>p(z)$，我们会放弃抽样结果，不纳入最终的统计。正是因为如此，抽样才会逼近复杂分布。
       • 而这里的接受是$“接受转移”$，如果$u>\alpha$，意味着$t+1$时刻状态和$t$时刻相同，并且应该纳入统计而不是放弃结果！！
     • 由于状态出现概率和转移概率满足细致平稳条件，因此状态出现概率即稳态分布概率。$长期转移$后，我们会得到想要的抽样分布效果。

   • 图解转移过程：
     

     • 要注意的是，当马氏链处在状态$i$时，我们并不知道如何选择下一个状态$j$。
     • 我们在这里采用抽样的方法，并$借助接受率，“强迫”转移过程逼近理想的遍历马氏链转移过程$。

   • 算法描述如下：

     • 初始化马氏链状态$X_0=x_0$
     • 从$t=0$开始，重复以下步骤：
       • 假设该$t$时刻状态为$X_t=x_t$；
       • 对$Q_{x_t,x}$抽样，随机得到一个状态$x_{next}$；
       • 在$[0,1]$上均匀抽样，得到一个$u$；
       • 若$u<\alpha (x_t,x_{next})=P(x_{next})Q_{x_{next},x_t}$，则接受转移$X_{t+1}=x_{next}$；
       • 否则不转移，即$X_{t+1}=x_t$。
2. 算法升级
   1. 存在问题
      如果$\alpha (x_t,x_{next})$太小，会导致转移很难发生，马氏链收敛过慢。
      我们回到细致平稳条件：
      $$P(i)Q_{ij}\alpha (i,j)=P(j)Q_{ji}\alpha (j,i)\text{\hspace{1em}(13.3.10)}$$
      我们希望在保持条件成立的前提下，使接受率尽量增大。
   2. 改进
      由于接受率的本质是概率，因此$\alpha (i,j)$和$\alpha (j,i)$中的较大者不应大于$1$。那么我们作下述改进即可：
      $$\begin{gathered} buf=\frac{P(j){\mathbf{Q}}_{ji}}{P(i){\mathbf{Q}}_{ij}} \\ \alpha (i,j)=min(buf,1)\text{\hspace{1em}(13.3.11)} \end{gathered}$$
      解释：
      • 如果$\alpha (i,j)$更大，那么应设为1。而第一项大于1，因此$\alpha (i,j)=1$，结束。
      • 如果$\alpha (j,i)$更大，那么为了等式成立，$\alpha (i,j)$必须等于buf。而$buf<1$，
   3. $改进算法$：

      • 初始化马氏链状态$X_0=x_0$

      • 假设该$t$时刻状态为$X_t=x_t$；

        • 对$Q_{x_t,x}$抽样，随机得到一个状态$x_{next}$；
        • 在$[0,1]$上均匀抽样，得到一个$u$；
        • 若$u<\alpha (x_t,x_{next})=min(\frac{P(next){\mathbf{Q}}_{next,t}}{P(t){\mathbf{Q}}_{t,next}},1)$，则接受转移$X_{t+1}=x_{next}$；
        • 否则不转移，即$X_{t+1}=x_t$。

        这样执行了$N$次，就可以得到 { X 1 , X 2 , ⋯   , X N } \{ X_{1},X_{2},\cdots,X_{N} \} {X1​,X2​,⋯,XN​}个样本点。

        这里的$P(t)=\frac{\hat{P}(t)}{t_p}$。其中$t_p$指的是归一化因子，几乎非常难以计算，所以一般是未知的。而$\hat{t}(Z)=likelihood\times prior$。所以这里的$P(t)$和$P(next)$就是$\hat{P}(t)$和$\hat{P}(next)$。由于归一化因子被抵消了，干脆就直接写成了$P(t)$和$P(next)$。

---

参考：

1. MCMC-1|机器学习推导系列（十五）
2. Metropolis-Hastings algorithm