We study general coordinate-wise MCMC schemes (such as Metropolis-within-Gibbs samplers), which are commonly used to fit Bayesian non-conjugate hierarchical models. We relate their convergence properties to the ones of the corresponding (potentially not implementable) Gibbs sampler through the notion of conditional conductance. This allows us to study the performances of popular Metropolis-within-Gibbs schemes for non-conjugate hierarchical models, in high-dimensional regimes where both number of datapoints and parameters increase. Given random data-generating assumptions, we establish dimension-free convergence results, which are in close accordance with numerical evidences. Applications to Bayesian models for binary regression with unknown hyperparameters and discretely observed diffusions are also discussed. Motivated by such statistical applications, auxiliary results of independent interest on approximate conductances and perturbation of Markov operators are provided.

Infectious disease transmission is often modelled by discrete-valued stochastic state-transition processes. Due to a lack of complete data, Bayesian inference for these models often relies on data-augmentation techniques. These techniques are often inefficient or time consuming to implement. We introduce a novel data-augmentation Markov chain Monte Carlo method for discrete-time individual-based epidemic models, which we call the Rippler algorithm. This method uses the transmission model in the proposal step of the Metropolis-Hastings algorithm, rather than in the accept-reject step. We test the Rippler algorithm on simulated data and apply it to data on extended-spectrum beta-lactamase (ESBL)-producing E. coli collected in Blantyre, Malawi. We compare the Rippler algorithm to two other commonly used Bayesian inference methods for partially observed epidemic data, and find that it has a good balance between mixing speed and computational complexity.

We consider adaptive increasingly rare Markov chain Monte Carlo (MCMC) algorithms, which are adaptive MCMC methods, where the adaptation concerning the "past” happens less and less frequently over time. Under a contraction assumption with respect to a Wasserstein-like function we deduce upper bounds of the convergence rate of Monte Carlo sums taking a renormalisation factor into account that is "almost” the one that appears in a law of the iterated logarithm. We demonstrate the applicability of our results by considering different settings, among which are those of simultaneous geometric and uniform ergodicity. All proofs are carried out on an augmented state space, including the classical non-augmented setting as a special case. In contrast to other adaptive MCMC limit theory, some technical assumptions, like diminishing adaptation, are not needed.

In recent years, differential privacy has been adopted by tech-companies and governmental agencies as the standard for measuring privacy in algorithms. In this article, we study differential privacy in Bayesian posterior sampling settings. We begin by considering differential privacy in the most common privatisation setting in which Laplace or Gaussian noise is injected into the output. In an effort to achieve better differential privacy, we consider adopting Huber's contamination model for use within privacy settings, and replace at random data points with samples from a heavy-tailed distribution (instead of injecting noise into the output). We derive bounds for the differential privacy level (ε,δ) of our approach, without requiring bounded observation and parameter spaces, a restriction commonly imposed in the literature. We further consider for our approach the effect of sample size on the privacy level and the rate at which (ε,δ) converges to zero. Asymptotically, our contamination approach is fully private with no information loss. We also provide examples of inference models for which our approach applies, with theoretical convergence rate analysis and simulation studies.

In Bayesian statistics, exploring high-dimensional multimodal posterior distributions poses major challenges for existing MCMC approaches. This paper introduces the Annealed Leap-Point Sampler (ALPS), which augments the target distribution state space with modified annealed (cooled) distributions, in contrast to traditional tempering approaches. The coldest state is chosen such that its annealed density is well-approximated locally by a Laplace approximation. This allows for automated setup of a scalable mode-leaping independence sampler. ALPS requires an exploration component to search for the mode locations, which can either be run adaptively in parallel to improve these mode-jumping proposals, or else as a pre-computation step. A theoretical analysis shows that for a d-dimensional problem the coolest temperature level required only needs to be linear in dimension, $\mathcal{O}(d)$, implying that the number of iterations needed for ALPS to converge is $\mathcal{O}(d)$ (typically leading to overall complexity $\mathcal{O}(d^3)$ when computational cost per iteration is taken into account). ALPS is illustrated on several complex, multimodal distributions that arise from real-world applications. This includes a seemingly-unrelated regression (SUR) model of longitudinal data from U.S. manufacturing firms, as well as a spectral density model that is used in analytical chemistry for identification of molecular biomarkers.

