Bayesian computation via empirical likelihood

Introduction

Bayesian statistical inference cannot easily operate when the likelihood function associated with the data is not entirely known, or cannot be computed in a manageable time, as is the case in most population genetic models (1, 2, 3). The fundamental reason for this difficulty with population genetics is that the statistical model associated with coalescent data needs to integrate over trees of high complexity. Similar computational issues with the likelihood function often occur in hidden Markov and other dynamic models (4). In those settings, traditional approximation tools based on stochastic simulation (5) are unavailable or unreliable. Indeed, the complexity of the latent structure defining the likelihood makes simulation of such structures too unstable to be trusted. Such settings call for alternative and often cruder approximations. The ABC methodology (1, 6) is a popular solution that bypasses the computation of the likelihood function (see (7) and (8) for surveys); (9) validate a conditional version of ABC that applies to hierarchical Bayes models in a wide generality.

The fast and polytomous development of the ABC algorithm is indicated by the rising literature in the domain, at both the methodological and the application levels. For instance, a whole new area of population genetic modelling (10, 8) has been explored thanks to the availability of such methods. However, both practitioners and theoreticians show a reluctance in adopting ABC, as some doubt about the validaty of the method (11, 12, 13). We propose in this paper to supplement the ABC approach with a generic and convergent likelihood approximation called the empirical likelihood that validates the new Bayesian computational technique as a convergent inferential method when the number of observations grows to infinity. The empirical likelihood perspective, introduced by (14), is a robust statistical approach that does not require the specification of the likelihood function. However, while it does not appear to have been used before in the settings that now rely on ABC, this data analysis method also is a broadly (albeit not universally) applicable and often fast approach which approximation differs from the one found in ABC algorithms, even though both are rooted in non-parametric statistics. Therefore, this methodology can be used both as a solution per se and as a benchmark against which to test the ABC output in many cases. This paper introduces the BC${}_{\text{el}}$ algorithm and illustrates its performances on selected representative examples, comparing the outcome with the true posterior density whenever available, and with an ABC approximation (15) otherwise.

Statistical Methods

The primary purpose of the ABC algorithm is to approximate simulation from the centerpiece of Bayesian inference, the posterior distribution $\pi(\boldsymbol{\theta}|\mathbf{y})\propto\pi(\boldsymbol{\theta})f(\mathbf{y}|\boldsymbol{\theta})$, when it cannot be numerically computed but when the distributions corresponding to both the prior $\pi$ and the likelihood $f$ can be simulated by manageable computer devices. The original (6) ABC algorithm is as follows: given a sample $\mathbf{y}$ of observations from the sample space, a sample of parameters $(\boldsymbol{\theta}_{1},\ldots,\boldsymbol{\theta}_{M})$ is produced by

The parameters of the ABC algorithm are the summary statistic $\eta$, the distance $\rho\{\cdot,\cdot\}$ and the tolerance level $\epsilon>0$. The basic justification of the ABC approximation is that, when using a sufficient statistic $\eta$, the distribution of the $\boldsymbol{\theta}_{i}$’s in the output of the algorithm converges to the genuine posterior distribution when $\epsilon$ goes to zero (16).

In practice, however, the statistic $\eta$ is non-sufficient and at best the approximation then converges to the genuine posterior $\pi(\boldsymbol{\theta}|\eta(\mathbf{y}))$ when $\epsilon$ goes to zero. This loss of information seems to be a necessary price to pay for the access to computable quantities. Furthermore, as argued below, it can be evaluated against the empirical likelihood approximation when the latter is available. Indeed, this approach does not require an information reduction through the choice of a tolerance zone or of a non-sufficient summary statistic.

2 Empirical likelihood

Owen (14) developed empirical likelihood techniques as a robust alternative to classical likelihood approaches. He demonstrated that, for some categories of statistical models, this approach inherited the convergence properties of standard likelihood at a much lower cost in assumptions about the model (as detailed in SI). While ABC algorithms do require a fully defined and often complex (hence debatable) statistical model, we argue that one should take advantage of the approximation device provided by the empirical likelihood to overcome most of the calibration difficulties encountered by ABC, at least as a convenient benchmark against which to test ABC solutions.

