|
关于马尔科夫蒙特卡罗模型，真的是重复看过很多次，始终不得其解，这次决心系统地学习一下，下面是我的理解过程，和大家分享一下，如果哪里讲得不对，欢迎评论指正。毕竟，MCMC模型，内容真的非常多，以它为名字的书籍，就有了以下两本。 马尔科夫链概述 马尔科夫链定义本身比较简单，它假设某一时刻状态转移的概率只依赖于它的前一个状态。 举个形象的比喻，假如每天的天气是一个状态的话，那个今天是不是晴天只依赖于昨天的天气，而和前天的天气没有任何关系。当然这么说可能有些武断，但是这样做可以大大简化模型的复杂度，因此马尔科夫链在很多时间序列模型中得到广泛的应用，比如循环神经网络RNN，隐式马尔科夫模型HMM等，当然MCMC也需要它。 如果用精确的数学定义来描述，则假设我们的序列状态是..., Xt−2, Xt−1, Xt, Xt+1, Xt+2, ...，那么我们的在时刻Xt+1的状态的条件概率仅仅依赖于时刻Xt，即： 既然某一时刻状态转移的概率只依赖于它的前一个状态，那么我们只要能求出系统中任意两个状态之间的转换概率，这个马尔科夫链的模型就定了。我们来看看下图这个马尔科夫链模型的具体的例子(来源于维基百科)。 这个马尔科夫链是表示股市模型的，共有三种状态：牛市（Bull market）, 熊市（Bear market）和横盘（Stagnant market）。每一个状态都以一定的概率转化到下一个状态。比如，牛市以 0.025 的概率转化到横盘的状态。这个状态概率转化图可以以矩阵的形式表示。如果我们定义矩阵阵P某一位置 P(i,j) 的值为 P(j|i)，即从状态 i 转化到状态 j 的概率，并定义牛市为状态 0，熊市为状态 1，横盘为状态2。这样我们得到了马尔科夫链模型的状态转移矩阵为： 讲了这么多，那么马尔科夫链模型的状态转移矩阵和我们蒙特卡罗方法需要的概率分布样本集有什么关系呢？这需要从马尔科夫链模型的状态转移矩阵的性质讲起。 马尔科夫链模型状态转移矩阵的性质 使用上面的状态转移矩阵为例。 假设我们当前股市的概率分布为：[0.3,0.4,0.3]，即30%概率的牛市，40%概率的熊盘与30%的横盘。然后这个状态作为序列概率分布的初始状态t0，将其带入这个状态转移矩阵计算 t1,t2,t3... 的状态。代码如下： …… …… 可以发现，从第 53 轮开始，我们的状态概率分布就不变了，一直保持在[0.625 0.3125 0.0625]，即62.5%的牛市，31.25%的熊市与6.25%的横盘。那么这个是巧合吗？ 我们现在换一个初始概率分布试一试，现在我们用 [0.7,0.1,0.2] 作为初始概率分布，然后再次执行程序。 可以看出，尽管这次我们采用了不同初始概率分布，但是，最终状态的概率分布还是会趋于同一个稳定的概率分布[0.625 0.3125 0.0625]。 也就是说我们的马尔科夫链模型的状态转移矩阵收敛到的稳定概率分布与我们的初始状态概率分布无关。这是一个非常好的性质，也就是说，如果我们得到了这个稳定概率分布对应的马尔科夫链模型的状态转移矩阵，则我们可以用任意的概率分布样本开始，带入马尔科夫链模型的状态转移矩阵，这样经过一些序列的转换，最终就可以得到符合对应稳定概率分布的样本。 这个性质不光对我们上面的状态转移矩阵有效，对于绝大多数的其他的马尔科夫链模型的状态转移矩阵也有效。同时不光是离散状态，连续状态时也成立。 同时，对于一个确定的状态转移矩阵P，它的 n 次幂 Pn 在当 n 大于一定的值的时候也可以发现是确定的，我们还是以上面的例子为例，计算代码如下： 我们可以发现，在 n≥6 以后，Pn 的值稳定不再变化，而且每一行都为[0.625 0.3125 0.0625]，这和我们前面的稳定分布是一致的。这个性质同样不光是离散状态，连续状态时也成立。 下面，我们用数学语言总结一下马尔科夫链的收敛性质： 如果一个非周期的马尔科夫链有状态转移矩阵P, 并且它的任何两个状态是连通的，那么有： 1) 2)
3)
4) π 是方程 πP=π 的唯一非负解，其中：
上面的性质中需要解释的有： 1) 非周期的马尔科夫链：这个主要是指马尔科夫链的状态转化不是循环的，如果是循环的则永远不会收敛。幸运的是我们遇到的马尔科夫链一般都是非周期性的。 2) 任何两个状态是连通的：这个指的是从任意一个状态可以通过有限步到达其他的任意一个状态，不会出现条件概率一直为 0 导致不可达的情况。 3) 马尔科夫链的状态数可以是有限的，也可以是无限的。因此可以用于连续概率分布和离散概率分布。 4) π通常称为马尔科夫链的平稳分布。 基于马尔科夫链采样 如果我们得到了某个平稳分布所对应的马尔科夫链状态转移矩阵，我们就很容易采用出这个平稳分布的样本集。 假设我们任意初始的概率分布是 π0(x)，经过第一轮马尔科夫链状态转移后的概率分布是 π1(x)，...，第 i 轮的概率分布是 πi(x)。假设经过n轮后马尔科夫链收敛到我们的平稳分布π(x)，即：
对于每个分布πi(x)，我们有：
现在我们可以开始采样了，首先，基于初始任意简单概率分布比如高斯分布 π0(x) 采样得到状态值 x0，基于条件概率分布 P(x|x0) 采样状态值 x1，一直进行下去，当状态转移进行到一定的次数时，比如到 n 次时，我们认为此时的采样集 (xn,xn+1,xn+2,...) 即是符合我们的平稳分布的对应样本集。 基于马尔科夫链的采样过程如下所示： 1) 输入马尔科夫链状态转移矩阵 P，设定状态转移次数阈值 n1，需要的样本个数 n2； 2) 从任意简单概率分布采样得到初始状态值 x0； 3) for t=0 to n1+n2−1: 从条件概率分布 P(x|xt) 中采样得到样本 xt+1 样本集 (xn1, xn1+1,...,xn1+n2−1) 即为我们需要的平稳分布对应的样本集。 如果假定我们可以得到我们需要采样样本的平稳分布所对应的马尔科夫链状态转移矩阵，那么我们就可以用马尔科夫链采样得到我们需要的样本集，进而进行蒙特卡罗模拟。但是一个重要的问题是，随意给定一个平稳分布 π，如何得到它所对应的马尔科夫链状态转移矩阵P呢？这是个大问题。我们绕了一圈似乎还是没有解决任意概率分布采样样本集的问题。 马尔科夫链的细致平稳条件 在解决从平稳分布 π，找到对应的马尔科夫链状态转移矩阵P之前，我们还需要先看看马尔科夫链的细致平稳条件。定义如下： 如果非周期马尔科夫链的状态转移矩阵 P 和概率分布 π(x) 对于所有的 i，j 满足： π(i)P(i,j)=π(j)P(j,i) 则称概率分布 π(x) 是状态转移矩阵 P 的平稳分布。 证明很简单,由细致平稳条件有：
将上式用矩阵表示即为： πP=π 即满足马尔可夫链的收敛性质。也就是说，只要我们找到了可以使概率分布 π(x) 满足细致平稳分布的矩阵 P 即可。这给了我们寻找从平稳分布 π，找到对应的马尔科夫链状态转移矩阵 P 的新思路。 不过不幸的是，仅仅从细致平稳条件还是很难找到合适的矩阵 P。比如我们的目标平稳分布是 π(x)，随机找一个马尔科夫链状态转移矩阵 Q，它是很难满足细致平稳条件的，即： π(i)Q(i,j)≠π(j)Q(j,i) 那么如何使这个等式满足呢？下面我们来看 MCMC 采样如何解决这个问题。 MCMC采样 由于一般情况下，目标平稳分布 π(x) 和某一个马尔科夫链状态转移矩阵 Q 不满足细致平稳条件，即 π(i)Q(i,j)≠π(j)Q(j,i) 我们可以对上式做一个改造，使细致平稳条件成立。方法是引入一个 α(i,j)，使上式可以取等号，即： π(i)Q(i,j)α(i,j)=π(j)Q(j,i)α(j,i) 问题是什么样的 α(i,j) 可以使等式成立呢？其实很简单，只要满足下两式即可： α(i,j)=π(j)Q(j,i) α(j,i)=π(i)Q(i,j) 这样，我们就得到了我们的分布 π(x) 对应的马尔科夫链状态转移矩阵 P，满足： P(i,j)=Q(i,j)α(i,j) 也就是说，我们的目标矩阵P可以通过任意一个马尔科夫链状态转移矩阵Q乘以 α(i,j) 得到。α(i,j) 我们有一般称之为接受率。取值在 [0,1] 之间，可以理解为一个概率值。即目标矩阵P可以通过任意一个马尔科夫链状态转移矩阵 Q 以一定的接受率获得。这个很像我们在上一篇文章《概率图模型-蒙特卡洛抽样》使用到的接受-拒绝采样，那里是以一个常用分布通过一定的接受-拒绝概率得到一个非常见分布，这里是以一个常见的马尔科夫链状态转移矩阵 Q 通过一定的接受-拒绝概率得到目标转移矩阵P,两者的解决问题思路是类似的。 好了，现在我们来总结下MCMC的采样过程。 1) 输入我们任意选定的马尔科夫链状态转移矩阵Q，平稳分布π(x)，设定状态转移次数阈值n1，需要的样本个数n2； 2) 从任意简单概率分布采样得到初始状态值x0； 3) for t=0 to n1+n2−1: a) 从条件概率分布Q(x|xt)中采样得到样本x∗ b) 从均匀分布采样u∼uniform[0,1] c) 如果u<α(xt,x∗)=π(x∗)Q(x∗,xt), 则接受转移xt→x∗，即xt+1=x∗ d) 否则不接受转移，t=max(t−1,0) 样本集(xn1,xn1+1,...,xn1+n2−1)即为我们需要的平稳分布对应的样本集。 上面这个过程基本上就是 MCMC 采样的完整采样理论了，但是这个采样算法还是比较难在实际中应用，为什么呢？问题在上面第三步的c步骤，接受率这儿。由于 α(xt,x∗) 可能非常的小，比如0.1，导致我们大部分的采样值都被拒绝转移，采样效率很低。有可能我们采样了上百万次马尔可夫链还没有收敛，也就是上面这个 n1 要非常非常的大，这让人难以接受，怎么办呢？这时就轮到我们的 M-H 采样出场了。 M-H采样 M-H 采样是 Metropolis-Hastings 采样的简称，这个算法首先由 Metropolis 提出，被 Hastings 改进，因此被称之为 Metropolis-Hastings 采样或 M-H 采样。M-H 采样解决了 MCMC 采样接受率过低的问题。 我们回到 MCMC 采样的细致平稳条件： π(i)Q(i,j)α(i,j)=π(j)Q(j,i)α(j,i) 我们采样效率低的原因是 α(i,j) 太小了，比如为 0.1，而 α(j,i) 为0.2。即： π(i)Q(i,j)×0.1=π(j)Q(j,i)×0.2 这时我们可以看到，如果两边同时扩大五倍，接受率提高到了 0.5，但是细致平稳条件却仍然是满足的，即： π(i)Q(i,j)×0.5=π(j)Q(j,i)×1 这样我们的接受率可以做如下改进，即： α(i,j)=min{π(j)Q(j,i)π(i)Q(i,j),1} 通过这个微小的改造，我们就得到了可以在实际应用中使用的M-H采样算法过程如下： 1) 输入我们任意选定的马尔科夫链状态转移矩阵Q，平稳分布π(x)，设定状态转移次数阈值n1，需要的样本个数n2； 2) 从任意简单概率分布采样得到初始状态值x0； 3) for t=0 to n1+n2−1: a) 从条件概率分布Q(x|xt)中采样得到样本x∗ b) 从均匀分布采样u∼uniform[0,1] c) 如果u<α(xt,x∗)=min{π(j)Q(j,i)π(i)Q(i,j),1}, 则接受转移xt→x∗，即xt+1=x∗ d) 否则不接受转移，t=max(t−1,0) 样本集(xn1,xn1+1,...,xn1+n2−1)即为我们需要的平稳分布对应的样本集。 很多时候，我们选择的马尔科夫链状态转移矩阵Q如果是对称的，即满足Q(i,j)=Q(j,i),这时我们的接受率可以进一步简化为： α(i,j)=min{π(j)π(i),1} M-H采样实例 在例子里，我们的目标平稳分布是一个均值3，标准差2的正态分布，而选择的马尔可夫链状态转移矩阵Q(i,j)的条件转移概率是以i为均值,方差1的正态分布在位置j的值。这个例子仅仅用来让大家加深对M-H采样过程的理解。毕竟一个普通的一维正态分布用不着去用M-H采样来获得样本。
可以看出，M-H 采样和真实分布之间，差异性非常小了。 M-H采样总结 M-H 采样已经可以很好的解决蒙特卡罗方法需要的任意概率分布的样本集的问题。但是 M-H 采样有两个缺点：一是需要计算接受率，在高维时计算量大。并且由于接受率的原因导致算法收敛时间变长。二是有些高维数据，特征的条件概率分布好求，但是特征的联合分布不好求。因此需要一个好的方法来改进 M-H 采样，这就是我们下面讲到的 Gibbs 采样。 重新寻找合适的细致平稳条件 如果非周期马尔科夫链的状态转移矩阵 P 和概率分布 π(x) 对于所有的 i,j 满足： π(i)P(i,j)=π(j)P(j,i) 则称概率分布 π(x) 是状态转移矩阵 P 的平稳分布。 在 M-H 采样中我们通过引入接受率使细致平稳条件满足。现在我们换一个思路。 从二维的数据分布开始，假设 π(x1,x2) 是一个二维联合数据分布，观察第一个特征维度相同的两个点：
容易发现下面两式成立：
由于两式的右边相等，因此我们有：
也就是：
观察上式再观察细致平稳条件的公式，我们发现在
这条直线上，如果用条件概率分布
作为马尔科夫链的状态转移概率，则任意两个点之间的转移满足细致平稳条件！ 这真是一个开心的发现，同样的道理，在在
这条直线上，如果用条件概率分布
作为马尔科夫链的状态转移概率，则任意两个点之间的转移也满足细致平稳条件。那是因为假如有一点
我们可以得到：
基于上面的发现，我们可以这样构造分布 π(x1,x2) 的马尔可夫链对应的状态转移矩阵P：
有了上面这个状态转移矩阵，我们很容易验证平面上的任意两点 E,F，满足细致平稳条件： π(E)P(E→F)=π(F)P(F→E) 二维Gibbs采样 利用上面找到的状态转移矩阵，我们就得到了二维Gibbs采样，这个采样需要两个维度之间的条件概率。具体过程如下： 1) 输入平稳分布 π(x1,x2)，设定状态转移次数阈值 n1，需要的样本个数 n2； 2) 随机初始化初始状态值 x1(0) 和 x2(0)； 3) for t=0 to n1+n2−1: a) 从条件概率分布 P(x2|x1(t)) 中采样得到样本 x2(t+1) b) 从条件概率分布 P(x1|x2(t+1)) 中采样得到样本x1(t+1) 样本集 {(x1(n1), x2(n1)), (x1(n1+1), x2(n1+1)), ... , (x1(n1+n2−1), x2(n1+n2−1))} 即为我们需要的平稳分布对应的样本集。 整个采样过程中，我们通过轮换坐标轴，采样的过程为： (x1(1),x2(1))→(x1(1),x2(2))→(x1(2),x2(2))→...→(x1(n1+n2−1),x2(n1+n2−1)) 用下图可以很直观的看出，采样是在两个坐标轴上不停的轮换的。当然，坐标轴轮换不是必须的，我们也可以每次随机选择一个坐标轴进行采样。不过常用的Gibbs采样的实现都是基于坐标轴轮换的。
多维Gibbs采样 上面的这个算法推广到多维的时候也是成立的。比如一个 n 维的概率分布 π(x1,x2,...xn)，我们可以通过在 n 个坐标轴上轮换采样，来得到新的样本。对于轮换到的任意一个坐标轴 xi 上的转移，马尔科夫链的状态转移概率为P(xi|x1,x2,...,xi−1,xi+1,...,xn)，即固定n−1个坐标轴，在某一个坐标轴上移动。 具体的算法过程如下： 1) 输入平稳分布π(x1,x2，...,xn)或者对应的所有特征的条件概率分布，设定状态转移次数阈值n1，需要的样本个数n2； 2) 随机初始化初始状态值(x1(0),x2(0),...,xn(0)) 3) for t=0 to n1+n2−1: a) 从条件概率分布 P(x1|x2(t),x3(t),...,xn(t)) 中采样得到样本x1(t+1) b) 从条件概率分布 P(x2|x1(t+1),x3(t),x4(t),...,xn(t))中采样得到样本x2(t+1) c) ... d) 从条件概率分布 P(xj|x1(t+1),x2(t+1),...,xj-1(t+1),xj+1(t),...,xn(t))中采样得到样本xj(t+1) e) ... f) 从条件概率分布 P(xn|x1(t+1),x2(t+1),...,xn-1(t+1))中采样得到样本xn(t+1) 样本集 {(x1(n1), x2(n1), ..., xn(n1)), ..., (x1(n1+n2−1), x2(n1+n2−1), ..., xn(n1+n2−1))} 即为我们需要的平稳分布对应的样本集。 整个采样过程和Lasso回归的坐标轴下降法算法非常类似，只不过Lasso回归是固定 n−1 个特征，对某一个特征求极值。而Gibbs采样是固定 n−1 个特征在某一个特征采样。 同样的，轮换坐标轴不是必须的，我们可以随机选择某一个坐标轴进行状态转移，只不过常用的Gibbs采样的实现都是基于坐标轴轮换的。 二维Gibbs采样实例 假设我们要采样的是一个二维正态分布 Norm(μ, σ)，其中： μ=(μ1,μ2)=(5,−1)
而采样过程中的需要的状态转移条件分布为：
具体的代码如下：
最后，绘图查看模拟出来的数据，是否我们需要的数据：
|
|
|
