Hamiltonian Monte Carlo Sampling Research Articles

Stan: A Probabilistic Programming Language.

Stan is a probabilistic programming language for specifying statistical models. A Stan program imperatively defines a log probability function over parameters conditioned on specified data and constants. As of version 2.14.0, Stan provides full Bayesian inference for continuous-variable models through Markov chain Monte Carlo methods such as the No-U-Turn sampler, an adaptive form of Hamiltonian Monte Carlo sampling. Penalized maximum likelihood estimates are calculated using optimization methods such as the limited memory Broyden-Fletcher-Goldfarb-Shanno algorithm. Stan is also a platform for computing log densities and their gradients and Hessians, which can be used in alternative algorithms such as variational Bayes, expectation propagation, and marginal inference using approximate integration. To this end, Stan is set up so that the densities, gradients, and Hessians, along with intermediate quantities of the algorithm such as acceptance probabilities, are easily accessible. Stan can be called from the command line using the cmdstan package, through R using the rstan package, and through Python using the pystan package. All three interfaces support sampling and optimization-based inference with diagnostics and posterior analysis. rstan and pystan also provide access to log probabilities, gradients, Hessians, parameter transforms, and specialized plotting.

Adaptive Hamiltonian MCMC sampling for robust visual tracking

Recent researches on visual tracking have shown significant improvement in accuracy by handling the large uncertainties induced by appearance variation and abrupt motion. Most studies concentrate on random walk based Markov chain Monte Carlo(MCMC) tracking methods which have shown inefficiency in sampling from complex and high-dimensional distributions. This paper proposes an adaptive Hamiltonian Monte Carlo sampling based tracking method within the Bayesian filtering framework. In order to suppress the random walk behavior in Gibbs sampling stage, the ordered over-relaxation method is used to draw the momentum item for the joint state variable. An adaptive step-size based scheme is used to simulate the Hamiltonian dynamics in order to reduce the simulation error and improve acceptance rate of the proposed samples. Furthermore, in designing the appearnce model, we introduce the locality sensitive histogram (LSH) to deal with appearance changes induced by illumination change. The proposed tracking method is compared with several state-of-the-art trackers using different quantitative measures: success rate and abruption capture rate. Extensive experimental results have shown its superiority to several other trackers.

Manual hierarchical clustering of regional geochemical data using a Bayesian finite mixture model

Interpretation of regional scale, multivariate geochemical data is aided by a statistical technique called “clustering.” We investigate a particular clustering procedure by applying it to geochemical data collected in the State of Colorado, United States of America. The clustering procedure partitions the field samples for the entire survey area into two clusters. The field samples in each cluster are partitioned again to create two subclusters, and so on. This manual procedure generates a hierarchy of clusters, and the different levels of the hierarchy show geochemical and geological processes occurring at different spatial scales. Although there are many different clustering methods, we use Bayesian finite mixture modeling with two probability distributions, which yields two clusters. The model parameters are estimated with Hamiltonian Monte Carlo sampling of the posterior probability density function, which usually has multiple modes. Each mode has its own set of model parameters; each set is checked to ensure that it is consistent both with the data and with independent geologic knowledge. The set of model parameters that is most consistent with the independent geologic knowledge is selected for detailed interpretation and partitioning of the field samples.

Hamiltonian Monte Carlo Sampling in Bayesian Empirical Likelihood Computation

SummaryWe consider Bayesian empirical likelihood estimation and develop an efficient Hamiltonian Monte Carlo method for sampling from the posterior distribution of the parameters of interest. The method proposed uses hitherto unknown properties of the gradient of the underlying log-empirical-likelihood function. We use results from convex analysis to show that these properties hold under minimal assumptions on the parameter space, prior density and the functions used in the estimating equations determining the empirical likelihood. Our method employs a finite number of estimating equations and observations but produces valid semiparametric inference for a large class of statistical models including mixed effects models, generalized linear models and hierarchical Bayes models. We overcome major challenges posed by complex, non-convex boundaries of the support routinely observed for empirical likelihood which prevent efficient implementation of traditional Markov chain Monte Carlo methods like random-walk Metropolis–Hastings sampling etc. with or without parallel tempering. A simulation study confirms that our method converges quickly and draws samples from the posterior support efficiently. We further illustrate its utility through an analysis of a discrete data set in small area estimation.

Fast Hamiltonian Monte Carlo Using GPU Computing

In recent years, the Hamiltonian Monte Carlo (HMC) algorithm has been found to work more efficiently compared to other popular Markov chain Monte Carlo (MCMC) methods (such as random walk Metropolis–Hastings) in generating samples from a high-dimensional probability distribution. HMC has proven more efficient in terms of mixing rates and effective sample size than previous MCMC techniques, but still may not be sufficiently fast for particularly large problems. The use of GPUs promises to push HMC even further greatly increasing the utility of the algorithm. By expressing the computationally intensive portions of HMC (the evaluations of the probability kernel and its gradient) in terms of linear or element-wise operations, HMC can be made highly amenable to the use of graphics processing units (GPUs). A multinomial regression example demonstrates the promise of GPU-based HMC sampling. Using GPU-based memory objects to perform the entire HMC simulation, most of the latency penalties associated with transferring data from main to GPU memory can be avoided. Thus, the proposed computational framework may appear conceptually very simple, but has the potential to be applied to a wide class of hierarchical models relying on HMC sampling. Models whose posterior density and corresponding gradients can be reduced to linear or element-wise operations are amenable to significant speed ups through the use of GPUs. Analyses of datasets that were previously intractable for fully Bayesian approaches due to the prohibitively high computational cost are now feasible using the proposed framework.

