HMC(Hamiltonian Monte Carlo抽样算法详细介绍)

[toc]

原文博客地址：http://blog.csdn.net/qy20115549/article/details/54561643

**未经本人，允许禁止转载。**

#Hamiltonian Monte Carlo简介

Hamiltonian Monte Carlo也叫Hybrid Monte Carlo，是一种快速抽样方法。在MCMC算法中随机游走的方式使得Markov链收敛于固定的分布p(x)  然效率不高。Hamiltonian or Hybrid Monte Carlo (HMC)这种MCMC算法应用的是物理系统中动力学的概念来计算Markov链中的未来状态，而不是概率分布。相比于随机游走，通过这种方式，能够更加高效的分析状态空间，从而达到更快的收敛。接下来，将一步一步介绍Hamiltonian dynamics(Hamiltonian动力学)中的基础性分析及相关概念。之后，介绍Hamiltonian动力学是如何用于MCMC抽样算法中Markov链建议函数的。

#Hamiltonian dynamics的物理含义
由于下面写的东西有公式，而敲公式太浪费时间，在此，把我的学习笔记，通过图片的方式，展现。如果有需要电子版的话，请发邮件。

![这里写图片描述](http://img.blog.csdn.net/20170115131146356?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcXkyMDExNTU0OQ==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)

物理学中的证明，大家也可以找一些关于Hamiltonian动力学方面的博客看看。

#Simulating Hamiltonian dynamics — the Leap Frog Method

![这里写图片描述](http://img.blog.csdn.net/20170115131405202?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcXkyMDExNTU0OQ==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)

![这里写图片描述](http://img.blog.csdn.net/20170115131519656?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcXkyMDExNTU0OQ==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)

##Example 1: Simulating Hamiltonian dynamics of an harmonic oscillator

![这里写图片描述](http://img.blog.csdn.net/20170115131748020?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcXkyMDExNTU0OQ==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)
![这里写图片描述](http://img.blog.csdn.net/20170115131839443?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcXkyMDExNTU0OQ==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)

![这里写图片描述](http://img.blog.csdn.net/20170115132035164?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcXkyMDExNTU0OQ==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)

以下是matlab程序：

```matlab
% EXAMPLE 1: SIMULATING HAMILTONIAN DYNAMICS
%            OF HARMONIC OSCILLATOR
% STEP SIZE
delta = 0.1;
 
% # LEAP FROG
L = 70;
 
% DEFINE KINETIC ENERGY FUNCTION
K = inline('p^2/2','p');
 
% DEFINE POTENTIAL ENERGY FUNCTION FOR SPRING (K =1)
U = inline('1/2*x^2','x');
 
% DEFINE GRADIENT OF POTENTIAL ENERGY
dU = inline('x','x');
 
% INITIAL CONDITIONS
x0 = -4; % POSTIION
p0 = 1;  % MOMENTUM
figure
 
%% SIMULATE HAMILTONIAN DYNAMICS WITH LEAPFROG METHOD
% FIRST HALF STEP FOR MOMENTUM
pStep = p0 - delta/2*dU(x0)';
 
% FIRST FULL STEP FOR POSITION/SAMPLE
xStep = x0 + delta*pStep;
 
% FULL STEPS
for jL = 1:L-1
    % UPDATE MOMENTUM
    pStep = pStep - delta*dU(xStep);
    % UPDATE POSITION
    xStep = xStep + delta*pStep;
    % UPDATE DISPLAYS
    subplot(211), cla
    hold on;
    xx = linspace(-6,xStep,1000);
    plot(xx,sin(6*linspace(0,2*pi,1000)),'k-');
    plot(xStep+.5,0,'bo','Linewidth',20)
    xlim([-6 6]);ylim([-1 1])
    hold off;
    title('Harmonic Oscillator')
    subplot(223), cla
    b = bar([U(xStep),K(pStep);0,U(xStep)+K(pStep)],'stacked');
    set(gca,'xTickLabel',{'U+K','H'})
    ylim([0 10]);
    title('Energy')
    subplot(224);
    plot(xStep,pStep,'ko','Linewidth',20);
        xlim([-6 6]); ylim([-6 6]); axis square
    xlabel('x'); ylabel('p');
    title('Phase Space')
    pause(.1)
end
% (LAST HALF STEP FOR MOMENTUM)
pStep = pStep - delta/2*dU(xStep);
```

#Hamiltonian dynamics and the target distribution 
<center>
![这里写图片描述](http://img.blog.csdn.net/20170115132542222?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcXkyMDExNTU0OQ==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)
![这里写图片描述](http://img.blog.csdn.net/20170115132721442?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcXkyMDExNTU0OQ==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)
</center>

#Hamiltonian Monte Carlo
<center>
![这里写图片描述](http://img.blog.csdn.net/20170115132856941?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcXkyMDExNTU0OQ==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)

![这里写图片描述](http://img.blog.csdn.net/20170115132954365?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcXkyMDExNTU0OQ==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)
</center>
或者更详细的伪代码：
<center>
![这里写图片描述](http://img.blog.csdn.net/20170115133031271?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcXkyMDExNTU0OQ==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)
</center>

##Example 2: Hamiltonian Monte for sampling a Bivariate Normal distribution
<center>
![这里写图片描述](http://img.blog.csdn.net/20170115133251383?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcXkyMDExNTU0OQ==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)
![这里写图片描述](http://img.blog.csdn.net/20170115133322009?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcXkyMDExNTU0OQ==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)
</center>

HMC的matlab代码为：

```matlab
% EXAMPLE 2: HYBRID MONTE CARLO SAMPLING -- BIVARIATE NORMAL
rand('seed',12345);
randn('seed',12345);
 
% STEP SIZE
delta = 0.3;
nSamples = 1000;
L = 20;
 
% DEFINE POTENTIAL ENERGY FUNCTION
U = inline('transp(x)*inv([1,.8;.8,1])*x','x');
 
% DEFINE GRADIENT OF POTENTIAL ENERGY
dU = inline('transp(x)*inv([1,.8;.8,1])','x');
 
% DEFINE KINETIC ENERGY FUNCTION
K = inline('sum((transp(p)*p))/2','p');
 
% INITIAL STATE
x = zeros(2,nSamples);
x0 = [0;6];
x(:,1) = x0;
 
t = 1;
while t < nSamples
    t = t + 1;
 
    % SAMPLE RANDOM MOMENTUM
    p0 = randn(2,1);
 
    %% SIMULATE HAMILTONIAN DYNAMICS
    % FIRST 1/2 STEP OF MOMENTUM
    pStar = p0 - delta/2*dU(x(:,t-1))';
 
    % FIRST FULL STEP FOR POSITION/SAMPLE
    xStar = x(:,t-1) + delta*pStar;
 
    % FULL STEPS
    for jL = 1:L-1
        % MOMENTUM
        pStar = pStar - delta*dU(xStar)';
        % POSITION/SAMPLE
        xStar = xStar + delta*pStar;
    end
 
    % LAST HALP STEP
    pStar = pStar - delta/2*dU(xStar)';
 
    % COULD NEGATE MOMENTUM HERE TO LEAVE
    % THE PROPOSAL DISTRIBUTION SYMMETRIC.
    % HOWEVER WE THROW THIS AWAY FOR NEXT
    % SAMPLE, SO IT DOESN'T MATTER
 
    % EVALUATE ENERGIES AT
    % START AND END OF TRAJECTORY
    U0 = U(x(:,t-1));
    UStar = U(xStar);
 
    K0 = K(p0);
    KStar = K(pStar);
 
    % ACCEPTANCE/REJECTION CRITERION
    alpha = min(1,exp((U0 + K0) - (UStar + KStar)));
 
    u = rand;
    if u < alpha
        x(:,t) = xStar;
    else
        x(:,t) = x(:,t-1);
    end
end
 
% DISPLAY
figure
scatter(x(1,:),x(2,:),'k.'); hold on;
plot(x(1,1:50),x(2,1:50),'ro-','Linewidth',2);
xlim([-6 6]); ylim([-6 6]);
legend({'Samples','1st 50 States'},'Location','Northwest')
title('Hamiltonian Monte Carlo')
```

#参考

(1)
<a href="https://theclevermachine.wordpress.com/2012/11/18/mcmc-hamiltonian-monte-carlo-a-k-a-hybrid-monte-carlo/" title="Title">
MCMC: Hamiltonian Monte Carlo (a.k.a. Hybrid Monte Carlo) </a>
(2)
Niederreiter H. Quasi‐Monte Carlo Methods[M]. John Wiley & Sons, Ltd, 2010.

#相关文章

Chen T, Fox E B, Guestrin C. Stochastic Gradient Hamiltonian Monte Carlo[C]//ICML. 2014: 1683-1691.

Neal R M. MCMC using Hamiltonian dynamics[J]. Handbook of Markov Chain Monte Carlo, 2011, 2: 113-162.

Shahbaba B, Lan S, Johnson W O, et al. Split hamiltonian monte carlo[J]. Statistics and Computing, 2014, 24(3): 339-349.

HMC(Hamiltonian Monte Carlo抽样算法详细介绍)

欢迎大家关注DataLearner官方微信，接受最新的AI技术推送

相关博客

最热博客