# Approximate Methods for Bayesian Computation

Radu V. Craiu  1 and Evgeny Levi  2 1Department of Statistical Sciences,
University of Toronto, Toronto, Canada, M5G 1Z5; email<EMAIL_ADDRESS>2Department of Statistical Sciences, University of Toronto, Toronto, Canada,
M5G 1Z5

###### Abstract

Rich data generating mechanisms are ubiquitous in this age of information and
require complex statistical models to draw meaningful inference. While
Bayesian analysis has seen enormous development in the last 30 years,
benefitting from the impetus given by the successful application of Markov
chain Monte Carlo (MCMC) sampling, the combination of big data and complex
models conspire to produce significant challenges for the traditional MCMC
algorithms. We review modern algorithmic developments addressing the latter
and compare their performance using numerical experiments.

###### keywords:

ABC, Bayesian synthetic likelihood, coresets, divide and conquer, Markov chain
Monte Carlo, subsampling

###### Contents

  1. 1 Introduction
  2. 2 Modern challenges for Bayesian computation: Massive Data
    1. 2.1 Divide and Conquer
    2. 2.2 Subsampling
      1. 2.2.1 Coresets
      2. 2.2.2 Random Subsampling
  3. 3 Modern Challenges for Bayesian Computation: Intractable Likelihoods
    1. 3.1 Approximate Bayesian Computation (ABC)
    2. 3.2 Bayesian Synthetic likelihood
  4. 4 Double Jeopardy
  5. 5 Numerical experiments
    1. 5.1 Description of the simulation settings
    2. 5.2 Logistic Regression
      1. 5.2.1 German Credit Data
    3. 5.3 Stochastic Volatility
  6. 6 Conclusion and future directions

##  1 Introduction

The data science revolution has led to multiple pressure points in statistics.
A statistical sample from a large population exhibits, in the 21st century,
very different characteristics than what one would have seen merely a few
years ago. The ubiquitous and almost continuous recording of many of our
activities has made it relatively easy to collect enormous amounts of
information that require analysis and interpretation. This drastic increase in
data volume imposes sober re-evaluations of most classical approaches to
statistical inference.

The computational side of a Bayesian statistician’s toolbox is perhaps most
challenged by these developments. The impetus of Bayesian statistics that has
been felt since the early 1990s has been given by the spectacular advances in
computation, especially those around Markov chain Monte Carlo (MCMC) sampling.
Thanks to methods in this class of algorithms, the statisticians have been
liberated to think freely about the Bayesian model components used for a given
problem, without worrying about the mathematical intractability of the
analysis.

Indeed, given a data set ${\mathbf{y}}$, most of the pairings of a sampling
density, $f({\mathbf{y}}|{\mathbf{\theta}})$, and a prior,
$p({\mathbf{\theta}})$, result in a posterior distribution

$\pi({\mathbf{\theta}}|{\mathbf{y}})={p({\mathbf{\theta}})f({\mathbf{y}}|{\mathbf{\theta}})\over\int
p({\mathbf{\theta}})f({\mathbf{y}}|{\mathbf{\theta}})d{\mathbf{\theta}}}$
((1))

that cannot be analyzed directly, usually because the denominator in ((1))
cannot be computed analytically. The latter fact impedes the calculation of
quantities of interest related to $\pi$, most of which can be expressed as

$I=\int
h({\mathbf{\theta}})\pi({\mathbf{\theta}}|{\mathbf{y}})d{\mathbf{\theta}},$
((2))

for some function $h$ that is determined by the question of interest. For
instance, if ${\mathbf{\theta}}$ is univariate and we let
$h({\mathbf{\theta}})={\mathbf{\theta}}^{r}$ in ((2)), then $I$ is equal the
$r$-th moment of $\pi$, or $h({\mathbf{\theta}})={\bf
1}_{(-\infty,t]}({\mathbf{\theta}})$ leads to the cumulative distribution
function (cdf) of $\pi$ at a point $t$.

The classical Monte Carlo method, devised by Von Neumann & Ulam (1951) at the
middle of the twentieth century, relies on sampling independently
$\\{{\mathbf{\theta}}_{1},\ldots,{\mathbf{\theta}}_{m}\\}$ from distribution
$\pi$ and approximating $I$ with

$\hat{I}={1\over m}\sum_{k=1}^{m}h({\mathbf{\theta}}_{k}).$ ((3))

However, the unknown constant in ((1)) creates a knowledge gap that an MCMC
algorithm closes by constructing and running a Harris-recurrent,
$\pi$-irreducible, aperiodic Markov chain whose stationary distribution is
exactly the posterior distribution $\pi({\mathbf{\theta}}|{\mathbf{y}})$. The
values taken by the chain make up the samples
${\mathbf{\theta}}_{1},\ldots,{\mathbf{\theta}}_{m}$. A couple of issues
emerge immediately. First, because $\pi$ is the chain’s stationary
distribution, the samples will be approximatively distributed with $\pi$ only
after the chain has entered its stationary regime. Second, due to the Markov
property, the samples are typically positively correlated (although see
Frigessi et al., 2000, Craiu & Meng, 2005, for exceptions) which reduces the
amount of information they contain about $\pi$. To see that, let us imagine
the extreme case in which the $m$ samples are perfectly correlated, in which
case they would provide very little information about $\pi$.

The success MCMC sampling had in boosting the use of Bayesian models is
largely due to the ease of implementation of some of its most popular
algorithms. For instance, the Metropolis-Hastings algorithm (henceforth, MH)
(Metropolis et al., 1953, Hastings, 1970) can be implemented using the
following recursive procedure.

Step 0

Initialize the Markov chain at ${\mathbf{\theta}}_{0}$ and choose a proposal
density $q(\cdot|\zeta)$ which may or may not depend on $\zeta$.

Step t

At the $t$-th step ($1\leq t\leq m-1$) do:

PR

Draw proposal $\omega_{t}$ from the density $q(\cdot|{\mathbf{\theta}}_{t})$;

AR

Set

