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

标签:#HMC##抽样# 时间:2017-01-16 14:01:22 作者:十七岁的雨季

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

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

这里写图片描述

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

Simulating Hamiltonian dynamics — the Leap Frog Method

这里写图片描述

这里写图片描述

Example 1: Simulating Hamiltonian dynamics of an harmonic oscillator

这里写图片描述
这里写图片描述

这里写图片描述

以下是matlab程序:

  1. % EXAMPLE 1: SIMULATING HAMILTONIAN DYNAMICS
  2. % OF HARMONIC OSCILLATOR
  3. % STEP SIZE
  4. delta = 0.1;
  5. % # LEAP FROG
  6. L = 70;
  7. % DEFINE KINETIC ENERGY FUNCTION
  8. K = inline('p^2/2','p');
  9. % DEFINE POTENTIAL ENERGY FUNCTION FOR SPRING (K =1)
  10. U = inline('1/2*x^2','x');
  11. % DEFINE GRADIENT OF POTENTIAL ENERGY
  12. dU = inline('x','x');
  13. % INITIAL CONDITIONS
  14. x0 = -4; % POSTIION
  15. p0 = 1; % MOMENTUM
  16. figure
  17. %% SIMULATE HAMILTONIAN DYNAMICS WITH LEAPFROG METHOD
  18. % FIRST HALF STEP FOR MOMENTUM
  19. pStep = p0 - delta/2*dU(x0)';
  20. % FIRST FULL STEP FOR POSITION/SAMPLE
  21. xStep = x0 + delta*pStep;
  22. % FULL STEPS
  23. for jL = 1:L-1
  24. % UPDATE MOMENTUM
  25. pStep = pStep - delta*dU(xStep);
  26. % UPDATE POSITION
  27. xStep = xStep + delta*pStep;
  28. % UPDATE DISPLAYS
  29. subplot(211), cla
  30. hold on;
  31. xx = linspace(-6,xStep,1000);
  32. plot(xx,sin(6*linspace(0,2*pi,1000)),'k-');
  33. plot(xStep+.5,0,'bo','Linewidth',20)
  34. xlim([-6 6]);ylim([-1 1])
  35. hold off;
  36. title('Harmonic Oscillator')
  37. subplot(223), cla
  38. b = bar([U(xStep),K(pStep);0,U(xStep)+K(pStep)],'stacked');
  39. set(gca,'xTickLabel',{'U+K','H'})
  40. ylim([0 10]);
  41. title('Energy')
  42. subplot(224);
  43. plot(xStep,pStep,'ko','Linewidth',20);
  44. xlim([-6 6]); ylim([-6 6]); axis square
  45. xlabel('x'); ylabel('p');
  46. title('Phase Space')
  47. pause(.1)
  48. end
  49. % (LAST HALF STEP FOR MOMENTUM)
  50. pStep = pStep - delta/2*dU(xStep);

Hamiltonian dynamics and the target distribution


这里写图片描述
这里写图片描述

Hamiltonian Monte Carlo


这里写图片描述

这里写图片描述

或者更详细的伪代码:

这里写图片描述

Example 2: Hamiltonian Monte for sampling a Bivariate Normal distribution


这里写图片描述
这里写图片描述

HMC的matlab代码为:

  1. % EXAMPLE 2: HYBRID MONTE CARLO SAMPLING -- BIVARIATE NORMAL
  2. rand('seed',12345);
  3. randn('seed',12345);
  4. % STEP SIZE
  5. delta = 0.3;
  6. nSamples = 1000;
  7. L = 20;
  8. % DEFINE POTENTIAL ENERGY FUNCTION
  9. U = inline('transp(x)*inv([1,.8;.8,1])*x','x');
  10. % DEFINE GRADIENT OF POTENTIAL ENERGY
  11. dU = inline('transp(x)*inv([1,.8;.8,1])','x');
  12. % DEFINE KINETIC ENERGY FUNCTION
  13. K = inline('sum((transp(p)*p))/2','p');
  14. % INITIAL STATE
  15. x = zeros(2,nSamples);
  16. x0 = [0;6];
  17. x(:,1) = x0;
  18. t = 1;
  19. while t < nSamples
  20. t = t + 1;
  21. % SAMPLE RANDOM MOMENTUM
  22. p0 = randn(2,1);
  23. %% SIMULATE HAMILTONIAN DYNAMICS
  24. % FIRST 1/2 STEP OF MOMENTUM
  25. pStar = p0 - delta/2*dU(x(:,t-1))';
  26. % FIRST FULL STEP FOR POSITION/SAMPLE
  27. xStar = x(:,t-1) + delta*pStar;
  28. % FULL STEPS
  29. for jL = 1:L-1
  30. % MOMENTUM
  31. pStar = pStar - delta*dU(xStar)';
  32. % POSITION/SAMPLE
  33. xStar = xStar + delta*pStar;
  34. end
  35. % LAST HALP STEP
  36. pStar = pStar - delta/2*dU(xStar)';
  37. % COULD NEGATE MOMENTUM HERE TO LEAVE
  38. % THE PROPOSAL DISTRIBUTION SYMMETRIC.
  39. % HOWEVER WE THROW THIS AWAY FOR NEXT
  40. % SAMPLE, SO IT DOESN'T MATTER
  41. % EVALUATE ENERGIES AT
  42. % START AND END OF TRAJECTORY
  43. U0 = U(x(:,t-1));
  44. UStar = U(xStar);
  45. K0 = K(p0);
  46. KStar = K(pStar);
  47. % ACCEPTANCE/REJECTION CRITERION
  48. alpha = min(1,exp((U0 + K0) - (UStar + KStar)));
  49. u = rand;
  50. if u < alpha
  51. x(:,t) = xStar;
  52. else
  53. x(:,t) = x(:,t-1);
  54. end
  55. end
  56. % DISPLAY
  57. figure
  58. scatter(x(1,:),x(2,:),'k.'); hold on;
  59. plot(x(1,1:50),x(2,1:50),'ro-','Linewidth',2);
  60. xlim([-6 6]); ylim([-6 6]);
  61. legend({'Samples','1st 50 States'},'Location','Northwest')
  62. title('Hamiltonian Monte Carlo')

参考

(1)

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.

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