随机模拟:马尔科夫链蒙特卡洛采样MCMC与EM算法「3」

本文是此系列的最后一篇拉！

最近学习了机器学习中的马尔科夫链蒙特卡洛(Markov Chain Monte Carlo, 简称MCMC) 相关的知识。

主要内容包括：

【1】蒙特卡洛原则，及其应用于采样的必要性（已经发布在头条）

【2】用于求解最大似然、近似推断、期望问题的经典采样算法：Metropolis-Hastings，Rejection，Importan，Metropolis和Gibbs算法。（本文属于此部分）

【3】马尔可夫链各个性质在蒙特卡洛采样问题中的应用，包括同质性，平移不变性

—————【3】—————

上一篇【2.3】中讨论了EM优化算法和采样的关系，介绍了MCMC类采样算法的基本思路，即不断改变假设分布来逼近目标分布，介绍了MCMC类算法中的Metropolis Hasting算法执行。

本篇将对该算法的有效性进行证明，并讨论MCMC中的其他算法，如Gibbs采样算法。在MHasting算法证明之前，先复习马尔科夫链的概念：

【马尔科夫链性质】

随机过程Xi，若满足下图中关系，则为马尔科夫链。此关系为，i时刻x的状态z仅与i-1时刻的状态相关，与再之前的状态构成基于条件i的条件独立性。

因此，这个T转移概率矩阵每一行的和都是1，可将链上任意一点隔断，其左侧与右侧点独立。

1、同质性 Homogeneous

如果链上每一个P(xi|xi-1)，对所有i都相同，则具有同质性。

2、平稳性和极限性 Stationary and limiting distributions

平稳性：p(x)是T上的平稳分布，不随i的变化改变

极限性：不管一开始p(x)是怎样的分布，最终收敛到一个稳定的分布，记为Π。（不自循环构成连结的马氏链都能收敛）

3、不可减性 Irreduciblity

对x的任何时刻的任何状态，都有一个大于0的概率转移到任何其他状态去。

4、非周期性 Aperiodicity

在链上如果出现现象：最少每隔d个状态，都存在Pii(n)>0,d>1则存在周期性。非周期性对应d=1。

5、遍历性 Ergodic

同时具有不可减性和非周期性的马氏链，称为具有遍历性。暗示我们可以得到平稳的极限p(x)。所有优秀的MCMC采样算法都必须满足遍历性，从而保证不管我们怎么取初始值，最终都能获得近似服从目标分布的采样。

以上是马氏链性质，接下来证明MHasting采样有效性。为证明，需先定义性质：

【细节平衡（可逆性）】

概率向量 pi=p(x)满足可逆性,a b 为不同x：pi(a)*Tab = pi(b)* Tba。称一个马氏链可逆若满足此性质。

可逆性可推出平稳性pi(a)=pi(b)。

接下来证明上一篇文章【2.3】中推导的MetropolisHasting采样算法的有效性。

【有效性证明】

该算法有效的含义为其所构成的马氏链具有可逆性。为方便查看，贴一下MH算法：

1、初始化x0

2、循环N-1次：对i=0，1，一直到N-1

从【0，1】均匀分布中随机采样小数u；

从假设分布q(x'|xi)中随机采样x'；

判断（记为A）：若u< A(x',x)=min {1, p_(x')*q(xi|x')/(p_(xi)*q(x'|xi))},x(i+1)=x'接受x'为新样本，

否则x(i+1)=x(i)。

以上为算法本身。从q中采样，再以A判断，则转移核为 T(x'|x)= q(x'|x) * A(x'|x),

观察A，不难发现若A(x',x)<=1,则A(x,x')=1,现在假设A(x',x)<1,则展开变形得到p_(x)T(x'|x)=p_(x')T(x|x'),得证可逆性。证明了MHasting算法推出稳定分布p_(x),而p_(x)被定义为未正态化的x的真实分布。

以上证明了MHasting算法有效性，下面给出两个MH算法的特例采样算法。

【M算法 Metropolis Algorithm】

M算法是MHasting算法的特例，其假设分布是随机游走的，有q(x|x')=q(x'|x),例如各向同性的高斯分布,因此在M算法中A=min{1，p_(x')/p_(x)}。A即接受概率。其他和MH完全相同。

【吉布斯采样算法 Gibbs Sampling】

假设N维度特征向量X，p(X)=p(x1,x2,-,xN),可以表示为对任意维度，可以从p(xj|x1,-,xj-1,xj+1,-,xN)获取样本。假设分布设置为下式：

xi为上一个样本，每个新样本的假设分布，为新样本状态j和前一个样本非j状态的条件概率，所有状态一起组成新样本的假设分布。

Gibbs采样算法可以被视为接受概率为1的MHasitng算法：

1、初始化X0={x10,x20,...,xn0}

2、假设已知一个样本Xi={x1i,x2i,...,xni}，对 Xi+1，分别对Xi+1的每一维度进行采样。即xi+1_j 从 p(xj|xi+1_1,...xi+1_j-1,xi_j+1,...xi_n0)。

执行这两步，就得到一个新样本，重复M次得到M个样本。

联系马尔科夫随机场的知识，可知这个状态间的条件概率也就是马尔科夫链垫的概念。在贝叶斯网络中，有向图，是节点xj的父节点，子节点和配偶节点。在无向图马尔可夫随机场中，是对所有直接邻点的条件概率。

在具体计算中，使用贝叶斯法则求所有后验概率，然后正态化是常用的计算方法。这样就能够得到某个状态的概率，和随机均匀抽样的u服从【0，1】相比，确定样本某维度的取值（u<P(j=T),则j=T，假设状态只取t/f）。若状态可能3种，将【0，1】按值三分，同样执行。

上图展示了j~{1，2}的采样过程，每次新维度的采样基于另一个已经确定的维度。

终于写完马尔科夫链上的蒙特卡洛采样算法了！EM算法，重要性采样、拒绝采样算法之前已经发布，如果需要请移步我的往期文章查看：）