$\theta_{t+1}=\left\\{\begin{array}[]{cc}\omega_{t}&\mbox{ with probability
}\alpha_{t}\\\ {\mathbf{\theta}}_{t}&\mbox{ with probability }1-\alpha_{t}\\\
\end{array}\right.$

where

$\alpha_{t}=\min\left\\{1,{\pi(\omega_{t}|{\mathbf{y}})q({\mathbf{\theta}}_{t}|\omega_{t})\over\pi(\theta_{t}|{\mathbf{y}})q(\omega_{t}|{\mathbf{\theta}}_{t})}\right\\}.$
((4))

Because of the form of the acceptance probability ((4)), its calculation is
not prevented by the unknown denominator in ((1)). Nevertheless, computation
of ((4)) hinges on the ability to calculate the sampling density
$f({\mathbf{y}}|{\mathbf{\theta}})$ for any parameter value
${\mathbf{\theta}}$ and to be able to do it $m$ times. The challenges posed to
Bayesian computation by the modern data and modelling environment have their
roots in this implicit assumption.

In very broad strokes, one can speak of two main challenges in modern Bayesian
computation. The first one concerns the computational price of calculating a
likelihood when the data is massive, say of order $N$ (think of $N$ as being
of the order of hundred of millions or even billions). Even in the tame case
of iid observations, in order to know ((4)) we will have to compute a
likelihood (or sampling density) that involves $N$ terms and this will have to
be repeated each time a new MCMC sample is produced. The cumulative cost is
unsustainable as a single MCMC iteration can take days. A second challenge
emerges when the model’s complexity keeps up with the data volume and yields
an intractable likelihood so that ((4)) simply cannot be computed
analytically. Finally, a meta-challenge appears in Bayesian analyses that
merge massive data with intractable models.

This paper discusses MCMC-adjacent methodology that is used to alleviate the
pressure posed on Bayesian computation by the above challenges. Space
constraints impedes the presentation of details and variants, but in all cases
we present the main ideas and refer the interested reader to the relevant
literature. In order to gauge their computational efficiency, we run numerical
experiments where, using publicly available software packages, the algorithms
are implemented on two statistical models.

In the next section we describe in more detail the challenges we just
described and Section 3 summarizes some of the solutions proposed to address
them. Numerical experiments meant to illustrate and compare different
algorithms are reported in Section 4. The paper ends with comments and
discussion of future directions for research.

##  2 Modern challenges for Bayesian computation: Massive Data

Consider data ${\mathbf{y}}$ collected on $N$ independent items so that
${\mathbf{y}}=\\{y_{1},\ldots,y_{N}\\}\in{\cal X}^{N}$ and denote by
$f({\mathbf{y}}|{\mathbf{\theta}})$ the sampling distribution which depends on
parameter ${\mathbf{\theta}}\in\Theta\subset{\bf R}^{d}$. At each iteration of
the MH sampler, one needs to compute
$f({\mathbf{y}}|\omega_{t})=\prod_{k=1}^{N}f(y_{k}|\omega_{t})$ where
$\omega_{t}$ is the proposal in ((4)). Modern applications often rely on data
that are large enough so that the repeated calculation of
$f({\mathbf{y}}|\omega_{t})$ is impractical, even impossible. It is also not
unusual for data size to be so large as to prohibit storage on a single
machine, so that computation of the likelihood involves also repeated
communication between multiple machines, thus adding significantly to the
computational burden.

Prompted by the obstacle of large data, computational Bayesians have designed
a number of approaches to alleviate the problem. Two general ideas are
currently standing out in terms of popularity and usage: divide-and-conquer
(DAC) strategies and subsampling with minimum loss of information.

###  2.1 Divide and Conquer

The DAC approach is based on partitioning the sample into a number of sub-
samples, called batches that are analyzed separately on a number of workers
(CPUs, GPUs, servers, etc). After the batch-specific estimates about the
parameter of interest are obtained, the results are combined so that the
analyst recovers a large part of, ideally all, the information that would have
been available if the whole sample were analyzed in the usual way, on a single
machine. While this idea seems applicable in a wide range of scenarios, there
are a couple of constraints that restrict its generality. First, the procedure
is computationally effective if it is designed to minimize, preferably
eliminate, communication between the workers before combining the batch-
specific results. Second, it is often difficult to produce an accurate
assessment of the resulting loss of information at the combining stage. Some
of the first proponents of DAC for MCMC sampling are Neiswanger et al. (2013),
Scott et al. (2016), and Wang & Dunson (2013). In their approach, the
subposterior distribution corresponding to the $j$th batch, is defined as

$\pi^{(j)}({\mathbf{\theta}}|{\mathbf{y}}^{(j)})\propto
f({\mathbf{y}}^{(j)}|{\mathbf{\theta}})[p({\mathbf{\theta}})]^{1/J}$ ((5))

where $f,p$ are as in ((1)), ${\mathbf{y}}^{(j)}$ is the data that was
assigned to batch $j$, $1\leq j\leq J$, and $J$ is the total number of
batches. With this choice, one immediately gets that
$\prod_{j=1}^{J}\pi^{(j)}\propto\pi({\mathbf{\theta}}|{\mathbf{y}})$. Both
Neiswanger et al. (2013) and Scott et al. (2016) consider ways to combine
samples from the subposteriors $\pi^{(j)}({\mathbf{\theta}})$, $1\leq j\leq
J$, in situations in which all posteriors, batch-specific and full data ones,
are Gaussian or can be approximated by mixtures of Gaussians. in this case,
one can demonstrate that a weighted average of samples from all the
$\pi^{(j)}$’s have density $\pi$. The use of the Weierstrass transform for
each posterior density, proposed in Wang & Dunson (2013), extends the range of
theoretical validity beyond Gaussian distributions. The authors also establish
error bounds between the approximation and the true posterior. Nemeth &
Sherlock (2018) use a Gaussian process (GP) approximation of each
subposterior. Once again, the Gaussian nature of the approximation makes
recombination possible and relatively straightforward. Limitations of the
method are strongly linked with those of GP-based estimation. For instance,
when the sub-posterior samplers are sluggish, large MCMC samples might be
needed which, in turn, make the calculation of the GP-based approximation very
expensive. The idea of using the values of the sub-posterior at each MCMC
sample is adopted also by Changye & Robert (2019) who propose to define the
subposteriors using
$\pi^{(j)}\propto\\{[p({\mathbf{\theta}})]^{1/J}f({\mathbf{y}}^{(j)}|{\mathbf{\theta}})\\}^{\lambda_{j}}$.
The scale factor $\lambda_{j}$ is used to control the uncertainty in the
subposterior. Alternative ways to define the sub-posteriors are produced by
Entezari et al. (2018) who use $\pi^{(j)}\propto
p({\mathbf{\theta}})[f({\mathbf{y}}^{(j)}|{\mathbf{\theta}})]^{J}$. The
intuitive idea is to ”match” the size of the original sample and the batch-
specific one. Their approach has been applied successfully to BART (Chipman et
al., 2010, Pratola, 2016) models.

###  2.2 Subsampling

Subsampling approaches are mostly developed under two assumptions. The first
one is that with massive data one expects a certain amount of redundancy, so
it is possible to the same likelihood when we eliminate a proportion of the
sample as long as the remaining observations are properly weighted. A simple
illustration is one in which $R$ observations are identical, so that $R-1$ of
them can be taken out of the likelihood calculation if the term corresponding
to the remaining one is raised to power $R$. The second idea is that one might
use only a small percentage of the sample to find accurate (e.g., unbiased)
approximations of the quantities needed to run an MCMC sampler. For instance,
in the case of a MH sampler, the pseudo-marginal approach of Andrieu & Roberts
(2009) demonstrates that the stationary distribution of the chain is the same
when the likelihoods involved in ((4)) are replaced with unbiased estimators.
The pseudo-marginal idea has largely impacted the methods based on subsampling
for MCMC. While some divide the latter into exact and approximate, we will
refrain from using similar taxonomy because, in our opinion, all subsampling
MCMC methods introduce some level of approximation into the computation of
interest.

Early efforts include those of Korattikara et al. (2014) and Bardenet et al.
(2014) who propose estimating the acceptance probability ((4)) using only a
random subset of the data. The latter authors demonstrate that, with
probability higher than a threshold set in place by the user, their method
yields estimates that are equal to the one produced by the full data
likelihood. However, one does not know in advance the size of the sample
needed at each iteration and thus must be able, in principle, to access most
of it at all times. A review of early subsampling methods can be found in
Bardenet et al. (2017).

####  2.2.1 Coresets

The process of establishing which sample points are redundant must have
theoretical backing, lest it leads to a very different posterior distribution
without any hope to control or assess the error incurred. The coreset approach
of Campbell & Broderick (2019) offers theoretical guarantees about the quality
of the approximation resulting from sample reduction. Consider the
loglikelihood obtained from $N$ iid observations

$l({\mathbf{\theta}}|{\mathbf{y}})=\sum_{i=1}^{N}l_{i}({\mathbf{\theta}}|y_{i}),$
((6))

where $l_{i}({\mathbf{\theta}}|y_{i})=\log f(y_{i}|{\mathbf{\theta}})$. The
aim of the coreset method is to find a set of weights $\\{w_{i}:\;1\leq i\leq
N\\}$, most of them zero, so that

$\|\Lambda({\mathbf{\theta}}|{\mathbf{W}},{\mathbf{y}})-l({\mathbf{\theta}}|{\mathbf{y}})\|\leq{\epsilon}\|l({\mathbf{\theta}}|{\mathbf{y}})\|,$
((7))

for all ${\mathbf{\theta}}\in\Theta$, where
${\mathbf{W}}=(W_{1},\ldots,W_{N})$ is the vector of weights, and
$\Lambda({\mathbf{\theta}}|{\mathbf{W}},{\mathbf{y}})=\sum_{i=1}^{N}W_{i}l_{i}({\mathbf{\theta}}|y_{i})$
is the weighted log-likelihood of the coreset. The weights found by Campbell &
Broderick (2019) are defined as

$W_{i}={\sigma\over\sigma_{i}}{M_{i}\over M},$ ((8))

where

$\sigma_{i}=\sup_{{\mathbf{\theta}}\in\Theta}\left\|{l_{i}({\mathbf{\theta}}|y_{i})\over
l({\mathbf{\theta}}|{\mathbf{y}})}\right\|$ ((9))

is called the sensitivity of the $i$-th observation,
$\sigma=\sum_{i=1}^{N}\sigma_{i}$, $M$ is the size of the coreset and
$(M_{1},\ldots,M_{N})\sim Multi\left(M,\left\\{{\sigma_{i}\over\sigma}:\;1\leq
i\leq N\right\\}\right)$ are multinomial draws. One can think of the
sensitivity in ((9)) as a measure of the influence of the $i$-th observation
on the whole likelihood as $\theta$ varies. As expected, the algorithm will
retain observations that correspond to relatively higher likelihood values,
but more importantly, it allows some evaluation of the error incurred when the
sample is reduced. The ideas that led to the weights in ((8)) illustrate the
general principles of the approach, but improvements are possible when one
considers other norms in ((7)) and ((9)). For instance, Campbell & Broderick
(2019) consider the $l_{i}$’s as vectors in a Hilbert space, link the norm to
the inner product in the space and include directionality in the selection of
the coresets. The latter allows replacing the simultaneous selection of the
coreset elements by a more intuitive procedure in which samples are
sequentially added to minimize the residual error. In Section 3 we implement
the coreset approach for logistic regression as presented in Huggins et al.
(2016). For this model, the coreset is build along the principles delineated
above and requires some specific tuning. The parameter space is taken to be an
Euclidian ball of radius $R$ which is a reasonable working assumption in the
case of a logistic regression with standardized covariates. The sensitivity
measure for each point is modified after $K$-clustering the entire sample. A
measure of spread within each cluster is used to construct upper bounds for
the sensitivity of each point in the sample. The intuition guiding this choice
is that clusters whose data vectors are tightly bundled together will be well
represented in the coreset by only a few points, while clusters with more
spread will need to contribute more points. Overall, the coreset construction
is intuitive and offers many possible directions for future research. The
biggest challenge is the evaluation of approximating error induced in the
posterior when replacing the full sample by the coreset, although some
promising initial results exist (Manousakas et al., 2020).

####  2.2.2 Random Subsampling

One can think of coreset subsampling as a static approach, in the sense that
the subsample is selected once and the Bayesian analysis is subsequently
conducted using the coreset in lieu of the original sample. A more dynamic
approach is considered by Quiroz et al. (2018) who propose to use all the data
for inference, just not at once. Their idea is to use a different subset of
the data each time the MCMC chain is updated. For instance, in the case of an
MH sampler, a different subset of individuals will contribute to the
likelihood needed in the calculation of ((4)), at each iteration. Following
the development of pseudo-marginal strategies, Andrieu & Vihola (2015) studied
the convergence properties of an MH or a random walk Metropolis sampler in
which the likelihood in ((4)) is replaced by an unbiased estimator. They have
shown that the efficiency of the MCMC sample increases when the variance of
the unbiased estimator decreases.

The use of subsampling within MCMC proposed by Quiroz et al. (2018) can be
applied quite generally and it is attractive because it addresses both the
construction of the unbiased estimator for the likelihood and the reduction of
its variance.

Given a random subsample of ${\mathbf{y}}$ of size $m$,
${\mathbf{y}}_{{\mathbf{u}}}=\\{y_{u_{1}},\ldots,y_{u_{m}}\\}$, where
${\mathbf{u}}=\\{u_{1},\ldots,u_{m}\\}$ are iid random variables uniformly
distributed over $\\{1,\ldots,N\\}$, the following estimator

$l_{m}({\mathbf{\theta}}|{\mathbf{y}}_{{\mathbf{u}}})={1\over
m}\sum_{k=1}^{m}l_{u_{k}}({\mathbf{\theta}}|y_{u_{k}}),$ ((10))

is unbiased for the average log-likelihood ${1\over
N}l({\mathbf{\theta}}|{\mathbf{y}})$. However, it usually has a large
variance, subjecting the pseudo-marginal chain that uses ((10)) instead of the
full-sample likelihood in ((4)) to an increased risk of poor mixing, since an
unusually high value of the likelihood at the current state of the chain will
make it unlikely to accept a proposal. A reduction in variance is desirable
and can be achieved via control variates (CV),
${\mathbf{q}}({\mathbf{\theta}})=\\{q_{1}({\mathbf{\theta}}),\ldots,q_{N}({\mathbf{\theta}})\\}$,
and via a modified estimator of ((10)),

$\tilde{l}_{m}({\mathbf{\theta}}|{\mathbf{y}}_{{\mathbf{u}}},{\mathbf{q}})=\sum_{i=1}^{N}q_{i}({\mathbf{\theta}})+{N\over
m}\sum_{k=1}^{m}(l_{u_{k}}({\mathbf{\theta}}|y_{u_{k}})-q_{u_{k}}({\mathbf{\theta}})).$
((11))

The notation implies that ${\mathbf{q}}$ might change with
${\mathbf{\theta}}$. Indeed, when the likelihood is unimodal, the construction
of the control variate follows Bardenet et al. (2017) who use for each
${\mathbf{\theta}}$, a Taylor series expansion of $l(\theta|{\mathbf{y}})$
around a fixed point, ${\mathbf{\theta}}^{*}$ which is a point centrally
located in $\Theta$ ( e.g. the maximum likelihood estimate) so that for all
$1\leq i\leq N$

$q_{i}({\mathbf{\theta}})=l_{i}({\mathbf{\theta}}^{*}|y_{i})+({\mathbf{\theta}}-{\mathbf{\theta}}^{*})^{T}{d\over
d{\mathbf{\theta}}}l_{i}(\theta^{*}|y_{i})+{1\over
2}({\mathbf{\theta}}-{\mathbf{\theta}}^{*})^{T}{d^{2}\over
d{\mathbf{\theta}}^{2}}l_{i}(\theta^{*}|y_{i})({\mathbf{\theta}}-{\mathbf{\theta}}^{*}).$
((12))

This control variate is called parameter expanded by Quiroz et al. (2018)
because it is obtained using an expansion in the parameter space. With this
modification, running a MH chain for, say, $M$ iterations requires the
evaluation of $N+mM$ item-specific likelihood terms,
$l_{i}({\mathbf{\theta}}|y_{i})$, for ((11)) and $MN$ for ((6)). This can
translate into significant reduction of computation effort when $m<<N$.

When the likelihood is multimodal or the Taylor approximation is poor when
${\mathbf{\theta}},{\mathbf{\theta}}^{*}$ are distanced, the authors discuss
an alternative construction that identifies a number of centroids
${\mathbf{y}}^{*}_{1},\ldots,{\mathbf{y}}^{*}_{r}$ via clustering of the data
and uses Taylor series expansions around each centroid to define the so-called
data expanded control variates. In Section 5 we implement the subsampling
method with parameter and data expanded control variates.

The reduction in variance requires a careful derivation in which the source of
variability is provided by the finite distribution of the random vector
$(u_{1},\ldots u_{m})$. The latter can be sampled at random at each iteration
or one can use the ideas in Deligiannidis et al. (2018) and allow dependence
between the $u$’s in consecutive iterations to further reduce the variance of
((11)). In Section 5x we implement the subsampling method with random or
correlated selection of indices, and parameter or data expanded control
variates.

Finally, we should also point out that while the estimators discussed are
unbiased for the loglikelihood, this does not translate into an unbiased
estimator for the likelihood itself. Therefore, an approximate correction term
is applied to reduce the bias but does not dissolve it, which means that the
pseudo-marginal theory cannot be applied mutatis mutandis in this case.
Therefore, the target distribution of the chain is perturbed and one must
assess the size of the error incurred. The authors produce a bound of the
perturbation error and provide empirical evidence that their bound is
conservative. Additional details and derivations can be found in Quiroz et al.
(2018).

##  3 Modern Challenges for Bayesian Computation: Intractable Likelihoods

So far we have looked at the pressure posed by the size of the sample on
Bayesian computation. However, there are other hurdles that accompany a
massive sample. Often, large data imply more information which, in order to be
used fully, requires a more complex model. As data become richer and modellers
are more ambitious, the likelihoods tend to get intractable, such as the ones
used in population genetics (Pritchard et al., 1999, Beaumont et al., 2002),
groundwater studies (Cui et al., 2018), hurricane surges (Plumlee et al.,
2021), or in climate change scenarios (Oyebamiji et al., 2015).

At first sight, it can be surprising that Bayesian inference can still be
conducted when the likelihood is intractable. The likelihood provides a
crucial analytical link between any parameter value and the probability of
observing a given data set. When such a link is not analytically tractable it
will have to be inferred from simulations. Central to the latter approach is
the ability to sample, given any value of the parameter, pseudo-data from the
model. To provide an intuition, imagine that infinite computational resources
are available. Then one can see that for any ${\mathbf{\theta}}\in\Theta$ it
is possible to simulate enough pseudo-data sets to approximate at any degree
of precision the distribution of the observed-data
$f({\mathbf{y}}_{0}|{\mathbf{\theta}})$, essentially filling the void left by
the intractability of the likelihood. However, computational resources are not
infinite so ingenious ways are needed to reduce computational burden. We
discuss here two algorithms, Approximate Bayesian computation (ABC) and
Bayesian Synthetic Likelihood (BSL), that have gained popularity in the
statistical and, more generally, the scientific communities.

###  3.1 Approximate Bayesian Computation (ABC)

Our discussion of ABC will be brief, given the recent and excellent reviews of
Robert (2014), Sisson et al. (2018a), and the comprehensive handbook of ABC
(Sisson et al., 2018b).

The ABC agorithm was initially proposed as an accept/reject sampler (Tavaré et
al., 1997). Given any ${\mathbf{\theta}}^{*}$ sampled from the prior
$p({\mathbf{\theta}})$, it assumed that is possible to generate pseudo-data
${\mathbf{y}}$ from $f({\mathbf{y}}|{\mathbf{\theta}}^{*})$. If the pseudo-
data and the original data are close enough, then the parameter
${\mathbf{\theta}}^{*}$ is an approximate draw from the posterior
$\pi({\mathbf{\theta}}|{\mathbf{y}}_{0})$. Let us frame next the “close
enough” and “approximate draw” in precise mathematical terms and provide some
justification for our choices.

Given ${\epsilon}>0$, a distance $d:{\bf R}^{p}\times{\bf
R}^{p}\rightarrow{\bf R}_{+}$ and summary statistic $S({\mathbf{y}})\in{\bf
R}^{p}$, the ABC algorithm has the following steps :

  * S1

Sample ${\mathbf{\theta}}^{*}\sim p({\mathbf{\theta}})$ and synthetic data
${\mathbf{y}}\sim f({\mathbf{y}}|{\mathbf{\theta}}^{*})$

  * S2

If $d(S({\mathbf{y}}),S({\mathbf{y}}_{0}))\leq{\epsilon}$ then accept
${\mathbf{\theta}}^{*}$ as a sample from the approximate posterior
$\pi_{\epsilon}(\theta|S({\mathbf{y}}_{0}))$, the marginal (in
${\mathbf{\theta}}$) of the joint distribution

$\pi_{{\epsilon}}(\theta,{\mathbf{y}}|S({\mathbf{y}}_{0}))\propto
p(\theta)f({\mathbf{y}}|\theta){\bf
1}_{\\{d(S({\mathbf{y}}),S({\mathbf{y}}_{0}))<{\epsilon}\\}}.$ ((13))

If it is possible to have ${\mathbf{y}}={\mathbf{y}}_{0}$ (for instance, if
${\mathbf{y}}_{0}$ is a discrete random variable with finite support) we can
choose $S({\mathbf{y}})={\mathbf{y}}$ and ${\epsilon}=0$, then the approximate
posterior is the true posterior, i.e.
$\pi_{\epsilon}({\mathbf{\theta}}|{\mathbf{y}}_{0})=\pi({\mathbf{\theta}}|{\mathbf{y}}_{0})$.
This is easier to see when both ${\mathbf{\theta}}$ and ${\mathbf{y}}$ take
discrete values. Then, one can easily see that

$\displaystyle{\rm Pr}({\mathbf{\theta}}={\mathbf{\theta}}_{0})$
$\displaystyle\propto$ $\displaystyle
p({\mathbf{\theta}}_{0})\Pr({\mathbf{y}}={\mathbf{y}}_{0}|{\mathbf{\theta}}={\mathbf{\theta}}_{0})\propto\pi({\mathbf{\theta}}_{0}|{\mathbf{y}}_{0})$
((14))

where ((14)) holds because of the algorithm’s construction with
$S({\mathbf{y}})={\mathbf{y}}$ and ${\epsilon}=0$. The above can be easily
extended to the case when $S$ is a sufficient statistics. In general, models
with the level of complexity that requires ABC, will not have a low-
dimensional sufficient statistics so the choice of $S$ is central to the
performance of the ABC algorithm (Fearnhead & Prangle, 2012, Marin et al.,
2014). The accept-reject form of the ABC sampler makes it inefficient when the
prior and posterior place most of their mass on different regions of $\Theta$.
Recognizing this, Marjoram et al. (2003) proposed an ABC-MCMC algorithm which
relies on building a Metropolis-Hastings (MH) transition kernel, with state
space $\\{(\theta,{\mathbf{y}})\in{\bf R}^{q}\times{\mathcal{X}}^{n}\\}$,
proposal distribution at iteration $t$, $q(\theta|\theta_{t})\times
f({\mathbf{y}}|\theta)$, and target

$\pi_{{\epsilon}}(\theta,{\mathbf{y}}|S({\mathbf{y}}_{0}))\propto
p(\theta)f({\mathbf{y}}|\theta){\bf
1}_{\\{d(S({\mathbf{y}}),S({\mathbf{y}}_{0}))<{\epsilon}\\}}$ ((15))

for which ((4)) can be computed exactly, because the intractable terms
involving the likelihood, $f({\mathbf{y}}|{\mathbf{\theta}})$, cancel out.
Alternatives to this include the pseudo-marginal approach of Lee et al.
(2012), or the sequential Monte Carlo implementation of Sisson et al. (2007),
Lee (2012), and Filippi et al. (2013).

###  3.2 Bayesian Synthetic likelihood

Indirect inference was developed in econometrics (Smith Jr, 1993, Gourieroux
et al., 1993) for complex data models which are intractable, but can be
sampled from. The central tenet is that a complex model of interest,
$f({\mathbf{y}}|{\mathbf{\theta}})$, can be well approximated using a
tractable sampling model $g({\mathbf{y}}|\phi)$ where
$\dim(\phi)>\dim({\mathbf{\theta}})$. In other words, the complex model can be
approximated by a simpler model whose parameter space is of larger dimension
and has a tractable likelihood. For instance, Bayesian estimation of
${\mathbf{\theta}}$ is possible if one estimates its functional connection
with $\phi$ (see, for instance, Gallant & McCulloch, 2009).

The BSL algorithm (Price et al., 2018) relies on the synthetic likelihood (SL)
approximation of Wood (2010) which falls squarely in the class of indirect
inference methods. The idea hinges on the assumption that the conditional
distribution $p(S({\mathbf{y}})|{\mathbf{\theta}})$ is well approximated by a
multivariate Gaussian ${\cal
N}(\mu({\mathbf{\theta}}),\Sigma({\mathbf{\theta}}))$ whenever
${\mathbf{y}}\sim f({\mathbf{y}}|{\mathbf{\theta}})$. The SL is defined as
$\mbox{SL}({\mathbf{\theta}})=n(S({\mathbf{y}});\mu({\mathbf{\theta}}),\Sigma({\mathbf{\theta}}))$
where $n(\cdot;\mu,\Sigma)$ is the density of a multivariate normal with mean
$\mu$ and variance $\Sigma$. One can estimate
$\mu({\mathbf{\theta}}),\Sigma({\mathbf{\theta}})$ numerically, for any
${\mathbf{\theta}}$. Given ${\mathbf{\theta}}$, it is enough to repeatedly
sample pseudo-data ${\mathbf{y}}_{j}\sim f({\mathbf{y}}|{\mathbf{\theta}})$
and compute $S({\mathbf{y}}_{j})$ for $1\leq j\leq K$ and then estimate
$\hat{\mu}({\mathbf{\theta}})={1\over K}\sum_{j=1}^{K}S({\mathbf{y}}_{j})$ and
$\hat{\Sigma}({\mathbf{\theta}})=\mbox{SamVar}(\\{S({\mathbf{y}}_{j}):\;1\leq
j\leq K\\})$ where SamVar is the sample variance of the computed statistics.
BSL is then based on the approximation
$\pi_{BSL}({\mathbf{\theta}}|S({\mathbf{y}}_{0}))\propto
p({\mathbf{\theta}})\mbox{SL}({\mathbf{\theta}}|S({\mathbf{y}}_{0}))$ which
can be explored via MCMC sampling using the following update rule at iteration
$t>0$:

PR

Generate ${\mathbf{\theta}}^{*}\sim q(\cdot|{\mathbf{\theta}}_{t})$, estimate
$\hat{\mu}_{{\mathbf{\theta}}^{*}},\hat{\Sigma}_{{\mathbf{\theta}}^{*}}$ from
$K$ pseudo-data $\\{{\mathbf{y}}_{j}\sim
f({\mathbf{y}}|{\mathbf{\theta}}^{*}):\;1\leq j\leq K\\}$, and compute
$SL({\mathbf{\theta}}^{*})={\cal
N}(S({\mathbf{y}}_{0});\hat{\mu}_{{\mathbf{\theta}}^{*}},\hat{\Sigma}_{{\mathbf{\theta}}^{*}})$.

AR

Set ${\mathbf{\theta}}_{t+1}={\mathbf{\theta}}^{*}$ with probability
$\alpha=\min\left(1,\frac{p({\mathbf{\theta}}^{*})SL({\mathbf{\theta}}^{*})q({\mathbf{\theta}}_{t}|{\mathbf{\theta}}^{*})}{p({\mathbf{\theta}}_{t})SL({\mathbf{\theta}}_{t})q({\mathbf{\theta}}^{*}|{\mathbf{\theta}}_{t})}\right)$
and ${\mathbf{\theta}}_{t+1}={\mathbf{\theta}}_{t}$ otherwise.

##  4 Double Jeopardy

The separate treatment of the challenges brought by the big data or
intractable models is artificial and we anticipate that, more and more, the
two challenges will have to be met simultaneously. Since the use of MCMC
within ABC or BSL procedures requires repeated generation of pseudo-data of
the same size and complexity as the observed ones, it incurs unmanageable
computational costs when the data are massive or the data generating procedure
is expensive.

Some of the methods described within the first challenge are amenable to being
used in combination with ABC or BSL. For instance, DAC strategies can be used
for an intractable model if each worker runs a separate ABC MCMC sampler for
each batch of data. The obvious caveat is the difficulty of ascertaining the
loss of information after the merging stage. Unfortunately, more generalizable
methods like those used for subsampling cannot be used within ABC or BSL.

A strategy customized to ABC and BSL samplers with large or complex data is
proposed by Levi & Craiu (2022). We describe here a variation of their
approach which combines finite adaptation ideas and pre-sampling of the
proposals. Assuming that a MH transition kernel is used to implement ABC MCMC
or BSL MCMC, the first $B$ samples are used to tune the proposal distribution.
For instance, if a Gaussian proposal is used then its covariance matrix can be
estimated using methods proposed by Haario et al. (2001), Roberts & Rosenthal
(2009) or, in the case of multimodal targets, by Craiu et al. (2009) or Pompe
et al. (2020). The computational effort is reduced because we rely on a set of
proposals that are generated in advance. This allows an embarassingly parallel
procedure that benefits from the use of multiple workers. The preprocessed
draws are collected in reference set
${\mathcal{Z}}=\\{(\xi_{h},s_{h}=(s_{h}^{(1)},\ldots,s_{h}^{(m)})^{T}):1\leq
h\leq H\\}$ where for each parameter value $\xi_{h}$ generated from the
proposal distribution, we sample $m$ pseudo-data
${\mathbf{w}}_{h}^{(1)},\ldots,{\mathbf{w}}_{h}^{(m)}\stackrel{{\scriptstyle
iid}}{{\sim}}f({\mathbf{w}}|\xi_{h})$ and set
$s_{h}^{(j)}=S({\mathbf{w}}_{h}^{(j)})$ for all $1\leq j\leq m$. Note that the
set ${\mathcal{Z}}$ is generated independently of the chain.

We illustrate here the use of ${\mathcal{Z}}$ to run the ABSL sampler. If the
chain’s proposal at $t$th iteration, ${\mathbf{\theta}}^{*}$, is identical to
one element, say $\xi_{h}\in{\mathcal{Z}}$, and $m$ is large, then we would
not need to generate ${\mathbf{y}}_{1},\ldots,{\mathbf{y}}_{m}\sim
f({\mathbf{y}}|{\mathbf{\theta}}^{*})$ since we already have the corresponding
pseudo-data statistics vectors $s_{h}$ which can be used to estimate
$\mu({\mathbf{\theta}}^{*})$, $\Sigma({\mathbf{\theta}}^{*})$ and thus
$\mbox{SL}({\mathbf{\theta}}^{*})$. While the intuition is attractive, it is
impractical to faithfully implement it. For instance, using a large value for
$m$ when creating ${\mathcal{Z}}$ might still be too costly and an exact match
with an element in the reference set is unattainable when the parameter space
is continuous. However, if ${\mathcal{Z}}$ contains enough $\xi$-values that
are close enough to ${\mathbf{\theta}}^{*}$, one can still use them for
estimating $\mbox{SL}({\mathbf{\theta}}^{*})$. Levi & Craiu (2022) build the
reference set with $m=1$ and propose the use of $K$-nearest neighbours (kNN)
estimators for $\mu({\mathbf{\theta}}^{*})$, $\Sigma({\mathbf{\theta}}^{*})$

$\begin{split}\tilde{\mu}({\mathbf{\theta}}^{*})&=\frac{\sum_{h=1}^{H}[W_{h}({\mathbf{\theta}}^{*}){1\over
m}\sum_{j=1}^{m}s_{h}^{(j)}]}{\sum_{h=1}^{H}W_{h}({\mathbf{\theta}}^{*})},\\\
\tilde{\Sigma}({\mathbf{\theta}}^{*})&=\frac{\sum_{h=1}^{H}[W_{h}({\mathbf{\theta}}^{*}){1\over
m}\sum_{j=1}^{m}(s_{h}^{(j)}-\hat{\mu}_{{\mathbf{\theta}}^{*}})(s_{h}^{(j)}-\hat{\mu}_{{\mathbf{\theta}}^{*}})^{T}]}{\sum_{h=1}^{H}W_{h}({\mathbf{\theta}}^{*})}.\end{split}$
((16))

where $W_{h}({\mathbf{\theta}}^{*})=1$ or
$W_{h}({\mathbf{\theta}}^{*})=1-\|\xi_{h}-{\mathbf{\theta}}^{*}\|/\|\xi^{*}-{\mathbf{\theta}}^{*}\|$
and $\xi^{*}=\max_{\xi\in{\mathcal{Z}}}\|\xi-{\mathbf{\theta}}^{*}\|$, i.e. is
the point in ${\mathcal{Z}}$ that is furthest away from
${\mathbf{\theta}}^{*}$. If $H$ is large, it is likely that most of its
elements will contribute little or not at all to the estimators ((16)).
Instead of summing over all the $H$ elements in ${\mathcal{Z}}$, it is then
advisable to use only the $K$ $\xi$’s that are closest to
${\mathbf{\theta}}^{*}$, where $K$ is user-defined and depends on the
available computational power. In our numerical experiments we have used
$W_{h}=1$ for all $1\leq h\leq K$ and $K=\lfloor\sqrt{H}\rfloor$. Clearly, the
estimators in ((16)) are consistent due to the properties of kNN estimators,
but are not unbiased, so pseudo-marginal arguments cannot be invoked to
justify the approach. Validity is demonstrated theoretically by showing that
the perturbation induced when using the modified transition kernel can be
controlled using the user-specified tuning parameters of the sampler (see
section 6 in Levi & Craiu, 2022, for details).

A similar approach is used by Levi & Craiu (2022) for the ABC MCMC chain that
targets the marginal posterior density of ${\mathbf{\theta}}$ resulting from
((15)), $\pi({\mathbf{\theta}}|S({\mathbf{y}}_{0}))\propto
p({\mathbf{\theta}})\Pr(d(S({\mathbf{y}}),S({\mathbf{y}}_{0}))<{\epsilon})|{\mathbf{\theta}})$.
Instead of using an unbiased estimator for
$\Pr(d(S({\mathbf{y}}),S({\mathbf{y}}_{0}))<{\epsilon})|{\mathbf{\theta}})$
which would require multiple pseudo-data generated from
$f({\mathbf{y}}|{\mathbf{\theta}})$, they construct the kNN-based estimator
from the collection ${\mathcal{Z}}$.

