S06P06-马尔科夫链蒙特卡罗法

马尔科夫链蒙特卡罗法(Markov Chain Monte Carlo,MCMC)是以马尔科夫链为概率模型的蒙特卡罗法，因此首先要明确什么是蒙特卡罗法

蒙特卡罗法

蒙特卡罗法也称为统计模拟方法，是通过从概率模型的随机抽样进行近似数值计算的方法。对于机器学习而言，就是基于数据对概率分布的特征进行推断，蒙特卡罗法在假设概率分布定义已知的情况下，通过抽样获得概率分布的随机样本，并通过随机样本对概率分布的特征进行分析
例如对如下图的函数$f(x)$，要求其阴影部分面积即在区间$[a,b]$上的积分，最准确的方法当然是直接计算积分，但是有时候函数表达式很复杂或者很难写出函数表达式时，就可以利用蒙特卡罗法

最简单的，首先在区间$[a,b]$上随机抽样$n$个样本点$x_1,x_2,\cdots,x_n$，然后用其均值来代表$f(x)$在区间$[a,b]$上所有值，则可计算
$$
\int_anf(x_i)}{n}
$$
这就是蒙特卡罗法的基本思想。但是上述方法隐含了一个假设，即$x$在$[a,b]$区间是均匀分布的，而绝大部分情况，并不是符合均匀分布的。如果我们用上面的方法，则模拟求出的结果很可能和真实值相差甚远，下面讨论符合任意分布时的求解方法

数学期望估计

假设随机变量$x\in\mathcal{X}$，其概率密度函数为$p(x)$，$g(x)$为定义在$\mathcal{X}$上的函数，用蒙特卡罗法求期望$E{p(x)}[g(x)]$
首先按照概率密度$p(x)$独立地抽取$n$个样本$x_1,x_2,\cdots,x_n$，抽样以后函数$g(x)$在该样本上的均值为
$$
\overline{g}={1\over n}\sum{i=1}^ng(xi)
$$
当$n$足够大时期望可以由样本均值近似估计，因此
$$
E{p(x)}[g(x)]\approx{1\over n}\sum_{i=1}^ng(x_i)
$$

积分计算

下面来回到求$f(x)$的积分，我们可以将$f(x)$分解为一个函数$g(x)$和一个概率密度函数$p(x)$的乘积，那么有
$$
\int\mathcal{X}f(x)dx=\int\mathcal{X}g(x)p(x)dx=E{p(x)}[g(x)]
$$
显然令$g(x)=\cfrac{f(x)}{p(x)}$即可，于是
$$
\int\mathcal{X}f(x)dx\approx{1\over n}\sum_{i=1}^n\frac{f(x_i)}{p(x_i)}
$$

抽样方法

上面解决了计算的问题，但是还有重要的抽样问题没有解决，即如何抽取一个符合$p(x)$分布的随机样本，蒙特卡罗法的核心就是随机抽样。一般的抽样方法有直接抽样法、接受-拒绝抽样法、重要性抽样法等，接受-拒绝抽样法、重要性抽样法适用于概率密度函数复杂不能直接抽样的情况

直接抽样法

直接抽样法又称为概率分布抽样法，就是直接根据概率分布函数来抽样。例如对于正态分布而言，可以写出其概率分布函数，然后根据均匀分布$U(0,1)$随机抽样$n$个随机数，再求出各随机数对应的概率分布函数的反函数即得到了$n$个随机样本，这$n$个随机样本显然是符合正态分布的
但是实际上很多分布非常复杂，很难写出概率分布函数，因此也很难得到符合该分布的随机样本，因此还有接受-拒绝抽样法、重要性抽样法等，下面主要介绍接受-拒绝采样法

接受-拒绝采样法

既然太复杂的分布在程序中没法直接采样，那么我们设定一个程序可抽样的分布比如高斯分布，然后按照一定的方法拒绝某些样本，达到接近原始分布的目的，这就是接受拒绝采样
假设一个概率密度函数为$p(z)$的分布不可以直接抽样，可以找一个可以直接抽样的分布，称为建议分布(proposal)，设其概率密度函数为$q(z)$，且满足$kq(z)\ge q(z)$，$k$为一个正的常数，如下图所示

