Archive for Metropolis-Hastings algorithm

mostly Monte Carlo in June

Posted in Statistics, University life with tags , , , , , , , , , , , , , , , , , , on May 30, 2026 by xi'an

The last episode of the academic year for our mostly Monte Carlo seminar, next week:

On Friday 05/06/26, from 3-5pm at PariSanté Campus

15h00: Sam Livingstoke (University College London)

Skew-symmetric numerical schemes for stochastic differential equations: strong convergence and multi-level extension
I will discuss recent work fusing together two strands of the applied mathematics and statistics literature, one concerned with developing flexible probability distributions for data that rely on a small number of parameters, and another concerned with developing numerical integration schemes to simulate stochastic processes.  The specific case that I will focus on uses the skew-symmetric family of probability distributions introduced by Adelchi Azzalini and co-authors to approximate the transition kernels of diffusion processes over small time steps, producing alternative numerical schemes to the classical Euler-Maruyama approach.  Applying the scheme to the overdamped Langevin diffusion leads to an unadjusted version of the Barker proposal Metropolis-Hastings algorithm.  In earlier work weak accuracy was established over finite and infinite time scales, crucially without needing a globally Lipschitz assumption on the drift of the stochastic differential equation.  I will review this and then discuss more recent work establishing strong convergence in the mean-squared sense using a novel coupling between the numerical and exact processes.  This also enables the development of a multi-level Monte Carlo scheme, which I will discuss the merits of with particular focus on the superlinear drift case, as compared to Euler and Tamed Euler alternatives.
This is joint work with Yuga Iguchi, Giorgos Vasdekis & Rui-Yang Zhang.
16h00: Dana Naderi (Université Paris Dauphine PSL)
Approximating evidence via bounded harmonic means

Efficient Bayesian model selection relies on the model evidence or marginal likelihood, whose computation often requires evaluating an intractable integral. The harmonic mean estimator (HME) has long been a standard method of approximating the evidence. While computationally simple, the version introduced by Newton and Raftery (1994) potentially suffers from infinite variance. To overcome this issue, Gelfand and Dey (1994) defined a standardized representation of the estimator based on an instrumental function and Robert and Wraith (2009) later proposed to use higher posterior density (HPD) indicators as instrumental functions. Following this approach, a practical method is proposed, based on an elliptical covering of the HPD region with non-overlapping ellipsoids. The resulting estimator, called the Elliptical Covering Marginal Likelihood Estimator (ECMLE), not only eliminates the infinite-variance issue of the original HME and allows exact volume computations, but is also able to be used in multimodal settings. Through several examples, we illustrate that ECMLE outperforms other recent methods such as THAMES and its improved version (Metodiev et al. 2025). Moreover, ECMLE demonstrates lower variance a key challenge that subsequent HME variants have sought to address-and provides more stable evidence approximations, even in challenging settings.

This is joint work with Kaniav Kamari, Dareen Wraith & myself (X).

proximal sampler

Posted in Books, pictures, R, Statistics, Travel, University life with tags , , , , , , , , , , , , on April 28, 2025 by xi'an

At the Columbia workshop last week, Andre Wibisono presented work related with a recent arXival on the exponentially fast convergence of both unadjusted Langevin and  proximal sampler algorithms under strong [definitely strong] log-concavity assumptions. The idea behind the proximal sampler is to target the demarginalised density

g(x,y) \propto \exp\{\log f(x) - ||x-y||^2/2\eta\}\quad\eta>0

by introducing an auxiliary Gaussian vector y, which preserves f(x) as the marginal distribution on the first component vector X. While the auxiliary Y is (obviously) conditionally Gaussian, the conditional of X is at least as challenging as simulating from f. Unless η is chosen small enough to regularize log g(x) into a strongly log-concave function, since

\log g(x,y) \le \log g(x^\star,y) -\beta||x-x^\star||^2

when x*=x*(y) is the maximiser of log g(x,y) (for a given value y) and β>0 is the appropriate log-concavity constant. This inequality means that an accept-reject can be implemented to simulate from the conditional of X given Y but it requires both the factor β and the derivation of x*(y), hence a pretty good understanding and a rather high regularity of the actual target f(x). Besides, the regularization term ||x-y||² means that y is approximately the previous value of the (sub)chain X, hence it creates a rappel force that slows down the exploration of the target.

Since the arXival does not contain numerical comparisons, I attempted one using the (2D) banana shaped distribution,

