Information geometry for approximate Bayesian computation

TL;DR

This paper applies information geometry and relative entropy to analyze and improve the Approximate Bayesian Computation (ABC) algorithm, providing theoretical insights and practical guidelines for its performance and sensitivity.

Contribution

It introduces an information-theoretic framework for analyzing ABC, including bias, error bounds, and sensitivity, which aids in optimizing the algorithm's efficiency.

Findings

Relative entropy analysis clarifies ABC behavior as a function of thresholds and data size.

The framework identifies important directions for reducing computational costs.

Simulation results validate theoretical error bounds and sensitivity analysis.

Abstract

The goal of this paper is to explore the basic Approximate Bayesian Computation (ABC) algorithm via the lens of information theory. ABC is a widely used algorithm in cases where the likelihood of the data is hard to work with or intractable, but one can simulate from it. We use relative entropy ideas to analyze the behavior of the algorithm as a function of the threshold parameter and of the size of the data. Relative entropy here is data driven as it depends on the values of the observed statistics. Relative entropy also allows us to explore the effect of the distance metric and sets up a mathematical framework for sensitivity analysis allowing to find important directions which could lead to lower computational cost of the algorithm for the same level of accuracy. In addition, we also investigate the bias of the estimators for generic observables as a function of both the threshold parameters and the size of the data. Our analysis provides error bounds on performance for positive tolerances and finite sample sizes. Simulation studies complement and illustrate the theoretical results.

Figures

Tables

Table 1. Table 1. Data for ball as acceptance region.

$(\text{tol},n)$	${\tau }^{\ast }=(\overline{X},S^2)$	$\epsilon$	$(\hat{\mu },\widehat{\text{sd}}\left(\mu \right))$	$(\widehat{{\sigma }^2},\widehat{\text{sd}}\left({\sigma }^2\right))$	${\hat{R}}_B$
$(0.05,100)$	$(0.167,1.061)$	$0.055$	$(0.161,0.109)$	$(1.035,0.075)$	$1484$
$(0.25,100)$	$(-0.022,0.965)$	$0.083$	$(-0.023,0.106)$	$(0.987,0.076)$	$608$
$(0.5,100)$	$(0.061,1.099)$	$0.111$	$(0.059,0.122)$	$(1.053,0.079)$	$384$
$(1,100)$	$(-0.181,0.866)$	$0.096$	$(-0.175,0.105)$	$(0.936,0.071)$	$389$
$(0.05,300)$	$(-0.112,1.003)$	$0.025$	$(-0.111,0.061)$	$(1.007,0.041)$	$6998$
$(0.25,300)$	$(0.022,0.974)$	$0.038$	$(0.024,0.061)$	$(0.989,0.041)$	$2902$
$(0.5,300)$	$(0.020,1.001)$	$0.046$	$(0.021,0.062)$	$(1.003,0.042)$	$1982$
$(1,300)$	$(-0.070,1.005)$	$0.054$	$(-0.069,0.063)$	$(1.007,0.043)$	$1490$
$(0.05,600)$	$(-0.014,1.025)$	$0.016$	$(-0.013,0.004)$	$(1.015,0.030)$	$17645$
$(0.25,600)$	$(0.023,0.933)$	$0.022$	$(0.023,0.042)$	$(0.968,0.028)$	$8820$
$(0.5,600)$	$(-0.009,0.986)$	$0.027$	$(-0.010,0.041)$	$(0.994,0.029)$	$5805$
$(1,600)$	$(0.016,0.972)$	$0.003$	$(0.014,0.004)$	$(0.986,0.028)$	$4172$
$(0.05,1000)$	$(-0.027,1.013)$	$0.011$	$(-0.027,0.031)$	$(1.017,0.023)$	$35639$
$(0.25,1000)$	$(0.016,0.971)$	$0.015$	$(0.017,0.032)$	$(0.986,0.022)$	$17730$
$(0.5,1000)$	$(-0.050,0.938)$	$0.0175$	$(-0.051,0.032)$	$(0.986,0.023)$	$13868$
$(1,1000)$	$(-0.012,0.995)$	$0.022$	$(-0.012,0.033)$	$(0.992,0.023)$	$8702$

Table 2. Table 2. Data for ellipse as acceptance region.

$(\text{tol},n)$	${\tau }^{\ast }=(\overline{X},S^2)$	$({\epsilon }_1,{\epsilon }_2)$	$(\hat{\mu },\widehat{\text{sd}}\left(\mu \right))$	$(\widehat{{\sigma }^2},\widehat{\text{sd}}\left({\sigma }^2\right))$	${\hat{R}}_E$
$(0.05,100)$	$(0.167,1.061)$	$(0.065,0.049)$	$(0.167,0.103)$	$(1.037,0.074)$	$1505$
$(0.25,100)$	$(-0.022,0.965)$	$(0.155,0.069)$	$(-0.027,0.121)$	$(0.991,0.077)$	$390$
$(0.5,100)$	$(0.061,1.099)$	$(0.173,0.094)$	$(0.063,0.133)$	$(1.061,0.082)$	$303$
$(1,100)$	$(-0.081,0.866)$	$(0.113,0.086)$	$(-0.178,0.106)$	$(0.937,0.072)$	$370$
$(0.05,300)$	$(-0.112,1.003)$	$(0.035,0.021)$	$(-0.113,0.061)$	$(1.004,0.042)$	$5660$
$(0.25,300)$	$(0.022,0.974)$	$(0.087,0.031)$	$(0.194,0.071)$	$(0.989,0.041)$	$1559$
$(0.5,300)$	$(0.020,1.001)$	$(0.106,0.038)$	$(0.024,0.077)$	$(0.998,0.040)$	$1102$
$(1,300)$	$(-0.070,1.005)$	$(0.091,0.046)$	$(-0.074,0.075)$	$(1.006,0.045)$	$1105$
$(0.05,600)$	$(-0.014,1.025)$	$(0.043,0.013)$	$(-0.015,0.046)$	$(1.014,0.028)$	$7617$
$(0.25,600)$	$(0.023,0.933)$	$(0.056,0.018)$	$(0.022,0.049)$	$(0.966,0.029)$	$4045$
$(0.5,600)$	$(-0.009,0.986)$	$(0.078,0.023)$	$(-0.006,0.056)$	$(0.991,0.029)$	$2679$
$(1,600)$	$(0.016,0.972)$	$(0.087,0.026)$	$(0.016,0.059)$	$(0.986,0.031)$	$1782$
$(0.05,1000)$	$(-0.027,1.033)$	$(0.028,0.009)$	$(-0.026,0.035)$	$(1.018,0.023)$	$18043$
$(0.25,1000)$	$(0.016,0.971)$	$(0.045,0.013)$	$(0.017,0.039)$	$(0.987,0.023)$	$7085$
$(0.5,1000)$	$(-0.050,0.938)$	$(0.036,0.015)$	$(-0.011,0.048)$	$(0.993,0.023)$	$6270$
$(1,1000)$	$(-0.012,0.985)$	$(0.068,0.018)$	$(-0.011,0.048)$	$(0.993,0.023)$	$3461$

Table 3. Table 3. Asymptotic and Monte Carlo comparison for the two different choices of acceptance region (ball and ellipse).

$(\text{tol},n)$	$(0.05,100)$	$(0.25,100)$	$(0.5,100)$	$(1,100)$	$(0.05,300)$	$(0.25,300)$	$(0.5,300)$	$(1,300)$
$\tilde{U}\left({\tau }^{\ast }\right)$	$0.986$	$0.763$	$0.884$	$1.041$	$0.947$	$0.648$	$0.651$	$0.861$
$\frac{{\hat{R}}_E}{{\hat{R}}_B}$	$1.014$	$0.641$	$0.790$	$0.949$	$0.809$	$0.537$	$0.556$	$0.741$
$(\text{tol},n)$	$(0.05,600)$	$(0.25,600)$	$(0.5,600)$	$(1,600)$	$(0.05,1000)$	$(0.25,1000)$	$(0.5,1000)$	$(1,1000)$
$\tilde{U}\left({\tau }^{\ast }\right)$	$0.531$	$0.592$	$0.540$	$0.574$	$0.592$	$0.532$	$0.574$	$0.511$
$\frac{{\hat{R}}_E}{{\hat{R}}_B}$	$0.431$	$0.462$	$0.461$	$0.436$	$0.506$	$0.401$	$0.452$	$0.397$

Equations

$$H(P|Q)={\int }_{\mathcal{X}}\log \frac{dP(x)}{dQ(x)}dP(x)$$

$$\text{maximize }\epsilon -\text{dependent acceptance region subject to }H(P|P^{\epsilon })({\tau }^{\ast })\le \text{tol},$$

$$H(P|P^{\epsilon })({\tau }^{\ast })$$

$$\text{maximize }\epsilon -\text{dependent acceptance region subject to }$$

$$\frac{1}{8}\frac{1}{(q+2{)}^2{f}^2({\tau }^{\ast })}{\mathbb{E}}_{f(\theta |{\tau }^{\ast })}\left[{\left(\sum _{i=1}^q{\epsilon }_i^2{w}_i({\tau }^{\ast },\theta )\right)}^2\right]\le \text{tol}.$$