先按照$q(z)$进行抽样，例如抽取到样本$z_0$，然后在均匀分布$U(0,1)$上随机抽样得到$u_0$，如果$u_0\gt\cfrac{p(z_0)}{kq(z_0)}$即落在图中阴影部分，那么拒绝该样本，如果$u_0\le\cfrac{p(z_0)}{kq(z_0)}$即落在图中红色曲线以下，那么接受该样本作为最终的采样样本。如此重复该过程就得到了一个符合$p(z)$分布的随机样本
对于每一个$z_0$，定义
$$
\alpha=\frac{p(z_0)}{kq(z_0)}
$$
称为接受率

有效性证明

证明该采样方法有效即证明采样出的随机样本符合$p(z)$分布，用$P(z)$表示$p(z)$的累积概率分布函数，也即证下式
$$
Pr(Z\le z|Z\text{ is accepted})=P(Z\le z)=P(z)
$$
证明：
因为$Z$和$U$相互独立，联合概率密度即边缘概率密度乘积，则样本$Z$被接受的概率
$$
\begin{aligned}
Pr(Z\text{ is accepted})&=P\left(U\le\frac{p(Z)}{kq(Z)}\right)\
&=\intz\int_0^{\frac{p(z)}{kq(z)}}p_U(u)q(z)dudz\
&=\int_zq(z)\int_0^{\frac{p(z)}{kq(z)}}p_U(u)dudz\
&=\int_zq(z)\frac{p(z)}{kq(z)}dz\
&={1\over k}
\end{aligned}
$$
$Z$和$U$的联合概率
$$
\begin{aligned}
Pr(Z\le z,Z\text{ is accepted})&=P\left(Z\le z,U\le\frac{p(Z)}{kq(Z)}\right)\
&=\int{-\infty}{\frac{p(t)}{kq(t)}}pU(u)q(t)dudt\
&=\int{-\infty}^zq(t)\frac{p(t)}{kq(t)}dt\
&={1\over k}\int_{-\infty}^zp(t)dt\
&={1\over k}P(Z\le z)
\end{aligned}
$$
根据贝叶斯公式
$$
\begin{aligned}
Pr(Z\le z|Z\text{ is accepted})&=\frac{Pr(Z\le z,Z\text{ is accepted})}{Pr(Z\text{ is accepted})}\
&=P(Z\le z)
\end{aligned}
$$
得证

方法优缺点

接受-拒绝采样法的优点是很容易实现，缺点是效率可能不高。因为当$p(z)$相对于$kq(z)$覆盖的面积很小时，容易导致拒绝比例非常高，抽样效率低下。而且，在高维空间中，一是很难找到合适的$k$值，二是高维空间即使$p(z)$和$kq(z)$很接近两者涵盖的面积差异也可能很大，与在低维空间的直觉不同

马尔科夫链

基本定义

马尔科夫链定义

设一个随机变量序列$X=(X1,X_2,\cdots,X_t,\cdots)$，每个随机变量的取值集合相同即$X_t\in S$，$S$称为状态空间。随机变量可以是离散的也可以是连续的。假设某一时刻$t(t\ge2)$的随机变量$X_t$只依赖于它前一时刻的随机变量$X{t-1}$，而与$X{t-2}$及以前的随机变量无关，这个性质称为马尔科夫性，即
$$
P(X_t|X{t-1},X{t-2},\cdots,X_1)=P(X_t|X{t-1})
$$
$P(Xt|X{t-1})$一般称为状态转移概率分布，即在$t-1$时刻$X_{t-1}$的状态确定的情况下，$t$时刻转移为$X_t$状态的概率
满足马尔科夫性的随机序列$X$称为马尔科夫链(Markov Chain)或者马尔科夫过程(Markov Process)

齐次马尔科夫性

若转移概率分布与时刻$t$无关，即
$$
P(X{t+s}|X{t-1+s})=P(Xt|X{t-1})
$$
称为时间齐次马尔科夫性，即每个时刻的状态转移概率分布都是相同的。满足时间齐次马尔可夫性的马尔科夫链称为齐次马尔可夫链

n阶马尔科夫链

