Weighted approximate Bayesian computation via Sanov’s theorem

We consider the problem of sample degeneracy in Approximate Bayesian Computation. It arises when proposed values of the parameters, once given as input to the generative model, rarely lead to simulations resembling the observed data and are hence discarded. Such “poor” parameter proposals do not contribute at all to the representation of the parameter’s posterior distribution. This leads to a very large number of required simulations and/or a waste of computational resources, as well as to distortions in the computed posterior distribution. To mitigate this problem, we propose an algorithm, referred to as the Large Deviations Weighted Approximate Bayesian Computation algorithm, where, via Sanov’s Theorem, strictly positive weights are computed for all proposed parameters, thus avoiding the rejection step altogether. In order to derive a computable asymptotic approximation from Sanov’s result, we adopt the information theoretic “method of types” formulation of the method of Large Deviations, thus restricting our attention to models for i.i.d. discrete random variables. Finally, we experimentally evaluate our method through a proof-of-concept implementation.


Introduction
Approximate Bayesian Computation (ABC) is a broad class of methods allowing Bayesian inference on parameters governing complex models. For such models, computing the likelihood, either analytically or numerically, is typically unfeasible. To overcome this critical problem, ABC dispenses with exact likelihood computation, and only requires the ability of simulating pseudo-data by sampling observations from a generative model, as detailed in Sect. 2.
In the literature, a variety of ABC methods have been proposed, see Sisson et al (2018, Ch. 4), and for recent reviews Lintusaari et al. (2017), Karabatsos et al. (2018). In the vast majority of these methods, the approximate likelihood function takes positive values only when the distance between the simulated and the observed data is lower than a predefined threshold. In other words, most ABC schemes involve-implicitly or explicitly-a rejection step, which often leads to discarding a huge number of proposals. This results in a waste of computational resources and/or in an inadequate sample size, that is, in sample degeneracy. Sample degeneracy may also cause serious distortions in the form of the approximate posterior distribution, at least when the number of iterations is not large enough. Indeed, accepting poor parameter proposals, i.e., those producing simulated data very rarely resembling the observed data, is a rare event. In the lack of accepted values, the posterior probability of such proposals will be approximated just as zero, in turn resulting in a distortion in the tails. This may be especially problematic for posterior distributions with long tails.
Our idea is to mitigate the problem of sample degeneracy by improving the approximation of the likelihood function. In particular, we speculate that taking into account the positive, however small, probability of rare events, i.e., poor proposals leading to simulated data resembling the observed data, allows avoiding the rejection step altogether and weighting all parameter proposals. To this end, we resort to the theory of Large Deviations (LDT). Our aim is to show how LDT provides a convenient way to define an approximate likelihood, as well as guarantees of its convergence to the true likelihood as the size of pseudo-datasets goes to infinity. In order to make the incorporation of LDT into ABC as smooth as possible, we rely on one of the less general formulations of Sanov's theorem. Accordingly, we only consider models for discrete i.i.d. random variables which, despite their apparent simplicity, will be shown to be of interest in several applications of ABC. This allows adopting a straightforward information theoretic formulation of LDT known as the method of types (Csiszár 1998;Cover and Thomas 2006). Related work In the literature, there have been many proposals aimed at improving the computational efficiency of basic ABC. Prangle (2016) proposed Lazy ABC, which saves computing time by abandoning simulations likely to lead to a poor match between the simulated and the observed data. To this end, at each iteration, the simulation is given up with a probability depending on the probability of its acceptance and on the expected required time for its completion. Unlike our method, Lazy ABC does not avoid rejection, but rather accelerates the process leading to discarding a proposal.
Another way to improve computational efficiency is to consider proposal distributions closer to the posterior on the parameter space, employing sophisticated sampling methods, such as MCMC (Marjoram et al. 2003), Population Monte Carlo (Beaumont et al. 2009) and Sequential Monte Carlo (Del Moral et al. 2012). In the same vein, Chiachio et al. (2014) proposed a sequential way of achieving computational efficiency by overcoming the difficulties in getting samples resembling the observed data. This latter was the first attempt to improve the acceptance rate by adopting a rare-event approach. In particular, Chiachio et al. (2014) combine the ABC scheme with a rare-event sampler that draws conditional samples from a nested sequence of subdomains. However, even this method cannot completely avoid rejections, and only partially mitigates the sample degeneracy problem.
In order to tackle the problem more systematically, clever proposal distributions should be combined with better approximations to the likelihood. Accordingly, Prangle et al. (2018) also resorted to a sequential approach, but explicitly considering a likelihood estimate that takes into account the probability of rare events. As a comparison, our method evaluates the probabilities of rare events based on theoretical results (LDT), rather than on Monte Carlo estimates of tail probabilities. Moreover, they focus on continuous data by showing that extensions to discrete data can be challenging and require application-specific solutions; in contrast, the method of types provides a natural way of dealing with discrete random variables by summarizing data via empirical distributions, thus avoiding the common practice of summarizing data by selecting ad hoc summary statistics.
Other methods have been proposed avoiding the selection of summary statistics and relying on empirical distributions. In particular, Park et al. (2016) rely on the maximum mean discrepancy between the embeddings of the simulated and the observed empirical distributions. They avoid rejection by weighting each parameter proposal by means of a kernel function defined on a non-compact support. Other interesting methods involve the Wasserstein distance (Bernton et al. 2019) or the Kullback-Leibler divergence (Jiang 2018) as measures of the discrepancy between observed and simulated data. In particular, Jiang (2018) approximates the likelihood by means of an estimator of the Kullback-Leibler divergence between the unknown distribution of the data given the true parameter, and given the parameter sampled at the current iteration. Exploiting the fact that the maximum likelihood estimator is the one minimizing that Kullback-Leibler divergence, they prove that their approximate posterior distribution converges to a restriction of the prior distribution on the region in which the above mentioned divergence is smaller than a predefined threshold. Although most of the above mentioned methods apply to continuous data, we note that ABC applications to discrete data appear frequently in population genetics, epidemiology, ecology and system biology (see Beaumont (2010) for an overview of the applications of ABC in these fields). In particular, in population genetics, discrete (possibly i.i.d.) data representing the genotyping at a few loci of different (unrelated) individuals have often been summarized through their empirical distributions (Marjoram et al. 2003;Buzbas and Rosenberg 2015, among others).
A very different way of bypassing the selection of summary statistics relies on the random forest method (Raynal et al. 2019). Here, regression random forests are trained by using a training-set composed of a large number of parameter proposals and pseudo-datasets sampled from the prior distribution and the generative model, respectively. Since all the summary statistics are involved as covariates, summary selection is avoided. The output of the algorithm is the predicted expected value of an arbitrary function of interest on the parameter space, conditional on the observed data.
Structure of the paper The rest of this paper is structured as follows. Section 2 contains background and preliminaries on ABC, focusing on the importance sampling scheme and the sample degeneracy issue. In Sect. 3 we introduce LDT by adopting the method of types. We also show how LDT allows poor parameter proposals to contribute to the representation of the approximate posterior distribution. Section 4 gives the LDW-ABC algorithm and compares it with R-ABC. In Sect. 5 we present the results of a toy example and an experiment conducted on a real world dataset. Section 6 contains some concluding remarks and ideas for future research. The Appendices contain the proofs, technical materials, and additional results from experiments.