Assume that the dataset $\mathbf{y}$ is composed of $n$ independent replicates $\mathbf{y}=(\mathbf{y}_{1},\ldots,\mathbf{y}_{n})$ of some random vector $Y$ with density $f$. Rather than defining the likelihood from the density $f$ as usual, the empirical likelihood method starts by defining parameters of interest, $\boldsymbol{\theta}$, as functionals of $f$, for instance as moments of $f$, and it then profiles a non-parametric likelihood. More precisely, given a set of constraints of the form

where the dimension of $h$ sets the number of constraints unequivocally defining $\boldsymbol{\theta}$, the empirical likelihood is defined as

While the convergence of the empirical likelihood is well-established (see SI and (18)), the Bayesian use of empirical likelihoods has been little examined in the past, apart from a Monte Carlo study in (19), and a probabilistic one in (20).

3 BCelel{}_{\text{el}}

The most natural use of the empirical likelihood approximation is to act as if this representation was an exact likelihood, as in (19). Incorporating this perspective into a basic sampler leads to the following algorithm: Algorithm 2: Basic BC${}_{\text{el}}$ sampler

The output of BC${}_{\text{el}}$ is a sample of size $M$ of parameters with associated weights, which operate as an importance sampling output (5). Thus, the performance of the algorithm can be evaluated through the effective sample size

which approximates the size of an iid sample with the same variance as the original sample. As shown in (21), this quantity is always between $1$ (corresponding to a very poor outcome) and $M$ (corresponding to an iid perfect outcome).

Any algorithm that samples from a posterior distribution (e.g., MCMC, Population Monte Carlo, SMC algorithms, see (5)) may instead use the empirical likelihood as a proxy to the exact likelihood. For instance, to speed up the computation in the population genetics model introduced below, we resorted to the adaptive multiple importance sampling (AMIS, (22)) which is easy to parallelize on a multi-core computer. While the original target distribution is $\pi(\boldsymbol{\theta})L(\boldsymbol{\theta}|\mathbf{y})$ and the AMIS algorithm uses several (multivariate) Student’s $t$ distributions, denoted $t_{3}(\cdot|\mathbf{m},\boldsymbol{\Sigma})$ (i.e., with three degrees of freedom, centered at mean $\mathbf{m}$ and with covariance matrix $\boldsymbol{\Sigma}$), as an importance sampling distribution, the algorithm can be adapted to the empirical likelihood in a straightforward manner:

Algorithm 3: BC${}_{\text{el}}$-AMIS sampler

The output is thus a weighted sample $\boldsymbol{\theta}_{t,i}$ of size $MT_{M}$.

In contrast with ABC, BC${}_{\text{el}}$ algorithms do not usually require simulations from the sampling model, given that (2) provides a converging and non-parametric approximation of the likelihood function. This feature thus induces very significative improvements in computing time when the production of pseudo-datasets is time consuming, since solving (1) is usually immediate. This is for instance the case in population genetics and the last section of the SI provides an illustration of a huge improvement in comparison with ABC in two experiments described below. However, the improvement in speed may vanish in cases when producing an iid structure connected with the constraint (1) requires simulations from the sampling model, as illustrated by a counter-example for point processes in the SI, BC${}_{\text{el}}$ and ABC then breaking even in terms of computing time. Even though the computing time required by BC${}_{\text{el}}$ is customarily negligible when compared with ABC (or does not induce any extra time as in the point process counter-example), we further caution against opposing both approaches solely based on computing times, since they differ in the approximations they provide to a genuine Bayesian analysis and thus should be used in conjunction.

Using empirical likelihoods means there is no calibration of the many tuning parameters of ABC algorithms; most significantly, the likelihood ratio acts as a natural distance and importance weights produce an implicit and self-defined quantile on the original sample simulated from the prior. Notwithstanding these appealing qualities, BC${}_{\text{el}}$ still requires calibration, in particular in the choices of the parameterization of the sampling distribution and of the corresponding constraints (1) defining the empirical likelihood. Some examples are discussed below. The BC${}_{\text{el}}$-AMIS sampler also implies computing values of the prior density, up to a constant, which may be an hindrance in peculiar cases.

4 Composite likelihood in population genetics