In the next section, we compare numerically the methods discussed so far using
a couple of examples.

##  5 Numerical experiments

In this section we present the performance of the discussed algorithms on two
models: logistic regression and stochastic volatility. We compare the accuracy
and computational efficiency of the described methods with a couple of
benchmark MCMC algorithms that are widely known to perform very well in these
cases.  Specifically, we measure the perfomance of the methods presented in
this paper against the Polya-Gamma (PG) sampler (Polson et al., 2013) for the
logistic regression, and the sequential ABC (ABC SMC) of Sisson et al. (2007),
Lee (2012) for the stochastic volatility model. The former is customized for
logistic regression and for the latter the length of ${\epsilon}$ sequence is
set at $15$.

###  5.1 Description of the simulation settings

The following variations of the algorithms described in previous sections are
implemented.

PG_DAC_J

DAC algorithm with PG sampler that follows the setup in Scott et al. (2016)
using ((5)). The samples from each batch are combined proportionally to the
inverse covariance matrices. $J$ denotes the number of batches.

RW_SS

Subsampling using Quiroz et al. (2018) with a Random Walk (RW) transition
kernel. There are four variations corresponding to pairing parameter or data
expansion with random or correlated index selection.

RW_SS_P_R_m

\- Parameter expansion and random index selection

RW_SS_D_R_K_m