Background on ABC
Let x ∈ X n be the observed data, which will be assumed to be drawn from a probability distribution in the family F = {P(·|θ) : θ ∈ Θ}.
In principle, given a prior distribution π(θ) on Θ, the aim of Bayesian inference is to provide information about the uncertainty on θ by deriving the posterior distribution π(θ|x) ∝ π(θ)P(x|θ) via Bayes' Theorem. When the likelihood function is intractable, ABC allows simulated inference providing a conversion of samples from the prior to samples from the posterior distribution, through comparisons between the observed data and the pseudo-datasets generated from a simulator. A simulator can be thought of as a probabilistic computer program taking as input a parameter value (or a vector thereof) θ ∈ Θ and returning a sample from the distribution P(·|θ). In general, no knowledge of the analytical form of the likelihood is necessary to write down such a program. More specifically, in the primal rejection sampling algorithm, whose origins can be traced back to Rubin (1984), Tavaré et al. (1997), Pritchard et al. (1999), the following actions are taken: 1. S ≥ 1 parameter values from the prior distribution π(·) are generated; 2. for each s ∈ {1, . . . , S}, given the parameter proposal θ (s) as input, the simulator generates a realization of a random variable Y ∈ X n distributed according to P(·|θ (s) ); 3. only parameter values leading to a pseudo-dataset equal to the observed data are accepted, thereby samples from the exact posterior are derived by conditioning on the event {Y = x}.
Introducing a twofold approximation scheme, as illustrated in Algorithm 1, might increase the efficiency of the algorithm outlined above. First, one introduces a summary statistic, s(·), which is a function from the sample space X n ⊆ R n to a lowerdimensional space S ⊂ R k , with k n. Second, exact matching of the simulated and the observed data is relaxed to similarity, expressed in terms of a predefined distance function d(·, ·) and tolerance threshold > 0.
Abbreviating s(x) by s x and s( y) by s y , the output of Algorithm 1 is a sample of pairs (θ (s) , s (s) y ) from the following approximate joint posterior distributioñ π(θ, s y |s x ) ∝ π(θ) P(s y |θ) 1{d(s y , s x ) ≤ } (1), that is, ignoring the simulated summary statistics, the output of the algorithm becomes a sample from the marginal posterior distribution Here Pr d(s Y , s x ) ≤ |θ is called the ABC approximate likelihood.
Remark 1 (Marginal samplers) Some ABC sampling schemes (Sisson et al. 2007;Marjoram et al. 2003, among others) allow directly sampling from the approximate marginal posterior distributionπ(θ|s x ). The key idea is thatπ(θ|s x ) can be estimated pointwise by by simulating M pseudo-datasets from P(·|θ (s) ) and computing s (i) y for i ∈ 1, . . . , M at each iteration s. As is apparent, the second term in (2) provides a Monte Carlo estimate of the ABC approximate likelihood.
Note that marginalizing the output of Algorithm 1 corresponds to the implementation of a marginal sampler with M = 1. In such a case, the indicator function represents a crude Monte Carlo estimate of the probability Pr d(s Y , s x ) ≤ |θ .
As pointed out by Sisson et al. (2018, Ch. 1), the use of the indicator function does not enable one to discriminate between whether the pseudo-dataset y coincides with the observed data and whether the pseudo-dataset just is close enough. This may lead to a waste of information. For this reason, the indicator function in (1) is often replaced by a kernel function:

Algorithm 2 IS-ABC
for s = 1, . . . , S do Draw θ (s) ∼ q Generate y ∼ P(·|θ (s) ) from the simulator Set the IS weight for θ (s) to ω s = K d(s y , s x ) · π(θ (s) ) q(θ (s) ) end for where κ · is a kernel function (e.g., triangular, Epanechnikov, Gaussian, etc.) defined on a compact support and decaying continuously from 1 to 0, see e.g. Beaumont et al. (2002). Now the ABC approximate likelihood becomes the convolution of the true model with the kernel K (Prangle et al. 2018): Note that this general setting encompasses also the case when R-ABC employs the uniform kernel as κ(·). The accuracy of the posterior distribution approximation depends both on how much information about the parameters is preserved by the summary statistics and on the magnitude of the threshold . In fact, as → 0, the approximate likelihoodL(θ ; x) converges to the true likelihood (Prangle et al. 2018, Appendix A) and, whenever sufficient summary statistics for θ have been chosen, the approximate posterior distributionπ(·|s x ) converges to the true posterior π(·|x) (Sisson et al. 2018, Ch. 1). On the other hand, as → ∞, the probability Pr d(s Y , s x ) ≤ |θ approaches 1 and samples are generated from the prior distribution. This establishes a trade-off between the statistical bias and the computational efficiency (Lintusaari et al. 2017): as the tolerance level decreases, the error of the approximation of the ABC posterior vs. the true posterior decreases at the cost of higher computational effort.

Importance sampling ABC and sample degeneracy
In the ABC literature, a great variety of methods to sample fromπ(θ, s y |s x ) have been proposed that go beyond the rejection scheme. Hereafter, we will adopt an importance sampling scheme (IS-ABC) which, as outlined by Karabatsos et al. (2018), encompasses R-ABC and most of the other ABC algorithms. Like the standard importance sampling, see Robert and Casella (2013, Ch. 3), IS-ABC consists of sampling pairs (θ (s) , s (s) y ) from an importance distribution and weighting each pair, avoiding the computation of the acceptance probabilities. More formally, let h : (Θ × S) → R be a function of interest and let E p [h(θ, s Y )] denote its expected value w.r.t. a probability distribution p over Θ × S. Suppose that we are interested in estimating Eπ [h(θ, s Y )], whereπ is our target distribution, i.e., the joint approximate posterior. In particular, by choosing h(·) to be a kernel function, this formulation also enables a kernel density estimation for the joint approximate posterior, a case which will be considered in Sect. 5.

