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;
}
}