\- Data expansion and random index selection

RW_SS_P_C_m

\- Parameter expansion and correlated index selection, the correlation $\rho$
is set at $\rho=0.9999$

RW_SS_D_C_K_m

\- Data expansion and correlated index selection, the correlation $\rho$ is
set at $\rho=0.9999$

Note that $m$ and $K$ indicate the number of observations that will be
evaluated with the actual log-likelihood and the number of clusters
respectively.

RW_CO_K_f

Coreset method for logistic regression proposed by Huggins et al. (2016). The
number of clusters and proportion of non-zero weights out of $N$ are specified
by $K$ and $f$, respectively. Note that the radius $R$ is calculated from the
average sum of squared distances within each cluster as suggested in Huggins
et al. (2016). The coreset is used with a random walk Metropolis (RWM)
sampling algorithm.

RW_ABC

ABC MCMC algorithm using a RWM transition kernel for target ((15)). Only one
pseudo data set is generated for each proposal ${\mathbf{\theta}}^{*}$.

RW_AABC

Approximate ABC MCMC algorithm proposed by Levi & Craiu (2022). Proposals from
the history of the chain are used to estimate the likelihood using the
k-nearest-neighbour approach with uniform weights. Only one pseudo data set is
generated at every iteration.

RW_BSL_m

