  Hamiltonian Monte Carlo (HMC) is a powerful tool for Bayesian statistical
inference due to its potential to rapidly explore high dimensional state space,
avoiding the random walk behavior typical of many Markov Chain Monte Carlo
samplers. The proper choice of the integrator of the Hamiltonian dynamics is
key to the efficiency of HMC. It is becoming increasingly clear that
multi-stage splitting integrators are a good alternative to the Verlet method,
traditionally used in HMC. Here we propose a principled way of finding optimal,
problem-specific integration schemes (in terms of the best conservation of
energy for harmonic forces/Gaussian targets) within the families of 2- and
3-stage splitting integrators. The method, which we call Adaptive Integration
Approach for statistics, or s-AIA, uses a multivariate Gaussian model and
simulation data obtained at the HMC burn-in stage to identify a system-specific
dimensional stability interval and assigns the most appropriate 2-/3-stage
integrator for any user-chosen simulation step size within that interval. s-AIA
has been implemented in the in-house software package HaiCS without introducing
computational overheads in the simulations. The efficiency of the s-AIA
integrators and their impact on the HMC accuracy, sampling performance and
convergence are discussed in comparison with known fixed-parameter multi-stage
splitting integrators (including Verlet). Numerical experiments on well-known
statistical models show that the adaptive schemes reach the best possible
performance within the family of 2-, 3-stage splitting schemes.
