Hamiltonian Monte Carlo算法理解

因为RBN，学习了一会儿MCMC，理解不太深，但在接触Bayesian GAN时碰到HMC又很晕，好在看了两篇博客，重新梳理了MCMC，对HMC有了点初步理解，整理如下。
两篇博客一篇是CSDN，一篇是WordPress上的。CSDN应该是翻译WordPress上的。

1. theclevermachine博客
   https://theclevermachine.wordpress.com/2012/11/18/mcmc-hamiltonian-monte-carlo-a-k-a-hybrid-monte-carlo/
2. CSDN翻译博客，http://blog.csdn.net/qy20115549/article/details/54561643

上述博客对细致平稳条件没有多加解释，因此理解有点困难。

1. Monte Carlo算法是什么？

这个总是容易忘记，Monte Carlo算法就是随机过程的模拟算法，模拟的输出就是随机过程的样本。通过对模拟得到的样本进行统计，我们就可以求积分，求最优解，当然求出来的都是可以实际应用的近似值（Approximation）。

1.1 求积分公式

如果需要求积分∫baf(x)dx，可以模拟一个概率分布为q(x)过程，采样得到N个样本，从而该积分可以用下式进行近似：

∫baf(x)dx=∫baf(x)q(x)q(x)dx≈1N∑if(xi)q(xi)

1.2 求最大值

如果需要求q(x)的最大值，可以将q(x)归一化为概率分布，模拟q(x)过程，采样得到N个样本，统计样本分布最密的区间就是最优值所在的区间。

2. 如何模拟？——简单的接受拒绝采样

对于任意的概率分布p(z)计算机很难直接模拟，计算机能够模拟的是均匀分布、正态分布等一些简单分布。接受拒绝采样算法就是从容易模拟的概率分布中得到样本，在根据一定的接受拒绝规则来再次选择样本。如下图：