BSL MCMC algorithm with a RWM transition kernel. The distribution of the
summary statistics is approximated by a Gaussian. The mean and covariance of
the distribution is estimated by generating $m$ pseudo data sets at each
proposed ${\mathbf{\theta}}^{*}$.

RW_ABSL

Approximate BSL MCMC algorithm proposed by Levi & Craiu (2022). Past results
are used to estimate the mean and covariance of the summary statistics
distribution using k-nearest-neighbour approach with uniform weights. Only one
pseudo data set is generated at every iteration. See more details in the
Supplementary Materials.

With the exception of PG and ABC SMC, all the approximate samplers rely on a
random walk Metropolis (RWM) kernel with Gaussian proposals to ensure
consistency and comparability. The RWM kernels used here benefit from a
finite-adaptation strategy, in which the covariance of the proposal is
modified using the method of Haario et al. (2001), during the first $B$
iterations that make up the burn-in period, and are kept fixed after that. The
ABC, AABC and ABC SMC samplers depend on the threshold ${\epsilon}$ and the
ingredients needed to compute the distance $d$ in ((15)), and are determined
following preliminary simulations, as detailed in Levi & Craiu (2022).
Additional details about each sampling design are provided in the
Supplementary Materials.

For standard MCMC samplers, their performance comparison is often reported in
terms of the effective sample size (ESS) per second of central processing unit
(CPU) time, denoted ESS/cpu. The ESS is interpreted as the number of
independent samples that would yield the same variance of the Monte Carlo
estimator. A higher ESS value indicates a more efficient MCMC sampling
algorithm, since it has been directly linked with the algorithm’s
computational uncertainty (e.g., Gong & Flegal, 2016, Vats et al., 2019). The
CPU time directly measures the computational cost in seconds so ESS/cpu can be
interpreted as a sampler’s speed of generating information about the target.

All the samplers discussed here will target a distribution different than the
posterior of interest. Thus, in order to fully compare these sampler, one must
consider the errors incurred because of this shift. Therefore, in addition to
metrics designed to measure the efficiency of a regular MCMC sampler, such as
ESS/cpu, we also use $R=50$ independent replicates to produce estimates of
Monte Carlo bias and variance. This led us to two measures of efficiency that
are used to convey the performance of each method: the root mean square error
(RMSE), and the ESS/cpu.

To fix the notions, let ${\mathbf{\theta}}^{rs}_{(t)}$ represent the posterior
samples from replicate $1\leq r\leq R$, iteration $B\leq t\leq M$ (only draws
obtained after burn-in are retained) and parameter component $1\leq s\leq d$.
Similarly, $\tilde{\mathbf{\theta}}^{rs}_{(t)}$ are posterior draws from the
benchmark chain (only draws obtained after the burn-in period are retained).
We also let ${\mathbf{\theta}}_{true}^{s}$ denote the true parameter value
that was used to generate the data. The following quantities are used for
comparing computational efficiency

$\begin{split}&\mbox{Bias}^{2}=\mbox{Mean}_{s}\left(\left(\mbox{Mean}_{tr}({\mathbf{\theta}}^{rs}_{(t)})-{\mathbf{\theta}}_{true}^{s}\right)^{2}\right),\\\
&\mbox{VAR}=\mbox{Mean}_{s}(\mbox{Var}_{r}(\mbox{Mean}_{t}({\mathbf{\theta}}^{rs}_{(t)}))),\\\
&\mbox{RMSE}=\sqrt{\mbox{Bias}^{2}+\mbox{VAR}},\end{split}$

where $\mbox{Mean}_{t}(a^{st})$ is defined as the average of $\\{a^{st}\\}$
over index $t$ and, similarly, $\mbox{Var}_{t}(a^{st})$ and
$\mbox{Cov}_{t}(a^{st})$ denote the sample variance and covariance,
respectively.

Using the coda library in R we compute ESS for each replicate and parameter’s
component $\mbox{ESS}^{rs}$. Letting $CPU^{r}$ denote the total CPU time used
for producing the MCMC samples in replicate $r$, we define ESS/cpu as:

$\begin{split}&\mbox{ESS/cpu}=Mean_{rs}(\mbox{ESS}^{rs}/CPU^{r}).\end{split}$

