31:T997,md`# Hamiltonian Markov Chain Monte Carlo Hamiltonian Monte Carlo (HMC) is an algorithm for drawing samples from a probability distribution for which we know its density and its gradient. In basic [random walk metropolis](https://observablehq.com/@herbps10/random-walk-metropolis), the sampler explores the parameter space by a random walk. The sampler makes a proposal by drawing from distribution centered on the current state. The performance of the sampler depends on the details of the proposal distribution. A proposal distribution that is too small will cause the sampler to explore the target distribution too slowly, while if the proposals are too far away from the current state then the likelihood of acceptance decreases and many proposals need to be made before one is accepted. HMC takes a different approach to proposing new states that allows the sampler to propose new states that are far away from the current state, but still have high probability of acceptance. HMC proposes new positions by simulating Hamiltonian trajectories over the surface defined by the negative log density. The algorithm uses leapfrog integration (see [visualization](https://beta.observablehq.com/@herbps10/leapfrog-integrator)) to solve the proposal trajectories.

1. Recall that momentum is the product of mass and velocity: ${tex.block` p = mv `} If we set the mass of the particle to ${tex`1`}, then momentum and velocity are equivalent. Whenever you see "momentum" in HMC, you can also think about it as "velocity".

HMC works by adding an auxiliary variable representing the momentum of a fictitious particle on the negative log-density surface¹. The particle starts each trajectory in a random direction and momentum. Leapfrog integration is used to estimate the path the particle takes over the negative log-density surface. After a specified number of steps, the current position of the particle is proposed as a new position for the sampler to take. There is an accept/reject step, which exists for different reasons than traditional MCMC; the accept/reject step is used to correct for numerical errors accumulated in the integration of the trajectory. Trajectories with high error are less likely to be accepted by the sampler. Inspired by the [The Markov-chain Monte Carlo Interactive Gallery](https://chi-feng.github.io/mcmc-demo/).

`32:Tfe8,{ const height = 500; const width = 500; const colors = ["rgb(32,142,183)", "rgb(117,231,221)", "rgb(25,79,70)", "rgb(90,224,110)"]; const svg = d3.select(DOM.svg(width, height)).attr('id', 'svg-plot'); const margin = { top: 50, right: 50, bottom: 50, left: 50 }; svg.append("svg:defs").append("svg:marker") .attr("id", "triangle") .attr("refX", 6) .attr("refY", 6) .attr("markerWidth", 30) .attr("markerHeight", 30) .attr("orient", "auto") .append("path") .attr("d", "M 0 0 12 6 0 12 3 6") .style("fill", "gray"); //const xlim = d3.extent(trajectory.p.map((d) => d[0])); const xlim = [-4, 4]; const x = d3.scaleLinear() .domain(xlim) .range([margin.left, width - margin.right]) //const ylim = d3.extent(trajectory.theta.map((d) => d[1])); const ylim = [-4, 4]; const y = d3.scaleLinear() .domain(ylim) .range([height - margin.bottom, margin.top]); const xAxis = g => g .attr("transform", `translate(0,${height - margin.bottom})`) .call(d3.axisBottom(x).tickSizeOuter(0)) .call(g => g.append("text") .attr("x", width - margin.right) .attr("y", -10) .attr("fill", "#000") .attr("font-weight", "normal") .attr("font-style", "italic") .attr("text-anchor", "end") .attr("font-size", "14pt") .text("x₁")) const yAxis = g => g .attr("transform", `translate(${margin.left},0)`) .call(d3.axisLeft(y)) .call(g => g.select(".tick:last-of-type text").clone() .attr("x", 4) .attr("text-anchor", "start") .attr("font-weight", "normal") .attr("font-style", "italic") .attr("font-size", "14pt") .text("x₂")) let lineGenerator = d3.line() .x((d) => x(d[0])) .y((d) => y(d[1])); const color = d3.scaleLog().domain(d3.extent(density_grid)).range(["white", "steelblue"]); // Plot density contours svg.selectAll('path') .data(contours) .enter() .append('path') .attr('d', d3.geoPath(d3.geoIdentity().scale((x.range()[1] - x.range()[0]) / Math.sqrt(density_grid.length)))) .attr('transform', `translate(${margin.left}, ${height - margin.top}) scale(1, -1)`) .attr('fill', (d) => color(d.value)) .attr('opacity', '0.1'); // Plot current trajectory svg.selectAll('.trajectory-path') .data([state.trajectory.p]) .enter() .append('path') .attr('class', 'trajectory-path') .attr('d', (d) => lineGenerator(d)) .attr('fill', 'transparent') .attr('stroke', 'steelblue') svg.selectAll('.trajectory-step') .data(state.trajectory.p) .enter() .append('circle') .attr('class', 'trajectory-step') .attr('transform', (d) => `translate(${x(d[0])}, ${y(d[1])})`) .attr('r', 2) .attr('fill', 'black'); // Plot samples if(drawLines == "true") { svg.selectAll('.traceplot') .data([state.samples]) .enter() .append('path') .attr('class', 'traceplot') .attr('d', (d) => lineGenerator(d)) .attr('fill', 'transparent') .attr('stroke', 'steelblue') } svg.selectAll('.sample') .data(state.samples) .enter() .append('circle') .attr('class', 'sample') .attr('transform', (d) => `translate(${x(d[0])}, ${y(d[1])})`) .attr('r', 2) .attr('fill', 'black'); const init_drag = () => { const x1 = d3.event.sourceEvent.pageX - document.getElementById('svg-plot').getBoundingClientRect().x; const y1 = d3.event.sourceEvent.pageY - document.getElementById('svg-plot').getBoundingClientRect().y; mutable p_init = [x.invert(x1), y.invert(y1)]; }; svg.append('g').selectAll('.init') .data([p_init]) .enter() .append('circle') .attr('class', 'init') .attr('transform', (d) => `translate(${x(d[0])}, ${y(d[1])})`) .attr('r', 4) .attr('fill', 'red') .style('cursor', 'grab') .call(d3.drag() .on('drag', init_drag) ); svg.append("g") .call(xAxis); svg.append("g") .call(yAxis); return svg.node() }33:T7fd,{ const height = 250; const colors = ["rgb(32,142,183)", "rgb(117,231,221)", "rgb(25,79,70)", "rgb(90,224,110)"]; const svg = d3.select(DOM.svg(width, height)); const margin = { top: 15, right: 50, bottom: 50, left: 50 }; const xlim = [0, state.samples.length]; const x = d3.scaleLinear() .domain(xlim) .range([margin.left, width - margin.right]) const ymin = d3.min(state.samples.map((x) => d3.min(x))); const ymax = d3.max(state.samples.map((x) => d3.max(x))); const y = d3.scaleLinear() .domain([ymin, ymax]) .range([height - margin.bottom, margin.top]); const xAxis = g => g .attr("transform", `translate(0,${height - margin.bottom})`) .call(d3.axisBottom(x).tickSizeOuter(0)) .call(g => g.append("text") .attr("x", width - margin.right) .attr("y", -4) .attr("fill", "#000") .attr("font-weight", "bold") .attr("text-anchor", "end") .text("Iteration")) const yAxis = g => g .attr("transform", `translate(${margin.left},0)`) .call(d3.axisLeft(y)); let lineGenerator = d3.line() .x((d) => x(d.x)) .y((d) => y(d.y)); svg.selectAll('.traceplot') .data([state.samples.map((d, i) => ({ x: i, y: d[0], variable: "x1"})), state.samples.map((d, i) => ({ x: i, y: d[1], variable: "x2"}))]) .enter() .append('path') .attr('class', 'traceplot') .attr('d', (d) => lineGenerator(d)) .attr('fill', 'transparent') .attr('stroke', (d) => d[0].variable == "x1" ? 'steelblue' : '#82b446') svg.append('text') .attr('x', (d) => width - margin.right + 15) .attr('y', (d) => y(state.samples[state.samples.length - 1][0])) .attr('fill', 'steelblue') .text("x₁"); svg.append('text') .attr('x', (d) => width - margin.right + 15) .attr('y', (d) => y(state.samples[state.samples.length - 1][1])) .attr('fill', '#82b446') .text("x₂"); svg.append("g") .call(xAxis); svg.append("g") .call(yAxis); return svg.node() }34:T741,state = { if(reset) { let state = { samples: [p_init], trajectory: null, divergences: 0, acceptances: 0 }; const kinetic = (x) => math.sum(x.map((d) => Math.pow(d, 2))) / 2; for(let iter = 1; iter < iterations; iter++) { const momentum = [rnorm(0, 1), rnorm(0, 1)]; const integration = leapfrog(steps, epsilon, log_density, log_gradient, state.samples[iter - 1], momentum); state.trajectory = integration; state.divergences += (integration.divergent ? 1 : 0); if(animate !== false && animateLeapfrog !== false) { for(let step = 1; step < steps; step++) { const slicedState = ({ samples: state.samples, trajectory: { K: integration.K.slice(0, step), P: integration.P.slice(0, step), p: integration.p.slice(0, step), q: integration.q.slice(0, step), }, divergences: state.divergences, acceptances: state.acceptances }); if(leapfrogDelay == 0) { yield slicedState; } else { yield Promises.delay(leapfrogDelay, slicedState) } } } const H_L = integration.K[integration.K.length - 1] + integration.P[integration.P.length - 1]; const H_0 = math.multiply(-1, log_density(state.samples[iter - 1])) + kinetic(momentum); const ratio = Math.min(Math.exp(H_L - H_0), 1); const u = Math.random(); if(u < ratio) { state.samples[iter] = integration.p[integration.p.length - 1]; state.acceptances += 1; } else { state.samples[iter] = state.samples[iter - 1]; } if(animate !== false && animateLeapfrog === false) { yield state; } } yield state; } }35:[] 2f:{"title":"Hamiltonian Markov Chain Monte Carlo","cells":[{"id":0,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":52,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":157,"value":"md`\n\n**Diagnostics**\n\nNumber of divergent transitions: ${state.divergences}. *Divergences are when the leapfrog integrator fails to converge, which happens when step size is set too high to allow the integrator to adequately move across regions of high curvature. You can cause divergent transitions by setting the correlation to a high value (0.95 or -0.95) and then increasing the step size.*\n\nAcceptance rate: ${Math.round(state.acceptances / state.samples.length * 10000) / 100}%. *A low acceptance rate indicates that the trajectories have high numerical integration error, which causes the sampler to reject them. You can increase the acceptance rate by decreasing the integrator step size.*\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":49,"value":"viewof iterations = slider({ title: \"Number of iterations\", min: 2, max: 2000, step: 2, value: 100 })","pinned":false,"hidden":false,"mode":"ojs"},{"id":47,"value":"viewof rho = slider({ title: \"Correlation\", min: -0.95, max: 0.95, step: 0.05, value: 0, description: \"Correlation in the underlying bivariate normal density.\" })","pinned":false,"hidden":false,"mode":"ojs"},{"id":45,"value":"viewof variance = slider({ title: \"Variance\", min: 0.2, max: 5, step: 0.1, value: 1, description: \"Variances of the bivariate normal density (both variables given the same variance.)\" })","pinned":false,"hidden":false,"mode":"ojs"},{"id":57,"value":"viewof epsilon = slider({ title: \"Step size\", min: 0.01, max: 1, value: 0.2, description: \"How far the integrator moves at each step\" })","pinned":false,"hidden":false,"mode":"ojs"},{"id":60,"value":"viewof steps = slider({ title: \"Trajectory length\", min: 1, max: 250, value: 10, step: 1, description: \"Number of steps to take in the trajectory\" })","pinned":false,"hidden":false,"mode":"ojs"},{"id":43,"value":"viewof animate = checkbox({options: [{ value: 'true', label: 'Animate' }], value: 'true' })","pinned":false,"hidden":false,"mode":"ojs"},{"id":137,"value":"viewof animateLeapfrog = checkbox({options: [{ value: 'true', label: 'Animate leapfrog trajectories' }], value: 'true' })","pinned":false,"hidden":false,"mode":"ojs"},{"id":141,"value":"viewof leapfrogDelay = slider({ title: \"Leapfrog tweening delay (ms)\", min: 0, max: 100, value: 50, step: 1 })","pinned":false,"hidden":false,"mode":"ojs"},{"id":41,"value":"viewof drawLines = checkbox({options: [{ value: 'true', label: 'Draw lines' }], value: \"true\" })","pinned":false,"hidden":false,"mode":"ojs"},{"id":39,"value":"viewof reset = button(\"Rerun\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":37,"value":"md`\n## Traceplot\nA traceplot of the samples illustrates the random walk of the sampler.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":35,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":176,"value":"md`\n** Learn more **\n\n*Papers/Book chapters*\n- [Hybrid monte carlo (Duane et al. 1987)](https://www.sciencedirect.com/science/article/pii/037026938791197X)\n- [MCMC Using Hamiltonian Dynamics (Neal 2011)](http://www.mcmchandbook.net/HandbookChapter5.pdf)\n- [A Conceptual Introduction to Hamiltonian Monte Carlo (Betancourt 2017)](https://arxiv.org/abs/1701.02434)\n\n*Videos*\n- [The intuition behind the Hamiltonian Monte Carlo algorithm](https://www.youtube.com/watch?v=a-wydhEuAm0): YouTube video by Ben Lambert which walks through the intuition behind the algorithm and presents an example\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":98,"value":"md`\n### HMC algorithm\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":21,"value":"$34","pinned":true,"hidden":false,"mode":"ojs"},{"id":33,"value":"density_grid = {\n const grid = [];\n let k = 0;\n for(let x = -4; x < 4; x += 0.1) {\n for(let y = -4; y < 4; y += 0.1, ++k) {\n grid[k] = log_density([y, x]);\n }\n }\n \n return grid;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":31,"value":"thresholds = d3.scaleLinear().domain(d3.extent(density_grid)).ticks(50);","pinned":false,"hidden":false,"mode":"ojs"},{"id":29,"value":"contours = d3.contours()\n .size([Math.sqrt(density_grid.length), Math.sqrt(density_grid.length)])\n .thresholds(thresholds)\n (density_grid);","pinned":false,"hidden":false,"mode":"ojs"},{"id":27,"value":"log_density = multivariate_normal_log_density([0, 0], covariance)","pinned":false,"hidden":false,"mode":"ojs"},{"id":54,"value":"log_gradient = multivariate_normal_log_gradient([0, 0], covariance)","pinned":false,"hidden":false,"mode":"ojs"},{"id":25,"value":"covariance = [[variance, rho * variance], [rho * variance, variance]]","pinned":false,"hidden":false,"mode":"ojs"},{"id":23,"value":"mutable p_init = [-2, -1]","pinned":false,"hidden":false,"mode":"ojs"},{"id":77,"value":"leapfrog(steps, epsilon, log_density, log_gradient, p_init, [-1, -1])","pinned":false,"hidden":false,"mode":"ojs"},{"id":19,"value":"// slightly modified, from https://stackoverflow.com/questions/25582882/javascript-math-random-normal-distribution-gaussian-bell-curve\nfunction rnorm(mean, sd) {\n let u = 0, v = 0;\n while(u === 0) u = Math.random(); //Converting [0,1) to (0,1)\n while(v === 0) v = Math.random();\n return (Math.sqrt( -2.0 * Math.log( u ) ) * Math.cos( 2.0 * Math.PI * v ) * sd) + mean;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":192,"value":"","pinned":false,"hidden":false,"mode":"ojs"},{"id":4,"value":"d3 = require('d3@5')","pinned":false,"hidden":false,"mode":"ojs"},{"id":64,"value":"math = require('mathjs@5.4.2/dist/math.js')","pinned":false,"hidden":false,"mode":"ojs"},{"id":9,"value":"import { leapfrog } from \"@herbps10/leapfrog-integrator\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":2,"value":"import { slider, checkbox, radio, button } from \"@jashkenas/inputs\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":17,"value":"import {multivariate_normal_log_density, multivariate_normal_log_gradient} from \"@herbps10/distributions-and-their-gradients\"","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W35","stdlib":"1","titling":"auto"}