1function [O, logW] = extend_orbit_backward(logmu, gradlogmu, theta_a, rho_a, logw_a, h, delta, L)
2 O = {};
3 logW = [];
4 theta = theta_a;
5 rho = rho_a;
6 logw = logw_a;
8 for iter = 1:L
9 ell_b = micro(logmu, gradlogmu, theta, -rho, h, delta);
10 ell = p_micro(ell_b);
11 [theta_1, rho_1] = leapfrog(logmu, gradlogmu, theta, -rho, h * (2^(-ell)), 2^ell);
12 rho_1 = -rho_1;
13 ell_f = micro(logmu, gradlogmu, theta_1, rho_1, h, delta);
15 if pmf_p_micro(ell, ell_f) == 0
16 logw = -Inf;
17 else
18 logw = logw + logmu(theta_1) - 0.5 * norm(rho_1)^2 - logmu(theta) + 0.5 * norm(rho)^2 + log(pmf_p_micro(ell, ell_f)) - log(pmf_p_micro(ell, ell_b));
19 end
21 theta = theta_1;
22 rho = rho_1;
23 O = [{theta, rho}; O];
24 logW = [logw; logW];
25 end
26end