Note that we consider the average over all parameters and replicates.
Generally, a sampler with a higher ESS/cpu is preferred because it yields a
higher amount of information per unit of time. The ABC SMC sampler produces
independent draws so its ESS is equal to the number of particles.

Finally, in order to frame the comparison in terms of unit-free measures, we
report the performance relative to the benchmark samplers. This means that
once we compute the RMSE for say method A, RMSEA, we report instead
$\mbox{RMSE}_{A}/\mbox{RMSE}_{Bench}$, where the denominator is the benchmark
sampler’s RMSE. Similarly, we also report the relative ESS/cpu performance.

###  5.2 Logistic Regression

This set of simulations contains the standard setting for the logistic
regression model. The design $N\times d$ matrix $X$ is generated by simulating
each variable independently from $Unif(0,1)$. The left most column is a column
of $1$s (intercept). For $i=1,\ldots,N$, $Y_{i}$ is Bernoulli with ${\rm
Pr}(Y_{i}=1)=\mbox{logistic}(x_{i}\cdot{\mathbf{\theta}}_{true})$. We
considered two values for the sample size $N$: $1,000$ and $10,000$; and two
sets of parameters:

  * •

$d=2$ with the true parameter of ${\mathbf{\theta}}_{true}=(-2,2)$

  * •

$d=10$ with the true parameter of
${\mathbf{\theta}}_{true}=(-2,2,-3,4,1,2,-3,-4,2,1)/3$

We set the prior distribution to be
$p({\mathbf{\theta}})\sim\mathcal{N}(0,4I_{d})$, where ${\mathbf{\theta}}\in
R^{d}$ and $I_{d}$ is $d\times d$ identity matrix.

All the samplers are run for $M=55,000$ iterations with burn-in set at
$B=15,000$.

For the DAC sampler we consider three values for the number of batches
$J=2,3,5$. We set $K=4$ for the Coreset method and compare four values for the
fraction $f=0.5,0.1,0.05,0.01$. Finally, the values of the tuning parameters
for the subsampling method were also variable. Specifically, the the number of
data clusters, $K\in\\{10,50\\}$, and the size of the subsample,
$m\in\\{20,100\\}$. Generally, $K$ and $m$ will depend on the sample size $N$
and parameter dimension $d$. The recommendation is to select larger values for
data expansion than parameter expansion. In addition, using correlated indices
require smaller values for these hyper-parameters. Figure 1 and 2 present the
simulation results for the scenarios with $N=10,000$ and $d\in\\{2,10\\}$. In
the supplemental material we include two additional scenarios $N=1,000,d=2$
and $N=10,000,d=2$. The height of the bars represents the value of the
relative measure and we add the dashed line at 1 to make it easier to separate
performance improvements from deteriorations.

Figure 1: Logistic model: Relative RMSE (top row) and ESS/cpu (bottom row)
when $N=10,000$ and $d=2$ for DAC-based samplers (left column), Coreset-based
samplers (center column), and Subsampling-based samplers (right column).
Figure 2: Logistic model: Relative RMSE (top row) and ESS/cpu (bottom row)
when $N=10,000$ and $d=10$ for DAC-based samplers (left column), Coreset-based
samplers (center column), and Subsampling-based samplers (right column).

From Figures 1 and 2 a few lessons emerge. Combining PG with DAC produces good
results, likely because the Gaussian approximation is accurate for such a
large sample. The ESS/cpu grows with the number of batches.

The performance of the coreset-based algorithm yield a relatively high RMSE.
To shed some light on this performance we can recover from Huggins et al.
(2016) the discrepancy, ${\epsilon}$, between the original likelihood and the
coreset one, as a function of the coreset size, mean sensitivity, and
parameter dimension for $\delta=0.10$.

Table 1 shows the average (over 50 replicates) discrepancy ${\epsilon}$ for
different values of the sample size $N$, parameter dimension $d$, and data
fraction $f$ divided by the average maximum value of the full data likelihood.

Table 1:  Coreset (logistic model): Relative average discrepancy ${\epsilon}$
for different values of sample size $N$, parameter dimension $d$, and data
fraction $f$. The numbers represent the average discrepancy divided by the
average maximum value of the full data likelihood.

| $N=1,000$ | $N=10,000$
---|---|---
Fraction | $d=2$ | $d=10$ | $d=2$ | $d=10$
$f=0.50$ | 3.217 | 3.602 | 0.685 | 1.160
$f=0.10$ | 7.845 | 9.717 | 3.570 | 3.154
$f=0.05$ | 11.992 | 13.972 | 6.138 | 4.461
$f=0.01$ | 28.595 | 31.663 | 16.383 | 10.091

It is not surprising that ${\epsilon}$ increases as the fraction (i.e, the
coreset size) decreases, but we also can see that the discrepancy is generally
quite large and this explains the poor performance of the sampler. The ESS/cpu
measure beats PG only when using 1% of the samples, but this comes at the
expense of a vastly inflated RMSE.

Overall, subsampling techniques show good results with very high ESS/cpu
without sacrificing the accuracy of the posterior, when $d=2$. The logistic
posterior tends to be unimodal so the parameter expansion methodology is more
suitable and clearly a larger concentration is achieved for $d=2$ than for
$d=10$. The deterioration of the performance is clearly visible for $d=10$
although the method still controls the RMSE at the PG level. Since the data do
not exhibit any clusters, it is not surprising that the data expansion
techniques are not competitive to the parameter expansion ones.

Note that the computational time for the calculation of the log-likelihood can
be significantly reduced using the vectorization trick available in R. This
method allows much faster calculation by executing operations on the entire
vectors of data instead of using a ’for’ loop that goes through all the $N$
records one by one. This technique enabled us to increase the sample size to
$100,000$. The comparison of the samplers using the vectorization-induced
speed-up can be found in the Supporting Materials.

####  5.2.1 German Credit Data

This concerns data with a sample size that is not exceedingly large, but the
dimension of the parameter is higher than we have considered so far.
Specifically, the german credit data consist of $1,000$ records and $49$
predictors including the intercept (see Biswas et al. (2019) for more
information). Most predictors are dummy variables taking only 0 and 1 values.
The target/response is binary, with 70% of them being cases, so the response
variable is quite balanced. Logistic regression is implemented to predict
$\Pr(Y=1)$ from the set of features. Before fitting the model, we transform
all the quantitative features by subtracting the minimum value and dividing by
the range so their values are in the $[0,1]$ interval.

All the samplers are run for $N=100,000$ iterations, burn-in is $B=50,000$ and
adaptation occurs every $500$ chain updates. The performance of the samplers
is presented in Figure 3. Note that the absolute value of the biases reported
in the top row panels are calculated with respect to the maximum likelihood
estimates, as the true parameter values are not known. We refer the reader to
the Supplementary Materials for additional metrics and details.

Figure 3: German Credit data: Relative $|\mbox{Bias}|$ (top row)and ESS/cpu
(bottom row) for DAC-based samplers (left column), Coreset-based samplers
(center column), and Subsampling-based samplers (right column).

The results are similar to the ones obtained in the previous subsection, but
some additional elements emerge. We can see that with 5 batches the DAC
approach is losing a bit in terms of bias and even more on the ESS/cpu side.
None of the subsampling-based methods (including coreset-based) can compete
with the PG sampler, likely because the signal to noise ratio is altered too
much when implementing any of these methods. We should also recognize that PG
is a Gibbs sampler which, unlike RWM samplers, will move at every iteration.
This makes a bigger difference when the parameter space has large dimensions
since then the RWM chain often gets stuck, especially if the posterior
exhibits strong dependence.

Based on these numerical experiments, we conclude that with a very large
sample size the first choice would be to use a DAC technique as long as the
Gaussian approximation is likely to be accurate. The latter assessment will
have to take into account the number of parameters and the nature of the model
and data. If the Gaussian approximation is unsuitable, the subsampling methods
can be used. The user will need to decide if the posterior is likely to be
concentrated, so that they can use a parameter expansion, or the data exhibits
multiple clusters, in which case a data expansion is needed. In the latter
case, an exploratory analysis is recommended to determine reasonable values
for the number of centroids, $K$. The size of the subsample $m$ is typically
decided based on the computational power available at the time of the analysis
- we recommend using the largest possible value that can be handled by the
system.

###  5.3 Stochastic Volatility

When analyzing stationary time series, it is frequently observed that there
are periods of high and low volatility, a phenomenon known as volatility
clustering (see for example Lux & Marchesi, 2000). One way to model such
behaviour is through a Stochastic Volatility (SV) model, where variances of
the observed time series depend on hidden states that themselves form a
stationary time series. We work with the following model which is indexed by
parameter
$\theta=({\mathbf{\theta}}_{1},{\mathbf{\theta}}_{2},{\mathbf{\theta}}_{3})$:

$\begin{split}&x_{1}\sim\mathcal{N}(0,1/(1-{\mathbf{\theta}}_{1}^{2}));{\hskip
8.61108pt}v_{i}\overset{iid}{\sim}\mathcal{N}(0,1);{\hskip
8.61108pt}w_{i}\overset{iid}{\sim}\mathcal{N}(0,1);{\hskip
8.61108pt}i=\\{1,\ldots,N\\},\\\
&x_{i}={\mathbf{\theta}}_{1}x_{i-1}+v_{i};{\hskip
8.61108pt}i=\\{2,\ldots,N\\},\\\
&y_{i}=\sqrt{\exp[{\mathbf{\theta}}_{2}+\exp({\mathbf{\theta}}_{3})x_{i}]}w_{i};{\hskip
8.61108pt}i=\\{1,\ldots,N\\}.\end{split}$ ((17))