ABC was first introduced by population geneticists (2, 10, 6) interested in statistical inference about the evolutionary history of species, as no likelihood-based approach existed apart from very rudimentary and hence unrealistic situations. This approach has been used in a number of biological studies (23, 24, 25), most of them including model choice. It is therefore crucial to obtain insights into the validity of such studies, particularly when they have economic, biological or ecological consequences (see, e.g., (26)). This can be achieved in part by running a comparison using BC${}_{\text{el}}$. Furthermore, given the major gain in computing time, due to the absence of replications of the data, BC${}_{\text{el}}$ can be applied to more complex biological models.

The main caveat when using the empirical likelihood in such settings is selecting a constraint (1) on the parameter of interest: in phylogeography, parameters like divergence dates, effective population sizes, mutation rates, etc., cannot be expressed as moments of the sampling distribution at a given locus. In particular, the data are not iid. However, when considering microsatellite loci with the stepwise mutation model (27) and evolutionary scenarios composed of divergence, we can derive the pairwise composite scores whose zero is the pairwise maximum likelihood estimator. Composite likelihoods have been proved consistent for estimating recombination rates, introducing an approximation of the dependency structure between nearby loci (28, 29, 30, 31). (See also (32) for composite likelihoods used in a likelihood-free setting.)

More specifically, we are approximating the intra-locus likelihood by a product over all pairs of genes in the sample at a given locus. Assuming that $y_{i}^{k}$ denotes the allele of the $i$-th gene in the sample at the $k$-th locus, and that $\phi$ is the vector of parameters, then the so-called pairwise likelihood of the data at the $k$-th locus, namely $\mathbf{y}^{k}$, is defined by

provide a constraint (1) in every way comparable to the score equations that give the maximum likelihood estimate and which is quite powerful for empirical likelihood derivations ((18), pp. 48–50). Hence the empirical likelihood of the full dataset $\mathbf{y}=(\mathbf{y}^{1},\ldots,\mathbf{y}^{K})$ given $\phi$ is computed with (2) under the (multidimensional) constraint that

where \rho(\theta)=\theta\big{/}\big{(}1+\theta+\sqrt{1+2\theta}\big{)}. If the two genes belong to individuals from demes having diverged at time $\tau$, then (33)

Results

2 Quantile distributions

Quantile distributions are defined by a closed-form quantile function $F^{-1}(p;\theta)$, and generally have no closed form for the density function. They are of great interest because of their flexibility and the ease with which they can be simulated by a simple inversion of the uniform distribution. A range of methods, including ABC approaches (10), have been proposed for estimation (see SI). We focus here on the four-parameter $g$-and-$k$ distribution, defined by its quantile function, denoted $Q(r;A,B,g,k)$ and equal to

where $z(r)$ is the $r$th standard normal quantile; the parameters $A,B,g$ and $k$ represent location, scale, skewness and kurtosis, respectively and $c$ measures the overall asymmetry (34, 35). We evaluated the BC${}_{\text{el}}$ algorithm for estimating this distribution using two values of $\theta=(A,B,g,k)$, two sets of priors and various combinations of $n,M$ and $p$, where $p$ is the number of percentiles used as constraints (see details in SI).

Figure 1 illustrates the true and fitted curves and a 95% credible region for the case with $n=100,M=5000$ and $p=3$. The corresponding posterior means (standard deviations) for the parameters $A,B,g,k$ were $3.08(0.14),1.12(0.23),1.79(0.25),0.41(0.12)$, respectively. The choice of sample size and number of constraints did not substantively affect the accuracy of parameter estimates, but the precision was noticeably improved for the larger sample size; see Figures S4, S5, and S6.

The accuracy and precision of the estimates were broadly comparable with the results obtained by (36) for the same distribution. Based on the whole experiment, the parameters $A$ and $B$ were well estimated in all cases, while the estimates of $g$ and $k$ were poorer for smaller values of $n$ and $M$. For small $n$ the estimates were more subject to the vagaries of sampling variation, whereas for small $M$ they were subject to the influence of a smaller number of very large importance weights. However, given the speed of BC${}_{\text{el}}$ compared with competing ABC algorithms, it is feasible to use even larger values of $M$ than considered in this experiment, since there is no requirement to simulate new datasets at each iteration. Moreover, this experiment is based on the very basic case of sampling from the prior; the results would be further improved by using an analogue of BC${}_{\text{el}}$-AMIS or alternative approaches similar to those proposed by (37) for ABC.

3 Dynamic models