以上定义的是一阶马尔科夫链，可以推广的n阶马尔科夫链，即$t$时刻的随机变量$Xt$只与其前$n$个时刻有关
$$
P(X_t|X{t-1},\cdots,X1)=P(X_t|X{t-n},\cdots,X_{t-1})
$$

离散状态马尔科夫链

一般情况下我们说的马尔科夫链都是指离散状态马尔科夫链。离散状态马尔科夫链$X=(X0,X_1,X_2,\cdots,X_t,\cdots)$的状态空间$S$为离散空间，令$X_t\in S=(s_1,s_2,\cdots,s_m)$
若马尔科夫链在$t-1$时刻处于状态$s_j$，$t$时刻转移到状态$s_i$，将转移概率记作
$$
p{ij}=P(Xt=s_i|X{t-1}=sj)
$$
满足$p{ij}\ge0,\quad\sum\limits{i=1}^mp{ij}=1$
一般也简化记作
$$
p{ij}=P(X_t=i|X{t-1}=j)
$$
那么有状态转移概率矩阵
$$
P=
\begin{bmatrix}
p{11}&p{12}&\cdots&p{1m}\
p{21}&p{22}&\cdots&p{2m}\
\vdots&\vdots&\ddots&\vdots\
p{m1}&p{m2}&\cdots&p{mm}
\end{bmatrix}
$$
$P$是一个随机矩阵(stochastic matrix)
考虑$X$在时刻$t(t=0,1,\cdots)$的概率分布，称为时刻$t$的状态分布
$$
\pi_t=\begin{bmatrix}\pi_t(1)\pi_t(2)\vdots\pi_t(m)\end{bmatrix}
$$
其中$\pi_t(i)=P(X_t=i),\quad i=1,2,\cdots,m$
那么特别地，在初始状态$X_0$有
$$
\pi_0=\begin{bmatrix}\pi_0(1)\pi_0(2)\vdots\pi_0(m)\end{bmatrix}
$$
通常我们设初始状态$\pi_0$向量中只有一个分量为1，其余分量均为0
对于时刻$t$的分布$\pi_t$可以计算
$$
\begin{aligned}
\pi_t(i)&=P(X_t=i)\
&=\sum{j=1}^mP(Xt=i,X{t-1}=j)\
&=\sum{j=1}^mP(X_t=i|X{t-1}=j)P(X{t-1}=j)\
&=\sum{j=1}^mp{ij}\pi{t-1}(j)\
&=\begin{pmatrix}p{11}&p{12}&\cdots&p{1m}\end{pmatrix}
\begin{pmatrix}\pi{t-1}(1)\pi{t-1}(2)\cdots\pi{t-1}(m)\end{pmatrix}\
&=\begin{pmatrix}p{11}&p{12}&\cdots&p{1m}\end{pmatrix}\pi{t-1}
\end{aligned}
$$
于是可得如下递推式
$$
\begin{aligned}
\pit&=\begin{bmatrix}\pi_t(1)\pi_t(2)\vdots\pi_t(m)\end{bmatrix}=\begin{bmatrix}\begin{pmatrix}p{11}&p{12}&\cdots&p{1m}\end{pmatrix}\pi{t-1}\begin{pmatrix}p{21}&p{22}&\cdots&p{2m}\end{pmatrix}\pi{t-1}\vdots\begin{pmatrix}p{m1}&p{m2}&\cdots&p{mm}\end{pmatrix}\pi{t-1}\end{bmatrix}\
&=\begin{bmatrix}
p{11}&p{12}&\cdots&p{1m}\
p{21}&p{22}&\cdots&p{2m}\
\vdots&\vdots&\ddots&\vdots\
p{m1}&p{m2}&\cdots&p{mm}
\end{bmatrix}\pi{t-1}\
&=P\pi{t-1}
\end{aligned}
$$
于是可递推得到
$$
\pit=P^t\pi_0
$$
这里将$P^t$称为$t$步转移矩阵
$$
P{ij}^t=P(X_t=i|X_0=j)
$$
表示初始时刻状态为$s_j$，在$t$时刻转移为$s_i$的$t$步转移概率，显然马尔科夫链的状态分布由初始分布和概率转移分布决定

