This article shows how coupled Markov chains that meet exactly after a random number of iterations can be used to generate unbiased estimators of the solutions of the Poisson equation. Through this connection, we rederive known unbiased estimators of expectations with respect to the stationary distribution of a Markov chain and provide conditions for the finiteness of their moments. We further construct unbiased estimators of the asymptotic variance of Markov chain ergodic averages, and provide conditions for the finiteness of the estimators' moments of any order. If their second moment is finite, the average of independent copies of such estimators converges to the asymptotic variance at the Monte Carlo rate, comparing favorably to known rates for batch means and spectral variance estimators. The results are illustrated with numerical experiments.

This document presents methods to remove the initialization or burn-in bias from Markov chain Monte Carlo (MCMC) estimates, with consequences on parallel computing, convergence diagnostics and performance assessment. The document is written as an introduction to these methods for MCMC users. Some theoretical results are mentioned, but the focus is on the methodology.

Importance sampling and independent Metropolis-Hastings are among the fundamental building blocks of Monte Carlo methods. Both require a proposal distribution that globally approximates the target distribution, and pointwise evaluation of the Radon-Nikodym derivative of the target distribution relative to the proposal, also called the weight function. We study the bias of importance sampling and independent Metropolis-Hastings, without assuming that the weight function is bounded. We show that the common random numbers coupling of independent Metropolis-Hastings is maximal. Using that coupling, we derive polynomial bounds on the total variation distance of the chain to its target distribution. We further consider bias removal techniques using couplings, and provide conditions under which the resulting unbiased estimators have finite moments, and under which their efficiency is comparable to that of importance sampling. Experiments illustrate unbiased estimators of the inverse of a normalizing constant, estimators of nested expectations, and combination of importance sampling with robust mean estimation methods.

Statisticians often use Monte Carlo methods to approximate probability distributions, primarily with Markov chain Monte Carlo and importance sampling. Sequential Monte Carlo samplers are a class of algorithms that combine both techniques to approximate distributions of interest and their normalizing constants. These samplers originate from particle filtering for state space models and have become general and scalable sampling techniques. This article describes sequential Monte Carlo samplers and their possible implementations, arguing that they remain under-used in statistics, despite their ability to perform sequential inference and to leverage parallel processing resources among other potential benefits.

Background In early 2020, the response to the SARS-CoV-2 pandemic focused on non-pharmaceutical interventions, some of which aimed to reduce transmission by changing mixing patterns between people. Aggregated location data from mobile phones are an important source of real-time information about human mobility on a population level, but the degree to which these mobility metrics capture the relevant contact patterns of individuals at risk of transmitting SARS-CoV-2 is not clear. In this study we describe changes in the relationship between mobile phone data and SARS-CoV-2 transmission in the USA. Methods In this population-based study, we collected epidemiological data on COVID-19 cases and deaths, as well as human mobility metrics collated by advertisement technology that was derived from global positioning systems, from 1396 counties across the USA that had at least 100 laboratory-confirmed cases of COVID-19. We grouped these counties into six ordinal categories, defined by the National Center for Health Statistics (NCHS) and graded from urban to rural, and quantified the changes in COVID-19 transmission using estimates of the effective reproduction number (R-t) between Jan 22 and July 9,2020, to investigate the relationship between aggregated mobility metrics and epidemic trajectory. For each county, we model the time series of R-t values with mobility proxies. Findings We show that the reproduction number is most strongly associated with mobility proxies for change in the travel into counties (0.757 [95% CI 0.689 to 0.857]), but this relationship primarily holds for counties in the three most urban categories as defined by the NCHS. This relationship weakens considerably after the initial 15 weeks of the epidemic (0.442 [-0.492 to -0.392]), consistent with the emergence of more complex local policies and behaviours, including masking. Interpretation Our study shows that the integration of mobility metrics into retrospective modelling efforts can be useful in identifying links between these metrics and R-t. Importantly, we highlight potential issues in the data generation process for transmission indicators derived from mobile phone data, representativeness, and equity of access, which must be addressed to improve the interpretability of these data in public health. Copyright (C) 2021 The Author(s). Published by Elsevier Ltd.

We consider Markov chain Monte Carlo (MCMC) algorithms for Bayesian high-dimensional regression with continuous shrinkage priors. A common challenge with these algorithms is the choice of the number of iterations to perform. This is critical when each iteration is expensive, as is the case when dealing with modern data sets, such as genome-wide association studies with thousands of rows and up to hundreds of thousands of columns. We develop coupling techniques tailored to the setting of high-dimensional regression with shrinkage priors, which enable practical, non-asymptotic diagnostics of convergence without relying on traceplots or long-run asymptotics. By establishing geometric drift and minorization conditions for the algorithm under consideration, we prove that the proposed couplings have finite expected meeting time. Focusing on a class of shrinkage priors which includes the ‘Horseshoe’, we empirically demonstrate the scalability of the proposed couplings. A highlight of our findings is that less than 1000 iterations can be enough for a Gibbs sampler to reach stationarity in a regression on 100,000 covariates. The numerical results also illustrate the impact of the prior on the computational efficiency of the coupling, and suggest the use of priors where the local precisions are Half- t distributed with degree of freedom larger than one.