$${\widehat{{\pi }_X[h]}}_{\epsilon ,K}=\frac{1}{K}\sum _{i=1}^{K}h({\theta }_i^{\epsilon }),$$

$$f_P(\theta )=f_{\Theta |T}(\theta |{\tau }^{\ast })=\frac{f_{T,\Theta }({\tau }^{\ast },\theta )}{f_T({\tau }^{\ast })}.$$

$$f_{P^{\epsilon }}(\theta )=f^{\epsilon }(\theta |{\tau }^{\ast })=\frac{\frac{1}{|D_{\epsilon }({\tau }^{\ast })|}{\int }_{D_{\epsilon }({\tau }^{\ast })}{f}_{T,\Theta }(\tau ,\theta )d\tau }{Z},$$

$$Z=\frac{1}{|D_{\epsilon }({\tau }^{\ast })|}{\int }_{D_{\epsilon }({\tau }^{\ast })}{f}_T(\tau )d\tau ,$$

$$D_{\epsilon }({\tau }^{\ast })$$

$$D_{\epsilon }({\tau }^{\ast })=\left\{\tau \in {\mathbb{R}}^q:\sum _{i=1}^q{\left|\frac{{\tau }_i-{\tau }_i^{\ast }}{{\epsilon }_i}\right|}^2\le 1\right\}.$$

$$H(P|P^{\epsilon })({\tau }^{\ast })$$

$$w_i({\tau }^{\ast },\theta )$$

$$=\frac{{\partial }_{{\tau }_i}^{2}f(\theta |{\tau }^{\ast })}{f(\theta |{\tau }^{\ast })}f({\tau }^{\ast })+2\frac{{\partial }_{{\tau }_i}f(\theta |{\tau }^{\ast })}{f(\theta |{\tau }^{\ast })}{\partial }_{{\tau }_i}f({\tau }^{\ast }).$$

$$H(P|P^{\epsilon })({\tau }^{\ast })$$

$$=-\int \log \left(1+\frac{f^{\epsilon }(\theta |{\tau }^{\ast })-f(\theta |{\tau }^{\ast })}{f(\theta |{\tau }^{\ast })}\right)f(\theta |{\tau }^{\ast })d\theta$$

$$=-\int \left(\frac{f^{\epsilon }(\theta |{\tau }^{\ast })-f(\theta |{\tau }^{\ast })}{f(\theta |{\tau }^{\ast })}-\frac{1}{2}{\left(\frac{f^{\epsilon }(\theta |{\tau }^{\ast })-f(\theta |{\tau }^{\ast })}{f(\theta |{\tau }^{\ast })}\right)}^2+\mathcal{O}(|f^{\epsilon }-f{|}^3)\right)f(\theta |{\tau }^{\ast })d\theta$$

$$=\frac{1}{2}\int \left({\left(\frac{f^{\epsilon }(\theta |{\tau }^{\ast })-f(\theta |{\tau }^{\ast })}{f(\theta |{\tau }^{\ast })}\right)}^2+\mathcal{O}(|f^{\epsilon }-f{|}^3)\right)f(\theta |{\tau }^{\ast })d\theta$$

$$=\frac{1}{2}\int \left({\left(\frac{{\epsilon }^T{\nabla }_{\epsilon }{f}^{\epsilon }(\theta |{\tau }^{\ast }){|}_{\epsilon =0}}{f(\theta |{\tau }^{\ast })}+\frac{1}{2}\frac{{\epsilon }^T{\nabla }_{\epsilon }^2{f}^{\epsilon }(\theta |{\tau }^{\ast }){|}_{\epsilon =0}\epsilon }{f(\theta |{\tau }^{\ast })}\right)}^2+\mathcal{O}(|\epsilon {|}^6)\right)f(\theta |{\tau }^{\ast })d\theta$$

$$=\frac{1}{2}\int \frac{{\epsilon }^T{\nabla }_{\epsilon }{f}^{\epsilon }(\theta |{\tau }^{\ast }){|}_{\epsilon =0}{\left({\nabla }_{\epsilon }{f}^{\epsilon }(\theta |{\tau }^{\ast }){|}_{\epsilon =0}\right)}^T\epsilon }{f(\theta |{\tau }^{\ast })}d\theta$$

$$+\frac{1}{2}\int \frac{{\epsilon }^T{\nabla }_{\epsilon }{f}^{\epsilon }(\theta |{\tau }^{\ast }){|}_{\epsilon =0}{\epsilon }^T{\nabla }_{\epsilon }^2{f}^{\epsilon }(\theta |{\tau }^{\ast }){|}_{\epsilon =0}\epsilon }{f(\theta |{\tau }^{\ast })}d\theta +$$

$$+\frac{1}{8}\int \frac{{\left|{\epsilon }^T{\nabla }_{\epsilon }^2{f}^{\epsilon }(\theta |{\tau }^{\ast }){|}_{\epsilon =0}\epsilon \right|}^2}{f(\theta |{\tau }^{\ast })}d\theta +\mathcal{O}(|\epsilon {|}^6)$$

$$=\frac{1}{2}{\mathbb{E}}_{f(\cdot |{\tau }^{\ast })}\left[\frac{{\epsilon }^T{\nabla }_{\epsilon }{f}^{\epsilon }(\theta |{\tau }^{\ast }){|}_{\epsilon =0}{\left({\nabla }_{\epsilon }{f}^{\epsilon }(\theta |{\tau }^{\ast }){|}_{\epsilon =0}\right)}^T\epsilon }{f^2(\theta |{\tau }^{\ast })}\right]+$$

$$+\frac{1}{2}{\mathbb{E}}_{f(\cdot |{\tau }^{\ast })}\left[\frac{{\epsilon }^T{\nabla }_{\epsilon }{f}^{\epsilon }(\theta |{\tau }^{\ast }){|}_{\epsilon =0}{\epsilon }^T{\nabla }_{\epsilon }^2{f}^{\epsilon }(\theta |{\tau }^{\ast }){|}_{\epsilon =0}\epsilon }{f^2(\theta |{\tau }^{\ast })}\right]+$$

$$+\frac{1}{8}{\mathbb{E}}_{f(\cdot |{\tau }^{\ast })}\left[\frac{{\left|{\epsilon }^T{\nabla }_{\epsilon }^2{f}^{\epsilon }(\theta |{\tau }^{\ast }){|}_{\epsilon =0}\epsilon \right|}^2}{f^2(\theta |{\tau }^{\ast })}\right]+\mathcal{O}(|\epsilon {|}^6).$$

$$H(P|P^{\epsilon })({\tau }^{\ast })$$

$${\mathbb{E}}_{f(\theta |T={\tau }^{\ast })}\left[w_i({\tau }^{\ast },\theta )\right]=0.$$

$${\mathbb{E}}_{f(\theta |T={\tau }^{\ast })}\left[\frac{{\left|{\epsilon }^T{\nabla }_{\epsilon }^2{f}^{\epsilon }(\theta |{\tau }^{\ast }){|}_{\epsilon =0}\epsilon \right|}^2}{f^2(\theta |{\tau }^{\ast })}\right]={\epsilon }^4\frac{1}{(q+2{)}^2{f}^2({\tau }^{\ast })}{\mathbb{E}}_{f(\theta |T={\tau }^{\ast })}\left[\frac{{\left({\Delta }_{\tau }f({\tau }^{\ast },\theta )-\frac{f({\tau }^{\ast },\theta )}{f({\tau }^{\ast })}{\Delta }_{\tau }f({\tau }^{\ast })\right)}^2}{f^2(\theta |{\tau }^{\ast })}\right]$$

$$f(x|\theta )=\zeta (x)\exp \{\sum _{i=1}^q{\eta }_i(\theta )S_i(x)-A(\theta )\}.$$

$$f(x_1,\cdots ,x_n|\theta )=\prod _{j=1}^n\zeta (x_j)\exp \{\sum _{i=1}^q{\eta }_i(\theta )Y_i(x)-nA(\theta )\},$$

$$f(y_1,\cdots ,y_q|\theta )=R(y_1,\cdots ,y_q)\exp \{\sum _{i=1}^q{\eta }_i(\theta )y_i-nA(\theta )\},$$

$$f({\tau }_1,\cdots ,{\tau }_q|\theta )=n^{q}R(n{\tau }_1,\cdots ,n{\tau }_q)\exp \{n\sum _{i=1}^q{\eta }_i(\theta ){\tau }_i-nA(\theta )\}.$$

$$f(\theta |\rho ,\kappa )=Z(\rho ,\kappa )\exp \{\sum _{i=1}^q{\eta }_i(\theta ){\rho }_i-\kappa A(\theta )\},$$

Full text

The author would like to thank the Department of Mathematics and Statistics of UMASS Amherst for the hospitality and excellent working environment during the fall of 2016, when most of this work was initiated. In particular, the author would like to thank Markos Katsoulakis and Luc Rey-Bellet for a number of useful and fruitful discussions and suggestions on this project. This work was partially supported by the National Science Foundation (NSF) CAREER award DMS 1550918.

1. Introduction

The goal of this paper is to explore Approximate Bayesian Computation (ABC) methods via the lens of information theoretic criteria offering a more geometric point of view on such methods. ABC is a very popular likelihood-free Bayesian inference method in many fields of science ranging from population genetics to genetic evolutions and material science, see for example [20, 19, 8, 9, 4, 23] to name just a few.