马尔科夫链的表示

马尔科夫链可用状态转移概率图来表示，节点表示状态，边表示状态之间的转移，边的数值即转移概率。例如假设天气符合马尔科夫性，如果今天是晴天，那么明天是晴天的概率是0.9，是雨天的概率是0.1；如果今天是雨天，那么明天是晴天的概率是0.5，是雨天的概率是0.5，那么可用如下状态转移概率图来表示

平稳分布

如果存在状态空间上的一个分布
$$
\pi=\begin{bmatrix}\pi(1)\pi(2)\vdots\pi(m)\end{bmatrix}
$$
使得
$$
\pi=P\pi
$$
则称$\pi$为马尔科夫链的平稳分布。显然如果平稳分布存在，那么以平稳分布为初始分布，之后任意时刻的状态分布都是平稳分布

平稳分布存在的充要条件

给定随机变量序列$X=(X0,X_1,X_2,\cdots,X_t,\cdots)$，状态空间为$S=(s_1,s_2,\cdots,s_m)$，转移矩阵$P=(p{ij})$，则分布$\pi=(\pi(1),\pi(2),\cdots,\pi(m))^T$为$X$的平稳分布的充要条件是$\pi$是下列方程组的解
$$
\begin{aligned}
&xi=\sum{j=1}^mp{ij}x_j,\quad i=1,2,\cdots,m\
&x_i\ge0,\quad i=1,2,\cdots,m\
&\sum{i=1}^mx_i=1
\end{aligned}
$$
一个马尔科夫链的平稳分布可能存在也可能不存在，可能只有一个也可能有多个

连续状态马尔科夫链