In order to tackle the problem of sampling from heavy tailed, high dimensional distributions via Markov Chain Monte Carlo (MCMC) methods, Yang, Latuszyński, and Roberts (2022) (arXiv:2205.12112) introduces the stereographic projection as a tool to compactify ℝ^d and transform the problem into sampling from a density on the unit sphere 𝕊^d. However, the improvement in algorithmic efficiency, as well as the computational cost of the implementation, are still significantly impacted by the parameters used in this transformation. To address this, we introduce adaptive versions three stereographic MCMC algorithms - the Stereographic Random Walk (SRW), the Stereographic Slice Sampler (SSS), and the Stereographic Bouncy Particle Sampler (SBPS) - which automatically update the parameters of the algorithms as the run progresses. The adaptive setup allows to better exploit the power of the stereographic projection, even when the target distribution is neither centered nor homogeneous. Unlike Hamiltonian Monte Carlo (HMC) and other off-the-shelf MCMC samplers, the resulting algorithms are robust to starting far from the mean in heavy-tailed, high-dimensional settings. To prove convergence properties, we develop a novel framework for the analysis of adaptive MCMC algorithms over collections of simultaneously uniformly ergodic Markov operators, which is applicable to continuous-time processes, such as SBPS. This framework allows us to obtain ℒ^2 and almost sure convergence results, and a CLT for our adaptive stereographic algorithms.

We consider a recently proposed class of MCMC methods which uses proximity maps instead of gradients to build proposal mechanisms which can be employed for both differentiable and non-differentiable targets. These methods have been shown to be stable for a wide class of targets, making them a valuable alternative to Metropolis-adjusted Langevin algorithms (MALA); and have found wide application in imaging contexts. The wider stability properties are obtained by building the Moreau-Yosida envelope for the target of interest, which depends on a parameter $\lambda$. In this work, we investigate the optimal scaling problem for this class of algorithms, which encompasses MALA, and provide practical guidelines for the implementation of these methods.