In ABC, the goal is to do Bayesian inference on a parameter $\theta\in\mathbb{R}^{m}$ that has some prior distribution, say $f(\theta)$, given that we observe data $X\in\mathbb{R}^{n}$. The density $f(x|\theta)$ is usually either not available or not in a tractable form. The data is usually summarized in form of appropriate statistics, $T(X):\mathbb{R}^{n}\mapsto\mathbb{R}^{q}$. The basic idea of ABC is to draw samples of $\theta$ from the prior distribution $f(\theta)$, then propagate the information forward producing simulated values of the sufficient statistic, say $T_{s}$ based on the distribution of $X|\theta$. Then if the simulated $T_{s}$ is close to the observed value $\tau^{*}=T_{o}=T(X)$ of the observed statistic, the sample $\theta$ is accepted, otherwise it is rejected. The procedure is then repeated a sufficiently large number of times allowing estimation of $\theta$ or of observables of interest.

ABC methods have been very popular, mainly due to their generality and simplicity in application. However, the generality and likelihood-free advantage that ABC has, comes with known issues on sensitivity in performance with respect not only to the distance metric $\|\cdot\|$ being chosen to compare simulated versus observed values, $T_{s}$ and $T_{o}$ respectively, but also on the choice of the threshold parameter $\varepsilon$ measuring whether the distance is large or not. In addition, the practical performance of ABC also depends on the choice of the observed statistics. One would like the statistics to be sufficient statistics, but that is usually hard to come up with in practice. It is also known that the basic ABC algorithm typically scales badly with the size of the data, say $n$, even though the use of sufficient statistics (when available) helps reduce the impact of big values of $n$. For these reasons, many clever modifications to the original algorithm have been proposed in the literature, including replacing the rejection mechanism by a kernel weighting the samples and doing localized regression [4], placing ABC in MCMC context [13], Gibbs sampling in ABC [23], using non-linear regression models in ABC [7] and others. Excellent tutorials on the different variants of ABC methods are [12, 21].

Despite significant developments in methodological approaches, theoretical developments are still in their infancy. Notable exceptions to this are the recent works [3, 6] where the behavior of the bias as a function of the threshold parameter for the Euclidean metric is explored and the expected cost of the basic version of the algorithm and limiting properties are studied.

The focus of our paper is different. We take a geometric point of view based on studying the relative entropy between the underlying and the simulated posterior distribution of $\theta|T(X)$. The relative entropy between a measure $P$ with respect to a measure $Q$ on a Euclidean space $\mathcal{X}$ is defined to be

[TABLE]

when the integral exists. The relative entropy is not a measure but it is a divergence in that it satisfies $H(P|Q)\geq 0$ and $H(P|Q)=0$ if and only if $P=Q$ almost everywhere. Information theoretic criteria have been successfully used in uncertainty quantification and sensitivity analysis, see for example [1, 2, 11, 15, 16].

This geometric point of view focuses on the entire distribution rather than on specific statistics. At the same time, it gives us the flexibility to quantify the effect of the distance metric and of the threshold parameter. For example, we can consider the threshold parameter to be a vector, i.e., $\varepsilon\in\mathbb{R}^{q}_{+}$ instead of a scalar. As we show, the latter allows us to study sensitivity of the relative entropy with respect to the different components of the vector of observed statistics. In practice this means that if certain directions are insensitive then one can afford having larger values of threshold parameters $\varepsilon_{i}$ in such directions. The latter can increase the acceptance region, allowing ABC to run more efficiently.

In this paper we study how the relative entropy of the underlying measure, say $P$, with respect to the simulated measure, say $P^{\varepsilon}$, scales as a function of $\varepsilon$ for $|\varepsilon|$ small enough. Then we can choose $\varepsilon$ such that the acceptance region is maximized for a fixed level of tolerance for the relative entropy.

The contribution of the paper is fourfold. Firstly, in Section 2 we compute the leading order term of the relative entropy, characterizing explicitly its dependence on the size of the vector $T$, $q$, and the threshold vector $\varepsilon\in\mathbb{R}^{q}_{+}$. Secondly, as discussed, we allow $\varepsilon\in\mathbb{R}^{q}_{+}$ to be multidimensional. We compute the weighted effect in each direction $i=1,\cdots,q$ expressed in terms of the posterior distribution $f(\theta|T(x))$ and make the computation explicit for the exponential family in Section 3. We demonstrate in specific examples that by exploiting the asymmetric effect of each direction, one is able to enlarge the acceptance region of ABC. In Section 5 we compute the theoretical mean rejection rate of ABC in both cases of $\varepsilon\in\mathbb{R}^{q}_{+}$ and of $\varepsilon\in\mathbb{R}^{1}_{+}$ and we compare the two. We find that the mean rejection rate is typically smaller in the case of $\varepsilon\in\mathbb{R}^{q}_{+}$.

As we shall see in Sections 6 and 7 allowing a threshold parameter to be different for each statistic $\tau_{i}$ can have measurable effects in the performance of the algorithm. For example, in Figure 1, we visually compare the acceptance regions for the ABC algorithm in the asymmetric case (i.e. acceptance region is an ellipse) and the symmetric case (i.e acceptance region is a ball) for the simulated example of Section 7. The improvement of performance is because the first direction (in this example, the mean) has less significance as compared to the second direction (in this example, the variance) and thus by taking advantage of this we can allow for a larger acceptance region. The overall tolerance error is the same in both cases.

Of course, allowing the threshold parameter to be potentially different for each different statistic $\tau_{i}$, as opposed to forcing them to be the same, will naturally lead to larger acceptance regions. As we discuss in more detail in Section 2 such choices (ellipse vs ball etc) correspond to special cases of distance metrics. Other choices of distance metrics are of course possible. In this paper we focus on explicitly characterizing the performance of ABC under the choice of the Euclidean metric being the distance metric via the lens of relative entropy. The relative entropy formulation gives us a mathematical framework to quantify the effect of the distance metric.

Third contribution of the paper is that we investigate in Section 4 the behavior of the bias of estimators as a function of both $\varepsilon\in\mathbb{R}^{q}_{+}$ and of the size of the data $n$. We compute the weighted effect from each direction $i=1,\cdots,q$ on the bias. We show that the weighted effects on both relative entropy and bias are based on the same weight functions.

Fourthly, in the case of exponential family and for sufficient statistics, we show that the typical order of the relative entropy is $n^{4}\sum_{i=1}^{q}\varepsilon_{i}^{4}/q$, whereas the typical order of the bias is $n^{2}\sum_{i=1}^{q}\varepsilon_{i}^{2}/q$. The latter naturally means that as $n$ increases, one should choose smaller values for $\varepsilon_{i}$. However, as we show in detail, this is a crude upper bound coming from the dominant term of the relative entropy and improvements are usually possible if one takes into account the specific effects of the weight functions. In particular, as we show in Section 6, the relative entropy in many cases can be of much smaller order in $n$ and the bias of certain observables of interest may even be independent of $n$.

We remark here that this inverse relation between $n$ and $|\varepsilon|$ is to be expected. After all, as $n\rightarrow\infty$ law of large numbers guarantees that sufficient statistics will converge to the true values of the corresponding parameters that they estimate. Hence, as $n$ gets larger, only samples that result in values for the simulated statistics that are closer to the values of the related observed statistics should be accepted. Essentially, this means that as $n$ increases, $|\varepsilon|$ should decrease as the variability in the estimation gets smaller. In this paper we quantify this inverse relation precisely.

Mathematically speaking, if we pick a level of tolerance, say tol, we study the following problem

[TABLE]

where $P$ is the posterior distribution with density $f(\theta|T(X)=\tau^{*})$ and $P^{\varepsilon}$ is the simulated posterior distribution with density $f^{\varepsilon}(\theta|T(X)=\tau^{*})$.

It is easy to see though that (1.1) is in general intractable. However, since in ABC the interest is in small values of $|\varepsilon|$ it makes sense to go to the limit as $|\varepsilon|\rightarrow 0$. We prove that for $|\varepsilon|$ small enough one has the expansion

[TABLE]

where the weight $w_{i}(\tau^{\ast},\theta)$ can be computed explicitly and in the case of exponential family one can also see the exact dependence on the number of data points $n$. Notice that the relative entropy now depends on the observed value of $\tau^{\ast}$, i.e., it is data driven.

The latter expansion for the constraint means that if one sets a tolerance level, say tol, and requires that $H(P|P^{\varepsilon})(\tau^{\ast})=\text{tol}$ and that $|\varepsilon|$ is sufficiently small, then one should choose $\varepsilon_{i}$ such that

[TABLE]

The latter expression makes the problem tractable and it also makes it clear that the value of $\varepsilon_{i}$ depends on the value of the weight factor $\mathbb{E}_{f(\theta|\tau^{\ast})}\left(w^{2}_{i}(\tau^{\ast},\theta)\right)$ in the $i^{\text{th}}$ direction and on the cross products $\mathbb{E}_{f(\theta|\tau^{\ast})}\left(w_{i}(\tau^{\ast},\theta)w_{j}(\tau^{\ast},\theta)\right)$ for $i\neq j$. This formulation allows us to exploit the asymmetric effect of different directions leading to potentially larger acceptance region for the basic version of the ABC algorithm. Notice also that (1.2) is a constrained optimization problem and in general is non-convex. We refer the interested reader to classical sources, such as [5, 10, 17] for more on numerical solution to constrained optimization as well as non-convex problems.