当状态空间$S$为连续时，就称为连续状态的马尔科夫链，与离散状态马尔科夫链类似，定义如下概率密度
$$
p(Xt=x|X{t-1}=x)
$$
称为转移核，即从状态$x^$转移到状态$x$概率密度，满足$\int{x\in S}k(x,x^)dx=1$
令$\pit(x)$为时刻$t$的状态分布，可以计算
$$
\begin{aligned}
\pi_t(x)&=\int{S}p(x,x\
&=\int{S}p(x)dx^\
&=\int{S}\pi{t-1}(x)dx^
\end{aligned}
$$
如果存在一个分布$\pi(x)$，使得
$$
\pi(x)=\int{S}\pi(x_)dx^
$$
那么称$\pi(x)$为该马尔科夫链的平稳分布

马尔科夫链的性质

下面均以离散状态马尔科夫链为例叙述，都可以自然推广到连续状态马尔科夫链
设马尔科夫链$X=(X_0,X_1,\cdots,X_t,\cdots)$，状态空间为$S$

1. 不可约
   对于任意状态$si,s_j\in S$，如果对于任意时刻$t\gt0$都满足
   $$
   P(X_t=i|X{t-1}=j)\gt0
   $$
   直观上讲就是指从任意某个状态都有可能转移到任意一个状态，则称此马尔科夫链是不可约的，否则就是可约的
2. 非周期
   对于任意状态$s_i\in S$，如果$t=0$时从$X_0=s_i$出发，在$t$时刻也到达状态$s_i$的所有时长${t:P(X_t=i|X_0=i)\gt0}$的最大公约数是1，则称马尔科夫链是非周期的(aperiodic)，否则是周期的(periodic)
   直观上讲就是指不会周期性的到达任意的某个状态
3. 正常返
   对于任意状态$si,s_j\in S$，定义概率$p{ij}^t$为从状态$X0=s_j$出发，在时刻$t$首次到达$X_t=s_i$的概率，即
   $$
   p{ij}^t=P(Xt=i,X_s\ne i,s=1,2,\cdots,t-1|X_0=j),\quad t=1,2,\cdots
   $$
   若$\lim\limits{t\to\infty}p_{ij}^t\gt0$，则称马尔科夫链是正常返的
   直观上讲就是说对于任意一个状态，就算前面所有的时刻都没有到达过该状态。在时间趋于无穷时，都有可能首次到达这个状态
4. 可逆马尔科夫链
   如果存在状态分布$\pi=(\pi(1),\pi(2),\cdots,\pi(i),\cdots)^T$，对于任意的状态$si,s_j\in S$，对任意时刻$t$都满足
   $$
   P(X_t=i|X{t-1}=j)\pi(j)=P(X{t-1}=j|X_t=i)\pi(i),\quad i,j=1,2,\cdots
   $$
   因为状态转移矩阵不变，那么即
   $$
   p{ij}\pi(j)=p{ji}\pi(i)
   $$
   则称此马尔科夫链为可逆马尔科夫链，上式称为细致平衡方程。直观上讲就是指可逆马尔可夫链在时间上无论是前向还是后向，概率分布是相同的，不仅可以用过去预测未来，也可以在未来推测过去
   根据细致平衡方程作如下计算
   $$
   \sum_jp{ij}\pi(j)=\sumjp{ji}\pi(i)=\pi(i)\sumjp{ji}=\pi(i)
   $$
   写成矩阵形式即
   $$
   P\pi=\pi
   $$
   很显然满足细致平衡方程的分布就是平稳分布

遍历定理

设马尔科夫链$X=(X0,X_1,\cdots,X_t,\cdots)$，状态空间为$S=(s_1,s_2,\cdots,s_i,\cdots)$。若该马尔科夫链是不可约、非周期且正常返的，则该马尔科夫链有唯一平稳分布$\pi=(\pi(1),\pi(2),\cdots,\pi(i),\cdots)^T$，而且转移概率的极限分布是马尔科夫链的平稳分布
$$
\lim{t\to\infty}P(Xt=i|X_0=j)=\pi(i),\quad i=1,2,\cdots;j=1,2,\cdots
$$
若$f(x)$是定义在状态空间的函数且满足$E\pi[|f(X)|]\lt\infty$，令
$$
\hat{f}t={1\over t}\sum{s=1}^tf(Xs)
$$
则有
$$
P(\hat{f}_t\to E\pi[f(X)])=1
$$
其中，$E\pi[|f(X)|]=\sum\limits{i}f(si)\pi(i)$为$f(X)$关于平稳分布$\pi$的期望
该定理的直观解释为：满足相应条件的马尔科夫链，无论其初始状态值为多少，在时间趋于无穷时，马尔科夫链的状态分布趋近于平稳分布，随机变量的函数的样本关于时刻的均值趋近于该函数的期望
理论上不知道迭代多长时刻状态分布才能接近平稳分布，在实际使用时一般取一个足够大的迭代次数$T$，认为迭代$T$次之后状态分布为平稳分布。这时候计算第$T+1$次迭代到第$N$次迭代的均值
$$
\hat{E}f={1\over{N-T}}\sum{t=T+1}^Nf(X_t)
$$
称为遍历均值

马尔科夫链蒙特卡罗法

在对一个目标概率分布做随机采样时，前面说过传统的蒙特卡罗法。现在介绍马尔科夫链蒙特卡罗法，这种方法更适用于随机变量多元、密度函数是非标准形式、随机变量各部分不独立等情况
例如下面这个问题：多元随机变量$x\in\mathcal{X}$，其概率密度函数为$p(x)$，$f(x)$为定义在$\mathcal{X}$上的函数，目标是获得符合$p(x)$分布的样本以及求$f(x)$的期望$E_{p(x)}[f(x)]$
如果应用马尔科夫链蒙特卡罗法解决这个问题，就是在状态空间$S$上定义一个满足遍历定理的马尔科夫链$X=(X_1,X_2,\cdots,X_t,\cdots)$，使其平稳分布为$p(x)$。那么根据上述遍历定理可得，在经过足够长时间$T$后，样本分布趋近平稳分布，样本函数的期望趋近遍历均值。于是取$T$时刻之后的样本即符合概率分布$p(x)$的样本
要注意的是马尔科夫链蒙特卡罗法得到的样本，相邻样本点是相关的。因此如果要得到相互独立的样本，可以在样本序列中再进行随机抽样，比如每隔一段时间取一次样本等方式

基本步骤

1. 在状态空间上构造一个满足遍历定理的马尔科夫链，使其平稳分布为$p(x)$
2. 从状态空间某一点出发，用构造的马尔科夫链生成样本序列$(x_0,x_1,\cdots,x_t,\cdots)$
3. 应用遍历定理，确定正整数$T,N(T\lt N)$，得到样本序列$(xT,\cdots,x_N)$，求得函数的期望
   $$
   \hat{E}f={1\over{N-T}}\sum{t=T+1}^Nf(x_t)
   $$

应用场景

马尔科夫链蒙特卡罗法在统计学习，尤其是贝叶斯学习中起着重要作用。它可以用在概率模型的学习和推理上
在贝叶斯学习中经常会有三种积分运算：规范化(normalization)、边缘化(marginalization)、数学期望(expectation)
设随机变量$y\in\mathcal{Y},x\in\mathcal{X},z\in\mathcal{Z}$
后验概率计算中需要规范化
$$
\int\mathcal{X}p(y|x)p(x)dx
$$
如果含有隐变量$z$，则需要边缘化
$$
p(x|y)=\int\mathcal{Z}p(x,z|y)dz
$$
计算函数关于后验概率分布的期望
$$
E{p(x|y)}[f(x)]=\int\mathcal{X}f(x)p(x|y)dx
$$
当观测数据和模型很复杂时，以上积分计算都很困难，马尔科夫链蒙特卡罗法就是一种解决方案

具体方法

根据上述讨论，可以看出马尔科夫链蒙特卡罗法的关键就在于得到一个平稳分布为目标概率分布的马尔科夫链。因此下面就以离散状态马尔科夫链介绍构造指定马尔科夫链的方法
假设状态空间为$S=(s_1,s_2,\cdots,s_i,\cdots)$，我们要抽样的目标分布为$\pi=(\pi(1),\pi(2),\cdots,\pi(i),\cdots)^T$

MCMC采样

考虑某个不可约马氏链的状态转移矩阵$Q$，其元素为$P(i|j)=q{ij}$，这个分布是容易采样的。在一般情况下，它是不满足细致平稳条件的，即
$$
\pi(i)q{ji}\ne\pi(j)q{ij}
$$
为了构造一个满足细致平稳条件的马氏链，可以引入一个$\alpha$使得下式成立
$$
\pi(i)q{ji}\alpha(i,j)=\pi(j)q{ij}\alpha(j,i)
$$
很显然，可以令
$$
\begin{gathered}
\alpha(i,j)=\pi(j)q{ij}\
\alpha(j,i)=\pi(i)q{ji}
\end{gathered}
$$
这样该马氏链就满足了细致平稳条件，设其状态转移矩阵为$P$，显然其元素为$p{ij}=q{ij}\alpha(j,i)$，平稳分布为目标分布$\pi$
也就是说，我们的目标矩阵$P$可以通过任意一个马尔科夫链状态转移矩阵$Q$乘以$\alpha(i,j)$得到。$\alpha(i,j)$我们一般称之为接受率，可以理解为一个概率值，即在原来的马氏链上以$q{ji}$的概率从状态$i$转移到状态$j$，然后我们以$\alpha(i,j)$的概率接受这个状态转移。这个很像接受-拒绝采样，接受-拒绝采样是以一个常用分布通过一定的接受-拒绝概率得到一个非常见分布，这里是以一个常见的马尔科夫链状态转移矩阵$Q$通过一定的接受-拒绝概率得到目标转移矩阵$P$，两者的解决问题思路是类似的
这个方法拓展到连续状态马尔科夫链同样有效，下面以连续的形式给出采样算法步骤

(1)初始化$X_0=x_0$
(2)从$t=0$开始重复以下步骤：

1. 假设$t$时刻状态为$X_t=x_t$；
2. 按建议状态转移核$k(x,x_t)$进行抽样，即随机得到一个状态$x^*$；
3. 从均匀分布$U(0,1)$上抽样得到$u$；
4. 若$u\le\alpha(xt,x)k(xt,x$，否则不接受转移$X_{t+1}=x_t$

(3)如此进行下去$t=0,1,\cdots$，在一定的迭代次数$T$后认为马尔科夫链收敛，得到符合目标分布$\pi$的样本$X{T+1},X{T+2},\cdots$

这个方法拓展到连续状态马尔科夫链同样是有效的。这个方法构造的马尔科夫链足够简单，容易实现，但是存在问题。在进行状态转移时，如果接受率$\alpha(x^*,x_t)$较小，那么容易造成大量的拒绝转移，于是遍历完马尔科夫链的所有状态要很大的迭代次数，收敛速度很慢，这样采样效率十分低下

M-H采样

M-H采样全称Metropolis-Hastings采样，它解决了普通MCMC采样在接受率低时采样效率低的问题
回到细致平稳条件
$$
\pi(i)q{ji}\alpha(i,j)=\pi(j)q{ij}\alpha(j,i)
$$
如果$\alpha(i,j),\alpha(j,i)$很小，假设$\alpha(i,j)=0.05,\alpha(j,i)=0.1$，即
$$
\pi(i)q{ji}\times0.05=\pi(j)q{ij}\times0.1
$$
显然，两边同时增加一定的倍数细致平稳条件仍然是满足的。比如两边同时扩大十倍
$$
\pi(i)q{ji}\times0.5=\pi(j)q{ij}\times1
$$
因此我们就可对接受率进行一些改造，将细致平稳条件$\alpha(i,j),\alpha(j,i)$两个数中的一个扩大为1，那么可有效避免接受率过小的情况。对细致平稳条件做个变换有
$$
\frac{\alpha(j,i)}{\alpha(i,j)}=\frac{\pi(i)q{ji}}{\pi(j)q{ij}}
$$
若$\cfrac{\pi(i)q{ji}}{\pi(j)q{ij}}\ge1$，则$\alpha(j,i)\ge\alpha(i,j)$，于是就将$\alpha(j,i)$扩大为1；若$\cfrac{\pi(i)q{ji}}{\pi(j)q{ij}}\lt1$，则将$\alpha(i,j)$扩大为1，于是有
$$
\alpha(j,i)=\min\left{1,\frac{\pi(i)q{ji}}{\pi(j)q{ij}}\right}
$$
这个方法拓展到连续状态马尔科夫链同样有效，下面以连续的形式给出M-H采样算法步骤

(1)初始化$X_0=x_0$
(2)从$t=0$开始重复以下步骤：

1. 假设$t$时刻状态为$X_t=x_t$；
2. 按建议状态转移核$k(x,xt)$进行抽样，即随机得到一个状态$x^*$，并计算接受率
   $$
   \alpha(x_t,x)k(xt,x,x_t)}\right}
   $$
