【转】Markov chain monte carlo

来源：http://www.xperseverance.net/blogs/2012/04/902/

本篇标题虽然是PRML，但是内容来自三本书，另外两本是《Introduction to Probability Models》（以下简称IPM）和《Markov Chain Monte Carlo In Practice》（以下简称MCMCP），都是经典教材，其中前一本内容讲的比较基础，例子很丰富，但是翻译版《应用随机过程-概率模型导论》个人不推荐，翻译较粗糙，只能用来对比中英文数学概念，呵呵。

引用清华大学教授@刘洋THU大牛在微博里的话：“PRML不适合初学者”。我觉得很在理，PRML虽然内容看起来诱人，但真的不太适合偶等初学的，书中有的地方吧，解释的很细腻，但是大部分内容都解释的不够基础和详细。所以光看一本书是不行的……

1. Markov chains

马尔科夫链的基本定义和性质就不废话了，只说一下，一个马尔科夫链可以用如下两项定义：

(1)初始状态概率p(z(0))：其中z的上标(i)表示马尔科夫链跳转的第i步，也就是时间i，z(0)表示在时间0这个马尔科夫链所在的初始状态，那么p(z(0))表示什么意思呢？我的理解大概是比如一个马尔科夫链有3个状态1、2、3，那么p(z(0))就是这个马尔科夫链一开始所在的初始状态的概率分布，比如{0.1,0.6,0.3}

(2)转移概率Tm(zm,z(m+1))≡p(z(m+1)|z(m))：这个PRML上叫做transition probability，MCMCP上叫transition kernel，就是转移矩阵啦。注意PRML上Tm(zm,z(m+1))这样的写法可不是表示联合概率，他是表示一个条件概率，表示从状态z(m)转移到z(m+1)的概率，这个感觉挺容易误会。这里只要记住是从前一个转移到后一个就行。一般随机过程都用Pij书，觉得PRML上用z(m)和z(m+1)更科学，应为i，j是状态，而这里的z是随机变量。

对于上述一个特定的随机变量z(m+1)，它的概率等于它所有上一步状态的条件概率取边缘：

p(z(m+1))=∑z(m)p(z(m+1)|z(m))p(z(m))

以前看随机过程都一直没有深入理解好这里，这个式子就是求第m+1步时的随机变量z在某一个状态时的概率。当m→∞时，它就成了极限概率，也就是通常随机过程书中写的π，事实上πj本来也就是这样求的：πj=limn→∞Pnij=∑i=0∞πiPij，这回要记好了这些个理解以免以后又忘记。

 

2. Time Reversible Markov Chains

时间可逆的马尔科夫链，这个在IPM里面解释的很好，PRML也太直接了点，偶等一点都不了解的人咋可能一下子就看懂嘛。由于以下内容来自IPM，就沿用IPM的符号系统了。看图片：

把一个马尔科夫链反过来看，它也是一个马尔科夫链，IPM上有这个的证明。设这样时间反过来的马尔科夫链的转移概率为Qij，则Qij=πjPjiπi，因为：

如果Qij=Pij，那么就说这个马尔科夫链时间可逆，这样的话把这个条件代入上式，得到πiPij=πjPji。对于时间可逆的马尔科夫链，有个充分必要条件：任意一条回路正向走和反向走的概率相等。比如只有三个状态组成一条回路，那么PikPkjPji=PijPjkPki

3. Markov Chain Monte Carlo

假设X是一个离散随机变量，而我们相对某个特殊函数h计算均值：

θ=E[h(X)]=∑j=1∞h(xj)P{X=xj}

这里有个疑问，IPM上说h(x_j)在计算上有困难，而MCMCP却说要产生独立的samples使他们符合想要的概率分布很困难，到底是哪一个原因有点令人费解。不过两者想要表示的意思应该都是相同的，就是直接用数学方法求θ=E[h(X)]很难，所以只能用统计模拟方法近期求θ，称为蒙特卡罗方法(Monte Carlo integration)。蒙特卡罗方法通过随机数产生一系列独立同分布的随机变量Xi，然后由强大数定理导出下式以估计θ：