In dynamic models, the difficulty with empirical likelihood stems from the dependence in the data $(y_{t})_{1\leq t\leq T}$. However, these models can be represented as transforms of unobserved iid sequences $(\epsilon_{t})_{1\leq t\leq T}$. The recovery of a converging empirical likelihood representation thus requires the reconstitution of the $\epsilon_{t}$’s as transforms of the data $\mathbf{y}$ and of the parameter $\theta$. Independence between the $\epsilon_{t}$’s is then at least as important as moment conditions. (This implies equivalent computing times for ABC and BC${}_{\text{el}}$.)

For instance, consider a simple dynamic model, namely the ARCH(1) model:

with a uniform prior over the simplex, i.e., $\alpha_{0},\alpha_{1}\geq 0$, $\alpha_{0}+\alpha_{1}\leq 1$. While this model can be handled by other means, since the likelihood function is available, we will compare here the behaviour of ABC and BC${}_{\text{el}}$ algorithms.

First, a natural empirical likelihood representation is based on the reconstituted $\epsilon_{t}$’s, defined as $y_{t}/\sigma_{t}$ when the $\sigma_{t}$’s are derived recursively. Figure 2 shows the result of estimating both parameters $\alpha_{0}$ and $\alpha_{1}$ when Algorithm ABC uses as summary statistics either the least square estimates of the parameters (derived from the series $(y_{t}^{2})$), which we label “optimal ABC” in connection with (38), or the mean of the series $\log(y_{t}^{2})$ supplemented by the two first autocorrelations of the series $(y_{t}^{2})$. The constraints in the empirical likelihood are either based on the three first moments of the reconstituted $\epsilon_{t}$’s or on the variance of those $\epsilon_{t}$’s complemented by both the correlations between the $y_{t-1}$’s and the $\epsilon_{t}$’s and between the $\epsilon_{t-1}$’s and the $\epsilon_{t}$’s. As seen from this experiment, BC${}_{\text{el}}$ does as well as the optimal ABC for the estimation of the parameters, but further brings a reduction in the variability of those estimates, thanks to the importance weights.

A much more complex dynamic model is the GARCH$(1,1)$ model of (39) that can be formalized as the observation of $y_{t}\sim\mathcal{N}(0,\sigma_{t}^{2})$ when

under the constraints $\alpha_{0},\alpha_{1},\beta_{1}>0$ and $\alpha_{1}+\beta_{1}<1$, that is, $y_{t}=\sigma_{t}\epsilon_{t}$. Given the constraints on the parameters, a natural prior is to choose an exponential distribution on $\alpha_{0}$, for instance an exponential $\mathcal{E}xp(1)$ distribution, and a Dirichlet $D_{3}(1,1,1)$ on $(\alpha_{1},\beta_{1},1-\alpha_{1}-\beta_{1})$. An ABC approach requires the choice of summary statistics, which are necessarily non-sufficient since the model is a state-space model. Following (38), we use the maximum likelihood estimator as summary statistics, relying on the R function garch for its derivation despite its lack of stability. Since (40) derived natural score constraints for the empirical likelihood associated with this model, we used their constraints to build our BC${}_{\text{el}}$ algorithm. Fig. 3 provides a comparison of both approaches with the MLE. It shows in particular that the ABC algorithm is unable to produce acceptable inference in this case, even in the most favorable case when it is initialized at a satisfactory maximum likelihood estimate (as shown by the bottom row). The BC${}_{\text{el}}$ algorithm is performing better, even though it fails to catch the correct range of $\beta_{1}$.

Another type of non-iid model relying on the superposition of an unknown number of gamma point processes and processed in (41) through a (non-Bayesian) alternative to ABC is discussed in the SI as an additional illustration of the possibilities of the empirical likelihood perspective for complex models, offering a free benchmark for evaluating the ABC outcome. Figure S7 shows a clear improvement brought by BC${}_{\text{el}}$ over the corresponding ABC outcome.

4 Population genetics

We compare our proposal with the reliable ABC-based estimates given by (3). We set up two toy experiments that are designed to defeat ABC, using pseudo observed data. The two evolutionary scenarios are given in Figure 4. In all experiments, we only consider microsatellite loci and assume that the effective population size is identical over all populations of the scenario.