3. 从均匀分布$U(0,1)$上抽样得到$u$；
4. 若$u\le\alpha(x*$，否则不接受转移$X_{t+1}=x_t$

(3)如此进行下去$t=0,1,\cdots$，在一定的迭代次数$T$后认为马尔科夫链收敛，得到符合目标分布$\pi$的样本$X{T+1},X{T+2},\cdots$

但是这个方法也有以下缺点：

1. 数据特征非常多，即高维情况下，M-H采样的接受率计算式计算会很耗时，因此会造成算法效率很低。而且因为还是存在拒绝转移的情况，有时候计算出了接受率又被拒绝，进一步导致效率变低
2. 在高维情况下，很难求出目标的各特征维度联合分布

Gibbs采样

在大数据时代，Gibbs采样解决了M-H算法的两个问题：第一是在采样过程中不存在拒绝，第二是不直接求联合概率分布，而是求各维度之间的条件概率分布
下面以连续状态马尔科夫链推导。先看二维的简单情况，即马尔科夫链的随机变量序列都是二维随机变量$X_t=(X_t{(2)}){(1)},X^{(2)})$。
考虑第一个维度取值相同的两个样本点$A(x{(2)}),B(x{(2)})$，有下式成立
$$
\begin{gathered}
p(x{(2)})p(y{(1)})=p(x{(2)}|x{(2)}|x^{(1)})\
p(x{(2)})p(x{(1)})=p(x{(2)}|x{(2)}|x^{(1)})\
\Rightarrow \pi(x{(2)})p(y{(1)})=\pi(x{(2)})p(x{(1)})
\end{gathered}
$$
也就是
$$
\pi(A)p(y{(1)})=\pi(B)p(x{(1)})
$$
同理取第二个维度相同的两个样本点$A(x{(2)}),C(y{(2)})$也成立
$$
\pi(A)p(y{(2)})=\pi(C)p(x{(2)})
$$
从这两个式子显然可以看出来，在$X{(1)}$这个维度相等时，用条件概率分布$p(X{(1)})$作为转移核满足细致平稳条件；在$X{(2)}$这个维度相等时，用条件概率分布$p(X{(2)})$作为转移核满足细致平稳条件
基于以上发现，我们可以构造一个马氏链，其转移核为
$$
k(y,x)=
\begin{cases}
p(y{(1)})&if&x{(1)}\
p(y{(2)})&if&x{(2)}\
0&else
\end{cases}
$$
根据上述推导显然这个转移核是满足细致平稳条件的，其平稳分布为目标分布$\pi(X)=\pi(X{(2)})$
于是我们可以利用不同维度之间的条件概率分布作为转移核构造马尔科夫链，这种方式称为Gibbs采样。一般情况下这个不同维度之间的条件概率分布是容易采样的，因此在状态转移时进行采样比较简单，从而提升采样效率，且不会出现拒绝样本
总结得到二维Gibbs采样步骤

