马尔可夫链蒙特卡罗算法（MCMC方法）

MCMC方法是什么

• 具体而言，假设我们要计算积分$$\mu ={\int }_{S}h(\theta )\pi (\theta |x)d\theta$$如果后验分布$\pi (\theta |x)$难以直接抽样，那么我们就可以构造一条马氏链$\theta$，使其状态空间为$S$且其平稳分布就是$\pi (\theta |x)$，将此马氏链运行一段时间使其收敛于平稳分布后，将产生一系列服从该平稳分布的随机样本${\theta }_0,{\theta }_1,{\theta }_2,...,{\theta }_{n-1}$，由该定理可得$$\mu \approx \frac{1}{n}\sum _{i=0}^{n-1}h({\theta }_i)$$这种以马尔可夫链生成随机样本，通过蒙特卡罗方法进行积分计算的方法就称为马尔可夫链蒙特卡罗算法，即$MCMC$方法

为什么需要MCMC方法

• 在之前的文章中有提到，接受-拒绝采样方法在高维随机变量的情形下很难找到合适的提议分布$q(x)$，而$MCMC$方法正是为了解决这个问题而生。$MCMC$方法通过马尔可夫链收敛后的平稳分布产生随机样本，随机变量的维度就是马氏链$X$的维度，因此无论随机变量的维数是多少，只要马尔可夫链能够收敛，$MCMC$方法都可以很好的工作。

MCMC方法中的若干术语

• 初始值： 初始值被用于初始化一条马氏链
• 预烧期： 在$MCMC$方法中用于保证马氏链收敛到平稳分布的运行时间，其迭代次数记为$B$，只要运行的时间足够长，去除预烧期对后验推断几乎是没有影响的
• 筛选间隔或抽样步长： 马氏链产生的样本并不是相互独立的，如果需要独立样本，可以通过监视产生样本的自相关图，然后选择筛选间隔或抽样步长$L$，使得步长$L$以后的样本自相关性很低，可近似为独立样本
• 迭代保持数： 指提供给贝叶斯分析的实际样本数，其计算公式为$$迭代保持数=迭代总次数-预烧期的迭代次数$$如果抽样步长$L>1$，则$$迭代保持数=\frac{迭代总次数-预烧期的迭代次数}{L}$$
• 蒙特卡罗误差： 度量了每一个估计因为随机模拟而导致的波动性，计算方法略

利用MCMC的输出结果描述目标后验分布

• MCMC输出结果给我们提供了一系列随机样本$${\theta }_1,{\theta }_2,{\theta }_3,...,{\theta }_n$$这组样本可以等价地看作从目标后验分布$\pi (\theta |x)$中生成的样本，因此，我们就可以通过这组样本计算后验分布的数字特征，包括但不限于以下所列：
  • 后验期望： $${\mu }^{\pi }(\theta )=\frac{{\sum }_{i=1}^n{\theta }_i}{n}$$
  • 后验方差： $$V^{\pi }(\theta )=\frac{{\sum }_{i=1}^n[{\theta }_i-{\mu }^{\pi }(\theta ){]}^2}{n}$$
  • 后验众数： 在生成的后验密度图中使得后验密度达到最大的点${\theta }_{mode}$(不等于 m a x { θ 1 , θ 2 , θ 3 , . . . , θ n } max\{\theta_1,\theta_2,\theta_3,...,\theta_n\} max{θ1​,θ2​,θ3​,...,θn​}！！)
• 进一步地，如果我们要应用$MCMC$方法的采样结果解决其它问题，那么原理也是相同，下面提供一个一般化的使用$MCMC$方法解决贝叶斯预测推断问题的例子：
  • 假设历史样本$X$服从一个参数为$\theta$的分布$f(x|\theta )$，$\theta$的先验分布为$\pi (\theta )$，现在要求随机变量$Z\sim g(z|\theta )$的后验预测密度
  • 按照贝叶斯预测推断的解题步骤，首先应该根据样本似然与先验分布求出后验分布，然后应用后验分布求出$Z$的后验预测密度。$\theta$的后验分布按照如下公式计算$$\pi (\theta |x)=\frac{f(x|\theta )\pi (\theta )}{{\int }_{\theta }f(x|\theta )\pi (\theta )d\theta }$$则$Z$的后验预测密度为$$m(z)={\int }_{\theta }g(z|\theta )\pi (\theta |x)d\theta$$假设这是一个难以计算的积分，那么应用$MCMC$方法，得到一组基于后验分布采样得到的随机样本${\theta }_1,{\theta }_2,{\theta }_3,...,{\theta }_n$，此时就有$$m(z)\approx \frac{{\sum }_{i=1}^{n}g(z|{\theta }_i)}{n}$$