target=function(x,sig,B,mu)-x[1]^2/2/sig-(x[2]+B*x[1]^2-mu)^2/2

with μ=σ=B=10. Comparing with a vanilla random walk Metropolis with three potential scales, chosen randomly at each iteration. Since I did not want to check whether or not the target was log-concave (and derive the corresponding β), I used the Normal distribution centred at proposal x*(y) of a Metropolis step, again with several scales. The following is the representation of the samples (sienna for MCMC, navy blue for proximal with β=50, dark green for β=5), with a lesser rate of tail exploration for the proximal samplers. It is thus unclear to me the theoretical characterisations of the method translate into practical efficiency beyond the most regular cases.

importance sampling and independent Metropolis–Hastings with unbounded weights

Posted in Books, Statistics with tags , , , , , , , , , , , , on December 12, 2024 by xi'an

George Deligiannidis, Pierre E. Jacob, El Mahdi Khribch, and Guanyang Wang just arXived a paper on the respective behaviours of importance sampling and independent Metropolis–Hastings (IMH) under the same proposal when the importance weight is unbounded but enjoys a p-th moment with p≥2. Both algorithms are sharing a lot, with importance sampling appearing as a rough Rao-Blackwellisation of Metropolis-Hastings, and its asymptotic variance being smaller than the asymptotic variance of Metropolis-Hastings. I was unable to check whether or not their conditions encompassed the highly interesting case when the integrand f is integrable under the target π, but not L²(π). (Theorem 2.3 does not seem to include this case.)

They consider a particular (!) version of Metropolis–Hastings (IMH) under the same proposal when the importance weight is unbounded but enjoys a p-th moment. Both algorithms are sharing a lot, with importance sampling appearing when N iid proposed values are drawn at once and accepted or rejected (again at once) with an acceptance ratio the average of the weights. Although this is already found in a 2010 paper by Christophe Andrieu and co-authors, and stem from an unbiased importance sampler, I was not aware of this version. My initial feeling (predictably) was pessimistic, but thinking about it, using the average weight brings into the sample simulations with small weights that would otherwise be discarded. Of course, a rejection proves N times more costly. But this is truly a form of Rao-Blackwellisation in the sense that it removes the weight variability to some extent (see p5) and it turns the outcome into an unbiased estimator. Despite the self-normalising behaviour! They also conclude that the rejection probability is at least c/√N  on average (Remark 4.1).

“We show that the bias of self-normalized importance sampling is of order N −1, and we obtain new bounds on the moments of the error in importance sampling. We then consider IMH, and show that the common random numbers coupling is optimal. Using this coupling, we show that the total variation distance between IMH at iteration t and π decays as tp-1.”

They also compare the biases in sampling importance resampling and independent Metropolis–Hastings, with the later getting the upper hand, but I do not see the justification in resampling when computing an integral. Since this does not a sample from the target, especially when the weights are unbounded, and adds to the variability of the estimator. They further propose a (telescopic) unbiased modification to the self-normalised importance sampling estimator, with an inefficiency twice as high. But a neat Rao-Blackwellisation trick brings it back to the same level!

biXarre, biXarre

Posted in Books, Statistics with tags , , , , , on May 2, 2024 by xi'an

combining normalizing flows and QMC

Posted in Books, Kids, Statistics with tags , , , , , , , , , , , , , on January 23, 2024 by xi'an

My PhD student Charly Andral [presented at the mostly Monte Carlo seminar and] arXived a new preprint yesterday, on training a normalizing flow network as an importance sampler (as in Gabrié et al.) or an independent Metropolis proposal, and exploiting its invertibility to call quasi-Monte Carlo low discrepancy sequences to boost its efficiency. (Training the flow is not covered by the paper.) This extends the recent study of He et al. (which was presented at MCM 2023 in Paris) to the normalising flow setting. In the current experiments, the randomized QMC samples are computed using the SciPy package (Roy et al. 2023), where the Sobol’ sequence is based on Joe and Kuo (2008) and on Matouˇsek (1998) for the scrambling, and where the Halton sequence is based on Owen (2017). (No pure QMC was harmed in the process!) The flows are constructed using the package FlowMC. As expected the QMC version brings a significant improvement in the quality of the Monte Carlo approximations, for equivalent computing times, with however a rapid decrease in the efficiency as the dimension of the targetted distribution increases. On the other hand, the architecture of the flow demonstrates little relevance. And the type of  RQMC sequence makes a difference, the advantage apparently going to a scrambled Sobol’ sequence.