(1)输入目标分布$\pi(X)=\pi(X{(2)})$，设定状态转移次数$T_1$，需要的样本个数$T_2$
(2)初始化$x_0=(x_0{(2)})^T$
(3)对于$t$从0到$T_1+T_2-1$进行如下步骤：

1. 从条件分布$p(X{(1)})$进行采样，得到第二个维度的样本$x_{t+1}^{(2)}$
2. 从条件分布$p(X{(2)})$进行采样，得到第一个维度的样本$x{t+1}^{(1)}$
   采样得到样本$x{t+1}=(x_{t+1}{(2)})$

最终${(x{T_1}),(x{T1+1}),\cdots,(x{T_1+T_2-1})}$即为符合目标分布$\pi(X)$的样本

从采样步骤可以看到在采样过程中，值不变的维度在不停的轮换，如此交替进行采样
以二维高斯分布为例，其Gibbs采样步骤反映在图像上即

按照同样的道理可以将上述方法推广到高维情况。假设目标分布为$k$维分布$\pi(X)=\pi(X{(j)},\cdots,X^{(k)})$，此时马氏链的转移核可以设为

1. 若当前状态为$x=(x{(2)},\cdots,xT$，在马氏链状态转移的过程中沿着某一个坐标进行转移，当沿着坐标$x{(j)}|x{(j-1)},x{(k)})$
2. 其他不沿某单个坐标转移的转移概率均设为0