We consider Markov chain Monte Carlo (MCMC) algorithms for Bayesian high-dimensional regression with continuous shrinkage priors. A common challenge with these algorithms is the choice of the number of iterations to perform. This is critical when each iteration is expensive, as is the case when dealing with modern data sets, such as genome-wide association studies with thousands of rows and up to hundred of thousands of columns. We develop coupling techniques tailored to the setting of high-dimensional regression with shrinkage priors, which enable practical, non-asymptotic diagnostics of convergence without relying on traceplots or long-run asymptotics. By establishing geometric drift and minorization conditions for the algorithm under consideration, we prove that the proposed couplings have finite expected meeting time. Focusing on a class of shrinkage priors which includes the “Horseshoe”, we empirically demonstrate the scalability of the proposed couplings. A highlight of our findings is that less than 1000 iterations can be enough for a Gibbs sampler to reach stationarity in a regression on 100, 000 covariates. The numerical results also illustrate the impact of the prior on the computational efficiency of the coupling, and suggest the use of priors where the local precisions are Half-t distributed with degree of freedom larger than one.

Sequential Monte Carlo (SMC) samplers form an attractive alternative to MCMC for Bayesian computation. However, their performance depends strongly on the Markov kernels used to rejuvenate particles. We discuss how to calibrate automatically (using the current particles) Hamiltonian Monte Carlo kernels within SMC. To do so, we build upon the adaptive SMC approach of Fearnhead and Taylor (2013), and we also suggest alternative methods. We illustrate the advantages of using HMC kernels within an SMC sampler via an extensive numerical study.

Couplings play a central role in the analysis of Markov chain Monte Carlo algorithms and appear increasingly often in the algorithms themselves, e.g. in convergence diagnostics, parallelization, and variance reduction techniques. Existing couplings of the Metropolis-Hastings algorithm handle the proposal and acceptance steps separately and fall short of the upper bound on one-step meeting probabilities given by the coupling inequality. This paper introduces maximal couplings which achieve this bound while retaining the practical advantages of current methods. We consider the properties of these couplings and examine their behavior on a selection of numerical examples.

Bayesian inference provides a framework to combine an arbitrary number of model components with shared parameters, allowing joint uncertainty estimation and the use of all available data sources. However, misspecification of any part of the model might propagate to all other parts and lead to unsatisfactory results. Cut distributions have been proposed as a remedy, where the information is prevented from flowing along certain directions. We consider cut distributions from an asymptotic perspective, find the equivalent of the Laplace approximation, and notice a lack of frequentist coverage for the associate credible regions. We propose algorithms based on the Posterior Bootstrap that deliver credible regions with the nominal frequentist asymptotic coverage. The algorithms involve numerical optimization programs that can be performed fully in parallel. The results and methods are illustrated in various settings, such as causal inference with propensity scores and epidemiological studies.

Agent-based models of disease transmission involve stochastic rules that specify how a number of individuals would infect one another, recover or be removed from the population. Common yet stringent assumptions stipulate interchangeability of agents and that all pairwise contact are equally likely. Under these assumptions, the population can be summarized by counting the number of susceptible and infected individuals, which greatly facilitates statistical inference. We consider the task of inference without such simplifying assumptions, in which case, the population cannot be summarized by low-dimensional counts. We design improved particle filters, where each particle corresponds to a specific configuration of the population of agents, that take either the next or all future observations into account when proposing population configurations. Using simulated data sets, we illustrate that orders of magnitude improvements are possible over bootstrap particle filters. We also provide theoretical support for the approximations employed to make the algorithms practical.

Global efforts to prevent the spread of the SARS-COV-2 pandemic in early 2020 focused on non-pharmaceutical interventions like social distancing; policies that aim to reduce transmission by changing mixing patterns between people. As countries have implemented these interventions, aggregated location data from mobile phones have become an important source of real-time information about human mobility and behavioral changes on a population level. Human activity measured using mobile phones reflects the aggregate behavior of a subset of people, and although metrics of mobility are related to contact patterns between people that spread the coronavirus, they do not provide a direct measure. In this study, we use results from a nowcasting approach from 1,396 counties across the US between January 22nd, 2020 and July 9th, 2020 to determine the effective reproductive number (R(t)) along an urban/rural gradient. For each county, we compare the time series of R(t) values with mobility proxies from mobile phone data from Camber Systems, an aggregator of mobility data from various providers in the United States. We show that the reproduction number is most strongly associated with mobility proxies for change in the travel into counties compared to baseline, but that the relationship weakens considerably after the initial 15 weeks of the epidemic, consistent with the emergence of a more complex ecosystem of local policies and behaviors including masking. Importantly, we highlight potential issues in the data generation process, representativeness and equity of access which must be addressed to allow for general use of these data in public health.