Now it is a standard fact that
where q(θ, s y ) is the importance distribution on Θ × S andω(θ, s y ) =π (θ, s y |s x ) q(θ, s y ) are the importance weights. In particular, in the ABC framework, the importance distribution can be q(θ, s y ) = q(θ )P(s y |θ) and, denoting by Z the normalizing constant for the joint posterior, the resulting importance weightsω(θ (s) , s (s) y ),ω s for short, arē By computing, at each iteration s, the following unnormalized weight an approximation for the constant Z is obtained: where the second equality is obtained by multiplying and dividing by q(θ, s y ). It follows that from the output of Algorithm 2 we can estimate (5) as where eachω s = ω s / S r =1 ω r is a normalized weight. Unlike the standard importance sampling scheme, IS-ABC implicitly involves a rejection step. In fact, usually the kernel density function, K (·), is such that a strictly positive weight is given to a pair (θ (s) , s y , s x ) ≤ . For example, looking at Algorithm 1, it is apparent that the primal rejection scheme is a special case of Algorithm 2 where the marginal importance distribution, q(θ ), is the prior distribution and the resulting importance weights are meaning that each pair is implicitly rejected or accepted depending on the value of ω s ∈ {0, 1}. In order to evaluate the efficiency of an importance sampling method, a widespread "rule of thumb" is to evaluate the Effective Sample Size (ESS), see Liu (2008, Ch 2). ESS represents the number of samples from the target distribution needed to get a Monte Carlo estimate with the same variance as the IS estimate in (7) given a budget of S iterations, and is defined by One of the major drawbacks of the importance sampling scheme is that the resulting Monte Carlo estimate in (7) is highly variable due to the problem of sample degeneracy, already mentioned previously, which in this context means that only a few of the proposed pairs (θ, s y ) have relatively high weights resulting in a small ESS. Generally speaking, sample degeneracy is caused by an importance distribution far from its target. In this case, parameter values from regions with a low target posterior density are very likely to be drawn under the importance distribution, so that they are often proposed and associated with very small weights.
In the ABC framework as well, an importance density q(θ ) far from the marginal targetπ(θ|s x ) can lead to sample degeneracy. In this setting, an additional issue is that the weights also depend on the distance d(s Y , s x ), hence on the random variable 1 s Y (Sisson et al. 2018). This implies that when a parameter θ * is proposed such that Pr(s Y = s x |θ * ) is close to zero, usually a very large number of zero-weighted pairs (θ * , s y ) will be generated before a distance smaller than will be observed, especially when is small.
In the next two sections we propose a method to define a kernel K (·) that improves the efficiency of IS-ABC in terms of ESS.

ABC and the theory of large deviations
Recall that in R-ABC, at each iteration s, the indicator function represents a crude estimate for the probability Pr(d(s Y , s x ) ≤ | θ (s) ) (see Remark 1). A possible approach to mitigate sample degeneracy is to provide a finer estimate for the ABC likelihood by evaluating that probability. In order to deal with rare events, we resort to LDT, which studies the exponential decay of the probability of rare events. We speculate that taking into account the positive probability of a large deviation event allows one to avoid rejection at all. This might provide a higher ESS, thus making the algorithm more efficient.
We will let and so on denote sequences of i.i.d. random variables, distributed according to an (intractable) probability distribution P θ ∈ F .

LDT via the method of types
Let x n be a sequence of n symbols drawn from X, say x n = (x 1 , . . . , x n ). The method of types moves the focus from the sequence x n itself to its type, defined as follows.
We let T n denote the set of n-types, that is, types with denominator n.
Note that the superscript n keeps track of the length of the sequence, which is also the denominator of the type. As is apparent, type is a function summarizing the information included in the observed sequence x n by mapping the n-dimensional observed sequence onto a |X|-dimensional summary statistic.
The following quantities play a crucial role in the method of types. Below, we stipulate that 0 · log 0 r = 0 and that r · log r 0 = +∞ if r > 0, where log denotes the logarithm to base 2. Given two probability distributions on X, P and Q, we consider -the entropy of P, defined as -the Kullback-Leibler divergence between P and Q, defined as With an abuse of notation, whenever the first argument of D(·||Q) is a set of probability distributions, say E, D(E||Q) stands for inf P∈E D(P||Q). When P * = argmin P∈E D(P||Q) exists, it is called the information projection of Q onto E.
be a sequence of i.i.d. random variables, distributed according to P θ = P(·|θ), for some θ ∈ Θ. In what follows, we let Pr(·|θ) be the probability measure on sequences induced by P θ . The joint probability of n i.i.d. extractions x n according to P θ , can be written as See Cover and Thomas (2006, Ch.11) for a proof. It follows from the Neyman-Fisher theorem (Cox and Hinkley 1979, Ch. 2.2) that types are always sufficient statistics for θ , whatever P θ .
Remark 2 (Types and ABC) While the number of sequences of length n is exponential in n, it is easy to show that the cardinality of T n is polynomial in n; in fact, |T n | ≤ (n + 1) |X| , see Cover and Thomas (2006, Ch.11). From an ABC perspective, it follows that using types as summary statistics could mitigate the computational problems related to the comparison between the observed dataset and the pseudo-dataset, especially for large n. Furthermore, summarizing data through their empirical distributions is a way of overcoming the difficulties in finding sufficient statistics when P θ is unknown (and Pr(·|θ) as well). Indeed, even when confined to discrete random variables, P(·| θ) is an unknown model, not necessarily a Multinomial model, see Sect. 5 for examples. With no knowledge of the analytical form of the likelihood, finding a sufficient summary statistic for θ , the vector of parameters given as an input to the simulator, is a central issue. In the literature there are several examples of models for conditionally independent discrete data in which the likelihood is analytically intractable and the required ABC method concerns empirical distributions. Examples are the ABC methods proposed by Joyce et al. (2012) and Buzbas and Rosenberg (2015) to make inference on the mutation and selection parameters governing the Fisher-Wright model (Fisher 1930). There, despite the conditional independence and the discretness of the observations, the likelihood function is difficult to evaluate since the normalizing constant depends on the parameters: for small values of the selection parameter, numerical solutions have been found by Genz and Joyce (2003), in other cases, likelihood-free methods are required.
Noting from (11) that the probability of the observed sequence decreases exponentially at a rate given by the Kullback-Leibler divergence between T x n and P θ , we can The Law of Large Numbers (LLN) states that as the length of a typical sequence goes to infinity, its type converges in probability to P θ , see Cover and Thomas (2006, Ch 11.2.1) for formulation of the LLN in terms of the method of types presented below.
On the other hand, observing a sequence whose type is far from P θ , called a nontypical sequence, is a rare event, and its probability obeys a fundamental result in LDT, Sanov's theorem; see Cover and Thomas (2006, Th.11.4.1).
where P * = argmin Suppose that E is composed of types of non-typical sequences. Then Sanov's theorem characterizes the exponential decrease rate of the probability of E. Taking into account this probability may provide a finer ABC approximation of the likelihood, as discussed in the next section.