Only data ${\mathbf{y}}=(y_{1},\ldots,y_{N})$ are observed, and
$(x_{1},\ldots,x_{N})$ are latent/hidden states. The parameter
${\mathbf{\theta}}_{1}\in(-1,1)$ controls the auto-correlation of hidden
states, while ${\mathbf{\theta}}_{2}$ and ${\mathbf{\theta}}_{3}$ are
unrestricted and relate to the hidden states influence on the variability of
the observed series. Given a hidden state, the distribution of the observed
variable is Gaussian. We introduce the following priors, independently for
each parameter:

$\begin{split}&{\mathbf{\theta}}_{1}\sim Unif[0,1],\\\
&{\mathbf{\theta}}_{2}\sim\mathcal{N}(0,1),\\\
&{\mathbf{\theta}}_{3}\sim\mathcal{N}(0,1).\end{split}$ ((18))

We set the true parameters to ${\mathbf{\theta}}_{true}=(0.95,-2,-1)$ and
consider three lengths of the time series $N=100$, $500$ and $1,000$. Note
that the model does not admit a closed form log-likelihood but allows
simulations of pseudo data sets. Therefore, for this model we only consider
simulation-based ABC samplers: AABC, BSL, ABSL, and the benchmark SMC ABC.The
summary statistics used for all the samplers is $S({\mathbf{y}})\in{\bf
R}^{6}$ and has the following components:

(C1)

Average of ${\mathbf{y}}^{2}$,

(C2)

Standard deviation of ${\mathbf{y}}^{2}$,

(C3)

Sum of the first 5 auto-correlations of ${\mathbf{y}}^{2}$,

(C4)

Sum of the first 5 auto-correlations of binary series $\\{{\bf
1}_{\\{y_{i}^{2}<\mbox{quantile}({\mathbf{y}}^{2},0.1)\\}}\\}_{i=1}^{N}$,

(C5)

Sum of the first 5 auto-correlations of binary series $\\{{\bf
1}_{\\{y_{i}^{2}<\mbox{quantile}({\mathbf{y}}^{2},0.5)\\}}\\}_{i=1}^{N}$,

(C6)

Sum of the first 5 auto-correlations of binary series $\\{{\bf
1}_{\\{y_{i}^{2}<\mbox{quantile}({\mathbf{y}}^{2},0.9)\\}}\\}_{i=1}^{N}$.

The $\mbox{quantile}({\mathbf{y}},\tau)$ is defined as the $\tau$-th quantile
of the sequence ${\mathbf{y}}$. We focus here on ${\mathbf{y}}^{2}$ and its
auto-correlations because the model parameters only affect its variability;
the auto-correlation of ${\mathbf{y}}$ is zero for any lag. The components
(C4)-(C6) have been considered because the auto-correlations of those binary
series, defined under different quantiles, are useful in characterizing a time
series (Schmitt et al., 2015, Dette et al., 2015). The ABC, AABC, BSL and ABSL
samplers are run for $M=55,000$ iterations. The burn-in period is of length
$B=15,000$, with adaptation taking place every other $200$ iterations.

Figures 4 and 5 present the simulation results when $N=500$ and, respectively,
$N=1000$. The ABC and BSL samplers exhibit loss in terms of both RMSE and
ESS/cpu when compared to the benchmark. The BSL is more costly since we
generate 20 pseudo-data sets at each iteration. Not surprisingly, using pre-
computation designs reduces the CPU time so we see a bump in efficiency for
AABC and ABSL. Less obvious is the reduction in RMSE which is due to the
increase the number of pseudo-data one can use while still saving
computational time and the higher acceptance rate. These findings mirror those
of Levi & Craiu (2022) and we refer the reader to that paper for more in-depth
explanations. In this example ABC-based samplers outperform the BSL ones. The
likely reason is that the Gaussian approximation on which the BSL relies is
not accurate for this choice of the summary statistic, $S({\mathbf{y}})$.

Figure 4: SV model: Relative RMSE (top row) and ESS/cpu (bottom row) when
$N=500$ for ABC-based samplers (left column) and BSL-based samplers (right
column).  Figure 5: SV model: Relative RMSE (top row) and ESS/cpu (bottom row)
when $N=1000$ for ABC-based samplers (left column) and BSL-based samplers
(right column).

Overall, we find reasons for cautious optimism in these numerical results.
They show that careful and controlled injection of noise in the transition
kernel can bring real practical benefits.

##  6 Conclusion and future directions

The Bayesian computational community finds itself at an inflexion point.
Traditional MCMC computation is no longer tenable for complex problems. The
new ideas and developments discussed here reduce significantly the
computational costs or bypass the intractability of the likelihood, but
introduce additional layers of approximation. The latter requires a careful
theoretical analysis to make sure that incurred errors are realistically
controllable via tuning parameters.

Complex models are often defined using high-dimensional parameters. MCMC
methods sample efficiently high-dimensional spaces as long as there are no
bottlenecks or regions of small probability that the chain has difficulty
traversing. Adaptive MCMC methods (Andrieu & Thoms, 2008, Hoffman et al.,
2014, Yang et al., 2019, Pompe et al., 2020) have been proven effective for
sampling in high-dimensional spaces with unfriendly geometries. Injecting
adaptive ideas into the world of sampling with intractable targets is hindered
by stringent conditions that need to be satisfied by an adaptive transition
kernel, e.g. the containment condition (Bai et al., 2011, Latuszyński &
Rosenthal, 2014). Some inroads have been made into eliminating the latter in
Craiu et al. (2015) and Rosenthal & Yang (2018) so we expect to see more
adaptive designs permeating in pseudodata-generation-type samplers.

Constraints on paper length and considerations of subject matter consistency
have prevented us from discussing methods that do not rely on MCMC sampling to
perform Bayesian inference such as variational Bayes (Blei et al., 2017) or
integrated nested Laplace approximation (Rue et al., 2017). These are active
research threads that continue to develop rapidly under the impetus provided
by the expansion of data science, explosive growth of machine learning methods
and other computationally demanding domains of information processing.
Creative intertwining of most of the ideas or methods mentioned in this paper
will likely continue well into the future, but we believe that entirely new
perspectives are also necessary in order to create the automatization of
computation that is required if widespread use of Bayesian methods is to be
seen in the 21st century.

## Acknowledgements

This research has been funded by NSERC of Canada. The authors thank Nancy Reid
for the invitation to write this article, Alicia Carriquiri for guidance in
defining the paper’s scope, and an anonymous referee for a number of
suggestions that have led to important improvements.