Non-Gaussian forecasts of weak lensing with and without priors

Assuming a Euclid-like weak lensing data set, we compare different methods of dealing with its inherent parameter degeneracies. Including priors into a data analysis can mask the information content of a given data set alone. However, since the information content of a data set is usually estimated with the Fisher matrix, priors are added in order to enforce an approximately Gaussian likelihood. Here, we compare priorless forecasts to more conventional forecasts that use priors. We find strongly non-Gaussian likelihoods for 2d-weak lensing if no priors are used, which we approximate with the DALI-expansion. Without priors, the Fisher matrix of the 2d-weak lensing likelihood includes unphysical values of $\Omega_m$ and $h$, since it does not capture the shape of the likelihood well. The Cramer-Rao inequality then does not need to apply. We find that DALI and Monte Carlo Markov Chains predict the presence of a dark energy with high significance, whereas a Fisher forecast of the same data set also allows decelerated expansion. We also find that a 2d-weak lensing analysis provides a sharp lower limit on the Hubble constant of $h > 0.4$, even if the equation of state of dark energy is jointly constrained by the data. This is not predicted by the Fisher matrix and usually masked in other works by a sharp prior on $h$. Additionally, we find that DALI estimates Figures of Merit in the presence of non-Gaussianities better than the Fisher matrix. We additionally demonstrate how DALI allows switching to a Hamiltonian Monte Carlo sampling of a highly curved likelihood with acceptance rates of $\approx 0.5$, an effective covering of the parameter space, and numerically effectively costless leapfrog steps. This shows how quick forecasts can be upgraded to accurate forecasts whenever needed. Results were gained with the public code from http://lnasellentin.github.io/DALI/

Estimation of realized stochastic volatility models using Hamiltonian Monte Carlo-Based methods

This study develops and compares performance of Hamiltonian Monte Carlo (HMC) and Riemann manifold Hamiltonian Monte Carlo (RMHMC) samplers with that of multi-move Metropolis-Hastings sampler to estimate stochastic volatility (SV) and realized SV models with asymmetry effect. In terms of inefficiency factor, empirical results show that the RMHMC sampler give the best performance for estimating parameters, followed by multi-move Metropolis-Hastings sampler. In particular, it is also shown that RMHMC sampler offers a greater advantage in the mixing property of latent volatility chains and in the computational time than HMC sampler.

Bayesian large-scale structure inference: initial conditions and the cosmic web

AbstractWe describe an innovative statistical approach for theab initiosimultaneous analysis of the formation history and morphology of the large-scale structure of the inhomogeneous Universe. Our algorithm explores the joint posterior distribution of the many millions of parameters involved via efficient Hamiltonian Markov Chain Monte Carlo sampling. We describe its application to the Sloan Digital Sky Survey data release 7 and an additional non-linear filtering step. We illustrate the use of our findings for cosmic web analysis: identification of structures via tidal shear analysis and inference of dark matter voids.

Riemann Manifold Langevin and Hamiltonian Monte Carlo Methods

SummaryThe paper proposes Metropolis adjusted Langevin and Hamiltonian Monte Carlo sampling methods defined on the Riemann manifold to resolve the shortcomings of existing Monte Carlo algorithms when sampling from target densities that may be high dimensional and exhibit strong correlations. The methods provide fully automated adaptation mechanisms that circumvent the costly pilot runs that are required to tune proposal densities for Metropolis–Hastings or indeed Hamiltonian Monte Carlo and Metropolis adjusted Langevin algorithms. This allows for highly efficient sampling even in very high dimensions where different scalings may be required for the transient and stationary phases of the Markov chain. The methodology proposed exploits the Riemann geometry of the parameter space of statistical models and thus automatically adapts to the local structure when simulating paths across this manifold, providing highly efficient convergence and exploration of the target density. The performance of these Riemann manifold Monte Carlo methods is rigorously assessed by performing inference on logistic regression models, log-Gaussian Cox point processes, stochastic volatility models and Bayesian estimation of dynamic systems described by non-linear differential equations. Substantial improvements in the time-normalized effective sample size are reported when compared with alternative sampling approaches. MATLAB code that is available from http://www.ucl.ac.uk/statistics/research/rmhmc allows replication of all the results reported.

LOCAL NON-GAUSSIANITY IN THE COSMIC MICROWAVE BACKGROUND THE BAYESIAN WAY

We introduce an exact Bayesian approach to search for non-Gaussianity of local type in Cosmic Microwave Background (CMB) radiation data. Using simulated CMB temperature maps, the newly developed technique is compared against the conventional frequentist bispectrum estimator. Starting from the joint probability distribution, we obtain analytic expressions for the conditional probabilities of the primordial perturbations given the data, and for the level of non-Gaussianity, f_nl, given the data and the perturbations. We propose Hamiltonian Monte Carlo sampling as a means to derive realizations of the primordial fluctuations from which we in turn sample f_nl. Although being computationally expensive, this approach allows us to exactly construct the full target posterior probability distribution. When compared to the frequentist estimator, applying the Bayesian method to Gaussian CMB maps provides consistent results. For the analysis of non-Gaussian maps, however, the error bars on f_nl do not show excess variance within the Bayesian framework. This finding is of particular relevance in the light of upcoming high precision CMB measurements obtained by the Planck satellite mission.

On the Continuity of the Generalized Inverse

On the Continuity of the Generalized Inverse