In the second experiment, we consider three populations, see Figure 4 (right): the last two populations diverged at time $\tau_{1}$ and their common ancestral population diverged from the first population at time $\tau_{2}$. The sample comprises thirty diploid individuals per population genotyped at a hundred independent loci. In contrast to the first experiment, all components of the composite scores are computed here by summing over all pairs of genes whatever the population to which they belong. The results given in Table 1 show that ABC and BC${}_{\text{el}}$ mainly agree on both parameters $\theta$ and $\tau_{1}$, but BC${}_{\text{el}}$ is slightly more accurate than ABC on $\tau_{2}$.

Table S1 gives a comparison of the computing times for both algorithms, showing the difference of magnitudes between them. This is due to the simulation of the the simulated datasets for ABC: While this difference should not be over-interpreted, it signals a potential for self-assessment and testing that is missing for ABC methods.

Discussion

When compared with ABC methods, the (often) significant time savings provided by BC${}_{\text{el}}$ due to the lack of pseudo-sample simulation may open wider ranges for processing models involving complex likelihoods. For instance, in population genetics, ABC is severely hindered by the time spent simulating a dataset when modelling isolation by distance in a continuously distributed population, or when studying a large set of SNP markers even on quite simple evolution scenarios. Moreover, when the dataset is composed of large sets of markers, the summary statistics proposed in ABC (in DIY-ABC, these are averages of some quantitative statistics over all loci) ignores some (statistical) information, while BC${}_{\text{el}}$ manages to recover most of it, more specifically to estimate divergence on large datasets. Improvements in accuracy of estimation and computational efficiency are also possible in other contexts as illustrated in the range of examples given above.

Even when BC${}_{\text{el}}$ requires the same computing time as ABC, it uses the outcome in a very different perspective and provides a benchmark likelihood that helps in evaluating the pertinence of the ABC approximation, as illustrated in the gamma point process of SI.

We acknowledge that a caveat of the empirical likelihood is that it requires a careful choice of the constraint (1). Those pivotal quantities have to be connected to the parameter in an identifying way, which may require complex manipulations as in the gamma process case or even be impossible. However, repeated experimentation is often available, as illustrated by the normal example and the population genetic experiments (where we computed the composite score on both a restricted set of pairs and all pairs of genes). Checking for the accuracy of the approximation means that a constraint in BC${}_{\text{el}}$should be tested on simulated datasets in controlled experiments where the true parameters are known, although much less than in ABC runs. Then we can test coverage of credibility intervals, and measure the error of various point estimates based on the output of the scheme.

Acknowledgments

The last two authors wish to thank Jean-Marie Cornuet for his help and availability. Their work has been partly supported by the Agence Nationale de la Recherche (ANR) through the 2009–2012 project Emile. The third author is grateful to Patrice Bertail, Chris Drovandi, Brunero Liseo, and Art Owen for useful discussions. Comments and suggestions from the whole PNAS editorial board greatly contributed to improve both the presentation and the scope of the paper.

References

Supplementary information (SI)

We also reproduce here an illuminating comment from Art Owen: “The interesting thing about Theorem 3.4 is what is not there. It includes no conditions to make $\hat{\theta}$ a good estimate of $\theta_{0}$, nor even conditions to ensure a unique value for $\theta_{0}$, nor even that any solution $\theta_{0}$ exists. Theorem 3.4 applies in the just determined, over-determined, and under-determined cases. When we can prove that our estimating equations uniquely define $\theta_{0}$, and provide a consistent estimator $\hat{\theta}$ of it, then confidence regions and tests follow almost automatically through Theorem 3.4.”.

Pairwise composite likelihoods in population genetics

We detail here the derivation of the composite likelihoods used for the version of the BC${}_{\text{el}}$ algorithm implemented in the case of the population genetics study.

First, we recall that we scale the time axis so that a pair of genes of the same deme coalesces at a random time with an exponential distribution with rate $1$. We now consider a given locus and two microsatellite genes from our sample that come from the same deme. We denote their respective allelic state by $x_{1}$ and $x_{2}$. Their most recent common ancestor (MRCA) dates back to a time $T$, where $T\sim\mathcal{E}(1)$. We assume that the mutation rate, namely $\theta/2$ does not vary along the whole history of our populations. Therefore, conditioned on $T$, the number of mutations between $x_{i}$ ($i=1,2$) and the MRCA is distributed according to a Poisson distribution with mean $\theta T/2$. Hence, conditional on $T$, the number $N_{0}$ of mutations between $x_{1}$ and $x_{2}$ is a Poisson variable with mean $\theta T$ and