## References

  * Andrieu & Roberts (2009) Andrieu C, Roberts GO. 2009. The pseudo-marginal approach for efficient Monte Carlo computations. The Annals of Statistics 37(2):697–725
  * Andrieu & Thoms (2008) Andrieu C, Thoms J. 2008. A tutorial on adaptive MCMC. Statist. Comput. 18:343–373
  * Andrieu & Vihola (2015) Andrieu C, Vihola M. 2015. Convergence properties of pseudo-marginal Markov chain Monte Carlo algorithms. The Annals of Applied Probability 25(2):1030–1077
  * Bai et al. (2011) Bai Y, Roberts GO, Rosenthal JS. 2011. On the containment condition for adaptive Markov chain Monte Carlo algorithms. Adv. Appl. Stat 21(1):1–54
  * Bardenet et al. (2014) Bardenet R, Doucet A, Holmes C. 2014. Towards scaling up Markov chain Monte Carlo: an adaptive subsampling approach, In International Conference on Machine Learning, pp. 405–413, PMLR
  * Bardenet et al. (2017) Bardenet R, Doucet A, Holmes C. 2017. On Markov chain Monte Carlo methods for tall data. The Journal of Machine Learning Research 18(1):1515–1557
  * Beaumont et al. (2002) Beaumont MA, Zhang W, Balding DJ. 2002. Approximate Bayesian computation in population genetics. Genetics 162(4):2025–2035
  * Biswas et al. (2019) Biswas N, Jacob PE, Vanetti P. 2019. Estimating convergence of Markov chains with L-lag couplings, In Advances in Neural Information Processing Systems, pp. 7389–7399
  * Blei et al. (2017) Blei DM, Kucukelbir A, McAuliffe JD. 2017. Variational inference: A review for statisticians. Journal of the American statistical Association 112(518):859–877
  * Campbell & Broderick (2019) Campbell T, Broderick T. 2019. Automated scalable Bayesian inference via hilbert coresets. The Journal of Machine Learning Research 20(1):551–588
  * Changye & Robert (2019) Changye W, Robert CP. 2019. Parallelising MCMC via random forests. arXiv preprint arXiv:1911.09698
  * Chipman et al. (2010) Chipman HA, George EI, McCulloch RE. 2010. Bart: Bayesian additive regression trees. The Annals of Applied Statistics 4(1):266–298
  * Craiu et al. (2015) Craiu RV, Gray L, Latuszyński K, Madras N, Roberts GO, Rosenthal JS. 2015\. Stability of adversarial Markov chains, with an application to adaptive MCMC algorithms. The Annals of Applied Probability 25(6):3592–3623
  * Craiu & Meng (2005) Craiu RV, Meng XL. 2005. Multiprocess parallel antithetic coupling for backward and forward Markov chain Monte Carlo. The Annals of Statistics 33(2):661–697
  * Craiu et al. (2009) Craiu RV, Rosenthal JS, Yang C. 2009. Learn from thy neighbor: Parallel-chain adaptive and regional MCMC. Journal of the American Statistical Association 104:1454–1466
  * Cui et al. (2018) Cui T, Peeters L, Pagendam D, Pickett T, Jin H, et al. 2018. Emulator-enabled approximate Bayesian computation (abc) and uncertainty analysis for computationally expensive groundwater models. Journal of hydrology 564:191–207
  * Deligiannidis et al. (2018) Deligiannidis G, Doucet A, Pitt MK. 2018. The correlated pseudomarginal method. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 80(5):839–870
  * Dette et al. (2015) Dette H, Hallin M, Kley T, Volgushev S, et al. 2015. Of copulas, quantiles, ranks and spectra: An $l_{1}$-approach to spectral analysis. Bernoulli 21(2):781–831
  * Entezari et al. (2018) Entezari R, Craiu RV, Rosenthal JS. 2018. Likelihood inflating sampling algorithm. Canadian Journal of Statistics 46(1):147–175
  * Fearnhead & Prangle (2012) Fearnhead P, Prangle D. 2012. Constructing summary statistics for approximate Bayesian computation: semi-automatic approximate Bayesian computation. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 74(3):419–474
  * Filippi et al. (2013) Filippi S, Barnes CP, Cornebise J, Stumpf MP. 2013. On optimality of kernels for approximate Bayesian computation using sequential Monte Carlo. Statistical applications in genetics and molecular biology 12(1):87–107
  * Frigessi et al. (2000) Frigessi A, Gåsemyr J, Rue H. 2000. Antithetic coupling of two Gibbs sampler chains. The Annals of Statistics 28:1128–1149
  * Gallant & McCulloch (2009) Gallant AR, McCulloch RE. 2009. On the determination of general scientific models with application to asset pricing. Journal of the American Statistical Association 104(485):117–131
  * Gong & Flegal (2016) Gong L, Flegal JM. 2016. A practical sequential stopping rule for high-dimensional Markov chain Monte Carlo. Journal of Computational and Graphical Statistics 25(3):684–700
  * Gourieroux et al. (1993) Gourieroux C, Monfort A, Renault E. 1993. Indirect inference. Journal of applied econometrics 8(S1):S85–S118
  * Haario et al. (2001) Haario H, Saksman E, Tamminen J. 2001. An adaptive Metropolis algorithm. Bernoulli 7:223–242
  * Hastings (1970) Hastings WK. 1970. Monte Carlo Sampling Methods Using Markov chains and their Applications. Biometrika 57(1):97–109
  * Hoffman et al. (2014) Hoffman MD, Gelman A, et al. 2014. The no-u-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. J. Mach. Learn. Res. 15(1):1593–1623
  * Huggins et al. (2016) Huggins J, Campbell T, Broderick T. 2016. Coresets for scalable Bayesian logistic regression. Advances in Neural Information Processing Systems 29:4080–4088
  * Korattikara et al. (2014) Korattikara A, Chen Y, Welling M. 2014. Austerity in MCMC land: Cutting the Metropolis-Hastings budget, In International Conference on Machine Learning, pp. 181–189, PMLR
  * Latuszyński & Rosenthal (2014) Latuszyński K, Rosenthal JS. 2014. The containment condition and adapfail algorithms. Journal of Applied Probability 51(4):1189–1195
  * Lee (2012) Lee A. 2012. On the choice of MCMC kernels for approximate Bayesian computation with SMC samplers, In Proceedings of the 2012 Winter Simulation Conference (WSC), pp. 1–12, IEEE
  * Lee et al. (2012) Lee A, Andrieu C, Doucet A. 2012. Discussion of constructing summary statistics for approximate Bayesian computation: semi-automatic approximate Bayesian computation. JR Stat. Soc. Ser. B Stat. Methodol 74(3):449–450
  * Levi & Craiu (2022) Levi E, Craiu RV. 2022. Finding our way in the dark: Approximate MCMC for approximate Bayesian methods. Bayesian Analysis 17(1):193 –221
  * Lux & Marchesi (2000) Lux T, Marchesi M. 2000. Volatility clustering in financial markets: a microsimulation of interacting agents. International journal of theoretical and applied finance 3(04):675–702
  * Manousakas et al. (2020) Manousakas D, Xu Z, Mascolo C, Campbell T. 2020. Bayesian pseudocoresets. Advances in Neural Information Processing Systems 33:14950–14960
  * Marin et al. (2014) Marin JM, Pillai NS, Robert CP, Rousseau J. 2014. Relevant statistics for Bayesian model choice. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 76(5):833–859
  * Marjoram et al. (2003) Marjoram P, Molitor J, Plagnol V, Tavaré S. 2003. Markov chain Monte Carlo without likelihoods. Proceedings of the National Academy of Sciences 100(26):15324–15328
  * Metropolis et al. (1953) Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E. 1953. Equation of state calculations by fast computing machines. The journal of chemical physics 21(6):1087–1092
  * Neiswanger et al. (2013) Neiswanger W, Wang C, Xing E. 2013. Asymptotically exact, embarrassingly parallel MCMC. arXiv preprint arXiv:1311.4780
  * Nemeth & Sherlock (2018) Nemeth C, Sherlock C. 2018. Merging MCMC subposteriors through gaussian-process approximations. Bayesian Analysis 13(2):507–530
  * Oyebamiji et al. (2015) Oyebamiji OK, Edwards NR, Holden PB, Garthwaite PH, Schaphoff S, Gerten D. 2015\. Emulating global climate change impacts on crop yields. Statistical Modelling 15(6):499–525
  * Plumlee et al. (2021) Plumlee M, Asher TG, Chang W, Bilskie MV. 2021. High-fidelity hurricane surge forecasting using emulation and sequential experiments. The Annals of Applied Statistics 15(1):460–480
  * Polson et al. (2013) Polson NG, Scott JG, Windle J. 2013. Bayesian inference for logistic models using Pólya–Gamma latent variables. Journal of the American statistical Association 108(504):1339–1349
  * Pompe et al. (2020) Pompe E, Holmes C, Latuszyński K. 2020. A framework for adaptive MCMC targeting multimodal distributions. The Annals of Statistics 48(5):2930–2952
  * Pratola (2016) Pratola MT. 2016. Efficient Metropolis–Hastings proposal mechanisms for Bayesian regression tree models. Bayesian analysis 11(3):885–911
  * Price et al. (2018) Price LF, Drovandi CC, Lee A, Nott DJ. 2018. Bayesian synthetic likelihood. Journal of Computational and Graphical Statistics 27(1):1–11
  * Pritchard et al. (1999) Pritchard JK, Seielstad MT, Perez-Lezaun A, Feldman MW. 1999. Population growth of human Y chromosomes: a study of Y chromosome microsatellites. Molecular biology and evolution 16(12):1791–1798
  * Quiroz et al. (2018) Quiroz M, Kohn R, Villani M, Tran MN. 2018. Speeding up MCMC by efficient data subsampling. Journal of the American Statistical Association
  * Robert (2014) Robert CP. 2014. Bayesian computational tools. Annual Review of Statistics and Its Application 1:153–177
  * Roberts & Rosenthal (2009) Roberts GO, Rosenthal JS. 2009. Examples of adaptive MCMC. J. Comput. Graph. Statist. 18:349–367
  * Rosenthal & Yang (2018) Rosenthal JS, Yang J. 2018. Ergodicity of combocontinuous adaptive MCMC algorithms. Methodology and Computing in Applied Probability 20(2):535–551
  * Rue et al. (2017) Rue H, Riebler A, Sorbye S, Illian J, Simpson D, Lindgren F. 2017. Bayesian computing with INLA: a review. Annual Review of Statistics and its Application 4:395–421
  * Schmitt et al. (2015) Schmitt TA, Schäfer R, Dette H, Guhr T. 2015. Quantile correlations: Uncovering temporal dependencies in financial time series. International Journal of Theoretical and Applied Finance 18(07):1550044
  * Scott et al. (2016) Scott SL, Blocker AW, Bonassi FV, Chipman HA, George EI, McCulloch RE. 2016. Bayes and big data: The consensus Monte Carlo algorithm. International Journal of Management Science and Engineering Management 11(2):78–88
  * Sisson et al. (2018a) Sisson S, Fan Y, Beaumont M. 2018a. Overview of approximate Bayesian computation. arXiv preprint arXiv:1802.09720
  * Sisson et al. (2018b) Sisson SA, Fan Y, Beaumont M. 2018b. Handbook of approximate Bayesian computation. Chapman and Hall/CRC
  * Sisson et al. (2007) Sisson SA, Fan Y, Tanaka MM. 2007. Sequential Monte Carlo without likelihoods. Proceedings of the National Academy of Sciences 104(6):1760–1765
  * Smith Jr (1993) Smith Jr AA. 1993. Estimating nonlinear time-series models using simulated vector autoregressions. Journal of Applied Econometrics 8(S1):S63–S84
  * Tavaré et al. (1997) Tavaré S, Balding DJ, Griffiths RC, Donnelly P. 1997. Inferring coalescence times from DNA sequence data. Genetics 145(2):505–518
  * Vats et al. (2019) Vats D, Flegal JM, Jones GL. 2019. Multivariate output analysis for Markov chain Monte Carlo. Biometrika 106(2):321–337
  * Von Neumann & Ulam (1951) Von Neumann J, Ulam S. 1951. Monte carlo method. National Bureau of Standards Applied Mathematics Series 12(1951):36
  * Wang & Dunson (2013) Wang X, Dunson DB. 2013. Parallelizing MCMC via Weierstrass sampler. arXiv preprint arXiv:1312.4605
  * Wood (2010) Wood SN. 2010. Statistical inference for noisy nonlinear ecological dynamic systems. Nature 466(7310):1102
  * Yang et al. (2019) Yang J, Levi E, Craiu RV, Rosenthal JS. 2019. Adaptive component-wise multiple-try Metropolis sampling. Journal of Computational and Graphical Statistics 28(2):276–289