I found that this package is significant slower than the hmc in matlab. The reason is that the interface combines the function and its derivative together, i.e., the function needs to return both loglikelihood and its derivative. But in fact, hmc only need to evaluate the loglikelihood twice during its sampling.

Hi Robert,

I found that in line 109 of the _hmc.pyx, there is an additional 0.5 scale in updating p (the full ), while as in standard HMC method (Neal's MCMC using Hamiltonian dynamics) there is no such term. May I ask if I have mis-understood anything?

Thanks,
Xuhui

1. Sokal's notes (page 16), citing an earlier paper use an procedure for choosing the window for summing the empirical ACF where the length of the window, M is chosen to be the smallest value such that M is at least c times the estimated autocorrelation time, where c is something like 4, 6, or 10.
2. Jenke also uses the same scheme, with c=6.
3. Evertz doesn't seem to like this procedure (page 59)
4. LaplacesDemon's in R uses the Geyer IPS estimator.
5. Geyer's ICS estimator.

I'm trying to sample a beta distribution with a = n + 44 and b = s + 102, where n and s are constants, scale constant C = 0.01 and starting value p = 0.2 for 1000 iterations.

My beta looks like

g(p) = Cp^(a-1)*(1-p)^(b-1).
log(g) = log(C) + (a-1)log(p) + (b-1)log(1-p)
d(log(g)) = 0 + (a-1)/p - (b-1)/(1-p)

import pyhmc

sample_size = 31
sum_spacings = 43
starting_value = 0.2
num_iterations = 1000
constant_c = 0.01

def log_posterior_grad(p, n, s):
    log_posterior_beta = numpy.log(constant_c) + (n + 44 - 1)*numpy.log(p) + (s + 102 -1 )*numpy.log(1-p)
    gradient = (n + 44 - 1)/p - (s + 102 -1 )/(1-p)
    return log_posterior_beta, gradient

samples = pyhmc.hmc(fun = log_posterior_grad, x0 = [starting_value ], args = (sample_size, sum_spacings, ), n_steps = 1000)

The dependencies "statsmodels" and "cython" should be added to the install_requires in the setup.py script.

In the readme and the docs, the number of samples specified in the simple example is written as n_samples=1e4. With a recent numpy update (not sure which) it won't automatically convert a float to an int when you get to the line in hmc.py: samples = np.zeros((n_samples, n_params)). Thus, you get an error if you run the sample script as-is. If you change it it to n_samples=10000 it works fine.