Hamiltonian Monte Carlo — Hmcsam • rmorie

The state is augmented with a momentum and the leapfrog integrator is run L times before a Metropolis accept/reject on exp(-Delta H). Leapfrog is exactly reversible and volume preserving, which is what makes the acceptance ratio simply exp(-Delta H); both properties are checked directly in the tests. Momenta and acceptance uniforms come from deterministic low-discrepancy streams so the chain is reproducible.

Usage

Hmcsam(log_p, grad_log_p, x0, step_size = 0.1, L = 10, n_iter = 200)

Arguments

log_p: Log target density of a numeric vector.
grad_log_p: Its gradient.
x0: Starting state.
step_size: Leapfrog step size.
L: Leapfrog steps per proposal.
n_iter: Number of proposals.

Value

List with estimate, mean, draws, accept_rate, mean_energy_error, n, method.

Details

Formula: p <- p + (eps/2) grad log p(x); x <- x + eps p; p <- p + (eps/2) grad log p(x); accept with exp(-Delta H).

References

Neal (2011), MCMC using Hamiltonian dynamics, in Brooks, Gelman, Jones and Meng (eds), Handbook of Markov Chain Monte Carlo, CRC Press, ch. 5. doi:10.1201/b10905

Examples

set.seed(1)
Hmcsam(function(x) -0.5 * sum(x^2), function(x) -x, x0 = c(0, 0),
       n_iter = 50)
#> $estimate
#> [1] -0.09463431
#> 
#> $mean
#> [1] -0.094634313 -0.009193392
#> 
#> $draws
#>                [,1]         [,2]
#>  [1,] -5.684265e-01  0.568426499
#>  [2,] -3.838915e-02  0.038389151
#>  [3,] -1.313609e+00  0.153301911
#>  [4,]  3.835831e-02 -0.664867549
#>  [5,] -1.118620e-01  0.933884620
#>  [6,]  5.682052e-03  0.016188314
#>  [7,] -8.481016e-01  0.347737299
#>  [8,]  6.528207e-01 -0.922993006
#>  [9,]  1.349495e-02  0.352798379
#> [10,]  4.953505e-01  0.124411835
#> [11,] -1.547714e+00  0.100191987
#> [12,] -2.251051e-01 -0.853977874
#> [13,] -3.555695e-01  0.575357526
#> [14,] -2.588488e-02 -0.064440254
#> [15,] -7.138899e-01  0.414736642
#> [16,]  1.026933e+00 -1.188461284
#> [17,]  1.049626e-01  0.058202208
#> [18,]  4.317799e-01 -0.134678975
#> [19,] -8.033238e-01  0.161303645
#> [20,]  4.743210e-01 -0.523488673
#> [21,]  2.230945e-01  1.532520371
#> [22,]  1.369648e-01  0.279610833
#> [23,] -8.641906e-01  0.436931207
#> [24,]  5.355479e-01 -1.005843822
#> [25,] -1.042473e-01  0.229321858
#> [26,]  3.743282e-01  0.007933151
#> [27,] -1.282529e+00  0.086933859
#> [28,] -1.562939e-02 -0.777183489
#> [29,] -1.914061e-01  0.731727600
#> [30,]  1.135685e-01  0.073920463
#> [31,] -5.709228e-01  0.547627957
#> [32,]  1.366635e+00 -1.379213222
#> [33,]  2.302018e-01 -0.112463606
#> [34,]  4.454745e-01 -0.277643359
#> [35,] -9.108343e-01  0.033053085
#> [36,]  3.323174e-01 -0.659026937
#> [37,]  9.678487e-02  1.128806004
#> [38,]  1.681486e-01  0.178883582
#> [39,] -6.816364e-01  0.490005467
#> [40,]  8.737149e-01 -0.737589653
#> [41,]  1.858085e-01  0.539882454
#> [42,]  6.482030e-01  0.275005542
#> [43,] -1.891775e+00  0.156741523
#> [44,] -4.426368e-01 -0.869003803
#> [45,] -4.988683e-01  0.516407742
#> [46,] -1.284297e-01 -0.123807523
#> [47,] -8.048359e-01  0.354390618
#> [48,]  8.857553e-01 -1.333933729
#> [49,] -7.313743e-05 -0.054716854
#> [50,]  3.479234e-01 -0.220969338
#> 
#> $accept_rate
#> [1] 1
#> 
#> $mean_energy_error
#> [1] 4.246953e-06
#> 
#> $n
#> [1] 50
#> 
#> $method
#> [1] "leapfrog integration of H = U + K with a Metropolis correction, Neal (2011) ch. 5"
#> 
#> attr(,"class")
#> [1] "morie_rich_result" "list"