Many machine learning methods employ Bayesian analyses where the posterior probability densities are often analytically intractable. Markov Chain Monte Carlo (MCMC) is a technique admitting sampling from most posterior densities. In theory, MCMC samplers can provide unbiased sampling of posterior densities. In practice, however, Markov chains may require a rather long time to converge to the target distribution or even trapped in a local mode.
This problem is especially acute for high dimensional multimodal target distributions which little knowledge of the target sampling landscape is available. MCMC samples could be drawn from just one of the modes. In high dimensional distributions that are multimodal, it is likely that no samples are drawn from major modes. This could lead to significant bias in the estimation of the mean and variance of the target distribution and results in false confidence on the quality of the inference. One may think that running multiple chains could alleviate this problem. For high dimensions, though, this may not be effective(Neal, 2012)
. The probability density function could be rather at in some subregions in a high dimensional space. This results in rather slow convergence and causes the strategy of running multiple chains almost computationally infeasible.
Recent efforts have been made to sample from multimodal distributions, such as Regenerative Darting Monte Carlo (RDMC, Ahn et. al. (2013)) and Wormhole Hamiltonian Monte Carlo (WHMC, Lan et. al. (2013)). Both methods update the sampler to include new modes as they are found during regeneration times, when the Markov chains probabilistically restarts. Markov chains probabilistically jump to other modes when they enter predefined regions determined by known modes. During regeneration times, a local search method, such as BFGS, is used to search for new modes and redefine the jump regions. It has been shown in the literature that WHMC converges rather fast if detailed and accurate knowledge about modes are provided. This requirement, however could be unrealistic in high dimensions such as 100 dimensions or higher due to the cost in discovering new modes.
In order to avoid rediscovering known modes in WHMC, Lan et. al. (2013) employed the so called residual potential in order to diminish the basins of attraction of known modes so that BFGS is more likely to find unknown modes when the residual potential is optimized. Even though the method has shown to be successful in specific cases, such as Gaussian mixture models, there are a few issues to this approach. First, a major drawback is that there is a possibility of fictitious modes as their method directly replaces the target distribution with the residual potential to be optimized. In fact, during our experiments we discovered that fictitious modes are indeed introduced, unstabilizing WHMC. Second, the starting points are still chosen randomly in the search space for the BFGS runs. This allows BFGS to rediscover already known modes and wasting computational resources. Another major drawback is that the independent sampler employed within WHMC requires the proposal distribution to be similar to the target distribution. Otherwise, the regeneration probability would be very small and the mode search practically never starts. This creates an awkward scenerio where one requires some knowledge of the target distribution for regeneration to be effective when no such knowledge is available apriori. Lastly, WHMC samples from the target distribution “on-the-fly”, where mode searching is conducted during sampling. This introduces an additional source of bias before major modes are found.
At first glance, modifying the objective function to subtract out contribution from known modes seems trivial. However, such a subtraction scheme tends to create regions with vanishing gradient if done perfectly, slowing the convergence of any gradient based search algorithms. In differential geometry, a geodesic is a generalization of a straight line on Riemmanian manifolds. The vanishing gradient problem can be avoided if one finds a manifold on which geodesics correspond to straight lines in the at regions of the search space and where the geodesics would pass through the maxima of the objective function. Recently,Fok et. al. (2016) showed that geodesics on manifolds conformally related to Euclidean spaces possesses such properties and a single geodesic can travel through multiple maxima. Such geodesics appear to be an excellent candidate for finding major modes.
We propose a novel mode searching method that outperforms the residual potential method by addressing the drawbacks experienced by WHMC. Our method first constructs a conformal geodesic that is able to visit the basins of many local maxima111The same trajectory was used as a component of a global optimization algorithm in Fok et. al. (2016). The novel contribution is our mode finding method as opposed to that with a residual potential.. The trajectory travels on the residual log objective function , where the modes have been subtracted out from the log objective function . This method addresses the first two concerns mentioned above with WHMC. First, a starting point is chosen with a bias towards undiscovered modes so that undiscovered modes are more likely to be found. Second, because the original function is maximized, there will be no fictitious maxima introduced. Our experiments show that our approach requires fewer BFGS calls to for mode searching and was stable in the experiment where WHMC was not able to converge. Lastly, we show theoretically that WHMC can be made more effcient by having fewer samples during the early runs where not many major modes are found. Therefore we simply force our sampler to search for new modes after a number of samples have been acquired so as to discover new modes as fast as possible in order to reduce the additional source of bias introduced by on-the- y sampling.
This article is organized as follows. Section 2 provides a brief review of MCMC literature. Our algorithm is proposed in Section 3. Numerical results are provided in Section 4. Section 5 concludes the paper with discussion and plan for future work.
2 Related Work
The first MCMC sampler developed is the Metropolis-Hastings algorithm (MH) (Metropolis et. al., 1953; Hastings, 1970). The Gibbs sampler was proposed more recently (Gelfand and Smith, 1990). It relies on a user-specified proposal function, one that is easy to sample from, to generate candidate points. Whether the candidate points are accepted is based on the acceptance probability. The Metroplis-Hastings algorithm is rather inefficient in today’s standards, however. First, its computational time scales with dimensionality as (Creutz, 1988) which makes it rather inefficient in high dimensions. As a consequence, the optimal acceptance rate is just 0.23 for random walk proposal functions222For a derivaton of the optimal acceptance rate, see Section 4 of (Neal, 2012) . It also mixes poorly in colinear regions of the target distribution.
A better alternative is the so called Hamiltonian Monte Carlo (HMC) (Neal, 2012). The candidate points are proposed by solving the Hamiltonian equations of motion under the potential energy , where is the target distribution and the kinetic energy , where is the momentum and is a constant real and symmetric mass matrix, usually set to be the identity. First, the momentum is sampled. Then the Hamiltonian equations of motion is solved using the leapfrog integrator and a candidate point is chosen given the trajectory length and the leapfrog step size. HMC is much more efficient than MH with optimal acceptance rate at 0.65. Also, its complexity scales with dimensions as which is less costly than MH in high dimensions. HMC has a few drawbacks. The leapfrog step size must be randomized to avoid periodic behaviors leading to poor mixing. Furthermore, the trajectory length and leapfrog step sizes need to be tuned from a prelimineary HMC run to obtain an optimal acceptance rate.
In recent years, MCMC samplers outperforming HMC has been proposed. One example is the Riemannian Manifold Hamiltonian Monte Carlo (RMHMC) (Girolami and Calderhead, 2011)
. In this algorithm the trajectories follow the Hamiltonian dynamics on a manifold defined by the Fisher-Rao metric tensor(Rao, 1945). In general, the Hamiltonian contains terms depending on and in an inseparable way. As a result the equations of motion are implicit and the solution requires a numerical intergrator such as fixed point iteration, which may be numerically unstable and computationally costly in some cases. To alleviate this, Riemannian MCMC with Lagrangian dynamics (LMC) (Lan et. al., 2012) was proposed. LMC uses a semi-explicit integrator, in which one of the updates is explicit.
Attempts have been made to tune HMC parameters adaptively. The No-U-Turn Sampler (NUTS) (Hoffman and Gelman, 2014) uses a recusrive algorithm to propose candidate points in a wide region of the target distribution and stops the trajectory as it makes a U-turn. The step sizes can also be adjusted adaptively in NUTS. The Adaptive Hamiltonian and Riemann Monte Carlo (AHRMC) sampler Wang et. al. (2013) uses a Bayesian optimization method to tune HMC parameters.
There exists samplers that specializes in sampling from multimodal distributions. Darting Monte Carlo (DMC) (Andricioaiei et. al., 2001; Sminchisescu and Welling, 2007) defines regions around locations of high posterior density and attempts to jump to another mode when entered. RDMC (Ahn et. al., 2013), a variant of DMC, was developed so that new regions can be defined when the Markov chains reaches their regeneration times. Rather than using the Fisher Information metric, Wormhole Hamiltonian Monte Carlo (WHMC) (Lan et. al., 2013) designs wormhole metric tensors in which the distances between modes on the corresponding manifold are reduced. Mode searches are performed with BFGS at regeneration times and the wormhole metric tensors are updated.
We show formally the source of bias introduced by on-the-fly sampling. Suppose we are given a set of samples , the probability density can be estimated as
where is a sample and is the kernel. Now, suppose that the target distribution is multimodal with modes and modes are known. Assume that is obtained on-the-fly, we can partition the samples into two sets , where denotes samples before all the modes are found and denotes samples obtained after the discovery of all modes. Similarly, and denotes the density estimate obtained from and , respectively and
their corresponding sample sizes. Since WHMC converges to a target distribution of known modes k, the true mean and variances (and higher statistical moments) are obtained fromas .
Without the loss of generality, consider an estimation of the mean using ,
Partitioning as above by splitting the sum in the function estimation into two corresponding terms, we get
The above expression has two terms. The first one is the bias introduced by on-the-fly sampling. In the limit , the first term above is negligible compared to the second and the estimated mean converges to the true mean. Another way to reduce the on-the-fly bias is to find all the major modes as soon as possible so that . This can be seen in the asymptotic limit () where the sum can be replaced by an integral . Note that in the asymptotic limit is in both and so the distinction between each set can be dropped. Finally, set and the denominator is canceled by the corresponding term in the numerator. One obtains and estimate of the true mean
The above steps only involves manipulation of the density estimate and therefore is valid for all higher statistical moments. From the above consideration one sees that the on-the-fly bias can be reduced by having an efficient mode finding algorithm such that major modes first so that , and that limiting the number of samples before all modes are found can speed up convergence. Even though it may be impractical to find all modes before sampling, we found that the proposed algorithm, when applied on-the-fly, still outperforms WHMC.
We perform comparisons between the tempered residual potential energy (TRPE) method of WHMC and SGEO-KDE using target distributions provided by Lan et. al. (2013). To justify the discussion in Section 3, the density around a mode is estimated by a Gaussian distribution with mean and covariance being the maximizer and inverse Hessian matrix at the mode in the first experiment involving Gaussian mixtures, so that the estimated density for mode k is very close to a Gaussian term in the target distribution. In the second experiment with the sensor network localization problem, KDEs were used. We found that SGEO-KDE converges much faster than TRPE (WHMC) in both experiments. Further, we found no ficititious modes using SGEO-KDE whereas the residual potential method ficititious modes are present which caused the WHMC sampler to be unstable.
The following diagnostic quantities are used to assess convergence, the relative error of mean (REM),
where denotes the components, denotes the MCMC estimate of the mean at time , and is the mean of the target distribution. Similarly, the relative error of covariance (RECOV) is defined as
where indices and donotes the components of the covariance matrix and its MCMC estimate .
Since BFGS is the computational bottleneck, especially in high dimensions. We show the average number of BFGS attempts as a measure of computation cost.
4.1 Mixture of Gaussians
We test SGEO-KDE on Gaussian mixture models used by Lan et. al. (2013) in this section, the estimates are obtained from 4 MCMC chains run in parallel. The HMC acceptance rate is tuned to be near optimal, in the range of 0.6 to 0.7. The wormhole in uence factor is set to be F = 0.1.
In addition to the convergence diagnostics, we also compute the average number of BFGS optimization needed to locate all the modes for the SGEO optimization algorithm and the method of tempered residual potential energy employed in WHMC. This quantity can be used as a measure for computational cost as BFGS optimization is the bottleneck for both algorithms.
Results for SGEO-KDE and TRPE (WHMC) for Gaussian mixture models with equal weights are shown in Table 1. The corresponding results for Guassian mixtures with unequal weights
where , are shown in Table 2. The results are taken after running the samplers for 800 seconds, and 2000 seconds for the = 100 case. Plots of convergence diagnostics over computational time are shown in Figure 1 and Figure 2 for the unequal weight and equal weight cases, respectively. From the results we see that SGEO-KDE outperforms TRPE (WHMC) in all diagnostic measures, with the exception of three cases in RECOV. In addition, the SGEO mode searching algorithm is shown to be superb in finding new modes, with only one failure out of a total 80 attempts in Tables 1 and 2. Furthermore, TRPE (WHMC) requires at least four times the number of BFGS attempts than SGEO-KDE. This shows that SGEO-KDE outperforms TRPE (WHMC) in terms of computational cost as the bottleneck is BFGS. This is confirmed by conisdering the REM measure from the most recent 1000 samples (third column of Figures 1 and 2). SGEO-KDE takes a much shorter time to converge to the true mean. The improvement is amplified as the number of dimensions increase. In the = 100 case, TRPE (WHMC) was not able to locate all the modes in 2000 seconds whereas the proposed method where all the modes were discovered at about 500 seconds. Furthermore, in higher dimensions (), the margin of improvement is larger when the Gaussian weights become unequal.
The results from HMC are included in the appendix. The HMC sampler is not able to sample from all the modes and the convergence is significantly worse than TRPE (WHMC) and SGEO-KDE.
|SGEO-KDE (all modes)||0.018371||0.95372||10|
|SGEO-KDE (on the fly)||0.083307||0.94318||10|
|SGEO-KDE (all modes)||0.030689||0.98234||10|
|SGEO-KDE (on the fly)||0.032558||0.98191||10|
|SGEO-KDE (all modes)||0.021922||0.99505||10|
|SGEO-KDE (on the fly)||0.17626||1.8093||10|
|SGEO-KDE (all modes)||0.024986||0.99685||10|
|SGEO-KDE (on the fly)||0.12135||2.9112||10|
|SGEO-KDE (all modes)||0.17218||0.95234||10.05|
|SGEO-KDE (on the fly)||0.14069||0.94555||10.05|
|SGEO-KDE (all modes)||0.14588||0.98155||10|
|SGEO-KDE (on the fly)||0.10588||0.9638||10|
|SGEO-KDE (all modes)||0.11964||0.99498||10|
|SGEO-KDE (on the fly)||0.1011||1.2841||10|
|SGEO-KDE (all modes)||0.07165||0.99667||10|
|SGEO-KDE (on the fly)||0.078958||2.5587||10|
4.2 Sensor Network Localization
Given a set of sensors scattered on a two dimensional plane, it is possible to infer the locations of each sensor given the distance measured between each sensors. The domain of the posterior distribution is a dimensional configuration space comprises of the locations of each sensor , where each
is a 2-dimensional coordinate vector denoting the position of sensor. The solution for the locations can be found by a maximum a- posterior (MAP) configuration estimate, and its uncentainties estimated by MCMC samples of the posterior. Due to measurement error, the posterior distribution can be multimodal, which each mode denoting a possible set of locations for the Ns sensors that maximizes the posterior locally. Obviously, the modes are not known aprori and hence this problem is an excellent application for SGEO-KDE. We run two chains in parallel for this experiment.
In general, each mode of the posterior is non-Gaussian. After finding a new mode, HMC samples are taken and used to obtain a Gaussian kernel estimate of the posterior around the new mode.Then the KDE estimate of each mode is subtracted out for SGEO. Obtaining the true mean and covariance matrix of the target distribution is intractible. Thus we perform a long SGEO-KDE run to obtain an estimate of the “true” mean and covariance of the target distribution.
If the prior distribution is taken to be uniform, the target distribution has the form given by (Ihler et. al., 2005)
where denotes whether sensors and detect each other. If so, denotes the distance measured between the sensors. The total number of unknown sensors is . In accordance to the literature, we choose and . The wormhole influence factor is , with the probability of jumping between modes being . The HMC parameters are set to be , and the mass matrix . The HMC acceptance rate is 0.60.
The results are shown in Table 3 and Figure 3. In Figure 3, it is shown that WHMC is not able to converge even after 1000 seconds in the plot on the right. For SGEO-KDE, convergence is achieved very quickly. The reason is that the tempered residual potential method finds fictitious modes that destabilize WHMC. To show this, we investigated whether another BFGS run starting from the recorded locations of the modes would lead to another point. Let be the -th recorded mode ( being the first mode discovered and denotes the last mode) from either the residual potential method or SGEO-KDE and be a local maximum of the target distribtion obtained by BFGS using as a starting point. Then we consider the displacement between and up to the -th mode as a measure for fictitious modes, that is
Figure 4 shows the cumulative displacement, for the residual potential method and SGEO-KDE over five runs. The results show that significant error tends to develop for residual potential after some modes have been discovered. For SGEO-KDE this trend is not observed and the cumulative error remains small as new modes are being discovere. This is consistent with the fact that the residual potential introduces ficititious modes as genuine ones are being found.
|SGEO-KDE (all modes)||7.01||0.14|
|SGEO-KDE (on the fly)||17.8||4.45|
5 Discussion and Conclusions
In principle, the mode searching algorithm described in this article can be applied to many MCMC methods without modification. We have also tried our algorithms on a Gaussian mixture model with 500 dimensions and obtained analogous results for mode searching. The next challenge is to extend our algorithm to much higher dimensions than described in this article.
Appendix A Review of MCMC algorithms
The HMC originates from Hamiltonian dynamics. The solution of the Hamiltonian equations of motion describes the trajectories of particles under a potential energy and the kinetic energy , where is the position and is the momentum of the particle. In HMC, the potential is chosen to be the negated log-density, i.e. , where is the target density. The kinetic energy is , where is a constant mass matrix, usually set to be the identity. Given a trajectory length , each candidate point is chosen as the last position of its trajectory.
Let be the Hamiltonian. The Hamiltonian equations of motion are
The above can be solved numerically by using the leapfrog integrator with step size parameter
Let a state be . Analogous to MH, the acceptance probability , where , , denotes the probability of being in state and denotes the transition kernel from state to . Since Hamiltonian dynamics is time reversible, we have . Then the acceptance probability is determined by the probability of being in states and . From statistical mechanics, the probability of finding a state is333We set the temperature parameter here. The distribution is also called the Boltzmann distribution used in Boltzmann machines.
here. The distribution is also called the Boltzmann distribution used in Boltzmann machines.
where is a normalization constant. The acceptance probability for HMC is
The algorithm for HMC is shown in A.1. The two input parameters are the leapfrog step size and the total number of leapfrog steps . In practice, these parameters need to be tuned so that the acceptance probability is close to the optimal value of 0.65.
Given input , generate
Obtain () with Leapfrog steps of size for steps.
Accept with probability , discard .
Despite the efficiency of HMC for target distribution with a single mode, the HMC trajectories can get trapped in only one of the modes when multiple are present. This causes biased sampling where only one of the modes are being sampled. Many methods exists to alleviate this problem. In the next subsection we review the state-of-the-art method for multimodal sampling using HMC, the Wormhole Monte Carlo (WHMC) sampler. Tables 4 and 5 shows the results with HMC on Gaussian mixture models.
When multiple modes exist in the target distribution, the samples and HMC trajectories can be confined within a local maximum. Given the locations of the modes, the WHMC employs the HMC sampler extended onto a Riemannian manifold, endowed with a metric that shortens the intermodal distances. Following Girolami and Calderhead (2011), the constant mass matrix in the kinetic energy is replaced by the positional dependent metric, . An immediate consequence of this is that, as we shall see, at least one of the leapfrog update equations become implicit and the updates require fixed point iterations.
Consider we have two known modes at and . The metric in WHMC comprises of two terms,
where is the wormhole metric defined by
with being a unit vector point from to and is a small and positive number . The mollifying function is defined by
where is the wormhole influence factor.
Notice that, using the triangle inequality,
if and only if is on the line connecting the two modes and that the wormhole influence factor controls the extent the wormhole metric influences the metric away from the intermodal line. Also, the wormhole metric shortens the distance between the two modes. Let , where so that in Euclidean space. Under , the 2-norm of is
That is to say that the wormhole metric shortens the distance between two modes by a factor of . For more than two modes, a network of wormholes connecting each mode is built.
To facilitate moving between modes, Lan et. al. (2013) introduced a vector field, , along the intermodel line, with norm proportional to the mollifying function and the projection of the velocity444The definition of the velocity in classical mechanics is and we are identifying the mass matrix with the metric here. onto the intermodal line
with the corresponding modification of the Hamilton equation
Lastly, Lan et. al. (2013) also introduced an auxiliary dimension to prevent interference between the wormholes and the HMC dynamics.
The generalized leapfrog integrator for WHMC is
We can see that the update equation for is implicit and fixed point iteration is needed for this update.
a.3 Updating Target Distribution During Regeneration Time
Regeneration times are periods where a Markov chain restarts itself probabilistically. States after regeneration is independent of states prior to the regeneration time. In thoery, another set of sampler parameters can be used for the sampler while leaving the MCMC estimates unaffected. In addition, Lan et. al. (2013) and Ahn et. al. (2013) performs BFGS at regeneration times to search for new modes and updates the proposal distribution on-the-fly.
Regeneration is checked at each step. The regeneartion probability is
is the proposal indepdence distribution. Note that, if and only if . Regeneration probability is dependent on the choice of the independence proposal function, and the choice of the constant . If the choice of the indenpendence sampler is poor, such as in the case where only a minority of the modes in a multimodal target distribution are known, then the regeneration probability can be too small that BFGS is never triggered to find new modes. We encountered such a scenario when testing the sensor network localization problem with WHMC.
Another problem with updating the target distribution on-the-fly is that sample estimates are biased by the samples from the previous target distribution. In the following sections, we show that the mean and covariance estimates of the target distribution are significantly more accurate if all modes of the target distribution are found first. To this end we employ a global optimization algorithm SGEO modified to discover all the modes before sampling.
The research of second and third authors are supported by National Science and Engineering Research Grants. We thank Shiwei Lan for providing the code for WHMC.
- Ahn et. al. (2013) S. Ahn, Y. Chen, and M. Welling. Distributed and adaptive darting Monte Carlo through regenerations. Journal of Machine Learning Research 31: W&CP, 31, 2013
- Andricioaiei et. al. (2001) I. Andricioaiei, J. Straub, and A. Voter. Smart darting Monte Carlo. 114(16), 2001.
- Brockwell and Kadone (2005) A.E. Brockwell and J.B. Kadane. Identification of regeneration times in mcmc simulation, with application to adaptive schemes. Journal of Computational and Graphical Statistics, 14(2):436–458, 2005.
- Brooks and Gelman (1998) S. P. Brooks and A. Gelman. General Methods for Monitoring Convergence of Iterative Simulations. Journal of Computational and Graphical Statistics 7(4):434–455, 1998.
- Creutz (1988) M. Creutz. Global Monte Carlo algorithms for many-fermion systems. Physical Review D, 38(4):1228–1238, 1988.
- Fok et. al. (2016) R. Fok, A. An and X. Wang. Geodesic and Contour Optimization Using Conformal Mapping. Journal of Global Optimization, 2016. doi:10.1007/s10898-016-0467-8
- Gelfand and Smith (1990) A. E. Gelfand and A. F. M. Smith. Sampling-Based Approaches to Calculating Marginal Densities. Journal of the American Statistical Association 85 (410): 398–409.
- Gilks et. al. (1998) Walter R. Gilks; Gareth O. Roberts; Sujit K. Sahu. Adaptive Markov Chain Monte Carlo through Regeneration. Journal of the American Statistical Association 93:443:1045–1054, 1998.
- Girolami and Calderhead (2011) M. Girolami and B. Calderhead. Riemann manifold Langevin and Hamiltonian Monte Carlo methods. J. R. Statist. Soc. B 73(2),:1–37, 2011.
- Hastings (1970) W. K. Hastings, Monte Carlo Sampling Methods Using Markov Chains and Their Applications. Biometrika 57(1):97–109, 1970.
- Hoffman and Gelman (2014) M. D. Hoffman and A. Gelman. Journal of Machine Learning Research 15:1351-1381, 2014
- Ihler et. al. (2005) A. T. Ihler, J. W. Fisher, R. L. Moses and A. S. Willsky. IEEE Journal on Selected Areas in Communication 73(2):809–819, 2005.
- Kristan et. al. (2011)
- Metropolis et. al. (1953) N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, and E. Teller. Equations of State Calculations by Fast Computing Machines. Journal of Chemical Physics 21 (6): 1087–1092, 1953.
- Mykland et. al. (1995) P. Mykland, L. Tierney, and B. Yu. Regeneration in markov chain samplers. Journal of the American Statistical Association, 90(429):233–241, 1995.
- Lan et. al. (2012) S. Lan, V. Stathopoulos, B. Shahbaba and M. Girolami. Lagrangian Dynamical Monte Carlo. arXiv preprint arXiv:1211.3759, 2012.
- Lan et. al. (2013) S. Lan, J. Streets and B. Shahbaba. Wormhole Hamiltonian Monte Carlo. arXiv preprint arXiv:1306.0063, 2013.
- Neal (2012) R.M. Neal, MCMC using Hamiltonian dynamics. arXiv preprint arXiv:1206.1901, 2012
- Rao (1945) C. R. Rao. Information and accuracy attainable in the estimation of statistical parameters. Bull. Calc. Math. Soc., 37, 81–91, 1945.
Sminchisescu and Welling (2007)
C. Sminchisescu and M.Welling. Generalized dartingmonte carlo. In Eleventh International Conference on Artificial Intelligence and Statistics (AISTATS2007), 2007. online proceedings.
- Wang et. al. (2013) Z. Wang, S. Mohamed and N. de Freitas, Adaptive Hamiltonian and Riemann Monte Carlo Samplers. International Conference on Machine Learning (ICML), 2013.