The paper is organized as follows. In Section 2 we present the general mathematical framework that we consider and the relative entropy computations. We explore the dependence of the relative entropy on the magnitude of $|\varepsilon|$ and characterize the leading order of the expansion as $|\varepsilon|\rightarrow 0$. In Section 3 we make the computations specific for exponential family. This allows to see the dependence on the size of the data $n$, but it also gives a very intuitive interpretation of the weight factors as the variance of the square of the natural parameters in the representation of the density. In Section 4 we calculate the asymptotics as $|\varepsilon|\rightarrow 0$ for the bias of the estimator for given observables characterizing the leading order term. We notice that the leading order term depends on the same weight functions $w_{i}(\tau^{*},\theta)$. In Section 5 we compare analytically the conditional (on the observed statistics) mean rejections rates of the algorithm that uses the ellipse as acceptance region versus the algorithm that uses the ball. Section 6 contains some examples to demonstrate that advantages can come from monitoring the weight factors and from considering different tolerance levels for different sufficient statistics. Section 7 contains the simulations studies that complement and illustrate the theoretical results. Conclusions and potential future work directions are presented in Section 8.

2. Information-theoretic criteria for ABC

Let us assume that the parameter of interest is $\theta\in\Theta\subseteq\mathbb{R}^{m}$ and that we observe a statistic of the data $X\in\mathbb{R}^{d}$, say $T(X)=\tau^{*}:\mathbb{R}^{d}\mapsto\mathbb{R}^{q}$. In this paper and in order to be able to carry over the mathematics, we assume that we are working with sufficient statistics.

Consider a sampling distribution $f_{X|\Theta}(x|\theta)$ and assume a prior density $f_{\Theta}(\theta)$. In the basic ABC algorithm we decide to accept or reject a sample of $\theta$ from the prior density $f_{\Theta}(\theta)$ based on whether the statistic formed by the sampled data based on $f_{X|\Theta}(x|\theta)$, say $T_{s}=T(X)$, is within a level of tolerance, say $\varepsilon$, from the observed value of the statistic, $\tau^{\ast}=T_{o}$.

Apart from the choice of observed statistics, the two main free parameters in ABC are (a): the distance metric $\|\cdot\|$ that we use to measure the distance of $T_{s}$ from $T_{o}$, and (b): the tolerance parameter $\varepsilon$.

The ultimate goal is to compute, for a given $h(\theta)$, quantities of the form $\pi_{X}[h]=\mathbb{E}\left[h(\Theta)|X\right]$. Notice that if $T(\cdot)$ is a sufficient statistic, then we actually have $\pi_{X}[h]=\mathbb{E}\left[h(\Theta)|T(X)\right]$.

Using the accept-reject mechanism just described, the estimator of the posterior distribution is

[TABLE]

where $K$ is the number of accepted samples.

Let us assume now that $T(X)$ and $\theta$ have joint distribution with a well defined everywhere positive density which we denote by $f_{T,\Theta}(\tau,\theta)$. Likewise, let $f_{T}(\tau)$ be the density of the sufficient statistics $T$. We allow for $\varepsilon$ to be a vector in $\mathbb{R}^{q}_{+}$.

Let $P$ be the distribution of $\Theta|T=\tau^{\ast}$ with density

[TABLE]

Let $D_{\varepsilon}(\tau^{\ast})$ be the appropriate acceptance set and let $|D_{\varepsilon}(\tau^{\ast})|$ be its volume. The simulated data is distributed according to the conditional distribution of $\theta$ given $\tau\in D_{\varepsilon}(\tau^{\ast})$, say $P^{\varepsilon}$. This means that in our case $P^{\varepsilon}$ will have density given by

[TABLE]

where

[TABLE]

is the necessary normalizing constant. We remark here that $f_{P^{\varepsilon}}(\theta)$ also reflects the accept/reject mechanism that is in the heart of the ABC algorithm. Namely, if the simulated $\theta$ is such that the corresponding statistic $\tau$ is in $D_{\varepsilon}(\tau^{\ast})$, then the specific simulated value of $\theta$ is accepted. Otherwise, it is rejected. We assume that $\lim_{\varepsilon\rightarrow 0}D_{\varepsilon}(\tau^{\ast})=\{\tau^{\ast}\}$ for all $\tau^{\ast}\in\mathbb{R}^{q}$. In this paper we focus on Euclidean distances and for $\varepsilon\in\mathbb{R}^{q}_{+}$ we set

[TABLE]

where $A(\varepsilon)$ is an $\varepsilon-$dependent positive definite matrix such that $\|\tau\|^{2}_{A(\varepsilon)}=\tau^{T}A^{-1}(\varepsilon)\tau$. We assume that for any value of the observed statistics $\tau^{\ast}$, $A(\varepsilon)$ is such that $\lim_{|\varepsilon|\rightarrow 0}|D_{\varepsilon}(\tau^{\ast})|=0$.

There are of course many ways to choose the matrix $A(\varepsilon)$. One possibility, which is also the usual practice, is to set $A(\varepsilon)=\varepsilon^{2}I$ with $\varepsilon\in\mathbb{R}^{1}_{+}$ in which case we get a ball for acceptance region

[TABLE]

Alternatively, we could pick $\varepsilon\in\mathbb{R}^{q}_{+}$, set $A(\varepsilon)=\text{diag}\left(\varepsilon_{1}^{2},\cdots,\varepsilon_{q}^{2}\right)$ and thus consider the ellipse

[TABLE]

Other choices for the acceptance region $D_{\varepsilon}(\tau^{\ast})$ and distance metrics are of course possible. In this paper we will focus on analyzing the standard choice of a ball (2.4) and comparing it to the case of the ellipse (2.5).

For notational convenience we will be writing $f^{\varepsilon}(\theta|\tau)$ for the perturbed model and $f(\theta|\tau)$ for the true model given $\tau$. We compare the distributions $P$ and $P^{\varepsilon}$ via relative entropy

[TABLE]

which here is naturally data-driven and depends on the observed value of the sufficient statistic $\tau^{\ast}$.

As we discussed in the introduction, we choose the threshold parameters $\varepsilon_{i}$ based on the optimization problem (1.1). Of course, such an optimization problem is difficult to solve as the relative entropy $H(P|P^{\varepsilon})(\tau^{\ast})$ is in general intractable. One can however consider the problem in the limit as $|\varepsilon|\rightarrow 0$. So our first task is to find a more specific expression for the first order approximation to $H(P|P^{\varepsilon})(\tau^{\ast})$ as $|\varepsilon|\rightarrow 0$.

Throughout the paper we will be assuming smoothness of the involved densities. In particular, we have the following assumption.

Condition 2.1**.**

We assume that $f^{\varepsilon}(\theta|\tau),f(\theta|\tau),f_{T}(\tau)$ are smooth, bounded away from zero and that $f^{\varepsilon}(\theta|\tau^{\ast})$ is at least two times continuously differentiable with respect to $\varepsilon$.

As $|\varepsilon|\downarrow 0$, it is easy to see that $f^{\varepsilon}(\theta|\tau^{\ast})$ converges to $f(\theta|\tau^{\ast})$. In particular we have the following proposition.

Proposition 2.2**.**

Assume that Condition 2.1 holds and that the acceptance region $D_{\varepsilon}(\tau^{\ast})$ is given by (2.5). As $|\varepsilon|\rightarrow 0$ we have that

[TABLE]

where the weight function $w_{i}(\tau^{\ast},\theta)$ is

[TABLE]

Proof.

Using Taylor series expansion on $f^{\varepsilon}(\theta|\tau^{\ast})$ with respect to $\varepsilon$, we have up to leading order for $|\varepsilon|$ sufficiently small

[TABLE]

To continue we need to study the contribution of the leading order term, which involves $\nabla_{\varepsilon}f^{\varepsilon}(\theta|\tau^{\ast})|_{\varepsilon=0}$ and $\nabla^{2}_{\varepsilon}f^{\varepsilon}(\theta|\tau^{\ast})|_{\varepsilon=0}$. We consider the case where $D_{\varepsilon}(\tau^{\ast})$ is an ellipse, as given by (2.5).

By Lemma A.1 in the appendix we have that $\nabla_{\varepsilon}f^{\varepsilon}(\theta|\tau^{\ast})|_{\varepsilon=0}=0$. In particular as $|\varepsilon|\downarrow 0$, we then obtain that

[TABLE]

Then, Lemma A.1 again, but for the second derivatives now, gives

[TABLE]

where the weight function $w_{i}(\tau^{\ast},\theta)$ is given by (2.7). This completes the proof of the proposition. ∎

For completeness purposes we mention that a simple calculation shows that for every $i=1,\cdots,q$

[TABLE]

Let us now comment on these results.

Remark 2.3*.*

1. (1)

Notice that the formula for the ellipse, i.e., (2.6)-(2.7), immediately reduces to the formula for the ball, when $\varepsilon_{1}=\cdots=\varepsilon_{q}=\varepsilon$. In particular, in the case of a ball we have the symmetric expression for the leading order term