第一步：红色p(z)难以直接采样，蓝色q(z)可以直接采样，先根据q(z)采样得到z0；注意要合理选择k使得kq(z)≥p(z)
第二步：再在[0 kq(z0)]采样一个u0，如果u0>p(z0),则拒绝z0。
因此，z0被接受概率为q(z0)(p(z0)/kq(z0)=p(z0)/k，所以在所有采样样本空间中被接受的样本分布为p(x)/k，在被接受样本空间中样本分布为p(x)，样本的接受概率为1/k，因此，如果k太大，则样本拒绝率太高。

3. MCMC模拟

马尔科夫链随机模拟的原理是基于马尔科夫链的平稳特性，对于满足一定条件的状态转移概率矩阵P=[pij]总对应一个马尔科夫平稳状态π。平稳状态和状态转移概率矩阵之间满足细致平稳条件，即两个状态相互转移的概率相等：

πipij=πjpji

对于需要模拟的概率分布π难以直接得到对应的状态转移概率矩阵，可以通过接受拒绝算法将任一状态转移概率矩阵根据细致平稳条件进行改造，从而使得马尔科夫链收敛到平稳状态π。
对于多维情况，则可应用Gibbs算法，从条件概率来构造状态转移概率矩阵。

3.1 Metropolis算法及Metropolis-Hasting算法

下面为Metropolis算法，其中α(i,j)，α(j,i)为接受拒绝概率，q(i,j)为状态i到状态j的初始转移概率。如果，α(i,j)=p(j)q(j,i)α(j,i)=p(i)q(i,j)，则下面细致平稳条件满足：

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

为了提高状态的接受概率，可以将等式两边的α同时放大，最高到1，这就是Metropolis-Hasting算法。

p(i)q(i,j)min[1,p(j)q(j,i)p(i)q(i,j)]=p(j)q(j,i)min[1,p(i)q(i,j)p(j)q(j,i)]

3.2 Gibbs算法

对于多维状态的情况任意两状态A=(X1,X2)B=(X1,X3)，由条件概率公式可得到细致平稳条件的状态转移概率矩阵：

p(A)p(X3|X1)=p(X1,X2)p(X3|X1)=p(X1)p(X2|X1)p(X3|X1)=p(X1,X3)p(X2|X1)=p(B)p(X2|X1)

4. Hamiltonian Monte Carlo算法

HMC算法应该也算是一种Metropolis-Hasting算法，通过Hamiltonian偏微分方程来实现初始状态转移，通过接受拒绝算法来满足细致平稳条件。与普通的Metropolis-Hasting算法的区别在于其初始状态转移是确定的，即概率为1（相比普通的q(i,j)≤1随机游走的方式提高了接受概率，从而提高了收敛速度）。只是由于Hamiltonian偏微分方程离散化计算误差导致产生了状态偏移，从而需要对偏移状态进行接受拒绝选择来达到细致平稳条件的满足。

4.1 Hamiltonian系统

系统方程：

H(q,p)=U(q)+K(p)∂qi∂t=∂H∂pi=∂K(p)∂pi∂pi∂t=∂H∂qi=−∂U(q)∂qi

其中U相当于系统势能，K相当于系统动能。系统能量是守恒的，即H不变。根据一个初始状态和时间能够计算目标状态，系统能量不变。
依据Hamiltonian系统构建如下Markov链有如下性质：

1）对于任一q的初始状态q0，按照分布1Zpexp(−K(p))采样得到p0
2） 从p0,q0根据Hamiltonian方程得到p1,q1，将q1赋给q0，回到第一步继续计算q1。
3）重复1，2步骤，得到q的采样序列。则q稳态的概率分布为1Zqexp(−U(q))

证明：
假设q的稳态分布为，1Zqexp(−f(q))，
所以Hamiltonian系统Markov链的状态转移概率等式为：

1Zpexp(−K(p1))1Zqexp(−f(q1))=1Zpexp(−K(p2))1Zqexp(−f(q2))

所以，K(p)+f(q)为常数，从而，

f(q)=U(q)+C

所以，q的稳态分布为1Zqexp(−U(q))

4.2 HMC算法

HMC的目标是模拟实现1Z1exp(−U(q))！
实现思路是：

1）首先构造一个容易模拟的正态分布1Zpexp(−K(p))
2）模拟Hamiltionian系统的Markov链实现对q的采样
3）根据前面Hamiltonian系统Markov链性质，q的稳态分布为1Zqexp(−U(q))

算法步骤：
任意取一个q的初始状态q1
重复如下计算：
1）依据正态分布采样得到p1和q1
2）以p1和q1为初值，根据Hamiltonian方程计算L步，得到(q2,p2)
3）计算(q2,p2)的接受拒绝概率：

α=min[1,exp(−H(q2,p2)+H(q1,p1))]

4）从[0,1]均匀分布抽样得到一随机数u，如果u≤α，则接受采样，即赋值q2→q1，回到第一步。
很容易证明该算法满足如下细致平稳条件等式：

1Zexp(−H(q1,p1))min[1,exp(−H(q2,p2)+H(q1,p1))]=1Zexp(−H(q2,p2))min[1,exp(−H(q1,p1)+H(q2,p2))]

之所以要进行接受拒绝操作是因为由于离散计算导致Hamiltonian系统Markov链转移概率发生变化。

4.3 HMC模拟采样双正态分布的python代码

仿照theclevermachine博客中的matlab代码，用python实现HMC模拟采样双正态分布，代码在csdn资源上：http://download.csdn.net/download/tkyjqh/10176310
发现matlab实现中有两个细节差异：
1）U函数没有除2
2）Hamiltonian计算并没有全部采用Leap-Frog算法。
在python实现时，试图改进。结果如下：
1）采用Leap-Frog算法，U函数不除2，接受概率1

2）采用Leap-Frog算法，U函数除2，接受概率1

3）采用Euller算法，U函数不除2，接受概率0.73

4）采用Euller算法，U函数除2，接受概率0.98

统计接受概率可以发现Leap-Frog算法基本上全接受，采用Euller算法，接受概率将降低，而且U是否除2将影响接受概率。