于是可以总结出高维Gibbs采样的步骤如下：

(1)输入目标分布$\pi(X)=\pi(X{(j)},\cdots,X^{(k)})$，设定状态转移次数$T_1$，需要的样本个数$T_2$
(2)初始化$x_0=(x_0{(2)},\cdots,x_0T$
(3)对于$t$从0到$T_1+T_2-1$进行如下步骤：

1. 从条件分布$p(X{(2)},xt{(k)})$进行采样，得到样本$x{t+1}^{(1)}$
2. 从条件分布$p(X{(1)},xt{(k)})$进行采样，得到样本$x{t+1}^{(2)}$
3. ……
4. 从条件分布$p(X{(1)},x{t+1}{(j-1)},x_t{(k)})$进行采样，得到样本$x{t+1}^{(j)}$
5. ……
6. 从条件分布$p(X{(1)},x{t+1}{(k-1)})$进行采样，得到样本$x{t+1}^{(k)}$
   得到样本$x{t+1}=(x{t+1}{(2)},\cdots,x_{t+1}T$

最终${(x{T_1}),(x{T1+1}),\cdots,(x{T_1+T_2-1})}$即为符合目标分布$\pi(X)$的样本

整个过程非常类似坐标下降法。实际上轮换坐标轴不是必须的，我们可以随机选择某一个坐标轴进行状态转移，但是常用的Gibbs采样的实现都是基于坐标轴轮换的，实现更简单
从M-H角度来看Gibbs算法，Gibbs算法相当于每一次在某单个维度进行状态跳转时，其接受率都为1，也就是不存在拒绝采样的部分，这也是Gibbs采样和M-H采样的区别之一