[TABLE]

We remark that in this case the effect of the individual weight is not being seen by the relative entropy. 2. (2)

Formula (2.6) makes clear that the typical dependence of the leading order term of the relative entropy on the size of the vector of sufficient statistics $q$ is of the order one.

Proposition 2.2 quantifies the effect of choosing $\varepsilon_{i}$ differently for different directions $i$. Depending on the weighted effects of $\mathbb{E}_{f(\theta|\tau^{\ast})}\left(w^{2}_{i}(\tau^{\ast},\theta)\right)$ and of the correlations $\mathbb{E}_{f(\theta|\tau^{\ast})}\left(w_{i}(\tau^{\ast},\theta)w_{j}(\tau^{\ast},\theta)\right)$ one may be able to obtain information in regards to which of the directions $i=1,\cdots,q$ are more sensitive. We will see specific examples where this is true in Sections 6 and 7.

3. The case of exponential family

In this section we focus on the case where $f(x|\theta)$ belongs in the exponential family and our purpose is twofold. Exponential family is a large class of distributions and its specific structure will allow us to get more explicit formulas and also see the effect that the size of the data, $n$, has on the behavior of the relative entropy.

We compute the weights $w_{i}(\tau^{\ast},\theta)$ given by (2.7) for an exponential family. Due to the specific structural form that the involved densities have, we are able to get detailed formulas for $w_{i}(\tau^{\ast},\theta)$. In addition, the size of the data, $n$, appears in a clear way. In particular, in the case of an exponential family, the weights $w_{i}(\tau^{\ast},\theta)$ are of order $n^{2}$ which also implies that the expansion of relative entropy in $\varepsilon$ is typically correct if the product $n|\varepsilon|$ is small. Note that the dependence on $n$ is also related to the fact that the relative entropy analysis here compares the distribution of the model $P$ with respect to the model $P^{\varepsilon}$.

It is important to remark here though that in ABC the size of the data $n$ is fixed and one can only change $\varepsilon$. We show that for fixed $n$ and sufficiently small $|\varepsilon|$, the leading order term of relative entropy will typically scale as $n^{4}\sum_{i=1}^{q}\varepsilon^{4}_{i}/q$. However, as we shall demonstrate in a specific example in Section 6.2, this is only a crude generic upper bound and improvements are in general possible.

Since we assume that we have an exponential family we can write that

[TABLE]

Hence the joint distribution for $n$ i.i.d data $X_{1},\cdots,X_{n}$ is given by

[TABLE]

where we have denoted the sufficient statistics by $Y_{i}(x)=\sum_{j=1}^{n}S_{i}(x_{j})$. It is then well known that the distribution of the vector $(Y_{1}(x),\cdots,Y_{q}(x))$ belongs in the exponential family too and the corresponding joint probability density function is given by

[TABLE]

for an appropriate function $R(y_{1},\cdots,y_{q})$, whose form depends on the specifics of the distribution at hand. Instead of the sufficient statistics $(Y_{1}(x),\cdots,Y_{q}(x))$ we would be interested in the vector $(T_{1}(x),\cdots,T_{q}(x))$ where for each $i=1,\cdots,q$, $T_{i}(x)=\frac{Y_{i}(x)}{n}$. This is done in order to single out the effect of the size $n$ of the data set. A simple transformation then shows that the density of $(T_{1}(x),\cdots,T_{q}(x))$ is given by

[TABLE]

Now, as it is known, the natural conjugate prior for $\Theta$ takes the form

[TABLE]

where $Z(\rho,\kappa)$ is the appropriate normalizing factor. Then a relatively straightforward computation shows that the posterior distribution of $\Theta$ takes the form

[TABLE]

Essentially, due to conjugancy, we move from the prior to the posterior by setting $\rho\mapsto\rho+n\tau$ and $\kappa\mapsto\kappa+n$. The marginal likelihood of the data for $\tau=(\tau_{1},\cdots,\tau_{q})$ becomes

[TABLE]

Next goal is to make Proposition 2.2 precise for the case of exponential family. We start with a special case to provide some intuition. The general case is presented in Subsection 3.1.

Lemma 3.1**.**

Let us assume that we are in the special case where $q=1$, that $n\varepsilon$ is sufficiently small and let us also assume that the function $R(n\tau)=\text{constant}$. Assume that the acceptance region $D_{\varepsilon}(\tau^{\ast})$ is given by (2.4) (or equivalently by (2.5) since $q=1$). Then, we have that

[TABLE]

Proof.

Let us denote $\tau_{\varepsilon}=\tau\varepsilon$. Based on the expression (3.1) we can write

[TABLE]

Hence, we obtain that the relative entropy satisfies

[TABLE]

Under the assumptions of the lemma we then have

[TABLE]

Next, we notice that

[TABLE]

Hence, we can write

[TABLE]

completing the proof of the lemma. ∎

Lemma 3.1 makes it clear that, at least in the case of exponential family, the typical order of the relative entropy is $(n\varepsilon)^{4}$ and that for a given small tolerance, say tol, one can choose

[TABLE]

We also see that there is a weight factor that takes the form $\frac{1}{72}\text{Var}_{f(\Theta|\tau^{\ast})}\left(\eta^{2}(\theta)\right)$. In Section 3.1 we generalize the computation to the $q\neq 1$ case and for $R(\cdot)\neq\text{constant}$ based on the expression given by Proposition 2.2. Of course the expression derived in Section 3.1 coincides with (3.3) with $q=1$ and for $R(\tau)=\text{constant}$, but it also shows that the individual weights may have non-trivial contribution, a point made clearer in Section 6.

3.1. Calculation of the weight factors in the general case.

The goal of this section is to derive the analogous statement to Lemma 3.1 in the general case $q\neq 1$ and $R(n\tau)\neq\text{constant}$. We start with some assumptions. For a given function $h:\mathbb{R}^{q+1}\mapsto\mathbb{R}$, let $\partial_{i}h$ be the partial derivative in the $i^{\text{th}}$ direction.

Condition 3.2**.**

For each $k\in\mathbb{N}$ and for all $i=1,\cdots,q$ we have that $\partial^{(k)}_{i}\log R(n\tau^{\ast})$ is bounded uniformly on $n$ as $n\rightarrow\infty$.

Condition 3.2 is trivially true for the normal distribution for example. We also require

Condition 3.3**.**

Define the function $g(\theta)=\sum_{i=1}^{q}\eta_{i}(\theta)\tau_{i}-A(\theta)$. Let us assume that the function g is smooth and has a maximum only at one interior non-degenerate critical point, say $\theta^{*}$, i.e., that $\nabla g(\theta^{*})=0$ and $\text{det}\nabla^{2}g(\theta^{*})<0$.

Lemma 3.4**.**

Assume that the acceptance region is given by by (2.5). Under Conditions 2.1, 3.2 and 3.3 we have as $|n\varepsilon|\rightarrow 0$,

[TABLE]

The proof of this lemma is at the end of this section. A few remarks follow.

Lemma 3.4 shows that for appropriately small $n|\varepsilon|$ and as long as the rest of the involved quantities are bounded uniformly on $n$, the leading order term of the relative entropy scales like $n^{4}\sum_{i=1}^{q}\varepsilon_{i}^{4}/q$. However, this is only a crude upper bound and is not using the individual properties of the weight factors $w_{i}$. The weight factors and the probability distribution under which the expectation is being calculated also depends on $n$, which implies that it is possible that for some $i$, the behavior is better than $n^{4}$. In particular, in Section 6 we show that this is the case indeed for a few simple examples.

Let us also remark that if the function $R(\tau)$ is constant, (see Section 6.2 for such a case), then (3.4) takes a simpler form. Indeed, we have the following corollary.

Corollary 3.5**.**

Under the setting of Lemma 3.4 and assuming that $R(\tau)$ is constant, we obtain

[TABLE]

In particular if also $q=1$ then we obtain the statement of Lemma 3.1.

Proof of Lemma 3.4.

We start by computing the weight factors $w_{i}(\tau^{\ast},\theta)$ given by (2.7). By taking derivatives on (3.1) we obtain

[TABLE]

Hence, by plugging these expressions into (2.7) we subsequently obtain

[TABLE]

Next, we notice that (3.1) implies that the $k^{th}$ derivative of $Z^{-1}(\rho+n\tau^{\ast},\kappa+n)=Z^{-1}(\rho_{1}+n\tau^{\ast}_{1},\cdots,\rho_{q}+n\tau^{\ast}_{q},\kappa+n)$ in the $i^{th}$ direction is for $i=1,\cdots,q$

[TABLE]

Hence, we can actually write

[TABLE]

and

[TABLE]

The latter imply that the weight function takes the form

[TABLE]

Therefore, as $|n\varepsilon|\rightarrow 0$, one expects that the following holds up to leading order term for the relative entropy

[TABLE]

Let us conclude the proof of the lemma by justifying the validity of the expansion (3.10) when $|n\varepsilon|$ is small enough. A simple calculation of the Taylor series expansion of the density $f^{\varepsilon}(\theta|\tau^{*})$ shows that the remainder terms depend on terms of the form $\mathbb{E}_{f(\theta|\tau^{\ast})}\eta^{k}_{i}(\theta)$ for $k\in\mathbb{N}$ and on $\partial^{(k)}_{i}\log R(n\tau^{\ast})$. In regards to $\partial^{(k)}_{i}\log R(n\tau^{\ast})$ we have assumed Condition 3.2.

