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

冷不防 2022-07-11 00:26 226阅读 0赞

*  Hamiltonian Monte Carlo简介
 *  Hamiltonian dynamics的物理含义
 *  Simulating Hamiltonian dynamics the Leap Frog Method
    
     *  Example 1 Simulating Hamiltonian dynamics of an harmonic oscillator
 *  Hamiltonian dynamics and the target distribution
 *  Hamiltonian Monte Carlo
    
     *  Example 2 Hamiltonian Monte for sampling a Bivariate Normal distribution
 *  参考
 *  相关文章

**本文作者：合肥工业大学 管理学院 钱洋 email：1563178220@qq.com 内容可能有不到之处，欢迎交流。**  
**未经本人，允许禁止转载。**

**下面是本文博客的另一个地址，该网站是师兄弄得一个专门做机器学习的网站，非常不错。**  
[http://www.datalearner.com/blog/1051484459699809][http_www.datalearner.com_blog_1051484459699809]

# 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的物理含义 #

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

![这里写图片描述][SouthEast]

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

# Simulating Hamiltonian dynamics — the Leap Frog Method #

![这里写图片描述][SouthEast 1]  
  
  
![这里写图片描述][SouthEast 2]

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

![这里写图片描述][SouthEast 3]  
![这里写图片描述][SouthEast 4]  
  
  
![这里写图片描述][SouthEast 5]  
  
以下是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 #

![这里写图片描述][SouthEast 6]  
![这里写图片描述][SouthEast 7]

# Hamiltonian Monte Carlo #

![这里写图片描述][SouthEast 8]

![这里写图片描述][SouthEast 9]  
  
或者更详细的伪代码：

![这里写图片描述][SouthEast 10]

## Example 2: Hamiltonian Monte for sampling a Bivariate Normal distribution ##

![这里写图片描述][SouthEast 11]  
![这里写图片描述][SouthEast 12]

HMC的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)  
[  
MCMC: Hamiltonian Monte Carlo (a.k.a. Hybrid Monte Carlo) ][MCMC_ Hamiltonian Monte Carlo _a.k.a. Hybrid Monte Carlo_]  
(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.

[http_www.datalearner.com_blog_1051484459699809]: http://www.datalearner.com/blog/1051484459699809
[SouthEast]: /images/20220711/993c1077ff9d4968ab88af2bd36a21dd.png
[SouthEast 1]: /images/20220711/af5230d3a9bb47db812f86df934d6ca0.png
[SouthEast 2]: /images/20220711/3b1019e74aa94ffd9db576c4d5ada2ba.png
[SouthEast 3]: /images/20220711/60da438174eb4f1b94964df408819bea.png
[SouthEast 4]: /images/20220711/63fae4c89b794c9e8c6741fde8e2dd01.png
[SouthEast 5]: /images/20220711/8fc0ecf50cab460bb8198aed63193f4d.png
[SouthEast 6]: /images/20220711/bed1dcaee0fd4e59979cb8480ce56af9.png
[SouthEast 7]: /images/20220711/1effaa9c5e024abeb9efb5d1e9647888.png
[SouthEast 8]: /images/20220711/0ae24e36899b46f487169e27353394ea.png
[SouthEast 9]: /images/20220711/e24cfdcce89c41ba8e3e12cb88a218b0.png
[SouthEast 10]: /images/20220711/9e09b05f1f8c40c2bfbf770da7e58b1f.png
[SouthEast 11]: /images/20220711/fa33275241f14cf599043987215a7804.png
[SouthEast 12]: /images/20220711/f0e2598a00f5469bb006ba7846dac162.png
[MCMC_ Hamiltonian Monte Carlo _a.k.a. Hybrid Monte Carlo_]: https://theclevermachine.wordpress.com/2012/11/18/mcmc-hamiltonian-monte-carlo-a-k-a-hybrid-monte-carlo/