1function [theta, rho] = leapfrog(logmu, gradlogmu, theta, rho, h, ell)
2 h2=h*h;
3 for k = 1:ell
4 g=gradlogmu(theta);
5 rho_half = rho + (h / 2) * g;
6 theta = theta + h * rho_half;
7 g=gradlogmu(theta);
8 rho = rho_half + (h / 2) * g;
9 end
10 % for k = 1:ell
11 % g=gradlogmu(theta);
12 % rho_half = rho + (h / 2) * g./(1+h2*abs(g));
13 % theta = theta + h * rho_half;
14 % g=gradlogmu(theta);
15 % rho = rho_half + (h / 2) * g./(1+h2*abs(g));
16 % end
17end