In regards to $\mathbb{E}_{f(\theta|\tau^{\ast})}\eta^{k}_{i}(\theta)$ we can be more explicit. We use classical Laplace asymptotics and Condition 3.3 to show that such quantities will be uniformly bounded in $n$ for any $k\in\mathbb{N}$. The latter and Condition 3.2 then imply that (3.4) is a valid expansion. Indeed, notice that we can write for $j=1,\cdots,q$ and $k\in\mathbb{N}$

[TABLE]

We have assumed that the function g is smooth and has a maximum only at one interior non-degenerate critical point, say $\theta^{*}$. That means $\nabla g(\theta^{*})=0$ and $\text{det}\nabla^{2}g(\theta^{*})<0$ (Condition 3.3). Then, classical Laplace asymptotics allows us to write as $n\rightarrow\infty$ up to leading order

[TABLE]

Since the rest of the terms in the expansion are of lower order with respect to $n$, we can conclude that (3.10) is a valid expansion of the relative entropy even when $n$ gets large as long as $|n\varepsilon|$ is small. If $n$ is held fixed, then (3.10) is a valid expansion of the relative entropy if $|\varepsilon|$ is sufficiently small. This concludes the proof of the lemma. ∎

4. Bias calculations

In this section we compute the bias of the estimator up to leading order with respect to $|\varepsilon|$. A similar computation has previously appeared in [3]. We present the result for the bias below, connecting it to the weight (2.7). In addition, the computations in [3] are for a ball as an acceptance region (2.4), whereas here we present the result for the more general case of an ellipse (2.5).

The bias of the estimator $\widehat{\pi_{X}[h]}_{\varepsilon,K}$ as defined by (2.1) is defined to be

[TABLE]

Algebraic computations similar to those in the case of a ball in [3] together with the computations of this section and the definition of the weights by (2.7), shows that

[TABLE]

where

[TABLE]

It is clear that the last display reduces to the form presented in [3] if $\varepsilon_{i}=\varepsilon$ for every $i=1,\cdots,q$ and after substitution of the form of the weight factors $w_{i}(\tau^{\ast},\theta)$ from (2.7).

Let us now make specific the bias computations of (4.1)-(4.2) for the case of an exponential family. By plugging in the expression for the weight $w_{i}(\tau^{*},\theta)$ from (3.9) in (4.2) we obtain

[TABLE]

Combined with (4.1), the latter expression shows that the typical order of the bias is $\mathcal{O}((n|\varepsilon|)^{2})$. However, as we will see in the specific examples of Section 6, by taking into account the weight functions, specific observables of interest do behave better in terms of the $n$-dependence for the bias, and the bias may even be independent of $n$.

5. Conditional mean rejection rate

The goal of this section is to compare the standard ABC algorithm that is using the ball as acceptance region (2.4) to the ABC algorithm that is using the ellipse as acceptance region (2.5). One possible way to do so, which can also be computed in practice, is the mean rejection rate per accepted particle for each one of the two algorithms. Let us define $R_{B}$ and $R_{E}$ to be the number of rejections before an accepted particle, conditional on the value of the observed statistic $\tau^{\ast}$, for ball and ellipse respectively. It is clear that both $R_{B}$ and $R_{E}$ are distributed according to geometric distribution

[TABLE]

where $p_{B}$ and $p_{E}$ being the acceptance probabilities for the case of acceptance region being a ball and an ellipse respectively. Writing $B_{\varepsilon}(\tau^{*})$ for the acceptance region for the ball, (2.4), and $D_{\varepsilon}(\tau^{*})$ for the acceptance region for the ellipse, (2.5), we have respectively

[TABLE]

where we recall that in the case of a ball $\varepsilon$ is positive scalar, whereas in the case of an ellipse $\varepsilon\in\mathbb{R}^{q}_{+}$.

Now, we are interested in approximating in the small $|\varepsilon|$ regime the quantity

[TABLE]

We write $U(\tau^{\ast})=\frac{\mathbb{E}\left[R_{E}|\tau^{\ast}\right]}{\mathbb{E}\left[R_{B}|\tau^{\ast}\right]}$ in order to emphasize that $U$ is a random quantity depending on the observed statistic, hence it is data driven itself.

We have the following lemma, with proof presented at the end of this section.

Lemma 5.1**.**

Asymptotically, as $|\varepsilon|\rightarrow 0$ we have that

[TABLE]

Let us make now Lemma 5.1 precise for the case of exponential family. As we show in Section 3, in the case of exponential family one sees clearly the dependence on the size of the data set, $n$, and as a matter of fact, in the generic situation, the relation (5.1) is then true if the products $n|\varepsilon|$ and $n\varepsilon$ are small.

Corollary 5.2**.**

In the case of exponential family, (5.1) reduces to

[TABLE]

where

[TABLE]

Interestingly (5.1) implies that if $\frac{\varepsilon_{i}}{\varepsilon}\gg 1$ for some directions $i=1,\cdots,q$ such that $\prod_{i=1}^{q}\frac{\varepsilon_{i}}{\varepsilon}>1$ then one can have benefits using the ellipse versus the ball since one can have that $U\ll 1$, which then implies that the mean rejection rate for the ellipse is smaller than the mean rejection rate for the ball, leading to potential computational gains. In Section 7 we demonstrate this in a specific simulation study.

Proof of Lemma 5.1.

We denote $\tau_{\varepsilon}=(\tau_{1}\varepsilon_{1},\cdots,\tau_{q}\varepsilon_{q})$. Let us start with $p_{E}$. Using Taylor series expansion up to second order, we have for $|\varepsilon|$ sufficiently small (with all other parameters such as $n$ being fixed)

[TABLE]

where in the last computation we use the symmetry of the domain and the relation

[TABLE]

It is clear that after setting $\varepsilon_{1}=\cdots=\varepsilon_{q}=\varepsilon$ we have

[TABLE]

Dividing the asymptotic expressions for $p_{E}$ and $p_{B}$ completes the proof of the lemma. ∎

Proof of Corollary 5.2.

Recalling now (3.6) and the identities (3.7)-(3.8) we may compute

[TABLE]

Recalling the definition of $q_{i}(n,\tau)$ from (5.2) we then have

[TABLE]

Plugging the expression for $\partial^{2}_{\tau_{i}}f_{T}(\tau)$ into (5.1) with $\tau=\tau^{\ast}$ and since by assumption $f_{T}(\tau^{\ast})\neq 0$ we conclude the proof of the Corollary. ∎

6. Examples

In this section, we present a few examples to illustrate the computations of this paper. In Section 6.1 we consider the case where we want to infer the rate of an exponential distribution. In Section 6.2 we consider the case where we want to infer the mean and variance of a normal distribution. In both cases we use the ABC framework and we demonstrate how the relative entropy computations can help in determining optimally values for the threshold parameter, based on relative entropy considerations. Of course, in these cases, we have explicit information for the likelihood and thus ABC methods are not necessary for inference. However, having access to the formulas allow us to make explicit and informative computations in regards to the behavior of the ABC algorithm.

6.1. Exponential distribution with unknown parameter

The goal of this section is to present a simple example where one can see how the computation of the weight factor can help in determining the tolerance parameter. In this example we will be inferring a scalar parameter, which implies that the $\varepsilon$ is scalar and the acceptance region is effectively a ball (2.4). We will see that even in this case, computation of the weight factor can be beneficial.

Let us assume that we have $n$ i.i.d $X_{1},\cdots,X_{n}\sim\text{Exp}(\theta)$ with unknown rate $\theta$. We then assume that the rate $\theta$ has gamma prior distribution, in particular we assume that $\theta\sim G(\alpha,\beta)$. Of course, the computations of Section 3 carry over here. One can pick as sufficient statistic for the parameter $\theta$ to be $T=\bar{X}$ where $\bar{X}=\frac{1}{n}\sum_{i=1}^{n}X_{i}$ is the empirical mean.

With this prior structure we have for the density posterior distribution

[TABLE]

which is the density of gamma $G(\alpha+n,\beta+n\tau)$ distribution. Let us set $\alpha_{n}=\alpha+n$ and $\beta_{n}=\beta+n\tau$. In the notation of Section 3 we have $\eta(\theta)=-\theta$ and $R(n\tau)=\frac{1}{\Gamma(n)}(n\tau)^{n-1}$.

Recall that the observed value of the statistic is $\tau^{*}=\bar{X}$. We compute

[TABLE]

[TABLE]

and

[TABLE]

Hence, we have the following approximation for the relative entropy

[TABLE]

It is easy to see now that if $n$ is large, then

[TABLE]

which then implies that for $n$ large such that $n\varepsilon$ is small enough

[TABLE]