θ=limn→∞∑i=1nh(Xi)n

接下来的主要任务是如何通过随机数模拟，产生一系列的独立同分布samples并且它们符合一个我们指定的概率分布？这就要用到Hastings-Metropolis算法。

先来说说马尔可夫链的极限概率。我们把一个随机过程的某个状态作为初始点，然后根据转移概率在这个过程上不断地sample，得到一个马尔可夫链（因为马尔可夫链的当前状态的概率只取决于前一状态，所以只要有个起点就可以不断sample），经过无数步sample之后，当前步所在状态的概率会逐渐独立于初始状态，并趋向于这个随机过程的极限概率。所以在无数步之后，我们得到的一系列samples就好像是从这个马尔可夫链独立同分布采样得到的结果。这个内容在MCMCP上解释的很清楚，而在IPM上说的就不是那么明白了。这个理解非常关键，马尔可夫蒙特卡罗方法之所以称为“马尔科夫…”，就是因为利用了下面这个性质：当用于采样的马尔可夫链的极限概率等于我们想要采样的指定概率分布的时候，在这个马尔可夫链上大量采样就好是从我们指定的概率分布上独立同分布采样的一样（这是MCMCP上说的）。PRML上有另一种说法：在这样的马尔科夫链上大量采样，相邻的采样间不是独立的，想要获得独立采样，只要每间隔m个结果保留采样就行了。

然后就是如何产生一个极限概率等于我们指定概率分布的马尔可夫链了，这部分IPM讲的很好。2012@4@10 23:20明日继续今日已晚。

4. Hastings-Metropolis算法

本算法的目的是要构建一个马尔科夫（设为A）链使这个链的极限概率等于某个指定的概率分布，并在这个马尔科夫链上采样得到X1,X2,…Xn。首先我们要构建这样一个目标马尔科夫链A，令Q是一个转移概率矩阵，它表示当Xn=i时，生成随机变量Y的概率为P{Y=j}=q(i,j)。

构建过程是这样的，假设现在随机过程处于状态i，开始游走到下一个状态，我们首先为Y随机挑选一个状态j，接下来不是直接令下一个状态Xn+1=j，而是Xn+1被赋值为j的概率是α(i,j)，在剩下的概率1−α(i,j)下Xn+1仍旧被赋值为i，即保持原状态。这样得到一个马尔科夫链A，其转移概率为：

Pi,j=q(i,j)α(i,j)

Pi,i=q(i,i)+∑k≠iq(i,k)(1−α(i,k))

注意这里要仔细区分Q和P，Q只是我们定义的一个状态之间的转换概率，它不是新构建的马尔科夫链的转移概率！P才是我们要构建的目标马尔科夫链A的transition probability。

然后可以通过指定的α(i,j)的值来把这马尔科夫链变成时间可逆的。根据时间可逆马尔科夫链的性质，即平稳概率π(j)（通常也是极限概率）满足：π(i)Pi,j=π(j)Pj,i，并推得：

π(i)q(i,j)α(i,j)=π(j)q(j,i)α(j,i)    式4.31

为了满足上面这个条件，我们令：

α(i,j)=min(π(j)q(j,i)π(i)q(i,j),1)   式4.32

这个式子设计成这样，是因为假如α(i,j)=π(j)q(j,i)π(i)q(i,j)小于1，那么α(j,i)=π(i)q(i,j)π(j)q(j,i)相当于把上下分子分母对调，就会大于1，这样由于min函数存在，α(j,i)=1，把这一对α(i,j)和α(j,i)的值代入4.31，就能够满足等式啦。这个设计挺巧妙。

最后Hastings-Metropolis算法就成了这个样子，算法里面的q(.|Xt)可以拥有任何形式，通常是人为选定的：