i.e., $N_{0}\sim\mathcal{G}e(\theta/(1+\theta))$, the geometric distribution with positive weight at . Finally, the difference between both genes is the accumulation over the $N_{0}$ mutations, i.e.,

where the $\epsilon_{k}$’s are iid Rademachers ($\pm 1$ with equal probability). Thus,

which proves that the pairwise likelihood is

with \rho(\theta)=\theta\big{/}\big{(}1+\theta+\sqrt{1+2\theta}\big{)}.

Two genes from different demes

We now consider two genes that come from two different demes that diverged from an ancestral deme at time $\tau$ in the past. We denote the allelic state of the two ancestors at time $\tau$ by $x_{1}^{0}$ and $x_{2}^{0}$, respectively. Then, $x_{1}-x_{2}=(x_{1}-x_{1}^{0})+(x_{1}^{0}-x_{2}^{0})+(x_{2}^{0}-x_{2})$, where $x_{1}^{0}-x_{2}^{0}$ follows a distribution whose Fourier transform is given by (4), while $(x_{j}-x_{j}^{0})$, $j=1,2$ are iid., whose distribution is given by the difference of two allelic states separated by a fixed time $\tau$. This distribution is derived in Equation (3) of (2):

where $I_{\delta}$ denotes the $\delta$th-order modified Bessel function of the first kind, given by $(n\geq 0)$

Using the independence between $(x_{1}-x_{1}^{0})$ and $(x_{1}-x_{1}^{0})$, we obtain

We then retrieve this distribution by computing Fourier transforms in the same vein as above. First, we note that the number $N_{1}$ of mutations between $x_{1}$ and its ancestor at time $\tau$ is a Poisson variate with parameter $\tau\theta/2$. And,

Finally, the distribution of $x_{2}-x_{1}$ is a (discrete) convolution product between the distributions given by (5) and (6) which yields

Quantile estimation

Examples of quantile distributions are the three-, four- and five-parameter Tukey’s lambda distributions and their generalizations and the Burr family of distributions; particular examples include the $g$-and-$h$ and $g$-and-$k$ distributions (3, 4, 5, 7).

Proposed methods for estimation of quantile distributions include maximum likelihood estimation using numerical approximations to the likelihood (7, 8, 9), moment matching (10, 11), location and scale-free shape functionals (12), percentile matching (6), quantile matching (13) and, more recently, ABC (14, 15, 16). Sequential Monte Carlo approaches for multivariate extensions of the $g$-and-$k$ have also been proposed (37).

There has been a number of ABC approaches proposed for this problem. For example, (15) adopted the ABC-MCMC algorithm of (18), in which draws of $\theta$ are based on a Metropolis algorithm with a Gaussian proposal distribution, and are accepted based on the rule $\rho(S(D),S(D\prime))<\epsilon)$, where D is the entire set of order statistics, $\rho$ is the Euclidean norm and $\epsilon$ is heuristically chosen after inspection of a histogram of $\rho(S,S\prime)$ obtained from a preliminary run using a very large value of $\epsilon$. This approach has recently been improved by (16) through more sophisticated MCMC approaches, the use of regression summary statistics for $D$ based on percentiles and their powers, and more automated choices of $\epsilon$. However, they still maintain a form of distance-based measure $\rho(S,S\prime)$ in accepting $\theta$.

Normal estimation

Figures S1–S3 evaluate the impact on the posterior distribution approximation of increasing the number of constraints in the empirical likelihood definition. Since this is a formal example, the true posterior distribution is available.

Quantile estimation