Establishing central limit theorems (CLTs) for ergodic averages of Markov chains is a fundamental problem in probability and its applications. Since the seminal work~\cite{MR834478}, a vast literature has emerged on the sufficient conditions for such CLTs. To counterbalance this, the present paper provides verifiable necessary conditions for CLTs of ergodic averages of Markov chains on general state spaces. Our theory is based on drift conditions, which also yield lower bounds on the rates of convergence to stationarity in various metrics. The validity of the ergodic CLT is of particular importance for sampling algorithms, where it underpins the error analysis of estimators in Bayesian statistics and machine learning. Although heavy-tailed sampling is of central importance in applications, the characterisation of the CLT and the convergence rates are theoretically poorly understood for almost all practically-used Markov chain Monte Carlo (MCMC) algorithms. In this setting our results provide sharp conditions on the validity of the ergodic CLT and establish convergence rates for large families of MCMC sampling algorithms for heavy-tailed targets. Our study includes a rather complete analyses for random walk Metropolis samplers (with finite- and infinite-variance proposals), Metropolis-adjusted and unadjusted Langevin algorithms and the stereographic projection sampler (as well as the independence sampler). By providing these sharp results via our practical drift conditions, our theory offers significant insights into the problems of algorithm selection and comparison for sampling heavy-tailed distributions (see short YouTube presentations~\cite{YouTube_talk} describing our \href{https://youtu.be/m2y7U4cEqy4}{\underline{theory}} and \href{https://youtu.be/w8I_oOweuko}{\underline{applications}}).

Enriching Brownian motion with regenerations from a fixed regeneration distribution $\mu$ at a particular regeneration rate $\kappa$ results in a Markov process that has a target distribution $\pi$ as its invariant distribution. For the purpose of Monte Carlo inference, implementing such a scheme requires firstly selection of regeneration distribution $\mu$, and secondly computation of a specific constant $C$. Both of these tasks can be very difficult in practice for good performance. We introduce a method for adapting the regeneration distribution, by adding point masses to it. This allows the process to be simulated with as few regenerations as possible and obviates the need to find said constant $C$. Moreover, the choice of fixed $\mu$ is replaced with the choice of the initial regeneration distribution, which is considerably less difficult. We establish convergence of this resulting self-reinforcing process and explore its effectiveness at sampling from a number of target distributions. The examples show that adapting the regeneration distribution guards against poor choices of fixed regeneration distribution and can reduce the error of Monte Carlo estimates of expectations of interest, especially when $\pi$ is skewed.

We develop the first exact Bayesian methodology for the problem of inference in discretely observed regime switching diffusions. Switching diffusion models extend ordinary diffusions by allowing for jumps in instantaneous drift and volatility. The jumps are driven by a latent, continuous time Markov switching process. We address the problem through an MCMC and an MCEM algorithm that target the exact posterior of diffusion parameters and the latent regime process. The algorithms are exact in the sense that they target the correct posterior distribution of the continuous model, so that the errors are due to Monte Carlo only. We illustrate the method on numerical examples, including an empirical analysis of the method's scalability in the length of the time series, and find that it is comparable in computational cost with discrete approximations while avoiding their shortcomings.

Equality-constrained models naturally arise in problems in which the measurements are taken at different levels of resolution. The challenge in this setting is that the models usually induce a joint distribution which is intractable. Resorting to instead sampling from the joint distribution by means of a Monte Carlo approach is also challenging. For example, a naive rejection sampler does not work when the probability mass of the constraint is zero. A typical example of such constrained problems is to learn energy consumption for a higher resolution level based on data at a lower resolution, for example, to decompose a daily reading into readings at a finer level. We introduce a novel Monte Carlo sampling algorithm based on Langevin diffusions and rejection sampling to solve the problem of sampling from equality-constrained models. Our method has the advantage of being exact for linear constraints and naturally deals with multimodal distributions on arbitrary constraints. We test our method on statistical disaggregation problems for electricity consumption datasets, and our approach provides better uncertainty estimation and accuracy in data imputation compared with other naive/unconstrained methods.

Piecewise Deterministic Markov Processes (PDMPs) such as the Bouncy Particle Sampler and the Zig-Zag Sampler, have gained attention as continuous-time counterparts of classical Markov chain Monte Carlo. We study their transient regime under convex potentials, namely how trajectories that start in low-probability regions move toward higher-probability sets. Using fluid-limit arguments with a decomposition of the generator into fast and slow parts, we obtain deterministic ordinary differential equation descriptions of early-stage behaviour. The fast dynamics alone are non-ergodic because once the event rate reaches zero it does not restart. The slow component reactivates the dynamics, so averaging remains valid when taken over short micro-cycles rather than with respect to an invariant law. Using the expected number of jump events as a cost proxy for gradient evaluations, we find that for Gaussian targets the transient cost of PDMP methods is comparable to that of random-walk Metropolis. For convex heavy-tailed families with subquadratic growth, PDMP methods can be more efficient when event simulation is implemented well. Forward Event-Chain and Coordinate Samplers can, under the same assumptions, reach the typical set with an order-one expected number of jumps. For the Zig-Zag Sampler we show that, under a diagonal-dominance condition, the transient choice of direction coincides with the solution of a box-constrained quadratic program; outside that regime we give a formal derivation and a piecewise-smooth update rule that clarifies the roles of the gradient and the Hessian. These results provide theoretical insight and practical guidance for the use of PDMP samplers in large-scale inference.

We investigate the increase in efficiency of simulated and parallel tempering MCMC algorithms when using non-reversible updates to give them "momentum". By making a connection to a certain simple discrete Markov chain, we show that, under appropriate assumptions, the non-reversible algorithms still exhibit diffusive behaviour, just on a different time scale. We use this to argue that the optimally scaled versions of the non-reversible algorithms are indeed more efficient than the optimally scaled versions of their traditional reversible counterparts, but only by a modest speed-up factor of about 42

This paper aims to provide differential privacy (DP) guarantees for Markov chain Monte Carlo (MCMC) algorithms. In a first part, we establish DP guarantees on samples output by MCMC algorithms as well as Monte Carlo estimators associated with these methods under assumptions on the convergence properties of the underlying Markov chain. In particular, our results highlight the critical condition of ensuring the target distribution is differentially private itself. In a second part, we specialise our analysis to the unadjusted Langevin algorithm and stochastic gradient Langevin dynamics and establish guarantees on their (Renyi) DP. To this end, we develop a novel methodology based on Girsanov's theorem combined with a perturbation trick to obtain bounds for an unbounded domain and in a non-convex setting. We establish: (i) uniform in n privacy guarantees when the state of the chain after n iterations is released, (ii) bounds on the privacy of the entire chain trajectory. These findings provide concrete guidelines for privacy-preserving MCMC.

Piecewise deterministic Markov processes (PDMPs) are a type of continuous-time Markov process that combine deterministic flows with jumps. Recently, PDMPs have garnered attention within the Monte Carlo community as a potential alternative to traditional Markov chain Monte Carlo (MCMC) methods. The Zig-Zag sampler and the Bouncy particle sampler are commonly used examples of the PDMP methodology which have also yielded impressive theoretical properties, but little is known about their robustness to extreme dependence or isotropy of the target density. It turns out that PDMPs may suffer from poor mixing due to anisotropy and this paper investigates this effect in detail in the stylised but important Gaussian case. To this end, we employ a multi-scale analysis framework in this paper. Our results show that when the Gaussian target distribution has two scales, of order $1$ and $\epsilon$, the computational cost of the Bouncy particle sampler is of order $\epsilon^{-1}$, and the computational cost of the Zig-Zag sampler is either $\epsilon^{-1}$ or $\epsilon^{-2}$, depending on the target distribution. In comparison, the cost of the traditional MCMC methods such as RWM or MALA is of order $\epsilon^{-2}$, at least when the dimensionality of the small component is more than $1$. Therefore, there is a robustness advantage to using PDMPs in this context.

Recent work has suggested using Monte Carlo methods based on piecewise deterministic Markov processes (PDMPs) to sample from target distributions of interest. PDMPs are non-reversible continuous-time processes endowed with momentum, and hence can mix better than standard reversible MCMC samplers. Furthermore, they can incorporate exact sub-sampling schemes which only require access to a single (randomly selected) data point at each iteration, yet without introducing bias to the algorithm's stationary distribution. However, the range of models for which PDMPs can be used, particularly with sub-sampling, is limited. We propose approximate simulation of PDMPs with sub-sampling for scalable sampling from posterior distributions. The approximation takes the form of an Euler approximation to the true PDMP dynamics, and involves using an estimate of the gradient of the log-posterior based on a data sub-sample. We thus call this class of algorithms stochastic-gradient PDMPs. Importantly, the trajectories of stochastic-gradient PDMPs are continuous and can leverage recent ideas for sampling from measures with continuous and atomic components. We show these methods are easy to implement, present results on their approximation error and demonstrate numerically that this class of algorithms has similar efficiency to, but is more robust than, stochastic gradient Langevin dynamics.

This article considers the challenge of designing football group draw mechanisms, which have a uniform distribution over all valid draw assignments, but are also entertaining, practical and transparent. Although this problem is trivial in completely symmetric problems, it becomes challenging when there are draw constraints that are not exchangeable across each of the competing teams, so that symmetry breaks down. We explain how to simulate the FIFA sequential draw method and compute the nonuniformity of its draws by comparison with a uniform rejection sampler. We then propose several practical methods of achieving the uniform distribution while still using balls and bowls in a way which is suitable for a televised draw. The solutions can also be carried out interactively. The general methodology we provide can readily be transported to different competition draws and is not restricted to football events.

Piecewise deterministic Markov processes provide scalable methods for sampling from the posterior distributions in big data settings by admitting principled sub-sampling strategies that do not bias the output. An important example is the Zig-Zag process of [Ann. Stats. 47 (2019) 1288 - 1320] where clever sub-sampling has been shown to produce an essentially independent sample at a cost that does not scale with the size of the data. However, sub-sampling also leads to slower convergence and poor mixing of the process, a behaviour which questions the promised scalability of the algorithm. We provide a large sample scaling analysis of the Zig-Zag process and its sub-sampling versions in settings of parametric Bayesian inference. In the transient phase of the algorithm, we show that the Zig-Zag trajectories are well approximated by the solution to a system of ODEs. These ODEs possess a drift in the direction of decreasing KL-divergence between the assumed model and the true distribution and are explicitly characterized in the paper. In the stationary phase, we give weak convergence results for different versions of the Zig-Zag process. Based on our results, we estimate that for large data sets of size n, using suitable control variates with sub-sampling in Zig-Zag, the algorithm costs O(1) to obtain an essentially independent sample; a computational speed-up of O(n) over the canonical version of Zig-Zag and other traditional MCMC methods

Particle filtering methods can be applied to estimation problems in discrete spaces on bounded domains, to sample from and marginalise over unknown hidden states. As in continuous settings, problems such as particle degradation can arise: proposed particles can be incompatible with the data, lying in low probability regions or outside the boundary constraints, and the discrete system could result in all particles having weights of zero. In this paper we introduce the Lifebelt Particle Filter (LBPF), a novel method for robust likelihood estimation in low-valued count problems. The LBPF combines a standard particle filter with one (or more) lifebelt particles which, by construction, lie within the boundaries of the discrete random variables, and therefore are compatible with the data. A mixture of resampled and non-resampled particles allows for the preservation of the lifebelt particle, which, together with the remaining particle swarm, provides samples from the filtering distribution, and can be used to generate unbiased estimates of the likelihood. The main benefit of the LBPF is that only one or few, wisely chosen, particles are sufficient to prevent particle collapse. Differently from other methods, there is no need to increase the number of particles, and therefore the computational effort, in regions of the parameter space that generate less likely hidden states. The LBPF can be used within a pseudo-marginal scheme to draw inferences on static parameters, theta, governing the system. We address here the estimation of a parameter governing probabilities of deaths and recoveries of hospitalised patients during an epidemic.

High-dimensional distributions, especially those with heavy tails, are notoriously difficult for off-the-shelf MCMC samplers: the combination of unbounded state spaces, diminishing gradient information, and local moves results in empirically observed ``stickiness'' and poor theoretical mixing properties -- lack of geometric ergodicity. In this paper, we introduce a new class of MCMC samplers that map the original high-dimensional problem in Euclidean space onto a sphere and remedy these notorious mixing problems. In particular, we develop random-walk Metropolis type algorithms as well as versions of the Bouncy Particle Sampler that are uniformly ergodic for a large class of light and heavy-tailed distributions and also empirically exhibit rapid convergence in high dimensions. In the best scenario, the proposed samplers can enjoy the ``blessings of dimensionality'' that the convergence is faster in higher dimensions.