implying that if there is a fixed tolerance level, say tol such that we require $H(P|P^{\varepsilon})(\tau^{\ast})\leq\text{tol}$, then one would need to choose $\varepsilon$ such that $n^{3}\varepsilon^{4}\frac{1}{18}\frac{1}{(\tau^{\ast})^{4}}\leq\text{tol}$. We remark here that the dependence on $n$ turned out to be better than $n^{4}$, which was due to the weight factor and that the value of the observed statistic affects the optimal choice of $\varepsilon$, in that the smaller it is, the smaller $\varepsilon$ should also become. The latter conclusion is also intuitive in the sense that the statistics $\tau^{\ast}=\bar{X}$ is the maximum likelihood estimator of $1/\theta$ and it measures the time that elapses on average between events. It is natural to expect that more frequent events would resort in smaller acceptance regions for ABC in order for the method to maintain its accuracy.

Let us conclude this section with the computation of the leading order of the bias. We choose as observable of interest the unknown rate $\theta$. We will show below that the leading order of the bias of this observable actually does not depend on $n$. For this purpose, we have for the specific example of exponential distribution

[TABLE]

Recalling the form of $\alpha_{n}$ and $\beta_{n}$ we then obtain that for large $n$,

[TABLE]

The latter implies that the bias does not depend on the size of the data set $n$ and it is of order $\varepsilon^{2}$. However, we also see that it is proportional to $\frac{1}{(\tau^{\ast})^{3}}$, which implies that if $\tau^{\ast}$ is small in value then the bias will be large for a fixed value of $\varepsilon$. The latter observation is in line with the qualitative dependence of the relative entropy with respect to the observed statistic $\tau^{*}$.

6.2. Normal distribution with unknown mean and unknown variance.

The goal of this section is to provide an explicit example where potential advantages appear when the acceptance region is an ellipse versus a ball, i.e., when one allows $\varepsilon_{i}$ to be different in each direction $i$.

We borrow the setup of the example from [21]. In particular, we consider a hierarchical binomial model, where we want to infer the probability of correct responses to a signal detection experiment for each subject in a given group, see [21]. Let $p_{i}$ be the probability of correct response of the i-th subject. Instead of modeling $p_{i}$ directly, we will model

[TABLE]

The logit function is useful as it transforms $p_{i}\in(0,1)$ to $\text{logit}(p_{i})\in(-\infty,\infty)$. In particular, set $X_{i}=\text{logit}(p_{i})$ and assume that we have $n$ i.i.d $X_{1},\cdots,X_{n}\sim N(\mu,\sigma^{2})$ with unknown mean $\mu$ and unknown variance $\sigma^{2}$.

In addition, we mention here that normal distribution has been used as a simulator tool in ABC not only due to its analytic structure but also due to its empirically observed robustness for irregular probability distributions in chaotic or nearly chaotic simulation dynamics, see [14, 24].

Let us denote by $\lambda=1/\sigma^{2}$ to be the precision parameter and assume the following normal-gamma prior for the pair $(\mu,\lambda)$

[TABLE]

Of course, normal distribution belongs to the exponential family, so the computations of Section 3 carry over here. One can pick as sufficient statistic for the parameter $\theta=(\mu,\lambda)$, $T=(\bar{X},S^{2})$ where $\bar{X}=\frac{1}{n}\sum_{i=1}^{n}X_{i}$ is the empirical mean and $S^{2}=\frac{1}{n-1}\sum_{i=1}^{n}(X_{i}-\bar{X})^{2}$ is the empirical variance. In this case we have $q=2$ for the range of the sufficient statistic $T$ and $m=2$ for the range of the parameter $\theta$. With this prior structure we have for the density posterior distribution

[TABLE]

where $Z(\alpha,\beta,\kappa)=\frac{\Gamma(\alpha)}{\beta^{\alpha}}\sqrt{\frac{2\pi}{\kappa}}$. Notice that if we define

[TABLE]

then we can re-express the formula for the posterior density as a product of a normal and gamma density as follows (with some abuse of notation)

[TABLE]

The density for the marginal likelihood takes the form

[TABLE]

Let us next rewrite the posterior density (6.2) in the canonical form. Writing $(n-1)S^{2}=n\frac{\sum_{i=1}^{n}X_{i}^{2}}{n}-n\bar{X}^{2}=n\overline{X^{2}}-n\bar{X}^{2}$, we see that when $\alpha=\frac{1+\kappa}{2}$ (6.2) takes the form (3.1) with the parameters

[TABLE]

and

[TABLE]

where now we consider (equivalently) as sufficient statistic to be the pair $(\bar{X},\overline{X^{2}})$.

Let us now compute the weight functions $w_{i}(\tau,\theta)$ given by (2.7) or equivalently $\text{Var}_{f(\Theta|\tau^{\ast})}\left(\eta^{2}_{i}(\theta)\right)$ and $\text{Cov}_{f(\Theta|\tau^{\ast})}\left(\eta^{2}_{1}(\theta),\eta^{2}_{2}(\theta)\right)$ as they appear in (3.5). This is the purpose of Corollary 6.1.

Corollary 6.1**.**

In the normal example studied in this section we have the following formulas

[TABLE]

[TABLE]

Lastly,

[TABLE]

Then, the following statement holds for the leading order term of the relative entropy.

Corollary 6.2**.**

We have that

[TABLE]

Proof of Corollaries 6.1 and 6.2 are deferred to the Appendix.

Now if the observed statistic $\tau^{\ast}=(\bar{X},\overline{X^{2}})$, or equivalently $\tau^{\ast}=(\bar{X},S^{2})$, is of order one, then it is easy to see that in terms of $n$ we have that $\beta_{n}=\mathcal{O}(nS^{2}/2)$, $\alpha_{n}=\mathcal{O}(n)$, $\kappa_{n}=\mathcal{O}(n)$ and $\mu_{n}=\mathcal{O}(\mu_{0}/n+\bar{X})$ as $n$ gets large. With some abuse of notation, and using Corollary 6.2, these suggest that for $n$ large and $\varepsilon$ sufficiently small such that $n|\varepsilon|$ is small we have up to leading order

[TABLE]

Expressions (6.7) and (6.8) make it clear that the weighted effect on each direction $i=1,2$ is different. It then follows that by exploiting the geometric structure of the relative entropy and the asymmetric effect of the sufficient statistic in different directions one may be able to afford to have larger acceptance region for the same overall tolerance level. Expression (6.8) also makes clear that the values of the observed sufficient statistics can have non-negligible consequences on the optimal choice of the threshold parameter. In Section 7, we confirm these points via a simulation study.

Remark 6.3*.*

We do point out that on this example we take advantage of the fact that, in a sense, mean and variance are independent directions. So allowing the possibility of choosing different $\varepsilon_{i}$ for each one of the directions $i=1,2$ is expected to be advantageous as opposed to forcing the choice $\varepsilon_{1}=\varepsilon_{2}=\varepsilon$. Notice though that their effect on the relative entropy is more complicated with the related cross-term being non-zero, see (6.6)-(6.8).

6.2.1. Mean rejection rate for the normal example

In this subsection we make precise the mean rejection rate computations of Section 5. In this case, we have that $q_{i}(n,\tau^{\ast})$ from (5.2) take the form

[TABLE]

and

[TABLE]

Thus, the leading order term of the ratio of the rejection rates in the limit as $n|\varepsilon|$ and $n\varepsilon$ go to zero, takes the form

[TABLE]

It is important to note here however that this is an approximation formula valid only in the limit and thus care is needed in its application and interpretation.

6.2.2. Bias computations for the normal example

Let us now discuss the behavior of the bias for the normal example. We will present the case of both mean and variance as observables. Let us start with the mean.

Let us set $h(\theta)=\mu$ and we want to compute the leading order of (4.1) for $h(\theta)=\mu$, i.e.,

[TABLE]

where in this case we have

[TABLE]

The result is summarized in the following Corollary.

Corollary 6.4**.**

For $i=1,2$, $C_{i}(\tau^{\ast})$ from (6.10) satisfy

[TABLE]

Therefore, we have

[TABLE]

Corollary 6.4 shows that if $\mu_{0}$ and $\bar{X}$ are small in value, then $C_{1}$, up to leading order, does not depend on $n$, but it grows inversely proportionally to $S^{2}$. The latter behavior is qualitatively the same as the dependence of the relative entropy with respect to $S^{2}$.

Let us then investigate the variance. Recalling that the precision $\lambda=1/\sigma^{2}$, let us now set $h(\theta)=\frac{1}{\lambda}$. We want to compute the leading order of (4.1) for $h(\theta)=\sigma^{2}=1/\lambda$, i.e.,

[TABLE]

where in this case we have

[TABLE]

The result is summarized in the following Corollary.

Corollary 6.5**.**

For $i=1,2$, $C_{i}(\tau^{\ast})$ from (6.12) satisfy

[TABLE]

Therefore, we have that for $n$ large enough and $\varepsilon$ small enough

[TABLE]

Proofs of Corollaries 6.4 and 6.5 are deferred to the appendix.

7. Simulation studies

In this section we present simulation studies to demonstrate some of the results of this paper. In particular, we show in the normal example of Section 6.2 that considering the ellipse instead of the circle as acceptance region, i.e., allowing $\varepsilon_{1}\neq\varepsilon_{2}$ has interesting advantages.