The BC${}_{\text{el}}$ experiment involved evaluation of Algorithm 1 for estimation of the parameters of the $g$-and-$k$ distribution using the two values of $\theta=(A,B,g,k)$, namely $\theta_{N}=(0,1,0,0)$, which corresponds to the standard normal distribution, and $\theta_{A}=(3,2,1,0.5)$; which was chosen by (36) as ‘an interesting, far-from-normal distribution’. The simulation experiment comprised multiple repetitions of BC${}_{\text{el}}$ using different combinations of sample size, $n=(100,500)$, number of iterations, $M=(1000,5000,10000)$, and number of constraints ($p=3,4,4,5,9$), corresponding to percentile sets $(0.2,0.5,0.8)$, $(0.2,0.4,0.6,0.8)$, $(0.1,0.4,0.6,0.9)$, $(0.1,0.25,0.5,0.75,0.9)$ and $(0.1,0.2,...,0.9)$. Two sets of priors were considered for $(A,B,g,k)$: $U^{4}$ (denoted as $P_{1}$) and $U(-5,5).U(0,5),U(-5,5),(-0.1,1)$ (denoted as $P_{2}$). Although the priors were set independently for each element of $\theta$, the four elements were drawn together at each iteration of the algorithm, so that the same importance weight $\omega_{i}$ was attached to the values $A^{(i)},B^{(i)},g^{(i)},k^{(i)}$ drawn in the $i$th iteration. The experiment was replicated ten times with different draws of samples of size $n$. Posterior means and standard deviations were computed for each parameter, and the overall goodness of fit to the true curve was assessed by comparing the true quantiles at $(0.05,0.10,...,0.95)$ with two measures: the estimated mean at each quantile (denoted by RSSm) and the average of the estimated quantile for each importance sample (RSSt).

Boxplots of the posterior means and standard deviations are shown in Figure S4 and Figure S5, respectively, for the four parameters, based on $\theta_{A}$, prior $P_{2}$ and the 20 replicates for $M=(5000,10000)$, for ten of the trials: $p=3,4,4,5,9$ for $n=100$ (trials 1-5) and $n=500$ (trials 6-10). Boxplots for the corresponding overall goodness of fit measures (RSSm, RSSt) are given in Figure S6.

Superposition of point processes

(41) discuss an alternative to ABC, using fractional design and linear interpolation. While their purpose is the non-Bayesian processing of models with intractable likelihood functions, they propose as their main example a model consisting in the superposition of $N$ renewal processes with waiting times $\tau_{ij}$ $(i=1,\ldots,M),\,j=1,\ldots)$ distributed as $\mathcal{G}(\alpha,\beta)$ variables, when $N$ is unknown. The renewal processes are thus

and the observations are made of the first $n$ values of the $\zeta_{ij}$’s,

This model offers an interesting testing ground for BC${}_{\text{el}}$ in that the data points $z_{t}$ are neither iid nor Markov. It is however possible to recover and exploit an iid structure in this case by first simulating a pseudo-dataset, $(z^{\star}_{1},\ldots,z^{\star}_{n})$, as in ABC settings, and then deriving a sequence of renewal processes indicators $(\nu_{1},\ldots,\nu_{n})$, as

These indicators are thus distributed from the prior distribution on the $\nu_{t}$’s and an iid sample of $\mathcal{G}(\alpha,\beta)$ variables can be derived from those indicators and the genuine data, leading to an associated empirical likelihood. As shown on Figure S7, when applied to a simulated dataset (as in (41)), the empirical likelihood approximation produces a better approximation than the corresponding ABC solution based on the same statistics as (41) (for exactly the same computational cost).

Time gains in population genetic models

In general, the speed of executing an ABC algorithm depends on many factors, including:

one’s ability to program an efficient simulator from the model distribution and to compute the selected summary statistics,

and the size of the Monte Carlo sample, denoted $M$,

and the speed of BC${}_{\text{el}}$ depends on:

the difficulty to optimize under the constraints (in the population genetics examples, this is not straightforward because of the Bessel functions and of the various series involved when the two individuals are not in the same deme),

and the size $M$ of the Monte Carlo sample.

While producing many simulations from the model distribution often is the stumbling block for ABC algorithms, the selection of the constraints and the time requirements of the optimization step are both highly variable and delicate to quantify. While all the experiments described in this paper induced no inflation in computing time and mostly significant reductions, we cannot exclude the possibility of $\text{BC}_{\text{el}}$ requiring more computing time than ABC.

Both population genetic experiments conducted in this paper analyse datasets with a large number of loci (one hundred). Thus ABC, which requires simulations of all loci to produce a simulated dataset, is quite time consuming and particularly so when the evolutionary scenario is more complex than the one in the first experiment. We compare here the computing times required by our implementation of the BC${}_{\text{el}}$-AMIS sampler and by DIYABC (3) on an Intel Xeon W3680 Plateform with GNU/Linux. Both methods were parallelized over five among the six cores of this CPU with the OpenMP API. Table 2 exhibits computation time averages on ten replicates of the estimation.

References