The Sliced-Wasserstein distance (SW) is being increasingly used in machine learning applications as an alternative to the Wasserstein distance and offers significant computational and statistical benefits. Since it is defined as an expectation over random projections, SW is commonly approximated by Monte Carlo. We adopt a new perspective to approximate SW by making use of the concentration of measure phenomenon: under mild assumptions, one-dimensional projections of a high-dimensional random vector are approximately Gaussian. Based on this observation, we develop a simple deterministic approximation for SW. Our method does not require sampling a number of random projections, and is therefore both accurate and easy to use compared to the usual Monte Carlo approximation. We derive nonasymptotical guarantees for our approach, and show that the approximation error goes to zero as the dimension increases, under a weak dependence condition on the data distribution. We validate our theoretical findings on synthetic datasets, and illustrate the proposed approximation on a generative modeling problem.

We consider a vector of $N$ independent binary variables, each with a different probability of success. The distribution of the vector conditional on its sum is known as the conditional Bernoulli distribution. Assuming that $N$ goes to infinity and that the sum is proportional to $N$, exact sampling costs order $N^2$, while a simple Markov chain Monte Carlo algorithm using 'swaps' has constant cost per iteration. We provide conditions under which this Markov chain converges in order $N \log N$ iterations. Our proof relies on couplings and an auxiliary Markov chain defined on a partition of the space into favorable and unfavorable pairs.

In state–space models, smoothing refers to the task of estimating a latent stochastic process given noisy measurements related to the process. We propose an unbiased estimator of smoothing expectations. The lack-of-bias property has methodological benefits: independent estimators can be generated in parallel, and CI can be constructed from the central limit theorem to quantify the approximation error. To design unbiased estimators, we combine a generic debiasing technique for Markov chains, with a Markov chain Monte Carlo algorithm for smoothing. The resulting procedure is widely applicable and we show in numerical experiments that the removal of the bias comes at a manageable increase in variance. We establish the validity of the proposed estimators under mild assumptions. Numerical experiments are provided on toy models, including a setting of highly informative observations, and for a realistic Lotka–Volterra model with an intractable transition density. Supplementary materials for this article are available online.

Markov chain Monte Carlo (MCMC) methods provide consistent of integrals as the number of iterations goes to infinity. MCMC estimators are generally biased after any fixed number of iterations. We propose to remove this bias by using couplings of Markov chains together with a telescopic sum argument of Glynn and Rhee (2014). The resulting unbiased estimators can be computed independently in parallel. We discuss practical couplings for popular MCMC algorithms. We establish the theoretical validity of the proposed estimators and study their efficiency relative to the underlying MCMC algorithms. Finally, we illustrate the performance and limitations of the method on toy examples, on an Ising model around its critical temperature, on a high-dimensional variable selection problem, and on an approximation of the cut distribution arising in Bayesian inference for models made of multiple modules.

Performing numerical integration when the integrand itself cannot be evaluated point-wise is a challenging task that arises in statistical analysis, notably in Bayesian inference for models with intractable likelihood functions. Markov chain Monte Carlo (MCMC) algorithms have been proposed for this setting, such as the pseudo-marginal method for latent variable models and the exchange algorithm for a class of undirected graphical models. As with any MCMC algorithm, the resulting estimators are justified asymptotically in the limit of the number of iterations, but exhibit a bias for any fixed number of iterations due to the Markov chains starting outside of stationarity. This "burn-in" bias is known to complicate the use of parallel processors for MCMC computations. We show how to use coupling techniques to generate unbiased estimators in finite time, building on recent advances for generic MCMC algorithms. We establish the theoretical validity of some of these procedures by extending existing results to cover the case of polynomially ergodic Markov chains. The efficiency of the proposed estimators is compared with that of standard MCMC estimators, with theoretical arguments and numerical experiments including state space models and Ising models.

We investigated whether communicative cues help observers to make sense of human interaction. We recorded EEG from an observer monitoring two individuals who were occasionally communicating with each other via either mutual eye contact and/or pointing gestures, and then jointly attending to the same object or attending to different objects that were placed on a table in front of them. The analyses were focussed on the processing of the interaction outcome (i.e. presence or absence of joint attention) and showed that its interpretation is a two-stage process, as reflected in the N300 and the N400 potentials. The N300 amplitude was reduced when the two individuals shared their focus of attention, which indicates the operation of a cognitive process that involves the relatively fast identification and evaluation of actor–object relationships. On the other hand, the N400 was insensitive to the sharing or distribution of the two individuals’ attentional focus. Interestingly, the N400 was reduced when the interaction outcome was preceded either by mutual eye contact or by a perceived pointing gesture. This shows that observation of communication “opens up” the mind to a wider range of action possibilities and thereby helps to interpret unusual outcomes of social interactions.