We present data for both the case of ball as acceptance region and for the case of an ellipse as acceptance region. In both cases, we fix a level of tolerance for the relative entropy, i.e., we require that $\varepsilon$ for the case of a ball or $\varepsilon_{1},\varepsilon_{2}$ for the case of an ellipse to be such that $H(P|P^{\varepsilon})(\tau^{\ast})\leq\text{tol}$. In the case of an ellipse, we choose $\varepsilon_{1}$ and $\varepsilon_{2}$ such that the area of the acceptance region, $\pi\varepsilon_{1}\varepsilon_{2}$, is maximized. In particular we solve the following optimization problem***In order to solve (7.1) we used R’s built-in function nloptr where the local and gradient-based optimization version of the Augmented Lagrangian algorithm was used.

[TABLE]

The measure of efficiency will be the average number of rejections per accepted particle, which is an indicative measure of computation time needed to complete the algorithm for a predefined level of tolerance error in terms of relative entropy. In addition, in order to showcase the effect of the size of the data, $n$, we record values for different values of $n$.

In Section 5 we derived asymptotic formulas for the mean rejection rates for both ellipse and ball. In the case of normal distribution, we derived that their ratio $U(\tau^{\ast})$ asymptotically, as $n|\varepsilon|$ and $n\varepsilon$ are small, satisfies (6.9).

We write $U(\tau^{\ast})$ for $U$ to emphasize that the ratio is a function of the observed statistic and thus a random quantity by itself. In addition, we remark here however that this approximation is accurate only in the limit as $n|\varepsilon|\downarrow 0$. To make clear that we are using the asymptotic formula we will write

[TABLE]

In all cases, the true values of the parameters under which the observed value of the statistic $\tau^{\ast}=(\bar{X},S^{2})$ is computed are $(\mu,\sigma^{2})=(0,1)$. The hyperparameters for the prior distributions are chosen to be $(\mu_{0},\kappa,\alpha,\beta)=(0,1,1,1)$. Also, in both cases the number of accepted sampled values, or particles, is fixed to $K=1000$. The goal is to make sure that the priors are appropriately chosen in order to concentrate the conclusions of the numerical studies on the effect of choosing a ball versus an ellipse as acceptance region.

Data for the ball as acceptance region, where $\varepsilon_{1}=\varepsilon_{2}=\varepsilon$ are presented in Table 1, whereas data for the ellipse as acceptance region, where $\varepsilon_{1}\neq\varepsilon_{2}$ are presented in Table 2. In the tables, we report the observed values of the sufficient statistic $\tau^{\ast}=(\bar{X},S^{2})$, the estimated values for $\mu$ and $\sigma^{2}$ along with their empirical standard deviations, as well as the average number of rejections per accepted particle, denoted by $\hat{R}_{B}$ and $\hat{R}_{E}$ for ball and ellipse respectively. Clearly, the more rejections per accepted particle the algorithm has, the more time it takes to run for the same level of accuracy.

Note that $\frac{\hat{R}_{E}}{\hat{R}_{B}}$ can be interpreted as the per unit time computational gain that the algorithm with the ellipse as acceptance region has compared to the algorithm that uses the ball as acceptance region (i.e. the choice of two different distance metrics). Clearly, the smaller $\frac{\hat{R}_{E}}{\hat{R}_{B}}$ than $1$ is, the better the algorithm that uses the ellipse versus the algorithm that uses the ball for acceptance region is.

In Table 3 we see a comparison of the asymptotic $\tilde{U}$ and empirical $\frac{\hat{R}_{E}}{\hat{R}_{B}}$ for the data presented in Tables 1 and 2. It is clear that the asymptotic formula $\tilde{U}$ overestimates the ratio of mean rejection rates due to the fact that it is only an asymptotically correct formula. However, it is evident from the simulation results that $\tilde{U}$ is of the correct order. In addition, for better presentation, we have plotted in Figure 2, the per unit time computational gain of the algorithm, $\frac{\hat{R}_{E}}{\hat{R}_{B}}$, for different values of tolerance levels taken from Table 3.

The conclusion is that for the same level of tolerance, allowing $\varepsilon_{1}\neq\varepsilon_{2}$ has computational advantages as it reduces the number of rejections per accepted particle. In addition, the tables also demonstrate that as $n$ increases, $\varepsilon_{i}$ decrease. This is of course to be expected due to the form that the weight functions take in (6.7). It is also noteworthy that the advantage of allowing $\varepsilon_{1}\neq\varepsilon_{2}$ seems to become more evident as $n$ increases even though the degree of the potential computational gain seems to diminish the larger the $n$ gets, see Figure 2 (at least for the simulation that is studied here).

8. Conclusions

In this paper we analyzed the popular ABC algorithm via the lens of relative entropy ideas. By computing the leading order term in the expansion of relative entropy as function of the threshold vector parameter, we can quantify the effect of the chosen distance metric. In this direction, we showed that advantages arise when one exploits potential asymmetries in the distribution of observed statistics. In particular, one can then allow for a larger acceptance region which then implies a smaller rejection rate, which in turn results in faster convergence for the same level of tolerance error. This is part of a larger question, i.e., the quantification of the effect of distance metrics on the behavior of ABC. In this paper, we have seen that analyzing the ABC through the lens of relative entropy gives us a way to do so.

In addition, we characterized precisely the effect of the number of data points on the performance of ABC by showing that in order to maintain the same level of tolerance, one should decrease the threshold parameter as the number of data points increases. In connection to that, Figure 2 seems to suggest that as the number of data points, $n$, increase, the gain in performance due to the choice of the distance metric seems to stabilize. It would be of interest to theoretically address this aspect of the algorithm and as a consequence quantify in more precise terms the effect that the choice of the distance metric has as the number of data points increase.

One important conclusion of our analysis is that the leading order term in the expansion of relative entropy depends on certain weight factors that then govern the appropriate choice of threshold parameter. In this paper, our goal is to illustrate the phenomenon and for the numerical examples that we presented we could compute the weight functions in closed form. This was done in order to focus on the more theoretical aspects of the algorithm. For more general models, closed form computation may not be possible and one would have to result to simulation. The weight functions can be potentially monitored online and computed on the fly using Monte Carlo methods. We leave the development of the latter computational methods for future work.

Appendix A Proof of intermediate results

Lemma A.1**.**

Assume Condition 2.1 and that the acceptance region $D_{\varepsilon}(\tau^{\ast})$ is given by (2.5). Let $f^{\varepsilon}(\theta|\tau^{\ast})$ be defined by (2.2). Then we have that $\nabla_{\varepsilon}f^{\varepsilon}(\theta|\tau^{\ast})|_{\varepsilon=0}=0$. In addition, we have that for all $i=1,\cdots,q$

[TABLE]

and for all $i,j=1,\cdots,q$ with $i\neq j$, $\partial^{2}_{\varepsilon_{i}\varepsilon_{j}}f^{\varepsilon}(\theta|\tau^{\ast})|_{\varepsilon=0}=0$.

Proof of Lemma A.1.

Let us denote $\tau_{\varepsilon}=(\tau_{1}\varepsilon_{1},\cdots,\tau_{q}\varepsilon_{q})$. Notice that by changing variables we can write

[TABLE]

Let $i=1,\cdots,q$. Taking now derivative with respect to $\varepsilon_{i}$, we obtain

[TABLE]

Due to the symmetry of the domain of integration around zero and because $s\mapsto s$ is an odd function, we obtain when we evaluate at $\varepsilon=0$ for all $i=1,\cdots,q$

[TABLE]

Hence in the case of the ellipse, i.e. (2.5), we actually obtain that the order $O(|\varepsilon|^{2})$ and $O(|\varepsilon|^{3})$ vanish. We subsequently obtain for the relative entropy

[TABLE]

Let us next compute the second derivative of the approximate posterior with respect to $\varepsilon$. We have for the diagonal terms of the Hessian $\nabla^{2}_{\varepsilon}f^{\varepsilon}(\theta|\tau^{\ast})$

[TABLE]

Evaluated at $\varepsilon=0$ we have

[TABLE]

Next using the symmetry of the domain of integration we compute

[TABLE]

Thus, we then have

[TABLE]

To complete the computations we also need the mixed-derivatives of the Hessian $\nabla^{2}_{\varepsilon}f^{\varepsilon}(\theta|\tau^{\ast})$. Hence, for $i,j=1,\cdots,q$ with $i\neq j$ we have

[TABLE]

Evaluated at $\varepsilon=0$ and using the symmetry of the domain and the fact that we are integrating odd polynomials we have

[TABLE]

Then, by (A.3) and (A.4)-(A.5) we conclude the proof of the lemma. ∎

Proof of Corollary 6.1.

We compute

[TABLE]

Similarly, we have

[TABLE]

In addition,

[TABLE]

concluding the proof of the Corollary. ∎

Proof of Corollary 6.2.

Using Lemma 3.4 and Corollary 6.1 we have that

[TABLE]

completing the proof of the Corollary. ∎

Proof of Corollary 6.4.

Recalling the form of $\eta_{1}(\theta)=\mu\lambda$ we have

[TABLE]

Similar computations and recalling that $\eta_{2}(\theta)=-\frac{\lambda}{2}$ give

[TABLE]

∎

Proof of Corollary 6.5.

Recalling the form of $\eta_{1}(\theta)=\mu\lambda$ we have

[TABLE]

Similarly, we have, for $n$ large

[TABLE]

Therefore, we get that the bias of $\sigma^{2}=\frac{1}{\lambda}$ behaves for $n$ large enough and $\varepsilon$ small enough as follows to leading order

[TABLE]

∎