马氏链收敛性诊断

• 样本路径图： 根据马氏链的迭代次数对生成的样本值作图，如果所有的值都在一个区域里，且没有明显的周期性和趋势性，则认为收敛性已经达到

MCMC采样方法

• 由于一般情况下，随机构造的一个马尔科夫链状态转移矩阵$Q$不满足细致平稳条件，即$${\pi }_i{Q}_{ij}\ne {\pi }_j{Q}_{ji}$$因此$MCMC$方法对上式做了一个改造，使得细致平稳条件恒成立，具体做法是引入了一个$\alpha$，即$${\pi }_i{Q}_{ij}\alpha (i,j)={\pi }_j{Q}_{ji}\alpha (j,i)$$
• 问题是什么样的$\alpha (i,j)$可以使等式成立呢？其实很简单，只要满足下两式即可：$$\{\begin{array}{c}\alpha (i,j)={\pi }_j{Q}_{ji}\\ \alpha (j,i)={\pi }_i{Q}_{ij}\end{array}$$这样我们就得到了转移矩阵$P$$$P_{ij}=Q_{ij}\alpha (i,j)$$也就是说，我们的目标矩阵$P$可以通过任意一个马尔科夫链状态转移矩阵$Q$乘以$\alpha (i,j)$得到。$\alpha (i,j)$我们一般称之为接受率，取值在$[0,1]$之间，可以理解为一个概率值，即目标矩阵$P$可以通过任意一个马尔科夫链状态转移矩阵$Q$以一定的接受率$\alpha$获得。
• 这个很像接受-拒绝采样，它是以一个提议分布的采样结果通过一定的接受-拒绝概率近似一个复杂分布的采样结果，这里是以一个常见的马尔科夫链状态转移矩阵$Q$通过一定的接受-拒绝概率得到目标转移矩阵$P$,两者的解决问题思路是类似的。
• 我们在MCMC采样方法中同样引入一个提议分布，用于近似转移矩阵的条件概率分布，即$$g(x^{\ast }|x)=Q(x^{\ast }|x)$$用于在马氏链的前一个状态为$x$的条件下，采样生成新的转移状态$x^{\ast }$的概率分布
• $MCMC$采样步骤：
  • 输入我们选定的提议分布$g(x|x_t)$，平稳分布$\pi (x)$，设定预烧期$B$，迭代保持数$T$
  • 从任意简单概率分布采样得到初始值$x_0$
  • for $t=0$ to $B+T-1$:
    1. 由提议分布$g(x|x_t)$中采样得到样本$x^{\ast }$
    2. 从均匀分布$U(0,1)$采样得到$u_0$
    3. 如果$u_0<\alpha (x_t,x^{\ast })=\pi (x^{\ast })g(x^{\ast }|x_t)$，则接受转移，即$x_{t+1}=x^{\ast }$，否则不接受转移，即$x_{t+1}=x_t$
  • 最终样本集 $(x_B,x_{B+1},...,x_{B+T-1})$便是我们所需要的从平稳分布中采集的样本。
  • 上面这个过程基本上就是$MCMC$采样的完整采样理论了，但是这个采样算法还是比较难在实际中应用，为什么呢？问题就出在上面的步骤3，由于$\alpha (x_t,x^{\ast })$可能非常的小，比如0.1，导致我们大部分的采样值都被拒绝转移，采样效率很低。有可能我们采样了上百万次马尔可夫链还没有收敛，怎么办呢？这时就轮到我们的$M-H$采样出场了。

Metropolis-Hastings算法

• 我们回到$MCMC$采样的细致平稳条件：$${\pi }_i{Q}_{ij}\alpha (i,j)={\pi }_j{Q}_{ji}\alpha (j,i)$$我们采样效率低的原因是$\alpha (i,j)$太小了，我们假设$\alpha (i,j)$为0.1，而$\alpha (j,i)$为0.2。即$${\pi }_i{Q}_{ij}\times 0.1={\pi }_j{Q}_{ji}\times 0.2$$这时我们可以看到，如果两边同时扩大五倍，接受率提高到了0.5，细致平稳条件却仍然是满足的，即$${\pi }_i{Q}_{ij}\times 0.5={\pi }_j{Q}_{ji}\times 1$$因此，我们可以通过类似的方式增大接受率，对接受率做如下改进，即 α ( i , j ) = m i n { π j Q j i π i Q i j , 1 } α(i,j)=min\{\frac{π_jQ_{ji}}{π_iQ_{ij}},1\} α(i,j)=min{πi​Qij​πj​Qji​​,1}
  通过这个微小的改造，我们就得到了可以在实际应用中使用的$Metropolis-Hastings$算法.
• 算法过程如下：
  • 输入我们选定的提议分布$g(x|x_t)$，平稳分布$\pi (x)$，设定预烧期$B$，迭代保持数$T$
  • 从任意简单概率分布采样得到初始值$x_0$
  • for $t=0$ to $B+T-1$:
    1. 由提议分布$g(x|x_t)$中采样得到样本$x^{\ast }$，每次迭代将会根据上一状态$x_t$对提议分布做出更新
    2. 从均匀分布$U(0,1)$采样得到$u_0$
    3. 如果$u_0<\alpha (x_t,x^{\ast })=\frac{\pi (x^{\ast })g(x_t|x^{\ast })}{\pi (x_t)g(x^{\ast }|x_t)}$（由于$u_0$本身就采样自$U(0,1)$，因此一定小于1），则接受转移，即$x_{t+1}=x^{\ast }$，否则不接受转移，即$x_{t+1}=x_t$
  • 最终样本集 $(x_B,x_{B+1},...,x_{B+T-1})$便是我们所需要的从平稳分布中采集的样本。

Metropolis抽样方法

• $Metropolis-Hastings$抽样方法是$Metropolis$抽样方法的推广，$Metropolis$抽样方法的不同之处主要还是在于接受率：
  • 在$Metropolis$抽样方法中，我们假设提议分布$g(x^{\ast }|x)$是对称的，即满足$g(x|x^{\ast })=g(x^{\ast }|x)$，这时我们的接受率可以进一步简化为 α ( i , j ) = m i n { π j π i , 1 } α(i,j)=min\{\frac{π_j}{π_i},1\} α(i,j)=min{πi​πj​​,1}

随机游动的Metropolis算法

• 随机游动的$Metropolis$算法又是$Metropolis$抽样方法的一个特例，随机游动的$Metropolis$算法的不同之处在于对马氏链生成新的状态的方式做出了改变
• 具体而言，随机游动的$Metropolis$算法的提议分布仍然满足对称性，即提议分布$g$满足$$g(x|x^{\ast })=g(x^{\ast }|x)$$但我们生成$x^{\ast }$的方式有所不同，我们引入了一个随机增量$Z$，$Z$从原始的提议分布$g(\cdot )$产生，然后令$x^{\ast }=x_t+Z$，这就要求提议分布满足$$g(x^{\ast }|x)=g(x^{\ast }-x|x)=g(Z)$$显然，提议分布$g$需要满足平移不变性
• 随机游动的$Metropolis$算法具体步骤如下：
  • 输入我们选定的提议分布$g(\cdot )$，平稳分布$\pi (x)$，设定预烧期$B$，迭代保持数$T$
  • 从任意简单概率分布采样得到初始值$x_0$
  • for $t=0$ to $B+T-1$:
    1. 由提议分布$g(\cdot )$中采样得到随机增量$Z$，令$x^{\ast }=x_t+Z$
    2. 从均匀分布$U(0,1)$采样得到$u_0$
    3. 如果$u_0<\alpha (x_t,x^{\ast })=\frac{\pi (x^{\ast })}{\pi (x_t)}$（由于$u_0$本身就采样自$U(0,1)$，因此一定小于1），则接受转移，即$x_{t+1}=x^{\ast }$，否则不接受转移，即$x_{t+1}=x_t$
  • 最终样本集 $(x_B,x_{B+1},...,x_{B+T-1})$便是我们所需要的从平稳分布中采集的样本。
• 随机游动的$Metropolis$算法对随机增量的方差比较敏感，当增量的方差太大时，大部分的候选点会被拒绝，而增量方差太小时，候选点就几乎都被接受。一种好的解决办法是监视拒绝概率，当拒绝概率在区间$[0.15,0.5]$之内才可以保证得到的马氏链具备较好的性质，下面是一个示例，假设随机增量$Z$从提议分布$N(0,{\sigma }^2)$中产生，当$\sigma$取值不同时，生成的马氏链如图
  

随机游动的Metropolis算法与Metropolis-Hastings算法在代码实现上的差异

• $Metropolis-Hastings$算法

  for i in range(1,m):
      xt=x[i-1]
      g=gamma(xt,1)
      y=g.rvs()
      num=f.pdf(y)*gamma(y,1).pdf(xt)
      den=f.pdf(xt)*gamma(xt,1).pdf(y)
      if u[i-1]<=num/den:
          x[i]=y
      else:
          x[i]=x[i-1]
          k=k+1
  

• 随机游动的Metropolis算法

  for i in range(1,m):
          xt=x[i-1]
          y=g.rvs()+xt
          num=f.pdf(y)
          den=f.pdf(xt)
          if u[i-1]<=num/den:
              x[i]=y
          else:
              x[i]=x[i-1]
              k+=1
  

Gibbs抽样方法

• $Gibbs$抽样是$Metropolis-Hastings$算法的另一种特殊情形，它特别适合于目标分布是多元的情形，它将从多元目标分布中抽样转化为了从一元目标分布抽样，固定n−1个特征在某一个特征采样，这便是$Gibbs$抽样的重要性所在

• 定义： 令$X=(X_1,X_2,X_3,...,X_k)$为$k$维随机变量，其联合分布$\pi (x)=\pi (x_1,x_2,x_3,...,x_k)$为目标抽样分布，定义$k-1$维随机变量$$X_{-j}=(X_1,...,X_{j-1},x_{j+1},...,X_k)$$并记$X_j|X_{-j}$的满条件分布的概率密度为$\pi (x)=\pi (x_j|x_{-j}),j=1,2,3,...,k$，则$Gibbs$抽样就是从这k个条件分布中产生候选点，以解决直接从$f$中进行抽样的困难，该算法步骤如下：

  • 输入平稳分布$\pi (x_1,x_2，...,x_n)$或者对应的所有特征的条件概率分布，设定预烧期为$B$，迭代保持数为$T$
  • 在$t=0$时，初始化$X_i^{(0)}(i=1,2,...,k)$（我们记$X_i^{(t)}$为$X_i$在第 $t$ 次迭代下的情形）
  • for $t=0$ to $B+T-1$:
    • 从条件概率分布$P(x_1|x_2^{(t)},x_3^{(t)},...,x_k^{(t)})$中采样得到样本$x_1^{t+1}$
    • 从条件概率分布$P(x_2|x_1^{(t+1)},x_3^{(t)},x_4^{(t)},...,x_k^{(t)})$中采样得到样本$x_2^{t+1}$
    • …
    • 从条件概率分布$P(x_k|x_1^{(t+1)},x_2^{(t+1)},...,x_{k-1}^{(t+1)})$中采样得到样本$x_k^{(t+1)}$
  • 样本集$\{(x_1^{(B)},x_2^{(B)},...,x_k^{(B)}),...,(x_1^{(B+T-1)},x_2^{(B+T-1)},...,x_k^{(B+T-1)})\}$即为我们需要的平稳分布对应的样本集。

• 一个示例：
  
  代码：

  import numpy as np 
  from scipy.stats import uniform,binom,beta
  import matplotlib.pyplot as plt 
  
  # 定义马氏链长度
  m=10000
  # 定义预烧期
  burn=1000
  # 初始化一条二维[xy]的马氏链
  x=np.zeros((m,2))
  n=50;a=1;b=1
  # 由于二项分布取值在0-1之间，故在U(0,1)上随机产生一个数值
  x[0,0]=uniform(0,1).rvs()
  # 由y|x的条件分布的期望值初始化x[0,1]
  x[0,1]=(x[0,0]+a)/(n+a+b)
  
  for i in range(1,m):
      x[i,]=x[i-1,]
      # print(x[i,] is x[i-1,])
      # 计算新的x
      x[i,0]=binom(n,x[i,1]).rvs()
      # 计算新的y
      x[i,1]=beta(x[i,0]+a,n-x[i,0]+b).rvs()
  
  x=x[burn:,:]
  

References

• https://www.cnblogs.com/pinard/p/6638955.html
• https://zhuanlan.zhihu.com/p/37121528