LDT in ABC
In this section we provide a formal explanation of what is meant by poor parameter proposals and how they can contribute to the representation of the approximate posterior distribution by means of LDT. We are interested in obtaining an approximation of the posterior distribution,π(θ|x n ), via R-ABC or an equivalent IS-ABC by assuming as given: (a) the marginal importance density q(θ ) to be the prior distribution on Θ; (b) > 0 as a threshold; (c) types as summary statistics; (d) the Kullback-Leibler divergence as distance function. For the sake of simplicity, from now on we will also assume T x n to be full support.
Given a budget of S iterations, both R-ABC and IS-ABC generate a sequence of pairs (T (s) (s) . We stress that the length of the simulated sequence, m, need not be equal to n, the length of the observed data sequence. Note also that, because of the independence assumption, choosing m = M · n with M ∈ N means that the algorithm simulates M pseudo-datasets at each iteration, like a marginal sampler (see Remark 1).
Looking at Algorithms 1 and 2, being the whole pair (θ (s) , T (s) y m ) accepted or rejected, one can define the joint acceptance region for these algorithms on the space Θ × T m . However, as the acceptance rule is based only on the simulated type, regardless of the proposed parameter value, the acceptance region can be projected onto the probability simplex Δ |X|−1 ⊃ T m .
Definition 2 (Acceptance region) Let Δ |X|−1 be the simplex of probability distributions over X and let T x n be the type of the observed sequence x n . The acceptance region B (T x n ), referred to as B for short, is defined for any ≥ 0, as Now we can define a poor parameter proposal as a parameter θ (s) such that T x n and the other types in the acceptance region are types of non-typical sequences w.r.t. P(· | θ (s) ).
Accordingly, sampling a poor parameter means that there is a large divergence between T x n and P(·|θ (s) ). On the other hand, with m large enough, T (s) y m is very likely to be close to P(·|θ (s) ), due to the Law of Large Numbers. Heuristically, this implies that the probability of simulating a sequence y m whose type is in the acceptance region is very small. Recalling that in R-ABC and in IS-ABC outlined at the beginning of this section a crude Monte Carlo estimate of the probability Pr(T Y m ∈ B | θ (s) ) is given by the indicator function 1{D(T (s) y m ||T x n ) ≤ }, the vast majority of the poor parameter proposals are discarded altogether. We propose to mitigate this problem by assigning strictly positive weights to each proposal θ (s) , even if T (s) y m is outside the acceptance region. To this end, we want to replace the indicator function with a finer estimate of the probability Pr(T Y m ∈ B | θ (s) ).
In principle, Sanov's theorem implies that, for m large enough, that probability can be approximated at each iteration by By replacing the indicator function in (1) with (14), the approximate posterior becomesπ Unfortunately, the computation of the probability in (14) is still not feasible when the model F = {P θ : θ ∈ Θ} is unknown, as we do not know how to compute D(B ||P θ (s) ). The following theorem provides an asymptotic approximation to circumvent the problem. A proof is provided in "Appendix A".
In essence, this result says that, as m increases and the type T y m converges to the distribution P θ that has generated y m , the information projection of T y m onto B converges to that of P θ onto B (see Fig. 1). From (14) and Theorem 3, for m large enough, 2 −m D(B ||T y m ) provides a feasible asymptotic estimate for the acceptance probability, Fig. 1 Acceptance region, B , types, T x n and T y m , and the probability distribution P θ that generated Replacing the indicator function in the ABC approximate posterior (1) with this estimate, we obtain the following new joint approximate posterior distribution:π

Weighted approximate Bayesian computation
The discussion in the previous section indicates that IS-ABC can be improved by resorting to a better approximation for the likelihood. In particular, the (implicit) rejection step can be avoided by evaluating the positive probability of rare events via Sanov's theorem. Indeed, an easy way of sampling from (17) is a large deviations version of IS-ABC, which we will call the weighted approximate Bayesian computation (LDW-ABC).
Starting from the definition of an acceptance region satisfying the hypothesis of Sanov's theorem, as in Definition 2, a sample from the approximate posterior distributionπ(θ, T y m |T x n ) can be obtained as described in Algorithm 3.

Algorithm 3 LDW-ABC
Looking at Algorithm 3, it is apparent that LDW-ABC is a specialization of the more general IS-ABC. More specifically, the sufficient summary statistics involved are the types, the distance function is the Kullback-Leibler divergence, and the kernel density function is defined as follows: At each iteration a positive weight is assigned to the proposed θ (s) . More precisely, the weight equals 0 only when D(B ||T y m ) = ∞. Each ω s is computed by approximating the divergence D(B ||T y m ) as described in "Appendix B". As a special case of the general IS-ABC, the output of Algorithm 3 is a weighted sample from the following approximate joint posterior distribution: which, by marginalizing out simulated types, becomes where T m denotes the set of the m-types. Hence, the likelihood approximated by Note that the quality of the approximation depends both on the threshold and on the size of pseudo-dataset, m. More precisely, the adjustment w.r.t. the likelihood approximate by R-ABC, 2 here denoted byL R ,m (θ ; T x n ) = T y m ∈B P θ (T y m ), depends on m and . In fact, from (18) and Definition 2, the approximate likelihood in (21) can be written as where the term 0 ≤ α ,m (θ ) ≤ 1 is the adjustment. The following lemma gives an upper bound for that adjustment α ,m (θ ), in two cases depending on P θ . The proof is in "Appendix A".
be the difference between the two likelihood functions approximated by LDW-ABC and R-ABC. Let B be the ABC acceptance region andB its interior. We have the following upper bounds, depending on θ , which hold for all m ≥ 1.
From Proposition 1 it follows that as m goes to infinity, α ,m (θ ) → 0 for almost all θ ∈ Θ. Therefore, the approximate likelihood from LDW-ABC achieves the approximate likelihood from R-ABC and preserves its asymptotic properties. Moreover, we speculate that LDW-ABC improves the efficiency by mitigating the sample degeneracy. An evaluation of ESS might be a way of appreciating the improvement induced by avoiding the implicit rejection.
Since the term ESS in (9) involves the evaluation of the variance of the normalized weights, var[ω], its exact computation is infeasible, as it depends on the unknown target normalizing constant. For this reason, we adopt the following estimate derived by Kong (1992) and Elvira et al. (2018): (with the proviso that ESS = 0 if all ω s 's are zero). Let ESS I S and ESS L D be, respectively, the value of ESS achieved by S iterations of IS-ABC and LDW-ABC by setting the same tuning parameters, distance function and importance density q(θ ). Explicitly, let us assume that the kernel function for IS-ABC is 1 within the acceptance region B and 0 outside. Heuristically, adding positive weights increases the numerator more than the denominator in (22), suggesting that a non null weight assigned by LDW-ABC to a parameter proposal rejected by IS-ABC is enough to have ESS L D > ESS I S . This is confirmed by the following simple result, whose proof is given in "Appendix A".

Proposition 2 (Empirical ESS) It holds that ESS L D ≥ ESS I S . Moreover this inequality is strict, provided that in at least one iteration of the algorithm there is generated a full support T y m falling outside B .
A tedious but straightforward analysis shows in fact that the event mentioned in the statement, upon which strict inequality holds, occurs with probability 1 as S −→ +∞. The above result will be empirically validated in the experiments of Sect. 5, thus providing further evidence that LDW-ABC achieves an improvement in terms of efficiency.
Below, we sum up the technical development so far with a discussion of the role of the parameters m and .
Remark 3 (On the role of the tuning parameters) Concerning the role of m, the size of pseudo-dataset, and , the tolerance, we can sum up the content of Propositions 1 and 2 as follows: 1. large m and small point to low ESS and low α ,m ; 2. small m and large point to high ESS and high α ,m .
If one regards ESS as a measure of efficiency, and α ,m as a measure of (lack of) accuracy w.r.t. the R-ABC likelihood, (but see also below), 1 and 2 above indicate how to trade off one for the other.
In particular, as Theorem 3 requires a relatively large m in order to get a good approximation for the posterior probability, 1 above says we can increase the tolerance to mitigate the resulting inefficiency. On the other hand, in cases where a small tolerance parameter is required, 2 above offers room to mitigate the resulting inefficiency by decreasing m.
Note, however, that when considering accuracy w.r.t. the target posterior density π(θ|T x n ), the adjustment α ,m cannot simply be regarded as a measure of imprecision: rather, it represents a compensation for those θ 's that would be assigned a too low probability by pure R-ABC. In this case, a sounder measure of precision can be obtained by directly comparing a kernel-estimated density (obtained with LDW-ABC weights) and the target posterior density, e.g., in terms of the mean integrated squared error (MISE). This measure is, however, impossible to evaluate analytically, since its calculation presupposes the knowledge of the target posterior density. From a more empirical point of view, further discussion of the consequences of different choices of and m on the performances of the posterior estimators is presented in Sect. 5, illustrated by a number of examples.

Experiments
In order to evaluate the performance of the proposed method, we have put a proofof-concept implementation of LDW-ABC at work on two examples. We compare the results obtained from LDW-ABC with those obtained from R-ABC. For both examples, there is a MCMC method for sampling from the exact posterior distribution, and the resulting posterior inference is taken as a reference for comparison.

Example 1: mixture of binomial distributions
Let X n = {X i } n i=1 be a sequence of i.i.d. discrete random variables distributed according to the following parametric finite mixture model: Here we assume a uniform prior distribution on the mixture weight λ and that (θ 1 , θ 2 ) are uniformly distributed on the set {(θ 1 , θ 2 ) : 0 ≤ θ 2 ≤ θ 1 ≤ 1} by imposing the following identifiability constraint: An analytical computation of the posterior distribution requires the evaluation of the likelihood The direct computation of (24) is infeasible, as even with a few hundred observations, it involves the expansion of the likelihood into 2 n terms. In the literature, there are several methods to deal with this problem, which allow sampling from the parameters' posterior distributions, see Marin et al. (2005). A widespread method is a Gibbs Sampling handling the finite mixtures issue as a missing data problem, see Diebolt and Robert (1994). Samples from the joint posterior distribution are obtained by means of a hierarchical model involving a vector of latent random variables, Z n = {Z i } n i=1 , where each Z i ∼ Bernoulli(1−λ) indicates to which component the i-th observation belongs: Here the generative model consists of simulating each of the n values from one of the two binomials according to the result of a Bernoulli(1-λ) experiment. The same generative model has been run by a plug-in of the true values of the parameters displayed in Table 1 (LHS) to obtain the observed data. We ran Algorithm 4 as detailed in Table 1 (RHS), and after burn-in and thinning we got 5000 values for each parameter regarded as drawn independently from the true posterior distributions. The posterior means and variances are displayed in Table 2.
In order to compare performance of LDW-ABC with that of R-ABC, the marginal importance distributions are set equal to the prior distributions.
We ran Algorithms 1 and 3 with S = 100,000 and with four different pairs (m, ). Since our speculation is that introducing the evaluation of the probability of rare events provides a better approximation in the tails of the distributions, we are interested in   comparing the shape of the posterior distributions obtained by LDW-ABC and by R-ABC taking the Gibbs posterior distributions as reference. Accordingly, besides the posterior estimates of the means and variances, we also reconstruct the posterior densities by means of a Gaussian kernel density estimation. The point estimates are compared via the MSE, and the kernel density estimates via the MISE. In particular, the corresponding estimates, MSE and MISE, are computed by averaging over 100 runs the squared errors and the integrated squared errors w.r.t. the output of the Gibbs sampler. The results are summarized in Table 3.
First, we note that both the MSE and the MISE achieved by LDW-ABC are always lower for LDW-ABC than for R-ABC. Hence, in our example, taking into account the probability of large deviation events has improved both the point estimates and the approximation of the posterior distributions. Moreover, as already pointed out in Sect. 4, LDW-ABC mitigates the sample degeneracy by achieving an ESS up to more than five times that achieved by R-ABC (see Table 4).
In order to evaluate the sample degeneracy, Table 4 (RHS) also displays the normalized perplexity, which equals 2 H (ω) /S, where H (ω) denotes the entropy of the normalized weights. Cappé et al. (2008) show that the normalized perplexity represents an estimate of 2 −D π(θ,T y m |T x n )||q(θ,T y m ) , meaning that when the perplexity is larger, the sample degeneracy is smaller.     The following comments are consistent with Remark 3. Concerning the role of the tuning parameters, m and , we note that by fixing a large m (e.g., 5000), as increases both ESS and the perplexity increase. Moreover, both MSE and MISE decrease. The same happens by reducing m with fixed to a small value (e.g., 0.005). This provides guidance on how to set the tuning parameters. In Fig. 2 three matrices of plots, one for each parameter, show the posterior densities: the size of the pseudo-dataset, m, equals 500 in the plots on the LHS of each panel, and 5, 000 on the RHS. The topmost plots show the approximate distributions with = 0.01, the others the distributions corresponding to = 0.005. According to Remark 3, we note that as m increases the blue lines (LDW-ABC) overlap the red ones (R-ABC). In principle, we would expect that both the algorithms achieve a better approximation of the posterior shapes with = 0.005 than = 0.01. However, in the case of R-ABC, we see a deviation from the true posterior distributions (dotted lines), when moving from the first to the second row of each matrix. The same deviation occurs for LDW-ABC, but only in the second column, when m = 5,000. This suggests that the quality of the R-ABC approximation is affected by a low value of the ESS which is in turn determined by a too exacting value of and m. In fact, when m = 500, LDW-ABC manages to mitigate the effect of a small , but it fails when a large value of m causes a too small ESS for the LDW-ABC as well. In the figures we also superimposed the ratioL ,m (θ ; T x n )/L R ,m (θ ; T x n ) = 1 + α ,m (θ ; T x n )/L R ,m (θ ; T x n ) evaluated pointwise and shown by the gray dashed lines. This quantity depends on the contribution of the adjustment w.r.t. the R-ABC likelihood and shows how the adjustment acts in modifying this latter when the R-ABC posterior density is underestimated (gray areas). Figure 3 shows the posterior cumulative density functions for θ 2 . The posterior cumulative density functions for the other two parameters are given in "Appendix C", Fig. 4. We also show the 90% credible intervals for the estimated cumulative density functions. The red areas are always larger than the blue areas, meaning that the estimates provided by R-ABC exhibit greater variability. This is more significant when = 0.005, due to the small acceptance probability.
To wrap up, as suggested by the MISE's, the posterior distributions approximated by LDW-ABC appear more faithful to the true shapes. Moreover, the ESS and the variability of the estimates are less sensitive to small values of .

Example 2: learning from anonymized data
The second example is aimed at comparing LDW-ABC and R-ABC at work on a real-world dataset. We consider a scenario in which the dataset contains microdata that have been anonymized in order to protect the privacy of the individuals involved. More specifically, our data are a subset of 5692 rows from the Adult dataset extracted by Barry Becker from the 1994 US Census database and available from the UCI machine learning repository (Kohavi and Becker 1996). The anonymizaton method we have adopted, Anatomy (Xiao and Tao 2006), is a group based anonymization scheme. Given a dataset consisting of a collection of rows, each one corresponding to an individual and containing his/her sensitive (e.g., disease, income) and nonsensitive attributes (e.g., gender, nationality, ZIP code), a group-based anonymization algorithm produces an obfuscated version of itself by partitioning the rows into groups. The idea is to process the set of rows in each group so that even knowing the nonsensitive attribute of an individual, one cannot identify his/her sensitive values. To reach this goal, Anatomy vertically and randomly permutes the nonsensitive features within each group, thus breaking the link between the sensitive and nonsensitive attributes.
An example of group based anonymization is in Table 5, adapted from Wong et al. (2011). The LHS table is the original table collecting medical data from eight individuals; here, Disease is considered as the only sensitive attribute. The RHS table is an example of an application of the Anatomy scheme: within each group, the nonsensitive part of the rows is vertically and randomly permuted, thus breaking the link between the sensitive and nonsensitive values. Generally speaking, the obfuscated table can be seen as the output of a generative mechanism: given the population parameters as input, the mechanism first generates a cleartext table by drawing a number of rows i.i.d., then applies the anonymization algorithm to this table and outputs the result. One is interested in the posterior distribution of the population parameters, given an observation of the obfuscated table. Clearly, the likelihood function involved in this mechanism is highly nontrivial, and also depends on the details of the anonymization algorithm. Below we make this precise by adopting the model proposed by Boreale et al. (2020).
Let a row of the original (cleartext) dataset be a pair (s, r ) ∈ S × R, for finite, nonempty sets S and R. Here, s and r represent the sensitive and nonsensitive attributes (or vectors of attributes). Given a multiset of n rows, d = {|(s 1 , r 1 ), . . . , (s n , r n )|}, Anatomy will first arrange d into a sequence of groups, x = g 1 , . . . , g k , the cleartext table. Each group in turn is a sequence of n i rows, g i = (s i,1 , r i,1 ), . . . , (s i,n i , r i,n i ). The obfuscated table is then obtained as the sequence x * = g * 1 , . . . , g * k , where the obfuscation of each group g i is a pair g * i = (σ i , ρ i ). Here, each σ i = s i,1 , . . . , s i,n i is the sequence of sensitive values occurring in g i ; each ρ i , called the generalized nonsensitive value, is the multiset of g i 's nonsensitive values: ρ i = {|r i,1 , . . . , r i,n i |}i.e., ρ i includes all those and only those values, with multiplicities, found in g i .
The model consists of the following random variables: Here, θ takes values on the set of full support probability distributions D over S × R and represents the joint probability distribution of the sensitive and nonsensitive attributes in the population. -X = G 1 , . . . , G k , which takes values in the set of cleartext tables X. Each group G i is in turn a sequence of n i ≥ 1 consecutive rows in X, G i = (S i,1 , R i,1 ), . . . , (S i,n i , R i,n i ) with S i, j ∼ θ S and R i, j ∼ θ R|s i, j . The number of groups k is not fixed, but depends on the anonymization scheme and on the specific tuples in d. -X * = G * 1 , . . . , G * k , which takes values in the set of obfuscated tables X * . We assume that X * solely depends on the table X and the underlying obfuscation algorithm, thus the above three random variables form a Markov chain: Here, our aim is to derive the posterior distribution for the population parameters θ by observing an instance of the anonymized table, x * . There is no tractable analytical expression for the likelihood L(θ ; x * ). In Kifer (2009) and Boreale et al. (2020) this problem is circumvented by defining an MCMC scheme for sampling from the joint posterior π(θ, x|x * ) and then discarding the cleartext tables (further details on the MCMC scheme are in Boreale et al. (2020, Section 5). Here we pursue an alternative solution based on ABC. Specifically, we simulate the anonymized tables y * according to the following generative model: 1. Generate a table of n i.i.d. rows (s, r ) ∈ S × R distributed according to the θ given as input; 2. Partition the table into k groups of dimensions {n i } k i=1 ; 3. Randomly permute the values of the nonsensitive attributes, r 1 , . . . , r n i , within each group g i .
Here, n is the number of rows and k is the number of groups in the observed anonymized table x * , while n i is the number of rows in g * i . In our example, the observed data x * consists of an obfuscated table composed of k = 1423 groups, with n i = 4 for each group g i . We take race (four possible values) as a nonsensitive attribute and workclass (four possible values) as the sensitive attribute. As in Boreale et al. (2020), we assume that θ S and the θ R|s 's are independently distributed according to non-informative Dirichlet prior distributions. The output of the ABC algorithms will be a sample from the approximate joint posterior distributionπ(θ R|S , T y * |T x * ) since the sensitive part is not changed by the anonymization algorithm and the posterior distribution for θ S exists in closed form. 3 Note that in order to satisfy the assumptions of Sanov's theorem, the m simulated rows must be independent and identically distributed. Recalling that the m pairs (s, r )'s are generated independently and that the permutation is completely at random, we conjecture that the i.i.d. assumption is satisfied. We have positively verified this assumption empirically via the permutation test based on the periodicity test statistic described in the National Institute of Standards and Technology Special Publication 800-90B (Turan et al. 2018).
As in the previous experiment, we consider the output of 100,000 MCMC runs as a reference, in order to compute the MSE's and MISE's, and thus comparing the accuracy of LDW-ABC and R-ABC. The posterior means derived via MCMC are displayed in "Appendix C" (Table 8).
By setting m = 100 and = 1, as far as point estimations are concerned, both LDW-ABC and R-ABC perform quite similarly. The results are displayed in "Appendix C" (Table 9). Nevertheless, the MSE's achieved by LDW-ABC are almost always smaller than the ones achieved by R-ABC. Concerning the approximations of the multivariate posterior distributions, looking at MISE we can conclude that LDW-ABC outperforms R-ABC (see Table 6). Moreover, by focusing on the improvement in efficiency, we note that the value of ESS for LDW-ABC is more than twice that for R-ABC.

Conclusions and future research
We have put forward an approach to address sample degeneracy in methods of approximate Bayesian computation (ABC). Our proposal consists in the definition of a convenient kernel function which, via the theory of large deviations, takes into account the probability of rare events-a poor parameter proposal generating pseudodata close to those observed. By adopting the information theoretic method of types, which involves summarizing data via their empirical distributions, we also by-pass the issue of selecting summary statistics. The proposed kernel function, being defined on a non-compact support, avoids any implicit or explicit rejection step, thus effectively increasing the effective sample size, as empirically verified in Sect. 5. Moreover, the resulting approximate ABC likelihood leads to a better approximation of the tails of the posterior distributions, that is, poor parameter proposals are assigned small but nonzero probability. We also provide formal guarantees of the convergence of our ABC approximate likelihood to the true likelihood.
The proposed method deals with the inefficiency in ABC algorithms by focusing on the kernel function. Although a variety of ABC sampling schemes addressing the same problem have been proposed, most of them (e.g., MCMC-ABC, SMC-ABC, PMC-ABC, SIS-ABC, etc.) handle the problem of finding a good marginal importance distribution, q(θ ), completely ignoring the choice of the kernel K (·). We speculate that these two approaches can be combined by adopting those sampling schemes rather than the involved IS-ABC. Finally, further research is called for in order to apply the proposed method to sequences of dependent random variables, such as Markov chains, and to continuous data. This will in turn require considering more sophisticated versions of the theory of large deviations, where the i.i.d. assumption is relaxed. Such extensions are required to make the algorithm applicable to more complex situations in which no other ways of sampling from the posterior distribution are available. However we want to emphasize that at the current stage, the method is already applicable in contexts in which the other available sampling methods (e.g., MCMC) can be computationally demanding. For example, an efficient ABC algorithm may be usefully applied to models involving high-dimensional latent variables even when an MCMC algorithm exists. In fact, ABC algorithms can be run in parallel, leading to a computational gain when many of the existing MCMC algorithms suffer from a slow mixing of the chain.
Funding Open access funding provided by Università degli Studi di Firenze within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A Proofs
In what follows, we will make use of a few basic notions and facts about the method of types and information projections, for which we refer the reader to (Csiszár et al. 2004, Chap. 1). The simplex of the distributions over X, given a subset Δ |X|−1 ⊆ R |X| , inherits the standard topology from R |X| . W.r.t. this topology, the function D(P||Q) is lower semi-continuous in the pair of arguments (P, Q), and continuous at (P, Q) whenever Q has full support, that is, whenever supp(Q) = {r ∈ X : Q(r ) > 0} = X. Convergence to Q in KL divergence, D(Q n ||Q) → 0, implies convergence in the standard topology, Q n → Q. As a function of P, D(P||Q) is strictly convex, and continuous whenever Q is full support. Hence for any convex and closed set E ⊆ Δ |X|−1 the information projection of Q onto E, P * = argmin P∈E D(P||Q), exists and is unique. The following is a fundamental result about information projections.
The support of E is defined as supp(E) = P∈E supp(P).

Proof of Theorem 3
Fix an infinite sequence τ ∈ X ∞ , τ = (y 1 , y 2 , y 3 , . . .). For each m ≥ 1, let T y m (τ ), T y m for short, denote the type of the first m symbols of τ , the sequence (y 1 , . . . , y m ). Assume τ is such that T y m → P θ as m → +∞. Note that, since P θ is full support, this implies that for all sufficiently large m, T y m is full support as well. Define P * and, for any such sufficiently large m, P * m as follows: Note that as T y m is full support and B = {P : D(P||T x n ) ≤ } is convex and closed, the projection P * m exists and is unique. Moreover, as T x n is by assumption full support, it is easily seen that supp(B ) = X: hence, by the first part of Theorem A.1, the projection P * m is full support as well. As B is closed and D(·||P θ ) is continuous, P * ∈ B . We can now apply the Pythagorean Inequality, considering P * m as a projection and P * ∈ E = B , and obtain As P θ is assumed to be full support, D(·||·) as a function of its second argument is continuous at P θ , hence lim m→∞ D(P * ||T y m ) = D(P * ||P θ ).
Assuming {P * m } converges, let P * * = lim m→∞ P * m , where clearly P * * ∈ B ; if {P * m } does not converge, we can equivalently take any convergent subsequence of it. Taking lim inf on both sides of (27), and exploiting (28) on the left-hand side, and lower semi-continuity on the right-hand side, we can write Summing up D(P * ||P θ ) ≥ D(P * ||P * * ) + D(P * * ||P θ ).
Recalling that P * is the information projection of P θ onto B , that P * * ∈ B and that D(·||·) is nonnegative, the only possibility for (29) to hold is that D(P * ||P * * ) = 0, which implies P * = P * * . In other words lim m→∞ P * m = P * .
Proof The fact that E is closed and convex ensures that the projection D(E||Q) exists and is finite. Consider the strictly descending chain of balls of radius δ n = 1/n centered at Q: By contradiction, assume that there exists 0 < γ < γ such that for each δ > 0, there is Q ∈ B δ (Q) such that D(E|| Q n ) < γ . In particular, we then have that for each n ≥ 1 there is Q n ∈ B δ n (Q) s.t. D(E|| Q n ) < γ .
We can therefore assume without loss of generality that lim n→∞ D(E||Q n ) < γ .
(if not, we can anyway extract from {D(E|| Q n )} a subsequence with the desired property). On the other hand, being lim n→∞ D(Q n ||Q) = 0, we have lim n→∞ Q n = Q. Being D(·||·) lower semi-continuous, we obtain But this contradicts (32).

Proof of Proposition 1
Let us consider the two cases separately, where the last inequality follows from a direct application of Sanov's Theorem.
-P θ ∈ B c . Choose any 0 < γ < γ = D(B c ||T x n ) (note that γ > 0) and apply Lemma A.1 with E = B and Q = P θ to obtain δ > 0 such that D(B ||Q ) ≥ γ for each Q ∈ B δ (P θ ). We can assume without loss of generality that δ ≤ γ . It follows that where (38) follows from (11) and the last step follows from an upper bound for the size of T m (see (Cover and Thomas 2006, Ch. 11, Th. 11.1.1)).

Proof of Proposition
for a function C that is > 0 in the domain of definition of ESS. Therefore, ∂ ∂ x i ESS is nonnegative when evaluated at any point (x 1 , . . . , x S ) in the domain of ESS with the following property: for each j = i s.t. x j > 0, one has 0 ≤ x i ≤ x j . If additionally at least one j = i exists s.t. x j > x i , then ∂ ∂ x i ESS is strictly positive. An execution of the IS algorithm consists of S ≥ 1 independent iterations of the main loop: let us denote by ω s and ρ s the unnormalized weights (6) generated using the LDW-ABC and IS-ABC kernel functions, respectively, at iteration s = 1, . . . , S, and by ω = (ω 1 , . . . , ω S ) and ρ = (ρ 1 , . . . , ρ S ) the resulting sequences. By definition, the set of indices s = 1, . . . , S can be partitioned into three subsets: the subset A where ρ s = ω s > 0, the subset B where ρ s = 0 and ω s > 0, and the subset C where ρ s = ω s = 0. Moreover, for each s ∈ A and s ∈ B, ω s > ω s . For notational simplicity, assume A = {1, . . . , h}, B = {h + 1, . . . , S } and C = {S + 1, . . . , S}, for some 0 ≤ h ≤ S ≤ S. Also assume, again only for notational simplicity, that ω h+1 ≥ ω h+2 ≥ · · · .

B Minimization of the KL divergence
In the proposed LDW-ABC, the minimization of the KL divergence between the acceptance region and the simulated type poses a computational difficulty. This is a constrained minimization problem on a space of dimension |X|. As |X| grows, this problem can rapidly become intractable. A practical work-around to this problem can be found by considering a suitable path from T y m to T x n , passing through P * . In Information Geometry, this path is represented by a linear interpolation on the logarithmic scale, the exponential geodesic (Nielsen 2018).
Definition 3 (Exponential geodesic) Let P 1 and P 2 be two probability distributions over X and let P ξ be the probability distribution such that for each r ∈ X log P ξ (r ) = ξ log P 1 (r ) + (1 − ξ) log P 2 (r ) + log c where ξ ∈ [0, 1] and c is a proper normalizing constant. The exponential geodesic between P 1 and P 2 is the following set of distributions γ e P 1 , P 2 = {P ξ : ξ ∈ [0, 1]}.
Our approach when minimizing the KL divergence between B (T x n ) and T y m is to focus on a path between the observed and the simulated type, that is the exponential geodesic γ e (T x n , T y m ). We search in this path the information projection P * , or an approximation of it. This reduces the dimension of the minimization problem from |X| to 1, that of the parameter ξ . Specifically, let P ξ * ∈ B (T x n ) be the element of γ e (T x n , T y ) defined as P ξ * (r ) = T x n (r ) ξ * · T y m (r ) 1−ξ * c * (r ∈ X) where ξ * = argmin ξ ∈[0,1] : T ξ ∈B D(P ξ ||T y m ).
We have empirically verified that D(P ξ * ||T y m ) approximates with very good accuracy D(B ||T y m ).
Hence, whatever |X|, D(B ||T y m ) is approximate by means of a minimization with respect to a single parameter, ξ . Table 7 summarizes the distribution of the distances approximation relative errors w.r.t. the true distance, over the S = 100,000 simulations in the experiment in Sect. 5.1, with m = 500 and = 0.005.   Table 8 shows the posterior means for each θ R|s estimated via MCMC. Such estimates are used as a benchmark in the computation of the MSE's shown in Table 9.