# Approximate Bayesian Computation via Classification

Yuexi Wang, Tetsuya Kaji and Veronika Ročková
Booth School of Business, University of Chicago The author gratefully
acknowledges the support from the James S. Kemper Faculty Fund at the Booth
School of Business and the National Science Foundation (Grant No. NSF
DMS-1944740)

###### Abstract

Approximate Bayesian Computation (ABC) enables statistical inference in
complex models whose likelihoods are difficult to calculate but easy to
simulate from. ABC constructs a kernel-type approximation to the posterior
distribution through an accept/reject mechanism which compares summary
statistics of real and simulated data. To obviate the need for summary
statistics, we directly compare empirical distributions with a Kullback-
Leibler (KL) divergence estimator obtained via classification. In particular,
we blend flexible machine learning classifiers within ABC to automate
fake/real data comparisons. We consider the traditional accept/reject kernel
as well as an exponential weighting scheme which does not require the ABC
acceptance threshold. Our theoretical results show that the rate at which our
ABC posterior distributions concentrate around the true parameter depends on
the estimation error of the classifier. We derive limiting posterior shape
results and find that, with a properly scaled exponential kernel, asymptotic
normality holds. We demonstrate the usefulness of our approach on simulated
examples as well as real data in the context of stock volatility estimation.

Keywords: Approximate Bayesian Computation; Classification; Likelihood-free
Inference;
Kullback-Leibler Divergence; Posterior Concentration

##  1 Introduction

Implicit in the deployment of many Bayesian methods is the assumption that the
likelihood function can be numerically evaluated. More often than not,
however, modern statistical models are prescribed by intractable likelihoods
or are defined through an underlying generating mechanism. Fortunately, it may
still be possible to simulate synthetic datasets from the model. The ability
to simulate from the likelihood has opened up new opportunities for simulating
from the posterior. For example, Approximate Bayesian Computation (ABC)
(Pritchard et al., 1999; Beaumont et al., 2002) emerged as a default
likelihood-free Bayesian inferential tool. It is an accept/reject posterior
sampling mechanism which obviates likelihood evaluations. Each iteration
proceeds by (1) simulating prior parameter guesses and fake data from the
likelihood, and then (2) accepting those parameter values whose fake data were
close to the observed data. A big challenge with ABC has been gauging the
similitude between observed and fake data.

Measures of similarity between data sets have traditionally been based on
summary statistics (see Blum et al. (2013) for an overview within the ABC
context). In other words, two datasets are considered similar if their summary
statistics are close. In the absence of expert knowledge, however,
constructing effective summary statistics can be challenging (Joyce and
Marjoram, 2008; Nunes and Balding, 2010; Blum et al., 2013) and one may need
to resort to automated strategies. One possibility is regressing parameter
values onto (functionals of) fake data in a pilot ABC run to train a flexible
mapping which can be substituted for summary statistics (Fearnhead and
Prangle, 2012; Jiang et al., 2017; Åkesson et al., 2020). Another possibility,
related to indirect inference, is to construct summary statistics from an
auxiliary model (Drovandi et al., 2011; Wood, 2010). One can also choose a
subset of candidate summary statistics that satisfy some optimality criterion
(Joyce and Marjoram, 2008; Nunes and Balding, 2010) or find an optimal
projection of a set of summary statistics onto a lower-dimensional map
(Fearnhead and Prangle, 2012). Alternatively, one can view the observed and
synthetic data sets as empirical distributions and directly use a distance
between them inside ABC (such as Kullback-Leibler (Jiang et al., 2018) or
Wasserstein (Bernton et al., 2019) or Maximum Mean (MM) discrepancy (Park et
al., 2016)). Our work fortifies this last point of view by focusing on the
Kullback-Leibler discrepancy estimated via classification.

The KL divergence is one of the most widely used discrepancy metrics. It
expresses the average information per observation to discriminate between two
probabilistic models (Kullback, 1958). In large deviations, for example, it
characterizes the exponential decay rate at which empirical measures converge
to their probabilities (see Sanov’s theorem in Den Hollander (2008)) and the
rate of decay of the probability of error in a binary hypothesis testing
problem (see Stein’s Lemma in Cover and Thomas (1991)). KL also naturally
connects to maximum likelihood estimation through its interpretation as the
expectation of the log-likelihood ratio. There exist many methods for
estimating the KL divergence. For example, Wang et al. (2009) proposed
nearest-neighbor techniques to obtain a mean-square consistent estimator. Wang
et al. (2005) proposed a histogram-based KL estimator based on partitioning of
the space into statistically equivalent intervals. Silva and Narayanan (2007)
and Silva and Narayanan (2010) went a step further and proposed using data-
driven partitions (including multivariate recursive partitioning) and
formulated sufficient consistency conditions. Alternatively, Nguyen et al.
(2007) proposed a variational approach by turning KL estimation into a
penalized convex risk minimization problem. Our work is different from the
approaches above as we propose a KL estimator based on classification.

We suggest embedding a machine learning classifier inside ABC to determine
whether or not fake and observed data are similar and, thereby, whether or not
the underlying parameter value should be kept in the ABC reference table. The
fundamental premise of this proposal is as follows: parameter values that
yield indistinguishable simulated datasets can be deemed close. Bayesian
inference via classification has been suggested before. Kaji and Rockova
(2021) developed a version of the Metropolis-Hastings algorithm, called MHC,
based on classification-based estimators of likelihood ratios. Gutmann et al.
(2018) proposed a classification strategy related to ours using a different
discrepancy metric. Our paper reframes the method of Gutmann et al. (2018) as
a genuine ABC algorithm with a KL divergence discriminator and provides
supporting theory which justifies its inferential potential.

In particular, we study statistical properties of the approximate posterior
which, in part, depend on the properties of the KL divergence estimator. We
consider the traditional accept/reject ABC version (with a uniform kernel) as
well as an exponential kernel variant which does not require the ABC tolerance
threshold. Similar to Frazier et al. (2018), we show that the choice of the
ABC acceptance threshold $\epsilon$ plays a critical role in the convergence
rate and in the limiting posterior shape. In practice, it is often not obvious
what the optimal threshold $\epsilon$ should be. Motivated by the connections
with the MHC algorithm of Kaji and Rockova (2021), we propose an exponential
kernel which yields ABC posteriors that correspond to the stationary
distribution of MHC. Our ABC kernel method can be thus regarded as a
parallelizable counterpart to the sequential MHC sampling, targeting the same
posterior approximation. The concentration and asymptotic shape behavior of
the ABC posterior, which can be derived from Kaji and Rockova (2021),
theoretically justify our exponential weighting scheme. In addition, the
exponential kernel can be motivated as an instantiation of General Bayesian
Inference (GBI) (Bissiri et al., 2016) with the KL estimator as the loss
function.

The remaining of the paper is structured as follows. In Section 2, we flesh
out the basic idea of ABC and introduce our framework with classification. In
Section 3, we investigate the posterior concentration and limiting shape
behaviors of the ABC posteriors. Section 4 shows performance on simulated
datasets and Section 5 further highlights the practical value of our approach
on real data. In Section 6, we conclude with a discussion.

Notation. We employ the operator notation for expectation, e.g., $P_{0}f=\int
fdP_{0}$. The $\epsilon$-bracketing number $N_{[]}(\epsilon,\mathcal{F},d)$ of
a set $\mathcal{F}$ with respect to a premetric $d$ is the minimal number of
$\epsilon$-brackets in $d$ needed to cover $\mathcal{F}$111A premetric on
$\mathcal{F}$ is a function $d:\mathcal{F}\times\mathcal{F}\to\mathbb{R}$ such
that $d(f,f)=0$ and $d(f,g)=d(g,f)\geq 0$.. The $\delta$-bracketing entropy
integral of $\mathcal{F}$ with respect to $d$ is
$J_{[]}(\delta,\mathcal{F},d)=\int_{0}^{\delta}\sqrt{1+\log
N_{[]}(\epsilon,\mathcal{F},d)}d\epsilon$. Next, $K(f,g)=\int f\log(f/g)d\mu$
denotes the Kullback-Leibler divergence between two density functions and
$V_{2}(f,g)=\int f\left|\log(f/g)\right|^{2}d\mu$. For real-valued sequences
$\\{a_{n}\\}_{n\geq 1}$ and $\\{b_{n}\\}_{n\geq 1}$, $a_{n}\lesssim b_{n}$
means that $a_{n}\leq C\,b_{n}$ for some generic constant $C>0$, $a_{n}\asymp
b_{n}$ means that $a_{n}\lesssim b_{n}\lesssim a_{n}$, and $a_{n}\gg b_{n}$
indicates a greater order of magnitude. For a sequence of random variables
$x_{n}$, $x_{n}=o_{P}(a_{n})$ if
$\lim_{n\to\infty}P(\left|x_{n}/a_{n}\right|\geq C)=0$ for every $C>0$, and
$x_{n}=O_{P}(a_{n})$ if for every $C>0$ there exists a finite $M>0$ and a
finite $N$ such that $P(\left|x_{n}/a_{n}\right|\geq M)\leq C$ for all $n>N$.
All limits are taken as $n\to\infty$. Take $\left\|\cdot\right\|$ to be the
Euclidean norm.

##  2 ABC without Summary Statistics

Our framework consists of a collection of i.i.d. observations
$\bm{X}=(X_{1},\ldots,X_{n})^{\prime}$ where each $X_{i}\in\mathcal{X}$ is
realized from a parametric model
$\\{P_{\theta_{0}}:\theta_{0}\in\Theta\subset\mathbb{R}^{d}\\}$. We assume
that $P_{\theta}$, for each $\theta\in\Theta$, admits a density $p_{\theta}$.
Hereafter, we will use the shorthand notation $p_{0}=p_{\theta_{0}}$ and
$P_{0}=P_{\theta_{0}}$. We are interested in Bayesian inference about
$\theta_{0}$ based on the posterior distribution

$\pi_{n}(\theta\mid\bm{X})\propto p_{\theta}^{(n)}(\bm{X})\pi(\theta)$ (1)

prescribed by the likelihood $p_{\theta}^{(n)}(\bm{X})$ and the prior density
$\pi(\theta)$. We are particularly interested in scenarios when the likelihood
function cannot be directly expressed/evaluated (such as discretely observed
diffusions (Sørensen, 2004) or generative models) but can be sampled from.

The now default ABC method for Bayesian likelihood-free inference constructs a
nested kernel-type approximation to the posterior distribution. The first
approximation occurs when the data is distilled into summary statistics to
obtain $\pi(\theta\,|\>S_{X})\propto\pi(S_{X}\,|\>\theta)\pi(\theta)$, where
$S_{X}=S(\bm{X})$ is a vector of summary statistics. The quality of this
approximation depends crucially on the informativeness of $S_{X}$. The actual
ABC approximation to the posterior (1) is then constructed via a kernel
function as $\pi_{ABC}(\theta\,|\>S_{X})=\int\pi(\theta,S\,|\>S_{X})dS$ with
$\pi(\theta,S)\propto T_{\epsilon}(\|S-S_{Y}\|)\pi(S\,|\>\theta)\pi(\theta)$,
where $T_{\epsilon}(\|u\|)=T(\|u\|/\epsilon)$ is a standard (smoothing) kernel
with a scale parameter $\epsilon>0$. The key two challenges with ABC are (1)
deriving low-dimensional summary statistics with a minimal loss of information
and (2) selecting the kernel and its the tolerance level $\epsilon$. To
remediate the reliance of ABC on summary statistics, we focus on viewing the
observed and fake data as empirical distributions and gauge the distance
between them. Regarding the choice of aggregation kernels, we consider the
traditional uniform kernel yielding an accept/reject algorithm and a smoothing
kernel free from $\epsilon$-tuning.

###  2.1 ABC with KL Divergence

Instead of summary statistics, we use the estimated KL divergence inside the
ABC algorithm. Our interest in the KL divergence as a discrepancy metric stems
partially from the following connection to the generalized Bayesian inference
(Bissiri et al., 2016). The posterior distribution (1) can be rewritten as a
generalized posterior
$\pi_{n}(\theta\mid\bm{X})\propto\pi(\theta)\exp(-n\times
KL(p_{0}^{(n)},p_{\theta}^{(n)}))$ where the parameter $\theta$ is linked to
data through the empirical Kullback-Liebler (KL) divergence
$KL(p_{0}^{(n)},p_{\theta}^{(n)})\equiv\frac{1}{n}\sum_{i=1}^{n}\log({p_{0}}/{p_{\theta}})(X_{i})$.
For when the KL divergence cannot be easily evaluated, we consider various
estimators in the next section. We denote a generic KL divergence estimator
obtained from observed data $\bm{X}\sim P_{0}^{(n)}$ and pseudo-data
$\tilde{\bm{X}}^{\theta}\sim P_{\theta}^{(n)}$ as
$\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})$. Adopting $\hat{K}(\cdot,\cdot)$ as
the ABC summary statistics, we consider a simple accept/reject ABC mechanism
detailed in Algorithm 1 below. While Jiang et al. (2017) used a nearest-
neighbor estimator of the KL divergence, we devise a different estimator based
on classification in Section 2.2.

For a pre-determined tolerance level $\epsilon>0$ repeat for $j=1,\dots,N$:

  1. 1.

Simulate $\theta_{j}$ from $\pi(\theta)$.

  2. 2.

Simulate
$\tilde{\bm{X}}^{\theta_{j}}=(\tilde{X}_{1}^{\theta_{j}},\ldots,\tilde{X}_{m}^{\theta_{j}})^{\prime}$
through i.i.d sampling from the model $p_{\theta_{j}}$.

  3. 3.

Accept $\theta_{j}$ when
$\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta_{j}})\leq\epsilon$.

Algorithm 1 KL-ABC with Accept-Reject

Algorithm 1 simulates pairs of parameter values and pseudo-data
$\\{\theta,\tilde{\bm{X}}^{\theta}\\}$ from the joint posterior density

$\hat{\pi}^{AR}\left(\theta,\tilde{\bm{X}}^{\theta}\mid\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon\right)=\frac{\pi(\theta)p_{\theta}^{(n)}(\tilde{\bm{X}}^{\theta})\mathbb{I}\big{(}\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon\big{)}}{\int\pi(\theta)p_{\theta}^{(n)}(\tilde{\bm{X}}^{\theta})\mathbb{I}\big{(}\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon\big{)}\mathrm{d}\tilde{\bm{X}}^{\theta}\mathrm{d}\theta},$
(2)

which margins towards the following approximate (Accept/Reject) posterior
density

$\displaystyle\hat{\pi}^{AR}_{\epsilon}(\theta\mid\bm{X})$
$\displaystyle=\int\hat{\pi}^{AR}\left(\theta,\tilde{\bm{X}}^{\theta}\mid\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon\right)\mathrm{d}\tilde{\bm{X}}^{\theta}$
$\displaystyle\equiv\frac{\pi(\theta)P_{\theta}^{(n)}\big{(}\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon\big{)}}{\int\pi(\theta)P_{\theta}^{(n)}\big{(}\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon\big{)}\mathrm{d}\theta}.$
(3)

The approximation (3) will be scrutinized theoretically later in Section 3.1.
In particular, we will later see that the convergence rate of (3) around
$\theta_{0}$ depends on the choice of $\epsilon$ (Frazier et al., 2018).
Algorithm 1 uses the uniform kernel which corresponds to the indicator
function $T_{\epsilon}(\|u\|)=\mathbb{I}(\|u\|\leq\epsilon)$. In practice, it
is difficult to balance out conflicting demands of smaller $\epsilon$
(yielding good approximability) and larger acceptance rates (yielding more
posterior samples). As a remedy, we propose a way to aggregate the ABC samples
through a scaled exponential kernel motivated by the connection between KL and
the log-likelihood ratio. This ABC variant requires no ad-hoc thresholding and
is summarized in Algorithm 2.

Repeat for $j=1,\dots,N$:

  1. 1.

Simulate $\theta_{j}$ from $\pi(\theta)$.

  2. 2.

Simulate
$\tilde{\bm{X}}^{\theta_{j}}=(\tilde{X}_{1}^{\theta_{j}},\ldots,\tilde{X}_{m}^{\theta_{j}})^{\prime}$
through i.i.d sampling from the model $p_{\theta_{j}}$.

  3. 3.

Assign $\theta_{j}$ a weight proportional to
$\exp\big{(}-n\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\big{)}$.

Algorithm 2 KL-ABC with Exponential Weighting

Algorithm 2 generates draws for the pair
$\\{\theta,\tilde{\bm{X}}^{\theta}\\}$ from a joint posterior density

$\hat{\pi}^{EK}(\theta,\tilde{\bm{X}}^{\theta}\mid\bm{X})=\frac{\pi(\theta)p_{\theta}^{(n)}(\tilde{\bm{X}}^{\theta})\exp\big{(}-n\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\big{)}}{\int\pi(\theta)p_{\theta}^{(n)}(\tilde{\bm{X}}^{\theta})\exp\big{(}-n\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\big{)}\mathrm{d}\tilde{\bm{X}}^{\theta}\mathrm{d}\theta},$
(4)

which leads to the approximated Bayesian posterior as

$\displaystyle\hat{\pi}^{EK}(\theta\mid\bm{X})$
$\displaystyle=\int\hat{\pi}^{EK}(\theta,\tilde{\bm{X}}^{\theta}\mid\bm{X})\mathrm{d}\tilde{\bm{X}}^{\theta}=\frac{\pi(\theta)P_{\theta}^{(n)}\exp\big{(}-n\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\big{)}}{\int\pi(\theta)P_{\theta}^{(n)}\exp\big{(}-n\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\big{)}\mathrm{d}\theta}.$
(5)

Smooth kernels have been used inside ABC before to rescale the acceptance
probability (e.g. Beaumont et al. (2002) employed the Epanechnikov kernel).
Wilkinson (2013) and Sisson et al. (2018) provide a thorough overview and
comparisons of the commonly used kernels. Our smoothed weights are directly
interpretable due to their linkage between the KL divergence and the log-
likelihood ratio. We later show in Section 3.1 that, with the scaled
exponential kernel, the ABC posterior corresponds to the stationary
distribution of the MHC algorithm of Kaji and Rockova (2021) and can be
regarded as a posterior under a misspecified model.

###### Remark 2.1

(Generating Fake Data) We assume
$\tilde{\bm{X}}^{\theta}=g_{\theta}(\tilde{\bm{X}})$, where
$\tilde{\bm{X}}\in\mathbb{R}^{m}$ are random variables arriving from
$\tilde{P}^{(m)}$ and where
$g_{\theta}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ is a deterministic
mapping. Generating random variable draws by passing $\tilde{\bm{X}}$ through
some mapping is commonly done in practice. For example, Gaussian random
variables $\tilde{\bm{X}}=\\{\tilde{X}^{\theta}_{i}\\}_{i=1}^{m}$ that follow
i.i.d. $N(\mu,\sigma^{2})$ distribution can be obtained by transforming
$\\{\tilde{X}_{i}\\}_{i=1}^{m}\overset{i.i.d.}{\sim}N(0,1)$ via
$\tilde{X}^{\theta}_{i}=\mu+\sigma\tilde{X}_{i}$. In other cases, one can use
uniform draws $\tilde{\bm{X}}$ and the inverse transform sampling.

###  2.2 Estimating KL Divergence via Classification

We adopt the ‘$-\log D$’ trick to estimate the KL divergence (Goodfellow et
al., 2014). More precisely, a flexible discriminator $D$ (such as a neural
network or logistic regression) is trained to maximize

$\mathbb{M}_{n,m}^{\theta}(D)=\mathbb{P}_{n}\log
D+\mathbb{P}_{m}^{\theta}\log(1-D),$ (6)

where we employ the operator notation for expectation, e.g.,
$\mathbb{P}_{n}f=\frac{1}{n}\sum_{i=1}^{n}f(X_{i})$ and
$\mathbb{P}_{m}^{\theta}f=\frac{1}{m}\sum_{i=1}^{m}f(\tilde{X}_{i}^{\theta})$.
This can also be regarded as a classification problem where we label
$\\{X_{i}\\}_{i=1}^{n}$ (‘real’ data) with 1 and
$\\{\tilde{X}_{i}^{\theta}\\}_{i=1}^{m}$ (‘fake’ data) with 0. The oracle
maximizer to (6) can be shown to be (Goodfellow et al., 2014, Proposition 1)

$D_{\theta}(X)=\frac{p_{0}(X)}{p_{0}(X)+p_{\theta}(X)}.$ (7)

The functional form of the oracle solution in (7) naturally suggests the
following KL estimator obtained from a trained discriminator
$\hat{D}^{\theta}_{n,m}$

$\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})=\mathbb{P}_{n}\log\frac{\hat{D}^{\theta}_{n,m}}{1-\hat{D}^{\theta}_{n,m}}=\frac{1}{n}\sum_{i=1}^{n}\log\frac{\hat{D}^{\theta}_{n,m}(X_{i})}{1-\hat{D}^{\theta}_{n,m}(X_{i})}.$
(8)

Next, we show that our classification-based KL estimator (8) converges to a
well-defined limit counterpart $K(p_{0},p_{\theta})$ under mild conditions. We
assume that the set of considered classifiers $\mathcal{D}$ resides in a sieve
$\mathcal{D}_{n}$ that expands with the sample size and its size is measured
by bracketing entropy $N_{[]}(\epsilon,\mathcal{F},d)$. Kaji et al. (2020)
prove the rate of convergence of such a classifier (under assumptions reviewed
later) with a Hellinger-type distance defined as

$d_{\theta}(D_{1},D_{2})=\sqrt{h_{\theta}(D_{1},D_{2})^{2}+h_{\theta}(1-D,1-D_{\theta})^{2}},$

where
$h_{\theta}(D_{1},D_{2})=\sqrt{(P_{0}+P_{\theta})(\sqrt{D}_{1}-\sqrt{D}_{2})^{2}}$.

###### Assumption 1

(Kaji et al., 2020, Assumption 3) Assume that $n/m$ converges and that an
estimator $\hat{D}^{\theta}_{n,m}$ exists that satisfies
$\mathbb{M}_{n,m}(\hat{D}^{\theta}_{n,m})\geq\mathbb{M}_{n,m}^{\theta}(D_{\theta})-O_{P}(\delta_{n}^{2})$
for a nonnegative sequence $\delta_{n}$. Moreover, assume that the bracketing
entropy integral satisfies
$J_{[]}(\delta_{n},\mathcal{D}^{\theta}_{n,\delta_{n}},d_{\theta})\lesssim\delta_{n}^{2}\sqrt{n}$
and that there exists $\alpha<2$ such that
$J_{[]}(\delta_{n},\mathcal{D}^{\theta}_{n,\delta_{n}},d_{\theta})/\delta^{\alpha}$
is increasing in $\delta$. Here
$\mathcal{D}^{\theta}_{n,\delta_{n}}=\\{D\in\mathcal{D}_{n}:d_{\theta}(D,D_{\theta})\leq\delta\\}$.

###### Lemma 2.1

(Kaji et al., 2020, Theorem 1’) Under Assumption 1,
$d_{\theta}(\hat{D}_{n,m}^{\theta},D_{\theta})=O_{P}^{*}(\delta_{n})$. 222We
use $P^{*}$ to denote outer expectation (definition see Section 1.2 of van der
Vaart and Wellner (1996)), here is the expectation of a “a smallest
measureable function g that dominates
$d_{\theta}(\hat{D}_{n,m}^{\theta},D_{\theta})$”.

In Assumption 1, we only assume that $\delta_{n}$ be a nonnegative sequence,
not necessarily vanishing. To ensure that the approximate posterior converges,
we require $\delta_{n}\to 0$ as $n\to\infty$. Kaji et al. (2020) show that if
the true likelihood ratio has a low-dimensional representation and an
appropriate neural network is used for the discriminator, the rate
$\delta_{n}$ depends only on the underlying dimension and it vanishes when
$n\to\infty$. This result can also be generalized to other network
configurations, including Schmidt-Hieber (2020) and Yarotsky (2017). Hence, we
implicitly assume $\delta_{n}=o(1)$ and note that $\delta_{n}$ will typically
be slower than $n^{-1/2}$.

To establish the rate of convergence of our approximated posterior, we have
the following assumption.

###### Assumption 2

(Kaji and Rockova, 2021, Assumption 2) There exists $\Lambda>0$ such that for
every $\theta\in\Theta$, $P_{0}(p_{0}/p_{\theta})$ and
$P_{0}(p_{0}/p_{\theta})^{2}$ are bounded by $\Lambda$ and

$\sup_{D\in\mathcal{D}^{\theta}_{n,\delta_{n}}}P_{0}\left(\frac{D_{\theta}}{D}\middle|\frac{D_{\theta}}{D}\geq\frac{25}{16}\right)<\Lambda,\quad\sup_{D\in\mathcal{D}^{\theta}_{n,\delta_{n}}}P_{0}\left(\frac{1-D_{\theta}}{1-D}\middle|\frac{1-D_{\theta}}{1-D}\geq\frac{25}{16}\right)<\Lambda,$

for $\delta_{n}$ in Assumption 1. The brackets in Assumption 1 can be taken so
that $P_{0}(\sqrt{u/l}-1)^{2}=O(d_{\theta}(u,l)^{2})$ and
$P_{0}(\sqrt{(1-l)/(1-u)}-1)^{2}=o(d_{\theta}(u,l))$.

The following theorem quantifies the rate of convergence of our estimator (8)
towards the empirical KL divergence
$\mathbb{P}_{n}\log\frac{p_{0}}{p_{\theta}}$.

###### Lemma 2.2

(Convergence Rate of Estimation Errors) Under Assumptions 1 and 2,

$\displaystyle\left|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-\mathbb{P}_{n}\log\frac{p_{0}}{p_{\theta}}\right|=O_{P^{*}}(\delta_{n}).$
(9)

Proof.  Since the estimation error can be rewritten as

$\displaystyle\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-\mathbb{P}_{n}\log\frac{p_{0}}{p_{\theta}}=-\mathbb{P}_{n}\Big{(}\log\frac{1-\hat{D}^{\theta}_{n,m}}{1-D_{\theta}}-\log\frac{\hat{D}^{\theta}_{n,m}}{D_{\theta}}\Big{)},$

it follows from Kaji and Rockova (2021, Theorem 4.1).

###### Remark 2.2

(Uniform Convergence Rate) The rate of estimation in Lemma 2.1 and Lemma 2.2
is characterized as point-wise. While for each $\theta\in\Theta$ the
estimation error is shrinking at the rate $\delta_{n}$, the multiplication
constant in front the rate could potentially depend on $\theta$. Assuming a
compact parameter support (which is not a-typical in deep learning models; see
e.g. Schmidt-Hieber (2020) or Polson and Rockova (2018) and Wang and Rockova
(2020)) and continuity of the multiplication constant, we can essentially
regard the rate as uniform. Hereafter, we thereby abuse the notation and
tacitly assume that $\delta_{n}$ is the worst rate over all
$\theta\in\Theta\subset\mathbb{R}^{d}$, i.e. the rate with the largest
multiplication constant.

Next, we investigate convergence around the actual KL divergence
$K(p_{0},p_{\theta})$. The next lemma will be utilized later in the proof of
Theorem 3.1. However, it is of independent interest as it shows how the speed
at which the joint error probability (accounting for randomness of both the
observed and fake data $(\bm{X},\tilde{\bm{X}})$) decays in terms of the
estimation error. Below, the probability $P$ corresponds to
$P_{0}^{(n)}\otimes\tilde{P}^{(m)}$, where $\tilde{P}^{(m)}$ is the measure
for $\tilde{\bm{X}}$ (see Remark 2.1).

###### Theorem 2.1

For a given $\theta\in\Theta$, we define for $u>0$ and $\delta_{n}>0$ as in
Lemma 2.1 and for an arbitrarily slowly increasing sequence $C_{n}>0$

$\rho_{n,\theta}(u;C_{n};\delta_{n})\equiv
P\biggl{(}\left|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K(p_{0},p_{\theta})\right|>2u,d_{\theta}(\hat{D}^{\theta}_{n,m},D^{\theta})\leq
C_{n}\delta_{n}\biggr{)}.$ (10)

Under Assumptions 1 and 2, we then have

$\rho_{n}(u;C_{n};\delta_{n})\equiv\sup_{\theta\in\Theta}\rho_{n,\theta}(u;C_{n};\delta_{n})=O\Big{(}\frac{C_{n}\delta_{n}}{u}+\frac{1}{nu^{2}}\Big{)}.$

The proof can be found in Section D.1. Note that the intersecting event above
has probability converging to one, according to Lemma 2.1.

###  2.3 Other KL Estimators

Our ABC results, presented later in Section 3, can be extended to other types
of KL estimators if similar estimation error results as in Theorem 2.1 can be
shown. One example is the k-Nearest Neighbor (kNN) estimator proposed in
Pérez-Cruz (2008). Wang et al. (2009) showed that this estimator is
asymptotically unbiased and mean-square consistent and they propose a data-
dependent choice of $k$ which can improve the convergence speed. Jiang et al.
(2018) assess data discrepancy inside ABC with the special case of 1-nearest
neighbor, which is defined as

$\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})=\frac{d}{n}\sum_{i=1}^{n}\log\frac{\min_{j}\left\|X_{i}-\tilde{X}^{\theta}_{j}\right\|}{\min_{j\neq
i}\left\|X_{i}-X_{j}\right\|}+\log\frac{m}{n-1}.$ (11)

where $d$ is the number of covariates in $\bm{X}$. Zhao and Lai (2020) provide
estimation error rate results for this kNN estimator.

Another route to estimate the KL divergence is via (data-dependent)
partitioning methods. Wang et al. (2005) proposed to estimate the Radon-
Nikodym derivative $\mathrm{d}P_{0}/\mathrm{d}P_{\theta}$ using frequency
counts on a statistically equivalent partition of $\mathbb{R}^{d}$. However,
the computational complexity of their method is exponential in $d$ and the
estimation accuracy deteriorates quickly as the dimension increases. Silva and
Narayanan (2010) further contributed to multivariate data-driven partition
schemes by using a Barron-type histogram-based density estimate. They provide
sufficient conditions on the partitions scheme to make the estimator strongly
consistent.

Lastly, Nguyen et al. (2007) adopted a variational approach to estimate KL by
reframing the estimation problem as a penalized convex risk minimization
problem, where they restrict the estimate to a bounded subset of a Reproducing
Kernel Hilbert Space (RKHS). Convergence rates are then obtained from
empirical process theory on nonparametric M-estimators (van de Geer, 2000). In
an independent contribution, Ghimire et al. (2020) used a discriminator in
RKHS to estimate KL using a similar approach to ours. They showed that the
estimator error bound is related to the complexity of the discriminator in
RKHS.

###  2.4 Generalizations to GAN-style Discrepancy Metrics

Beyond the KL divergence, our classification framework allows us to consider
other discrepancy metrics. Alternatively to (8), we could instead estimate the
reversed KL divergence

$\hat{K}_{\text{reverse}}(\tilde{\bm{X}}^{\theta},\bm{X})=\frac{1}{m}\sum_{i=1}^{m}\log\frac{1-\hat{D}^{\theta}_{n,m}}{\hat{D}^{\theta}_{n,m}}(\tilde{X}_{i}^{\theta})$

which converges to $K(p_{\theta},p_{0})$ and which is still uniquely minimized
at $p_{\theta}=p_{0}$. One can show that the estimation error of this reversed
KL estimator is still $O_{P^{*}}(\delta_{n})$ by following the same techniques
used in Lemma 2.2. The reversed KL divergence is widely used in variational
inference (Jordan et al., 1999; Wainwright and Jordan, 2008). Forward and
reversed KL’s could perform differently when the function class inside the
variational approach is not rich enough. The reversed KL is zero-forcing/mode-
seeking, while the forward KL is mass-covering/mean-seeking (Bishop, 2006).
Another related metric, deployed by Gutmann et al. (2018), is the
classification accuracy (CA) defined as

$\text{CA}(\bm{X},\tilde{\bm{X}}^{\theta})=\frac{1}{n+m}\Big{(}\sum_{i=1}^{n}\hat{D}^{\theta}_{n,m}(X_{i})+\sum_{i=1}^{m}\big{(}1-\hat{D}^{\theta}_{n,m}(\tilde{X}^{\theta}_{i})\big{)}\Big{)}.$
(12)

Since $\hat{D}(\cdot)$ and $\log\frac{\hat{D}}{1-\hat{D}}(\cdot)$ are linked
by a logistic transformation which is Lipschitz-continous, CA can be roughly
regarded as a weighted average of the forward KL divergence and the reversed
KL divergence. All three metrics are associated with Jensen-Shannon (JS)
divergence, each having a slightly different emphasis.

There are many other divergences that can be estimated using different
configurations of (6). Given a set $\mathcal{F}$ of functions from
$\mathcal{X}$ to $\mathbb{R}$, we can define

$d_{\mathcal{F}}(p_{0},p_{\theta})=\sup_{f\in\mathcal{F}}\mathbb{P}_{n}f(X_{i})-\mathbb{P}_{m}^{\theta}f(\tilde{X}^{\theta}_{i}),$
(13)

as the Integral Probability Metric (IPM) (Müller, 1997) associated with the
function class $\mathcal{F}$. Arjovsky et al. (2017) show that different
choices of the function class $\mathcal{F}$ can dramatically change the
topology, as we discuss below, and the regularity of
$d_{\mathcal{F}}(p_{0},p_{\theta})$ as a loss function. Several widely used
discrepancy metrics can be obtained from IPMs. For example, the 1-Wasserstein
distance or the Earth-Mover distance (Villani, 2009) are obtained when
$\mathcal{F}$ is the set of 1-Lipschitz functions. Arjovsky et al. (2017)
proposed Wasserstein GANs, which approximate the 1-Wasserstein distance with
gradient-clipped neural networks. Another example is the Total Variation (TV)
distance, obtained when $\mathcal{F}$ is the set of all measurable functions
bounded between -1 and 1. Next, the Maximum Mean (MM) discrepancy is obtained
when $\mathcal{F}=\\{f\in\mathcal{H}:\left\|f\right\|_{\infty}\leq 1\\}$ for
$\mathcal{H}$ some Reproducing Kernel Hilbert Space (RKHS) associated with a
given kernel $k:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$. Our
classification-based GAN framework thus provides an alternative route towards
implementing ABCs with the Wasserstein distance (Bernton et al., 2019) as well
as MM (Park et al., 2016).

One important aspect of these various metrics is that they induce topologies
of different strengths. A divergence makes it easier for a sequence of
distributions to converge when it induces a weaker topology. Arjovsky et al.
(2017) show that KL induces the strongest one, followed by JS and TV, with
Wasserstein distance being the weakest. However, the computation costs of the
Wasserstein distance are high, i.e. $O((n+m)^{3}\log(n+m))$ to calculate it
exactly (Burkard et al., 2012). The GAN version approximation also requires
weight clipping to enforce the Lipschitz constraint which can be problematic
to implement. Our KL estimator, on the other hand, is more scalable and places
no constraints on the discriminator. The neural network discriminator can
converge faster than traditional nonparametric methods if some low-
dimensionality assumption is satisfied (Kaji et al., 2020; Bauer and Kohler,
2019). In addition, the KL estimator has a natural interpretation through its
link to the log-likelihood ratio. The KL divergence appears in many
theoretical results in Bayesian non-parametrics. For example, Schwartz’s
Theorem (Schwartz, 1965) requires that, for posterior consistency, the prior
should assign positive probability to any KL neighborhood of $p_{0}$. By
adopting the KL estimator within our framework, we can relate ABC to classical
posterior convergence rate results.

##  3 Frequentist Analysis of the ABC Posterior

One way to assess the quality of the posterior distribution is through the
speed at which it contracts around the truth $\theta_{0}$ as
$n\rightarrow\infty$. While the ABC posterior is ultimately an approximation,
it might still concentrate about $\theta_{0}$ at a reasonable rate.
Concentration rates are typically quantified in terms of a prior concentration
(measured in terms of a combination of the KL divergence and the KL variation)
and the entropy of the model. We have a similar prior mass condition (see
(8.4) in Section 8.2 in Ghosal and Van der Vaart (2017)). We denote the KL-
neighborhood of $p_{0}$ by

$B_{2}(p_{0},\epsilon)=\\{{\theta}:K(p_{0},p_{\theta})<\epsilon^{2}\\}.$ (14)

###### Assumption 3

(Prior Mass) There exist some constants $\kappa>0$ and $\xi>0$ such that for
every $0<\epsilon<\xi$ and some constant $C>0$, the prior probability
satisfies $\Pi[B_{2}(p_{0},\epsilon)]\geq C\epsilon^{\kappa}$.

Next, we assume that the parameter $\theta$ is identifiable in the sense that
the KL divergence is locally compatible with the Euclidean norm. This
assumption is adopted from Assumption 3(ii) of Frazier et al. (2018).

###### Assumption 4

(Identification) The density function $p_{\theta}$ is continuous in $\theta$
and for every $\theta$ in some open neighborhood of $\theta_{0}$ satisfies

$\left\|\theta-\theta_{0}\right\|\leq L\times K(p_{0},p_{\theta})^{\alpha}$

for some $L>0$ and $\alpha>0$.

Similarly as in (5.1) in Kleijn and van der Vaart (2006), Assumption 4 ensures
posterior concentration around $\theta_{0}$ when $K(p_{0},p_{\theta})\to 0$.
This holds for many distributions. For example, for the exponential
distribution with a rate parameter $\theta$, we have
$K(p_{0},p_{\theta})=\frac{\theta}{\theta_{0}}-\log(\frac{\theta}{\theta_{0}})-1$.
Since $\log(1+x)=x-\frac{x^{2}}{2}+o(x^{2})$ when $x\to 0$, we have
$K(p_{0},p_{\theta})\geq\frac{1}{2}\theta_{0}^{-2}(\theta-\theta_{0})^{2}$.
For multivariate normal distribution with a known variance $\Sigma$ and an
unknown location $\mu$, we have
$K(p_{0},p_{\theta})=\frac{1}{2}(\mu-\mu_{0})\Sigma^{-1}(\mu-\mu_{0})\geq\frac{1}{2}\rho(\Sigma)^{-1}\left\|\mu-\mu_{0}\right\|^{2}$,
where $\rho(\Sigma)$ is the spectral radius, i.e., the largest eigenvalue, of
a matrix $\Sigma$.

###  3.1 Posterior Concentration Rate

First, we focus on the uniform kernel
$T_{\epsilon}(x)=\mathbb{I}(\left|x\right|\leq\epsilon)$ used in Algorithm 1.
Recall (from Remark 2.1) that
$\tilde{\bm{X}}^{\theta}=g_{\theta}(\tilde{\bm{X}})$ for some suitable mapping
$g_{\theta}(\tilde{\bm{X}})$ where $\tilde{\bm{X}}\sim\tilde{P}^{(m)}$. The
ABC joint posterior (2) is a weighted aggregation of uniform kernels, i.e

$\hat{\pi}^{AR}\left(\theta,\tilde{\bm{X}}\mid\hat{K}[\bm{X},g_{\theta}(\tilde{\bm{X}})]\leq\epsilon\right)=\frac{\tilde{\pi}(\tilde{\bm{X}})\pi(\theta)\mathbb{I}\big{(}\hat{K}[\bm{X},g_{\theta}(\tilde{\bm{X}})]\leq\epsilon\big{)}}{\int\tilde{\pi}(\tilde{\bm{X}})\pi(\theta)\mathbb{I}\big{(}\hat{K}[\bm{X},g_{\theta}(\tilde{\bm{X}})]\leq\epsilon\big{)}\mathrm{d}\tilde{\bm{X}}\mathrm{d}\theta},$
(15)

which yields the following ABC posterior distribution

$\displaystyle\hat{\Pi}^{AR}_{\epsilon_{n}}(A\mid\bm{X})$
$\displaystyle=\frac{\int_{\theta\in
A}\pi(\theta)\tilde{P}^{(m)}\big{(}\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n}\big{)}\mathrm{d}\theta}{\int_{\Theta}\pi(\theta)\tilde{P}^{(m)}\big{(}\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n}\big{)}\mathrm{d}\theta}\quad\text{for
a Borel-measurable $A\subset\Theta$}.$ (16)

The following theorem (a modification of Theorem 1 in Frazier et al. (2018))
quantifies the concentration rate in terms of the tolerance threshold
$\epsilon_{n}$ as well as the rate at which the classification-based KL
estimator can estimate $\mathbb{P}_{n}\log(p_{0}/p_{\theta})$ (as formulated
in (9)).

###### Theorem 3.1

Let Assumptions 1, 2 and 3 hold and take $\delta_{n}$ as in (9) in Lemma 2.2.
Then, as $n\to\infty$ and with $\epsilon_{n}=o(1)$ such that
$n\epsilon_{n}^{2}\rightarrow\infty$ and $C_{n}\delta_{n}=o(\epsilon_{n})$ for
some arbitrarily slowly increasing sequence $C_{n}>0$ we have

$P_{0}^{(n)}\Pi[K(p_{0},p_{\theta})>\lambda_{n}\mid\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n}]=o(1),$
(17)

where
$\lambda_{n}=\epsilon_{n}+M_{n}C_{n}\delta_{n}\epsilon_{n}^{-\kappa}+\sqrt{M_{n}}n^{-1/2}\epsilon_{n}^{-\kappa/2}$
for some arbitrarily slowly increasing sequence $M_{n}>0$. Moreover, if
Assumption 4 also holds, as $n\to\infty$, we have

$P_{0}^{(n)}\Pi[\left\|\theta-\theta_{0}\right\|>L\lambda_{n}^{\alpha}\mid\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n}]=o(1).$
(18)

The proof of the theorem is provided in Section D.2. Thus, the convergence
rate of our ABC posterior depends on three components: the accept-reject
threshold $\epsilon_{n}$, the estimation error of the KL estimator
$\delta_{n}$ and the rate of discrepancy $n^{-1/2}$ between the empirical and
true KL divergence. Since $\delta_{n}$ will typically be greater than the
parametric rate $n^{-1/2}$, the overall convergence rate is then driven by
$\lambda_{n}=\epsilon_{n}+\widetilde{M}_{n}\delta_{n}\epsilon_{n}^{-\kappa}$,
where $\widetilde{M}_{n}$ is an arbitrarily slowly increasing sequence.

In practice, it is unclear how to properly choose $\epsilon_{n}$. In Algorithm
2, we proposed to weight the draws using a scaled exponential kernel
$\exp(-n\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta}))$. We denote the ABC posterior
under the exponential kernel as $\hat{\Pi}^{EK}(\cdot\mid\bm{X})$ where

$\displaystyle\hat{\Pi}^{EK}(A\mid\bm{X})$
$\displaystyle=\frac{\int_{A}\tilde{P}^{(m)}\exp\big{(}-n\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\big{)}\pi(\theta)\mathrm{d}\theta}{\int_{\Theta}\tilde{P}^{(m)}\exp\big{(}-n\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\big{)}\pi(\theta)\mathrm{d}\theta}.$
(19)

To gain more insights into the ABC posterior behavior under the exponential
kernel, we take a closer look at the “likelihood function” above

$\displaystyle\tilde{P}^{(m)}\exp\big{(}-n\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\big{)}$
$\displaystyle=\frac{p_{\theta}^{(n)}}{p_{0}^{(n)}}\tilde{P}^{(m)}\mathrm{e}^{u_{\theta}},$

where
$u_{\theta}(\bm{X},\tilde{\bm{X}}^{\theta})=-n\times\Big{(}\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-\mathbb{P}_{n}\log\frac{p_{0}}{p_{\theta}}\Big{)}.$
From the equations above, we can write

$\hat{\pi}^{EK}(\theta\mid\bm{X})\propto
p_{\theta}^{(n)}(\bm{X})\times\mathrm{e}^{\hat{u}_{\theta}(\bm{X})}\times\pi(\theta)\quad\text{with}\quad\hat{u}_{\theta}(\bm{X})=\log\int\mathrm{e}^{u_{\theta}(\bm{X},\tilde{\bm{X}}^{\theta})}\mathrm{d}\tilde{P}^{(m)}(\tilde{\bm{X}}).$
(20)

When $\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})$ is the classification-based
estimator, $\hat{u}_{\theta}(\bm{X})$ can be related to the random generator
setting of the Metropolis-Hastings MHC algorithm in Kaji and Rockova (2021)
which has (20) as its stationary distribution. Similarly as in Appendix
Section 5 of Kaji and Rockova (2021), we can regard the posterior
approximation in (20) as a posterior $\hat{\pi}^{EK}(\theta\mid\bm{X})\propto
q_{\theta}^{(n)}\tilde{\pi}(\theta)$ under a misspecified likelihood

$q_{\theta}^{(n)}=\frac{p_{\theta}^{(n)}(\bm{X})\mathrm{e}^{\hat{u}_{\theta}(\bm{X})}}{C_{\theta}}\quad\text{
where }\quad
C_{\theta}=\int_{\mathcal{X}}p_{\theta}^{(n)}(\bm{X})\mathrm{e}^{\hat{u}_{\theta}(\bm{X})}\mathrm{d}\bm{X}$
(21)

and a modified prior $\tilde{\pi}(\theta)\propto\pi(\theta)C_{\theta}$. Since
the likelihood is misspecified, the ABC posterior concentrates around a
projection point $\theta^{*}$ defined as

$\theta^{*}=\arg\min_{\theta\in\Theta}-P_{0}^{(n)}\log[q_{\theta}^{(n)}/p_{0}^{(n)}],$
(22)

which corresponds to the mis-specified model that is closest to $P_{0}^{(n)}$
in the KL sense (Kleijn and van der Vaart, 2006). Kaji and Rockova (2021)
study posterior concentration of (20). We recall their Theorem 4.5 in Appendix
A. Unlike in Theorem 3.1, the posterior concentration rate here depends both
on the estimation error $\delta_{n}$ and the actual concentration rate of the
true posterior, not the acceptance threshold.

###### Remark 3.1

(Vanishing Bias) To better understand the severity of the centering bias of
the misspecified model, we note

$\displaystyle-P_{0}^{(n)}\log[q_{\theta^{*}}^{(n)}/p_{0}^{(n)}]$
$\displaystyle\leq-
P_{0}^{(n)}\log[q_{\theta_{0}}^{(n)}/p_{0}^{(n)}]=P_{0}^{(n)}\log\frac{p_{0}^{(n)}}{p_{0}^{(n)}\mathrm{e}^{\hat{u}_{\theta_{0}}(\bm{X})}/C_{\theta_{0}}}$
$\displaystyle=\log
C_{\theta_{0}}-P_{0}^{(n)}{\hat{u}_{\theta_{0}}(\bm{X})}=\log
P_{0}^{(n)}\mathrm{e}^{\hat{u}_{\theta_{0}}(\bm{X})}-P_{0}^{(n)}{\hat{u}_{\theta_{0}}(\bm{X})}.$

This is essentially the Jensen gap. If we have this Jensen gap vanishing when
$n\to\infty$, then we can conclude that the centering bias is also vanishing,
and the ABC posterior in (19) will eventually concentrate at the right
location.

###### Remark 3.2

(Connection to Generalized Bayesian Inference) The exponential kernel has an
alternative interpretation as a generalized posterior within the Generalized
Bayesian Inference (GBI) framework (Bissiri et al., 2016). Schmon et al.
(2020) show that one can view the ABC posterior constructed from the scaled
exponential kernel as

$p_{l}(\theta\mid\bm{X})\propto\int\exp(-w\cdot
l(\bm{X},\tilde{\bm{X}}^{\theta}))\pi(\theta)\mathrm{d}\tilde{P}^{(m)}(\tilde{\bm{X}}),$

where $\l(\cdot,\cdot)$ is a “loss function” which can be any distance between
the summary statistics or a data discrepancy metric used in ABC, and where the
weight $w>0$ can be chosen to adjust the influence of the loss on the
posterior. The exponential kernel can thus be applied more generally to other
summary statistics or data discrepancy measures.

###  3.2 Shape of the Limiting ABC Posterior Distribution

We now analyze the limiting shape of
$\hat{\Pi}^{AR}_{\epsilon_{n}}(\cdot\mid\bm{X})$ defined in (16). We focus on
the case when $\epsilon_{n}\gg\delta_{n}^{1/(\kappa+1)}$, where $\kappa$ was
defined in Assumption 3, since the posterior is not guaranteed to converge
when the decision threshold $\epsilon_{n}$ is smaller than the estimation
error $\delta_{n}$ of the KL estimator.

###### Assumption 5

Assume that for every $\varepsilon>0$, we have
$\inf_{\|\theta-\theta_{0}\|>\varepsilon}K(p_{0},p_{\theta})>0$. In addition,
assume that $\log p_{\theta}$ is twice differentiable with respect to $\theta$
and that, for every $\theta$ in some neighborhood of $\theta_{0}$, the
remainder of the second order Taylor expansion of
$K(p_{0},p_{\theta})=P_{0}\log\frac{p_{0}}{p_{\theta}}(x)$ around $\theta_{0}$
is comparatively small relative to the second-order term, i.e.

$\displaystyle K(p_{0},p_{\theta})$
$\displaystyle=\nabla_{\theta=\theta_{0}}K(p_{0},p_{\theta})(\theta-\theta_{0})+\frac{1}{2}(\theta-\theta_{0})^{\prime}\nabla^{2}_{\theta=\theta_{0}}K(p_{0},p_{\theta})(\theta-\theta_{0})\\{1+o(1)\\}$
$\displaystyle=\frac{1}{2}(\theta-\theta_{0})^{\prime}I(\theta_{0})(\theta-\theta_{0})\\{1+o(1)\\},$

where
$I(\theta)=\nabla^{2}_{\theta}K(p_{0},p_{\theta})=P_{0}[(\nabla_{\theta}\log
p_{\theta})^{2}]$ is the Fisher information matrix.

When Assumption 5 is satisfied, the condition in Assumption 4 is immediately
satisfied as well and the identification of $\theta_{0}$ is guaranteed. In the
open-neighborhood of $\theta_{0}$, since we have

$(\theta-\theta_{0})^{\prime}I(\theta_{0})(\theta-\theta_{0})=2K(p_{0},p_{\theta})\\{1+o(1)\\},$

and $I(\theta_{0})$ is positive definite, the convergence in the bilinear form
$(\theta-\theta_{0})^{\prime}I(\theta_{0})(\theta-\theta_{0})$ ensures the
convergence of the parameters.

###### Theorem 3.2

Assume that the prior function $\pi(\cdot)$ is continuous around $\theta_{0}$.
Then, under Assumptions 1, 2, 3 and 5, if
$\lim_{n}\delta_{n}/\epsilon_{n}^{\kappa+1}\to 0$, the average posterior
distribution of $\epsilon_{n}^{-1/2}(\theta-\theta_{0})$ converges to the
uniform distribution over the ellipse $\\{w:w^{\prime}I(\theta_{0})w\leq 2\\}$
where $I(\theta)$ is the Fisher information matrix defined in Assumption 5. In
particular, as $n\to 0$, we have

$P_{0}^{(n)}\int
f\Big{(}\epsilon_{n}^{-1/2}(\theta-\theta_{0})\Big{)}\hat{\Pi}^{AR}_{\epsilon_{n}}(\theta\mid\bm{X})\to\int_{u^{\prime}I(\theta_{0})u\leq
2}f(u)\mathrm{d}u/\int_{u^{\prime}I(\theta_{0})u\leq 2}\mathrm{d}u$

for every continuous and bounded function
$f(\cdot):\mathcal{X}\rightarrow\mathbb{R}$.

The proof is provided in Section D.3.

###### Remark 3.3

Theorem 3.2 is adapted from the case (i) in Theorem 2 of Frazier et al.
(2018). We only consider situations when
$\epsilon_{n}\gg\delta_{n}^{1/(\kappa+1)}$ with the prior shrinkage parameter
$\kappa$ defined in Assumption 3. In other words, we assume that the ABC
decision threshold $\epsilon_{n}$ is dominating both the estimation error
$\delta_{n}$ and the asymptotic error $n^{-1/2}$ and, thereby, determines the
posterior concentration rate. It is not entirely obvious how the posterior
would behave when the threshold $\epsilon_{n}$ shrinks faster than the
estimation error $\delta_{n}$.

For the asymptotic behavior of the ABC posterior (20) induced by the
exponential kernel, we resort to BvM characterizations under misspefication in
LAN models (Kleijn and van der Vaart, 2012). When the posterior (20)
concentrates around $\theta^{*}$ in (22) at the rate $\epsilon_{n}^{*}$, one
can show under a suitable LAN condition that the approximate posterior
converges to a sequence of normal distributions in total variation at the rate
$\epsilon_{n}^{*}$. The centering and the asymptotic covariance matrix both
depend on $\theta^{*}$. The formal statement and the proof is in Kaji and
Rockova (2021).

##  4 Simulations

In this section, we illustrate our approach and make comparisons with other
likelihood-free inference techniques. Within our KL-ABC framework, we include
two types of KL estimators. One is obtained with the logit discriminator
score, which we refer to as KL estimation via classification (KLC), and the
other one is estimated via the kNN method (kNN) with $k=1$ (Jiang et al.,
2018). For both estimators, we aggregate ABC samples with the accept-reject
kernel as in Algorithm 1 and the exponential kernel as in Algorithm 2. The
latter will be denoted with a suffix ‘exp’, e.g. KLC-exp or kNN-exp. The
discriminator used for each dataset will be specified later. We provide more
discussions on discriminator calibrations in Appendix C. Besides the two
examples presented here, we have another analysis on the g-and-k distribution
in Appendix B.

The ABC discrepancy metrics we choose for comparisons are (1) the
classification accuracy (CA) (Gutmann et al., 2018) defined as (12); (2) the
Wasserstein (WA) distance under the Euclidean metric (Bernton et al., 2019)
defined as
$\text{WA}(\bm{X},\tilde{\bm{X}}^{\theta})=\min_{\gamma}[\sum_{i=1}^{n}\sum_{j=1}^{m}\gamma_{ij}\left\|X_{i}-\tilde{X}^{\theta}_{j}\right\|^{2}]^{1/2}$
s.t. $\gamma^{\prime}{\bf 1}_{m}={\bf 1}_{n},\gamma^{\prime}{\bf 1}_{n}={\bf
1}_{m}$ with $0\leq\gamma_{ij}\leq 1$; (3) $\ell_{2}$-distance between summary
statistics (SS) and we use the semi-automatic (SA) method (Fearnhead and
Prangle, 2012) if no candidate summary statistics are given; (4) auxiliary
likelihood $AL=\frac{1}{m}\ln
p_{A}(\tilde{\bm{X}}^{\theta}\mid\hat{\phi}(\tilde{\bm{X}}^{\theta}))-\frac{1}{m}\ln
p_{A}(\tilde{\bm{X}}^{\theta}\mid\hat{\phi}(\bm{X}))$ proposed by Drovandi et
al. (2011), where $p_{A}(x\mid\phi)$ is a $d$-dimensional Gaussian
distribution with $\phi$ being the sample mean and covariance. For the
classification accuracy (CA), we use the same discriminator as the one in KL
estimation. In each experiment, unless otherwise noted, we set the tolerance
threshold $\epsilon$ adaptively such that $1\,000$ of $100\,000$ (i.e. the top
1%) proposed ABC samples are accepted.

###  4.1 M/G/1-queuing Model

Because queuing models are usually easy to simulate from, but have no
tractable likelihoods, they have been frequently used as test cases in the ABC
literature, see e.g. Fearnhead and Prangle (2012) and Bernton et al. (2019).
Here, we choose the same setup as in Jiang et al. (2017). Each datum is a
5-dimensional vector consisting of the first five inter-departure times
$x_{i}=(x_{i1},x_{i2},x_{i3},x_{i4},x_{i5})^{\prime}$. In the model, the
service times $u_{ik}$ follow a uniform distribution
$U[\theta_{1},\theta_{2}]$, and the arrival times $w_{ik}$ are exponentially
distributed with the rate $\theta_{3}$. We only observe the interdeparture
times $x_{i}$, given by the process
$x_{ik}=u_{ik}+\max(0,\sum_{j=1}^{k}w_{ij}-\sum_{j=1}^{k-1}x_{ij})$. We
perform ABC on $n=500$ observed samples which are generated from the true
parameter $\theta_{0}=(1,5,0.2)$. The prior on
$(\theta_{1},\theta_{2}-\theta_{1},\theta_{3})$ is uniform on
$[0,10]^{2}\times[0,0.5]$. 333We place the uniform prior on
$\theta_{2}-\theta_{1}$ instead of $\theta_{2}$, since $\theta_{2}$ must be
larger than $\theta_{1}$. This is used in Jiang et al. (2017).

We train the discriminator as a $\ell_{1}$-penalized logistic regressor onto
degree-2 polynomials of the data, including quadratic and interaction terms.
The penalty parameter $\lambda$ is selected via 5-fold cross-validation. We
implement the discriminator with the R package glmnet.

Figure 1: ABC reconstructed posteriors under M/G/1-queuing model with
$\theta_{0}=(1,5,0.2)^{\prime}$. The vertical black lines mark the true
values. The blue curves in kNN and KLC boxes mark a smoothed density
calculated from the exponential kernel.

The shape of the ABC posteriors is given in Figure 1. The plot reveals that,
among the three parameters, $\theta_{1}$ is the hardest one to estimate where
all methods, except for the semi-automatic ABC (SA), gave relatively flat
posterior estimates. Regarding $\theta_{2}$, our methods, SA and W2 seem to
center well around the truth and W2 gives the tightest estimation. Regarding
$\theta_{3}$, it is obvious that SA does not provide enough information. Next,
we repeat the experiments on 10 different datasets, and summarize the average
squared estimation errors and the width of the 95% credible intervals in Table
1. Overall, we see that although SA captures $\theta_{1}$ and $\theta_{2}$
well, it fails to identify the correct location of $\theta_{3}$. Our KLC
method outperforms the kNN estimator in estimating all three parameters. The
two variants of KLC and CA give relatively good and robust estimates for all
three parameters. These methods tend to be similar in both the point estimate
and credible intervals, which is not entirely unexpected since they are
derived from the same discriminators.

| $\theta_{1}=1$ | $\theta_{2}=5$ | $\theta_{3}=0.2$
---|---|---|---
Method | $(\hat{\theta}_{1}-\theta_{1})^{2}$ | 95% CI width | $(\hat{\theta}_{2}-\theta_{2})^{2}$ | 95% CI width | $(\hat{\theta}_{3}-\theta_{3})^{2}$ | 95%CI width
KLC | 0.197 | 3.116 | 0.217 | 4.599 | 0.308 | 0.064
KLC-exp | 0.169 | 2.851 | 0.312 | 3.708 | 0.234 | 0.030
kNN | 0.525 | 3.135 | 0.106 | 3.986 | 3.659 | 0.094
kNN-exp | 1.057 | 2.664 (0.8) | 0.431 |  3.331 (0.9) | 2.634 | 0.072
AL | 0.259 | 3.423 | 0.057 | 5.651 | 0.575 | 0.044
CA | 0.191 | 3.107 | 0.235 | 4.602 | 0.280 | 0.064
SA | 0.180 | 2.457 | 0.355 | 3.514 | 45.297 | 0.446
W2 | 0.595 | 3.631 | 0.039 | 4.871 | 3.846 | 0.052 (0.8)

Table 1: ABC performance on the M/G/1 queuing model over 10 repetitions with
top 1% ABC samples selected. Most of the 95% CIs have full coverage with the
rest having their coverage marked in brackets. The bold fonts mark the best
model in each metric.

###  4.2 Lotka-Volterra Model

The Lotka-Volterra (LV) predator-prey model (Wilkinson, 2018) describes
population evolutions in ecosystems where predators interact with prey. It is
one of the classical stochastic kinetic network model examples. The state of
the population is prescribed deterministically via a system of ordinary
differential equations (ODEs). Inference for such models is challenging
because the transition density is intractable. However, simulation from the
model is possible, which makes it a natural candidate for ABC methods.

The model monitors population sizes of predators $X_{t}$ and prey $Y_{t}$ over
time $t$. The changes in states are determined by four parameters
$\theta=(\theta_{1},\ldots,\theta_{4})^{\prime}$ controlling: (1) the rate
$r_{1}^{t}=\theta_{1}X_{t}Y_{t}$ of a predator being born; (2) the rate
$r_{2}^{t}=\theta_{2}X_{t}$ of a predator dying; (3) the rate
$r_{3}^{t}=\theta_{3}Y_{t}$ of a prey being born; (4) the rate
$r_{4}^{t}=\theta_{4}X_{t}Y_{t}$ of a prey dying. Given the initial population
sizes $(X_{0},Y_{0})$ at time $t=0$, the dynamics can be simulated using the
Gilliespie algorithm (Gillespie, 1977). The algorithm samples times to an
event from an exponential distribution (with a rate $\sum_{j=1}^{4}r_{j}^{t}$)
and picks one of the four reactions with probabilities proportional to their
individual rates $r_{j}^{t}$.

We use the same simulation setup as Kaji and Rockova (2021). Each simulation
is started at $X_{0}=50$ and $Y_{0}=100$ and state observations are recorded
every 0.1 time units for a period of 20 time units, resulting in a series of
$T=201$ observations each. The real data ($n=20$ time series) are generated
with true values $\theta_{0}=(0.01,0.5,1,0.01)^{\prime}$. The predator-prey
interaction dynamic is very sensitive to parameter changes. For example,
Figure 4 of Kaji and Rockova (2021) shows that a slight perturbation in
$\theta_{2}$ leads to significant changes in the population renewal cycle. We
rely on the ability of the discriminator to tell such different patterns
apart. The sensitivity of the model to minor parameter changes is confirmed
with heat-map plots of the estimated KL divergence as a function of
$(\theta_{1},\theta_{4})^{\prime}$ in Figure 2. Figure 2(a) provides a plot of
the estimated KL over the region $[0,0.1]^{2}$ where, apparently, the majority
of the region is flat and uninformative with a sharp spike around the true
values at $\theta_{1}=\theta_{4}=0.01$. We thus narrow the investigation down
to a smaller region $[0,0.02]^{2}$ in Figure 2(b). Again, the curvature in the
estimated KL around the truth is quite steep. This may pose some issues for
Metropolis-Hasting algorithms, since the majority of the prior region is
uninformative and improper initialization could lead to extremely slow
convergence.

(a) $(\theta_{1},\theta_{4})\in(0,0.1)^{2}$

(b) $(\theta_{1},\theta_{4})\in(0,0.02)^{2}$

Figure 2: Estimated $\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})$ (via
Classification) over a grid of $(\theta_{1},\theta_{4})$ values. The other two
parameters are fixed as $\theta_{2}=0.5,\theta_{3}=1$.

Previous ABC analyses of this model suggested various summary statistics
including the mean, log-variance, autocorrelation (at lag 1 and 2) of each
series as well as their cross-correlation (Papamakarios and Murray, 2016). For
the discriminator of our method, we choose the $\ell_{1}$-penalized (LASSO)
logistic regression classifier with $m=n$ and with a penalty $\lambda$
selected via 5-fold cross-validation (as implemented in the R package glmnet).

Similar to Kaji and Rockova (2021), we use an informative prior $\theta\in
U(\Xi)$ with a restricted domain
$\Xi=[0,0.1]\times[0,1]\times[0,2]\times[0,0.1]$ so that the computation is
more economic and efficient. A typical snapshot of the ABC posteriors is
plotted in Figure 3. From the plot, we can see that the difficulty in
estimating different parameters varies a lot and $\theta_{3}$ is the most
challenging among all. While other methods give relatively flat posteriors for
$\theta_{3}$, our method (KLC) combined with the scaled exponential kernel
identifies the correct location of the parameter with a much tighter
posterior. The average performance of ABC methods under this model is
summarized in Table 2. We can see that KLC with the exponential kernel gives
the tightest CI most of the time, while maintaining relatively small
estimation errors.

Figure 3: ABC posterior densities for the Lotka-Volterra model with
$\theta_{0}=(0.01,0.5,1,0.01)^{\prime}$. The black vertical lines mark the
true parameter values. The blue curves in kNN and KLC boxes represent the
smoothed density calculated from the exponential kernel. The ABC posteriors
are plotted from the top 1% out of $10^{5}$ samples.

| $\theta_{1}=0.01$ | $\theta_{2}=0.5$ | $\theta_{3}=1$ | $\theta_{4}=0.01$
---|---|---|---|---
Method | $(\hat{\theta}_{1}-\theta_{1})^{2}$ | 95% CI width | $(\hat{\theta}_{2}-\theta_{2})^{2}$ | 95% CI width | $(\hat{\theta}_{3}-\theta_{3})^{2}$ | 95% CI width | $(\hat{\theta}_{4}-\theta_{4})^{2}$ | 95% CI width
(scale) | $(10^{-5})$ |  |  |  | $(10^{-2})$ |  | $(10^{-5})$ |
KLC | 42.949 | 0.0447 | 0.0169 | 0.626 | 0.119 | 0.945 | 17.159 | 0.0472
KLC-exp | 0.066 | 0.0035 | 0.0025 | 0.285 | 0.405 | 0.643 | 0.021 | 0.0063
kNN | 0.563 | 0.0146 | 0.0011 | 0.602 | 0.054 | 0.934 | 0.203 | 0.0133
kNN-exp | 0.049 | 0.0054 (0.9) | 0.0004 | 0.243 (0.9) | 0.346 | 0.599 (0.9) | 0.066 | 0.0062 (0.9)
CA | 36.477 | 0.0438 | 0.0160 | 0.629 | 0.118 | 0.945 | 15.684 | 0.0462
Stats | 1.639 | 0.0178 | 0.0226 | 0.721 | 1.019 | 0.901 | 3.767 | 0.0166
W2 | 1.265 | 0.0167 | 0.0048 | 0.599 | 0.011 | 0.917 | 4.080 | 0.0205

Table 2: ABC performance evaluations on the Lotka-Volterra Model, averaged
over 10 repetitions, with the top 1% ABC samples selected. Most of the 95% CIs
have full coverage with the rest having their coverage marked in brackets. The
bold fonts mark the best model in each metric.

##  5 Empirical Analysis

We further demonstrate our approach on the nontrivial problem of estimating
stock volatility using merely daily observations on high, low and closing
prices. All of these price observations are typically available to investors.
We use a similar data generating process as in Magdon-Ismail and Atiya (2003),
assuming that the assets follow a Brownian motion with a constant drift and
volatility. In particular, suppose that the log-price processes
$X_{j}(t),i=1,\ldots,d$, are correlated Brownian motions, that is
$E[X_{i}(s)X_{j}(t)]=\sigma_{ij}\min\\{s,t\\}$, and that the joint movement of
the log-price processes $\bm{X}(t)=(X_{1}(t),\ldots,X_{d}(t))^{\prime}$
follows a multivariate Brownian motion as

$d\bm{X}(t)=\bm{\mu}dt+\Sigma d\bm{W}(t),$ (23)

where $\bm{\mu}=(\mu_{1},\ldots,\mu_{d})^{\prime}$ and
$\Sigma=[\sigma_{ij}]_{1\leq i,j\leq d}$ denote the drift and the volatility
of the log processes, respectively. We write

$H_{j}=\max_{0\leq t\leq 1}X_{j}(t),\qquad L_{j}=\min_{0\leq t\leq
1}X_{j}(t),\qquad S_{j}=X_{j}(1),$

for the high, low and final log price, respectively, over a fixed time
interval $[0,1]$. We want to estimate the drift $\bm{\mu}$ and the volatility
matrix $\Sigma$ merely from observing these three prices over a period of
time.

We impose a normal-inverse-Wishart prior $(\bm{\mu},\Sigma)\sim
NIW(\bm{\mu}_{0},\lambda,\Phi,\nu)$. This distribution can be sampled from in
two steps: (1) sample $\Sigma$ from an inverse Wishart distribution
$\Sigma\mid\Phi,\nu\sim W^{-1}(\Phi,\nu)$; (2) sample $\bm{\mu}$ from a
multivariate normal distribution $\bm{\mu}\mid\bm{\mu}_{0},\lambda,\Sigma\sim
N(\bm{\mu}_{0},\frac{1}{\lambda}\Sigma)$. Since $\Sigma$ is a semi-positive
definite matrix, we model the parameters through its Cholesky root
$\Sigma^{1/2}$. Without loss of generality, we only consider the case
$W(0)=(0,\ldots,0)^{\prime}$, since the closing prices on the previous day or
the opening prices of today are usually known.

###  5.1 Synthetic Data

To compare various likelihood-free estimators, we first generate synthetic
data for 1 000 trading days. For each day (of length $t=1$), we simulate the
Brownian motion using 500 time steps to obtain the high, low and closing price
data for each particular window. We first restrict our attention to the case
of just two assets, which leaves us with $5$ parameters to estimate:
$\mu_{1},\mu_{2}$ and the upper triangular root of $\Sigma$, denoted with
$L={\scriptsize\left[\begin{array}[]{cc}l_{11}&0\\\
l_{12}&l_{22}\end{array}\right]}$. We first illustrate how the covariance
parameter $\sigma_{12}$ impacts the co-movement of asset prices. Holding
$\mu_{1}=\mu_{2}=0$ and $\sigma_{11}=\sigma_{22}=1$ fixed, we plot time series
realizations of the closing prices for three particular choices of
$\sigma_{12}$ in Figure 4. The patterns are as expected where the prices tend
to co-fluctuate when $\sigma_{12}$ is closer to one. The success of our method
depends on how well the discriminator can tell apart these trajectories.

(a) $\sigma_{12}=0.05$ (b) $\sigma_{12}=0.5$ (c) $\sigma_{12}=0.95$

Figure 4: Closing prices time series for three choices of $\sigma_{12}$
($\mu_{1}=\mu_{2}=0,\sigma_{11}=\sigma_{22}=1$).

For our simulation, we set $\mu=(0,0)^{\prime}$ and
$\sigma_{11}=\sigma_{22}=1,\sigma_{12}=0.5$. We choose relatively non-
informative prior hyper-parameters
$\bm{\mu}_{0}=(0,0)^{\prime},\lambda=1,\Phi=I_{d}$ and $\nu=d$. Posterior
distributions reconstructed with different ABC methods are given in Figure 5
and the averaged performance over 10 repetitions is summarized in Table 3. It
is clear that our exponential kernel methodsplace more mass around the true
location of the parameters, and that the shape of kNN-exp is slightly less
regular than KLC-exp. Although KLC-exp induces a larger bias in estimating the
drift $(\mu_{1},\mu_{2})$, it does a better job at capturing the correct
location of the volatilities.

(a) ABC Figure 5: Posterior densities on simulated volatility data ($d=2$)

|  | KLC | KLC-exp | kNN | kNN-exp | CA | SS | WA
---|---|---|---|---|---|---|---|---
$\mu_{1}=0$ | MSE$(\times 10^{-4})$ | 2.004 | 2.113 | 0.638 | 0.159 | 1.497 | 1.216 | 0.170
95% CI width | 1.760 | 0.331 | 1.480 | 0.442 | 1.660 | 1.734 | 0.831
$\mu_{2}=0$ | MSE$(\times 10^{-4})$ | 0.734 | 7.439 | 0.363 | 0.259 | 0.142 | 0.021 | 0.277
95% CI width | 1.628 | 0.431 | 1.435 | 0.538 | 1.564 | 1.569 | 0.813
$\sigma_{11}=1$ | MSE$(\times 10^{-3})$ | 0.311 | 0.111 | 0.167 | 0.460 | 0.002 | 1.404 | 2.634
95% CI width | 0.963 | 0.371 | 0.954 | 0.329 | 0.896 | 0.763 | 1.140
$\sigma_{12}=0.5$ | MSE$(\times 10^{-3})$ | 1.035 | 0.103 | 0.427 | 0.107 | 0.861 | 0.109 | 1.528
95% CI width | 1.193 | 0.355 | 1.092 | 0.324 | 1.122 | 0.628 | 0.988
$\sigma_{22}=1$ | MSE$(\times 10^{-3})$ | 0.416 | 0.304 | 0.124 | 0.018 | 0.007 | 3.687 | 4.244
95% CI width | 0.959 | 0.252 | 0.974 | 0.314 | 0.920 | 0.783 | 1.157

Table 3: Performance on the stock volatility estimation example, averaged over
10 repetitions, with top 1% selected. All 95% CIs have full coverage of the
true parameters. The bold fonts mark the best model in each row.

###  5.2 Real Data Analysis

Following the example in Rogers and Zhou (2008), we examine a small data set
of stock prices focusing on two stocks: Boeing (BA) and Proctor & Gamble (PG).
The prices were obtained from NYSE (Yahoo Finance), starting from 3rd January
2011 and consisting of $1\,000$ trading days. Since the off-market trades
follow a different mechanism than the market trade, we only model price
changes from the opening price to the closing price each day where the log
price differences $X(t)$ are all computed based on the opening prices of that
day.

We observe that the fluctuations in the prices are much smaller than in our
simulated time series and we thereby choose the hyperparameters based on the
mean and covariance of the closing prices. Figure 6 gives the ABC posterior
distributions estimated from different methods and the corresponding summaries
of the distributions are provided in Table 4. We again observe that our
methods with the exponential kernels return much narrower posterior
distributions. The estimates on the volatilities are quite close except for
the ABC with summary statistics (SS). The drifts of both BA and PG are not
significantly different from zero.

Figure 6: Posterior densities estimated from log-prices of BA and PG

|  | KLC | KLC-exp | kNN | kNN-exp | CA | SS | W2
---|---|---|---|---|---|---|---|---
$\mu_{1}$ $(\times 10^{-3})$ | $\bar{\mu}_{1}$ | -0.774 | 1.354 | 0.250 | 0.779 | 0.342 | -0.071 | -0.126
$l$ | -10.380 | -1.981 | -6.816 | -1.145 | -7.099 | -9.532 | -3.987
$u$ | 8.202 | 3.514 | 7.229 | 3.385 | 7.395 | 10.121 | 3.736
$\mu_{2}$ $(\times 10^{-3})$ | $\bar{\mu}_{2}$ | 0.009 | -0.039 | 0.491 | 0.822 | 0.445 | 0.450 | 0.453
$l$ | -6.591 | -1.718 | -4.329 | -1.406 | -4.604 | -9.831 | -3.491
$u$ | 6.531 | 4.563 | 5.503 | 2.166 | 5.538 | 10.171 | 4.450
$\sigma_{11}$$(\times 10^{-4})$ | $\bar{\sigma}_{11}$ | 1.406 | 1.142 | 1.049 | 1.054 | 1.101 | 0.738 | 1.130
$l$ | 0.824 | 0.076 | 0.633 | 0.844 | 0.634 | 0.155 | 0.626
$u$ | 2.362 | 1.444 | 1.640 | 1.202 | 1.819 | 1.976 | 1.938
$\sigma_{12}$$(\times 10^{-5})$ | $\bar{\sigma}_{12}$ | 2.811 | 2.446 | 1.767 | 2.084 | 1.906 | -0.376 | 2.052
$l$ | -2.379 | -0.37 | -2.359 | 0.494 | -2.390 | -6.660 | -3.137
$u$ | 8.226 | 3.629 | 6.037 | 4.298 | 6.442 | 6.065 | 7.253
$\sigma_{22}$$(\times 10^{-4})$ | $\bar{\sigma}_{22}$ | 0.692 | 0.526 | 0.492 | 0.482 | 0.547 | 0.552 | 0.495
$l$ | 0.407 | 0.089 | 0.277 | 0.365 | 0.315 | 0.133 | 0.173
$u$ | 1.105 | 0.698 | 0.802 | 0.559 | 0.883 | 1.285 | 0.984

Table 4: Posterior estimates on analysis of BA and PG. For each parameter, we
report three summary statistics, the posterior means, the lower limit of the
95% CI intervals ($l$) and the upper limit of the 95% CI intervals ($u$).

##  6 Discussion

This paper develops an ABC variant using a classification-based KL estimator
as a discrepancy measure. By deploying a flexible classifier, the empirical KL
divergence can be estimated with a vanishing error. In addition, inspired by
the connection between the KL divergence and the log-likelihood ratio, we
propose a scaled exponential kernel to aggregate ABC samples. This smoothing
variant avoids the need for choosing the ad hoc threshold $\epsilon_{n}$ and
fully utilizes information returned from all ABC samples. Under mild
conditions, we show that the posterior concentration rate of the accept-reject
ABC depends on the estimation error $\delta_{n}$ and the accept-reject
threshold $\epsilon_{n}$, while the rate of the smooth version depends on the
estimation error $\delta_{n}$ and the contraction rate of the actual posterior
distribution. Our methodology can also be related to many other likelihood-
free inference methods, including ABC with Classification Accuracy (Gutmann et
al., 2018), Wasserstein distance ABC (Bernton et al., 2019), and Generalized
Posteriors (Schmon et al., 2020). Along with our theoretical investigations,
we demonstrate competitive performance of our methods on benchmark examples.
Our theoretical analysis provides theoretical justifications for the method of
Gutmann et al. (2018).

## References

  * Åkesson et al. (2020) Åkesson, M., Singh, P., Wrede, F., and Hellander, A. (2020).  Convolutional neural networks as summary statistics for approximate Bayesian computation.  arXiv.
  * Arjovsky et al. (2017) Arjovsky, M., Chintala, S., and Bottou, L. (2017).  Wasserstein generative adversarial networks.  In International conference on machine learning, pages 214–223. PMLR.
  * Bauer and Kohler (2019) Bauer, B. and Kohler, M. (2019).  On deep learning as a remedy for the curse of dimensionality in nonparametric regression.  The Annals of Statistics, 47(4):2261–2285.
  * Beaumont et al. (2002) Beaumont, M. A., Zhang, W., and Balding, D. J. (2002).  Approximate Bayesian computation in population genetics.  Genetics, 162(4):2025–2035.
  * Bernton et al. (2019) Bernton, E., Jacob, P. E., Gerber, M., and Robert, C. P. (2019).  Approximate Bayesian computation with the Wasserstein distance.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 81(2):235–269.
  * Bishop (2006) Bishop, C. M. (2006).  Pattern recognition and Machine Learning.  Springer-Verlag, Berlin, Heidelberg.
  * Bissiri et al. (2016) Bissiri, P. G., Holmes, C. C., and Walker, S. G. (2016).  A general framework for updating belief distributions.  Journal of the Royal Statistical Society. Series B, Statistical methodology, 78(5):1103.
  * Blum et al. (2013) Blum, M. G., Nunes, M. A., Prangle, D., and Sisson, S. A. (2013).  A comparative review of dimension reduction methods in approximate Bayesian computation.  Statistical Science, 28(2):189–208.
  * Burkard et al. (2012) Burkard, R., Dell’Amico, M., and Martello, S. (2012).  Assignment problems: revised reprint.  SIAM.
  * Cover and Thomas (1991) Cover, T. M. and Thomas, J. A. (1991).  Entropy, relative entropy and mutual information.  Elements of information theory, 2(1):12–13.
  * Den Hollander (2008) Den Hollander, F. (2008).  Large deviations, volume 14.  American Mathematical Soc.
  * Drovandi et al. (2011) Drovandi, C. C., Pettitt, A. N., and Faddy, M. J. (2011).  Approximate Bayesian computation using indirect inference.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 60(3):317–337.
  * Fearnhead and Prangle (2012) Fearnhead, P. and Prangle, D. (2012).  Constructing summary statistics for approximate bayesian computation: Semi-automatic approximate bayesian computation.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 74(3):419–474.
  * Frazier et al. (2018) Frazier, D. T., Martin, G. M., Robert, C. P., and Rousseau, J. (2018).  Asymptotic properties of approximate bayesian computation.  Biometrika, 105(3):593–607.
  * Ghimire et al. (2020) Ghimire, S., Gyawali, P. K., and Wang, L. (2020).  Reliable estimation of kullback-leibler divergence by controlling discriminator complexity in the reproducing kernel hilbert space.  arXiv.
  * Ghosal and Van der Vaart (2017) Ghosal, S. and Van der Vaart, A. (2017).  Fundamentals of nonparametric Bayesian inference, volume 44.  Cambridge University Press.
  * Gillespie (1977) Gillespie, D. T. (1977).  Exact stochastic simulation of coupled chemical reactions.  The journal of physical chemistry, 81(25):2340–2361.
  * Goodfellow et al. (2014) Goodfellow, I. J., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A. C., and Bengio, Y. (2014).  Generative adversarial nets.  In Advances in Neural Information Processing Systems, pages 2672–2680.
  * Gutmann et al. (2018) Gutmann, M. U., Dutta, R., Kaski, S., and Corander, J. (2018).  Likelihood-free inference via classification.  Statistics and Computing, 28(2):411–425.
  * Jiang et al. (2018) Jiang, B., Wu, T.-Y., and Wong, W. H. (2018).  Approximate Bayesian computation with Kullback-Leibler divergence as data discrepancy.  In International Conference on Artificial Intelligence and Statistics, pages 1711–1721. PMLR.
  * Jiang et al. (2017) Jiang, B., Wu, T.-y., Zheng, C., and Wong, W. H. (2017).  Learning summary statistic for approximate Bayesian computation via deep neural network.  Statistica Sinica, pages 1595–1618.
  * Jordan et al. (1999) Jordan, M. I., Ghahramani, Z., Jaakkola, T. S., and Saul, L. K. (1999).  An introduction to variational methods for graphical models.  Machine learning, 37(2):183–233.
  * Joyce and Marjoram (2008) Joyce, P. and Marjoram, P. (2008).  Approximately sufficient statistics and Bayesian computation.  Statistical applications in genetics and molecular biology, 7(1).
  * Kaji et al. (2020) Kaji, T., Manresa, E., and Pouliot, G. (2020).  An adversarial approach to structural estimation.  arXiv.
  * Kaji and Rockova (2021) Kaji, T. and Rockova, V. (2021).  Metropolis-hastings via classification.  arXiv.
  * Kleijn and van der Vaart (2006) Kleijn, B. J. and van der Vaart, A. W. (2006).  Misspecification in infinite-dimensional bayesian statistics.  The Annals of Statistics, 34(2):837–877.
  * Kleijn and van der Vaart (2012) Kleijn, B. J. and van der Vaart, A. W. (2012).  The bernstein-von-mises theorem under misspecification.  Electronic Journal of Statistics, 6:354–381.
  * Kullback (1958) Kullback, S. (1958).  Information theory and statistics.  Wiley, New York.
  * Magdon-Ismail and Atiya (2003) Magdon-Ismail, M. and Atiya, A. F. (2003).  A maximum likelihood approach to volatility estimation for a Brownian motion using high, low and close price data.  Quantitative Finance, 3(5):376.
  * Müller (1997) Müller, A. (1997).  Integral probability metrics and their generating classes of functions.  Advances in Applied Probability, pages 429–443.
  * Nguyen et al. (2007) Nguyen, X., Wainwright, M. J., and Jordan, M. I. (2007).  Estimating divergence functionals and the likelihood ratio by penalized convex risk minimization.  In Proceedings of the 20th International Conference on Neural Information Processing Systems, pages 1089–1096.
  * Nunes and Balding (2010) Nunes, M. A. and Balding, D. J. (2010).  On optimal selection of summary statistics for approximate Bayesian computation.  Statistical applications in genetics and molecular biology, 9(1).
  * Papamakarios and Murray (2016) Papamakarios, G. and Murray, I. (2016).  Fast $\varepsilon$-free inference of simulation models with bayesian conditional density estimation.  In Proceedings of the 30th International Conference on Neural Information Processing Systems, pages 1036–1044.
  * Park et al. (2016) Park, M., Jitkrittum, W., and Sejdinovic, D. (2016).  K2-ABC: Approximate Bayesian computation with kernel embeddings.  In Artificial Intelligence and Statistics, pages 398–407. PMLR.
  * Pérez-Cruz (2008) Pérez-Cruz, F. (2008).  Kullback-Leibler divergence estimation of continuous distributions.  In 2008 IEEE international symposium on information theory, pages 1666–1670. IEEE.
  * Polson and Rockova (2018) Polson, N. G. and Rockova, V. (2018).  Posterior concentration for sparse deep learning.  In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pages 938–949.
  * Pritchard et al. (1999) Pritchard, J. K., Seielstad, M. T., Perez-Lezaun, A., and Feldman, M. W. (1999).  Population growth of human Y chromosomes: a study of Y chromosome microsatellites.  Molecular biology and evolution, 16(12):1791–1798.
  * Rogers and Zhou (2008) Rogers, L. C. and Zhou, F. (2008).  Estimating correlation from high, low, opening and closing prices.  The Annals of Applied Probability, 18(2):813–823.
  * Schmidt-Hieber (2020) Schmidt-Hieber, J. (2020).  Nonparametric regression using deep neural networks with relu activation function.  The Annals of Statistics, 48(4):1875–1897.
  * Schmon et al. (2020) Schmon, S. M., Cannon, P. W., and Knoblauch, J. (2020).  Generalized posteriors in approximate Bayesian computation.  In Third Symposium on Advances in Approximate Bayesian Inference.
  * Schwartz (1965) Schwartz, L. (1965).  On Bayes procedures.  Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete, 4(1):10–26.
  * Silva and Narayanan (2007) Silva, J. and Narayanan, S. (2007).  Universal consistency of data-driven partitions for divergence estimation.  In 2007 IEEE International Symposium on Information Theory, pages 2021–2025. IEEE.
  * Silva and Narayanan (2010) Silva, J. and Narayanan, S. S. (2010).  Information divergence estimation based on data-dependent partitions.  Journal of Statistical Planning and Inference, 140(11):3180–3198.
  * Sisson et al. (2018) Sisson, S. A., Fan, Y., and Beaumont, M. A. (2018).  Overview of abc.  In Handbook of approximate Bayesian computation, pages 3–54. Chapman and Hall/CRC.
  * Sørensen (2004) Sørensen, H. (2004).  Parametric inference for diffusion processes observed at discrete points in time: a survey.  International Statistical Review, 72(3):337–354.
  * van de Geer (2000) van de Geer, S. A. (2000).  Empirical Processes in M-estimation, volume 6.  Cambridge university press.
  * van der Vaart and Wellner (1996) van der Vaart, A. W. and Wellner, J. A. (1996).  Weak convergence and empirical processes: with applications to statistics.  Springer.
  * Villani (2009) Villani, C. (2009).  Optimal transport: old and new, volume 338.  Springer.
  * Wainwright and Jordan (2008) Wainwright, M. J. and Jordan, M. I. (2008).  Graphical models, exponential families, and variational inference.  Now Publishers Inc.
  * Wang et al. (2005) Wang, Q., Kulkarni, S. R., and Verdú, S. (2005).  Divergence estimation of continuous distributions based on data-dependent partitions.  IEEE Transactions on Information Theory, 51(9):3064–3074.
  * Wang et al. (2009) Wang, Q., Kulkarni, S. R., and Verdú, S. (2009).  Divergence estimation for multidimensional densities via $k$-nearest-neighbor distances.  IEEE Transactions on Information Theory, 55(5):2392–2405.
  * Wang and Rockova (2020) Wang, Y. and Rockova, V. (2020).  Uncertainty quantification for sparse deep learning.  In International Conference on Artificial Intelligence and Statistics, pages 298–308. PMLR.
  * Wilkinson (2018) Wilkinson, D. J. (2018).  Stochastic modelling for systems biology.  CRC press.
  * Wilkinson (2013) Wilkinson, R. D. (2013).  Approximate bayesian computation (abc) gives exact results under the assumption of model error.  Statistical applications in genetics and molecular biology, 12(2):129–141.
  * Wood (2010) Wood, S. N. (2010).  Statistical inference for noisy nonlinear ecological dynamic systems.  Nature, 466(7310):1102–1104.
  * Yarotsky (2017) Yarotsky, D. (2017).  Error bounds for approximations with deep relu networks.  Neural Networks, 94:103–114.
  * Zhao and Lai (2020) Zhao, P. and Lai, L. (2020).  Analysis of k nearest neighbor KL divergence estimation for continuous distributions.  In 2020 IEEE International Symposium on Information Theory (ISIT), pages 2562–2567. IEEE.

##  Appendix A Frequentist’s analysis on the exponential kernel

We thus study the concentration in terms of a KL neighborhood around
$Q^{(n)}_{\theta^{*}}$ defined as

$B(\epsilon,Q^{(n)}_{\theta^{*}};P_{0}^{(n)})=\\{Q^{(n)}_{\theta^{*}}\in\mathcal{Q}^{(n)}:\tilde{K}(\theta^{*},\theta)\leq
n\epsilon^{2},\tilde{V}(\theta^{*},\theta)\leq n\epsilon^{2}\\},$ (24)

where $\tilde{K}(\theta^{*},\theta)\equiv
P_{0}^{(n)}\log\frac{q_{\theta^{*}}^{(n)}}{q_{\theta}^{(n)}}$ and
$\tilde{V}(\theta^{*},\theta)\equiv
P_{0}^{(n)}\left|\log\frac{q_{\theta^{*}}^{(n)}}{q_{\theta}^{(n)}}-\tilde{K}(\theta^{*},\theta)\right|^{2}$.

The following corollary is directly adopted from Theorem 5.1 of Kaji and
Rockova (2021).

###### Corollary A.1

Denote with $\widetilde{Q}_{\theta}^{(n)}$ a measure defined through
$\mathrm{d}\widetilde{Q}_{\theta}^{(n)}=({p_{0}^{(n)}}/{q_{\theta^{*}}^{(n)}})\mathrm{d}P_{\theta}^{(n)}$
and let $d(\cdot,\cdot)$ be a semi-metric on $\mathcal{P}^{(n)}$. Suppose that
there exists a sequence $\epsilon_{n}>0$ satisfying $\epsilon_{n}\to 0$ and
$n\epsilon_{n}^{2}\to\infty$ such that for every $\epsilon>\epsilon_{n}$ there
exists a test $\phi_{n}$ (depending on $\epsilon$) such that for every
$J\in\mathbb{N}_{0}$

$P_{0}^{(n)}\phi_{n}\lesssim e^{-n\epsilon^{2}/4}\quad\text{ and
}\sup_{Q_{\theta}^{(n)}:d(Q_{\theta}^{(n)},P_{\theta^{*}}^{(n)})>J\epsilon}\widetilde{Q}_{\theta}^{(n)}(1-\phi_{n})\leq
e^{-nJ^{2}\epsilon^{2}/4}.$ (25)

Let $B(\epsilon,Q^{(n)}_{\theta^{*}};P_{0}^{(n)})$ be as in (24) and let
$\tilde{\Pi}_{n}(\theta)$ be a prior distribution with a density
$\tilde{\pi}(\theta)\propto C_{\theta}\pi(\theta)$ with $C_{\theta}$ as in
(21). Assume that there exists a constant $L>0$ such that, for all $n$ and
$j\in\mathbb{N}$,

$\frac{\tilde{\Pi}_{n}\left(\theta\in\Theta:j\epsilon_{n}<d(Q_{\theta}^{(n)},P_{\theta^{*}}^{(n)})\leq(j+1)\epsilon_{n}\right)}{\tilde{\Pi}_{n}\left(B(\epsilon_{n},Q^{(n)}_{\theta^{*}};P_{0}^{(n)})\right)}\leq
e^{n\epsilon_{n}^{2}j^{2}/8}.$ (26)

Then for every sufficiently large constant $M$, as $n\to\infty$,

$P_{0}^{(n)}\Pi^{EK}\left(Q_{\theta}^{(n)}:d(Q_{\theta}^{(n)},P_{\theta^{*}}^{(n)})\geq
M\epsilon_{n}\mid X^{(n)}\right)\to 0.$ (27)

Next, we want to show the shape of the posterior is actually asymptotically
gaussian around $\theta^{*}$. The following corollary follows from Theorem 2.1
of Kleijn and van der Vaart (2012) and Lemma 8.1 of Kaji and Rockova (2021).

###### Corollary A.2

(Bernstein von-Mises) Assume that the posterior (21) concentrates around
$\theta^{*}$ at the rate $\epsilon_{n}^{*}$ and that for every compact
$K\in\mathbb{R}^{d}$

$\sup_{h\in
K}\left|\log\frac{q_{\theta^{*}+\epsilon_{n}^{*}h}^{(n)}(\bm{X})}{q_{\theta^{*}}^{(n)}(\bm{X})}-h^{\prime}\tilde{V}_{\theta^{*}}\tilde{\delta}_{n,\theta^{*}}-\frac{1}{2}h^{\prime}\tilde{V}_{\theta^{*}}h\right|\to
0$ (28)

for some random vector $\tilde{\delta}_{n,\theta^{*}}$ and a non-singular
matrix $\tilde{V}_{\theta^{*}}$. Then the approximated posterior
$\Pi^{EK}(\cdot)$ converges to a sequence of normal distributions in total
variation at the rate $\epsilon_{n}^{*}$, i.e.

$\sup_{B}\left|\Pi^{EK}\left({\epsilon_{n}^{*}}^{-1}(\theta-\theta^{*}\right)\in
B\mid\bm{X})-N_{\tilde{\delta}_{n,\theta^{*}},\tilde{V}_{\theta^{*}}^{-1}(B)}\right|\to
0\quad\text{ in }P_{0}^{(n)}-probability.$ (29)

##  Appendix B Simulation: Univariate g-and-k Distributions

Another classical example in the ABC literature is the univariate $g$-and-$k$
distribution. It is defined implicitly by its inverse distribution function

$F^{-1}(x)=A+B\Big{[}1+0.8\frac{1-\mathrm{e}^{-gz_{x}}}{1+\mathrm{e}^{-gz_{x}}}\Big{]}(1+z_{x}^{2})^{k}z_{x},$

where $z_{x}$ is the $x$-th quantile of the standard normal distribution, and
parameters $A,B,g,k$ are related to location, scale, skewness and kurtosis.
The probability density function has no analytical form but can be numerically
calculated with high precision since it only involves one-dimensional
inversions and differentiations of the quantile function. Therefore, Bayesian
inference can be carried out.

We generate $n=500$ observations from the model by letting
$A=3,B=1,g=2,k=0.5$. Among the four parameters, $g$ is the hardest to
identify, as observed in previous analyses. From a pilot run, narrow the range
of our uniform priors to $A\sim U([2,4]),B\sim U([0,0.25]),g\sim U([0,10])$
and $k\sim U([0,2.5])$.

We estimate the KL divergence by a neural network with two hidden layers with
10 nodes each. The first layer is activated with the rectified linear unit
(ReLU) function, while the second layer is equipped with a hyperbolic tangent
function (tanh). In addition to observed data $\bm{X}$, we also include
higher-order terms $\bm{X}^{2}$ and $\bm{X}^{3}$ in the input.444We are trying
to manually create more features for this univariate data. We also explore
even higher order terms and find the including up to $\bm{X}^{3}$ works the
best in practice. We compare the posteriors obtained from different ABC
methods in Figure 7. We can see that, indeed, the parameter $g$ is the most
difficult one to estimate. While the majority of methods return a flat
posterior, KLC and kNN with the exponential kernel still place most of the
posterior mass around the true value. AL and SA misidentify the location of
one or two parameters. While CA and KLC give relatively similar ranking on the
ABC sampled parameters, KLC with the exponential kernel can utilize all the
samples and re-adjust the shape according to weights. In addition, KLC-exp
returns more spiky posteriors than kNN-exp. This is also reflected in the
performance summaries in Table 5 which show that KLC-exp has a narrower CI
than kNN-exp on average in terms of estimating $A,B$ and $k$. We observe that
kNN could be advantageous than KLC in estimating $g$-and-$k$ distributions,
since the data is univariate and the use of classifiers may not be necessary.

Figure 7: Posterior densities for the $g$-and-$k$ distributions
$(A=3,B=1,g=2,k=0.5)$. The black vertical lines mark the true parameter
values. The blue curves in kNN and KLC boxes represent the weighted density
calculated from the exponential kernel. The ABC posteriors were plotted with
the top 1% of $10^{5}$ samples.

| $A=3$ | $B=1$ | $g=2$ | $k=0.5$
---|---|---|---|---
Method | $(\hat{A}-A)^{2}$ | 95% CI width | $(\hat{B}-B)^{2}$ | 95% CI width | $(\hat{g}-g)^{2}$ | 95% CI width | $(\hat{k}-k)^{2}$ | 95 %CI width
KLC | 0.0143 | 0.535 | 0.0369 | 1.528 | 12.840 | 8.356 | 0.00341 | 1.044
KLC-exp | 0.0156 | 0.455 | 0.0339 | 1.229 | 12.516 | 8.191 | 0.00800 | 0.784
kNN | 0.0024 | 0.843 | 0.0486 | 1.770 | 2.695 | 7.558 | 0.00421 | 1.086
kNN-exp | 0.0032 | 0.610 | 0.0410 | 1.333 | 1.673 | 4.109 | 0.00480 | 0.805
AL | 0.0176 | 1.325 | 0.2211 | 1.593 | 8.704 | 9.480 | 0.04451 | 0.835
CA | 0.0125 | 0.503 | 0.0466 | 1.568 | 13.077 | 8.162 | 0.00574 | 1.080
SA | 0.1191 | 0.424 (0) | 0.0598 | 2.379 | 8.800 | 9.478 | 0.56399 | 2.379
WA | 0.0145 | 0.627 | 0.0240 | 1.034 | 11.621 | 8.566 | 0.01060 | 0.776

Table 5: Performance summaries on $g$-and-$k$ distributions over $10$
repetitions with top 1% ABC samples selected. Most of the 95% CIs have full
coverage with the rest having their coverage marked after the CI width. The
bold fonts mark the best model in each metric.

##  Appendix C Discriminator calibration

We provide more intuition for how the choice of the discriminator, including
the $m/n$ ratio and the number of sets of latent variables needed to
approximate the KL divergence, affects our results. We also discuss possible
modifications we make to the KL estimation to eliminate bias. We use the
M/G/1-queuing model as an example.

First, we want to investigate how the choice of the discriminator and the
$m/n$ ratio could impact the estimation of the KL divergence. We consider
three discriminators: (1) logistic regression on degree-2 polynomial terms of
$\bm{X}$ (LRD); (2) a neural network with one layer of 10 ReLu nodes and two
layers of 10 tanh nodes with input $\bm{X}$ (NND1); (3) a neural network with
one layer of 10 ReLu nodes and one layer of 10 tanh nodes with input as the
degree-2 polynomial terms of $\bm{X}$ (NND2). For each discriminator, we do a
grid search over $m/n\in\\{1,2,3,5\\}$. We plot the estimated KL divergence
(rescaled) in Figure 8. From the plot we can see it is quite obvious that LRD
best identifies the location of the true values, and the fluctuations/variance
in its estimations are the smallest. It also worth noting that the
difficulties in estimating different parameters are significantly different
here. All discriminators are able to learn the location of $\theta_{3}$ with
little variance, but only LRD is able to learn $\theta_{1}$ and $\theta_{2}$
well. The $m/n$ ratio does not seem to have significant impact for LRD, but we
observe for NND1, NND2 and some other models that the variance in the
estimation tends to get smaller when $m/n>1$ but not overly large to avoid
bias. We continue our experiments with LRD and $m/n=3$.

Figure 8: Estimated KL divergence (rescaled) under different discriminators
(nlatent=10). We plot the marginal changes in estimated KL divergence with
respect to different parameters and different $m/n$ ratios here. The estimated
values are rescaled by subtracting the minimum in each scenario so the changes
are more comparable.

Next, we want to know how many different sets of latent variables are needed
to control the variance in estimation. We focus on the best discriminator we
have previously found (i.e. LRD with $m/n=3$). We include the cases where
$nlatent\in\\{1,5,10\\}$ here. The results are shown in Figure 9. It is
expected that with more sets of latent variables to average over, the variance
in estimation gets smaller. The main challenge for the M/G/1 model lies in
estimating $\theta_{1}$ correctly and stably. With only one or five sets of
latent variable, it seems hopeless to identify $\theta_{1}$. To trade off
variance and computation costs, we set $\text{nlatent}=10$ in our experiments.

Figure 9: Estimated KL divergence (rescaled) using LRD and $m/n=3$. We plot
the marginal changes in estimated KL divergence with respect to different
parameters and different number of sets of latent variables. The estimated
values are rescaled by subtracting the minimum in each scenario so the changes
are more comparable.

##  Appendix D Proofs

###  D.1 Proof of Theorem 2.1

Denote $K_{n}=\mathbb{P}_{n}\log\frac{p_{0}}{p_{\theta}}$. Using Chebyshev’s
inequality and Assumption 2 as

$\displaystyle P_{0}^{(n)}[\left|K_{n}-K(p_{0},p_{\theta})\right|>u]$
$\displaystyle=P_{0}^{(n)}\Big{(}\left|(\mathbb{P}_{n}-P_{0})\log\frac{p_{0}}{p_{\theta}}\right|>u\Big{)}\leq\frac{1}{u^{2}}P_{0}^{(n)}\Big{[}\left|(\mathbb{P}_{n}-P_{0})\log\frac{p_{0}}{p_{\theta}}\right|^{2}\Big{]}$
$\displaystyle=\frac{1}{u^{2}}P_{0}^{(n)}\Big{[}\left|\mathbb{P}_{n}\big{(}\log\frac{p_{0}}{p_{\theta}}-P_{0}\log\frac{p_{0}}{p_{\theta}}\big{)}\right|^{2}\Big{]}=\frac{1}{nu^{2}}P_{0}\Big{[}\left|\log\frac{p_{0}}{p_{\theta}}-P_{0}\log\frac{p_{0}}{p_{\theta}}\right|^{2}\Big{]}$
$\displaystyle\leq\frac{16(2+\Lambda)h(p_{0},p_{\theta})}{nu^{2}},$

where the last inequality follows from Lemma 2.1 (iii) of Kaji and Rockova
(2021). Next, note

$\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K_{n}=-\mathbb{P}_{n}\Big{(}\log\frac{1-\hat{D}^{\theta}_{n,m}}{1-D_{\theta}}-\log\frac{\hat{D}^{\theta}_{n,m}}{D_{\theta}}\Big{)}.$

Recall the outer expectation $P^{*}$ in Lemma 2.1 and the set of classifiers
$D^{\theta}_{n,\delta}$ defined in Assumption 1, we can bound

$\displaystyle
P\biggl{(}\left|\mathbb{P}_{n}\Big{(}\log\frac{1-\hat{D}^{\theta}_{n,m}}{1-D_{\theta}}-\log\frac{\hat{D}^{\theta}_{n,m}}{D_{\theta}}\Big{)}\right|>u,d_{\theta}(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\biggr{)}$ $\displaystyle\qquad\qquad\leq
P^{*}\Big{(}\sup_{D\in\mathcal{D}^{\theta}_{C_{n}\delta_{n}}}\left|\mathbb{P}_{n}\Big{(}\log\frac{1-D}{1-D_{\theta}}-\log\frac{D}{D_{\theta}}\Big{)}\right|>u\Big{)}$
$\displaystyle\qquad\qquad\leq\frac{1}{u}\mathbb{E}^{*}\sup_{D\in\mathcal{D}^{\theta}_{C_{n}\delta_{n}}}\left|\mathbb{P}_{n}\Big{(}\log\frac{1-D}{1-D_{\theta}}-\log\frac{D}{D_{\theta}}\Big{)}\right|$

by Markov’s inequality. The proof of Theorem 4.1 of Kaji and Rockova (2021)
shows that the expectation is $O(C_{n}\delta_{n})$. Using the triangle
inequality and the Bonferroni inequality, since $h(p_{0},p_{\theta})\leq 2$,
we can then write

$\displaystyle
P\biggl{(}\left|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K(p_{0},p_{\theta})\right|>2u,d_{\theta}(\hat{D}^{\theta}_{n,m},D^{\theta})\leq
C_{n}\delta_{n}\biggr{)}$ $\displaystyle\quad\leq
P\biggl{(}\left|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K_{n}\right|+\left|K_{n}-K(p_{0},p_{\theta})\right|>2u,d_{\theta}(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\biggr{)}$ $\displaystyle\quad\leq
P\biggl{(}\left|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K_{n}\right|>u,d_{\theta}(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\biggr{)}+P_{0}^{(n)}[\left|K_{n}-K(p_{0},p_{\theta})\right|>u]$
$\displaystyle\quad\leq
O\Big{(}\frac{C_{n}\delta_{u}}{u}\Big{)}+\frac{32(2+\Lambda)}{nu^{2}}.$

###  D.2 Proof of Theorem 3.1

Throughout, we continue to assume that
$\tilde{\bm{X}}^{\theta}=g_{\theta}(\tilde{\bm{X}})$ and we denote with
$P=P_{0}^{(n)}\otimes\tilde{P}^{(m)}$ the joint measure for
$(\bm{X},\tilde{\bm{X}})$. Below, we will be using the notation $\Pi(\cdot)$
to denote the generic probability, i.e. for $(\theta,\tilde{\bm{X}})$ or for
the conditional probability $\theta$ given $\tilde{\bm{X}}$. Later, we will
define a high-probability event $\Omega_{n}(C,\varepsilon_{n})$ such that
$P_{0}^{(n)}[\Omega_{n}(C,\varepsilon_{n})^{c}]=o(1)$ for some $C\in(0,1)$.
Given $\delta_{n}>0$ from our assumptions, we can write for every
$\lambda_{n}>0$ and $\epsilon_{n}>0$ and for every arbitrarily slowly
increasing sequence $C_{n}>0$

$P_{0}^{(n)}\Pi\left(K(p_{0},p_{\theta})>\lambda_{n}\middle|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n}\right)\leq\Pi_{1}+o(1)+\\\
P_{0}^{(n)}\Pi\left(K(p_{0},p_{\theta})>\lambda_{n}\middle|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n},d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\right)\mathbb{I}[\Omega_{n}(C,\varepsilon_{n})],$ (30)

where, using Lemma 2.1 and the fact that the rate $\delta_{n}$ is uniform (see
Remark 2.2),

$\Pi_{1}\equiv
P_{0}^{(n)}\Pi(d(\hat{D}^{\theta}_{n,m},D_{\theta})>C_{n}\delta_{n})\leq\sup_{\theta\in\Theta}P(d(\hat{D}^{\theta}_{n,m},D_{\theta})>C_{n}\delta_{n})=o(1).$

Consider the joint event, for some $\delta^{\prime}>0$,

$A_{\epsilon_{n}}(\delta^{\prime})=\\{(\tilde{\bm{X}},\theta):\hat{K}[\bm{X},g_{\theta}(\tilde{\bm{X}})]\leq\epsilon_{n}\\}\cap\\{K(p_{0},p_{\theta})>\delta^{\prime}\\}.$

For every $(\tilde{\bm{X}},\theta)\in A_{\epsilon_{n}}(\delta^{\prime})$ we
have

$\displaystyle
K(p_{0},p_{\theta})\leq\hat{K}(\bm{X},g_{\theta}(\tilde{\bm{X}}))+\left|\hat{K}(\bm{X},g_{\theta}(\tilde{\bm{X}}))-K(p_{0},p_{\theta})\right|\leq\epsilon_{n}+\left|\hat{K}(\bm{X},g_{\theta}(\tilde{\bm{X}}))-K(p_{0},p_{\theta})\right|.$

Hence $(\tilde{\bm{X}},\theta)\in A_{\epsilon_{n}}(\delta^{\prime})$ implies
that

$\left|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K(p_{0},p_{\theta})\right|>\delta^{\prime}-\epsilon_{n},$

and choosing $\delta^{\prime}\geq\epsilon_{n}+t_{\epsilon}$ leads to

$\Pi[A_{\epsilon_{n}}(\delta^{\prime})]\leq\int_{\Theta}\tilde{P}^{(m)}\left[\left|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K(p_{0},p_{\theta})\right|>t_{\epsilon}\right]d\Pi(\theta).$

Using (30), we now focus on the conditional probability, given
$d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq C_{n}\delta_{n}$,

$\Pi\left(K(p_{0},p_{\theta})>\epsilon_{n}+t_{\epsilon}\middle|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n},d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\right)\\\
\leq\frac{\int_{\Theta}\tilde{P}^{(m)}\left[\left|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K(p_{0},p_{\theta})\right|>t_{\epsilon}\mid
d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\right]d\Pi(\theta)}{\int_{\Theta}\tilde{P}^{(m)}\left[\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n}\mid
d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq C_{n}\delta_{n}\right]d\Pi(\theta)}.$
(31)

We now find a lower bound for the denominator. Recall the KL neighborhood
$B_{2}(p_{0},\epsilon_{n})$ defined in (14). Since

$\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq
K(p_{0},p_{\theta})+\left|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K(p_{0},p_{\theta})\right|\leq\epsilon_{n}/2+K(p_{0},p_{\theta}),$

provided that
$\left|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K(p_{0},p_{\theta})\right|\leq\epsilon_{n}/2$.
The denominator can be then bounded by

$\displaystyle\int_{\Theta}\tilde{P}^{(m)}[\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n}\mid
d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq C_{n}\delta_{n}]d\Pi(\theta)$
$\displaystyle\qquad\geq\int_{B_{2}(p_{0},\epsilon_{n}/{2})}\tilde{P}^{(m)}\left[\left|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K(p_{0},p_{\theta})\right|\leq\epsilon_{n}/2\mid
d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq C_{n}\delta_{n}\right]d\Pi(\theta)$
$\displaystyle\qquad\geq\Pi[{B_{2}(p_{0},\epsilon_{n}/{2})}]-\int_{B_{2}(p_{0},\epsilon_{n}/{2})}\tilde{P}^{(m)}\left[\left|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K(p_{0},p_{\theta})\right|>\epsilon_{n}/2\mid
d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq C_{n}\delta_{n}\right]d\Pi(\theta).$

Denoting

$Z(\bm{X})\equiv\int_{B_{2}(p_{0},\epsilon_{n}/{2})}\tilde{P}^{(m)}\left[\left|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K(p_{0},p_{\theta})\right|>\epsilon_{n}/2\mid
d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq C_{n}\delta_{n}\right]d\Pi(\theta)$

we can write, for every $C>0$, using Fubini’s theorem and Markov’s inequality

$\displaystyle P_{0}^{(n)}\left(Z(\bm{X})>C\right)$
$\displaystyle\leq\frac{1}{C}\int_{B_{2}(p_{0},\epsilon_{n}/{2})}P\left[\left|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K(p_{0},p_{\theta})\right|>\epsilon_{n}/2\mid
d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq C_{n}\delta_{n}\right]d\Pi(\theta)$
$\displaystyle=\frac{\Pi(B_{2}(p_{0},\epsilon_{n}/{2}))}{C}\frac{\sup\limits_{\theta\in\Theta}P\left[\left|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K(p_{0},p_{\theta})\right|>\epsilon_{n}/2,d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\right]}{\sup\limits_{\theta\in\Theta}P\left[d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\right]}$
$\displaystyle=\frac{\Pi(B_{2}(p_{0},\epsilon_{n}/{2}))}{C(1+o(1))}\sup\limits_{\theta\in\Theta}\rho_{n,\theta}(\varepsilon_{n}/2;C_{n};\delta_{n}),$

where we have used Theorem 2.1 and the fact that
$P(d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq C_{n}\delta_{n})=1+o(1)$ for every
$\theta\in\Theta$ from Lemma 2.1. We now define an event, for some $0<C<1$ and
$\epsilon_{n}>0$,

$\Omega_{n}(C,\epsilon_{n})=\left\\{\bm{X}:Z(\bm{X})\leq
C\times\Pi(B_{2}(p_{0},\epsilon_{n}/{2}))\right\\}.$

Using Theorem 2.1, we have

$\rho_{n,\theta}(\epsilon_{n}/2;C_{n};\delta_{n})=O\Big{(}\frac{C_{n}\delta_{n}}{\epsilon_{n}}+\frac{1}{n\epsilon_{n}^{2}}\Big{)}\quad\text{
for every}\quad\theta\in\Theta.$

Choosing $\epsilon_{n}>0$ such that $\epsilon_{n}=o(1)$ and
$n\epsilon_{n}^{2}\rightarrow\infty$ and $C_{n}\delta_{n}=o(\epsilon_{n})$ we
have $P_{0}^{(n)}[\Omega_{n}(C,\epsilon_{n})^{c}]=o(1)$ for every $C\in(0,1)$.
On the event $\Omega_{n}(C,\epsilon_{n})$, for some $0<C<1$, we can lower-
bound the denominator with

$\int_{\Theta}\tilde{P}^{(m)}[\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n}\mid
d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}]d\Pi(\theta)>(1-C)\times\Pi(B_{2}(p_{0},\epsilon_{n}/{2}).$

Using this bound and applying Fubini’s theorem, we can further write

$P_{0}^{(n)}\Pi\left(K(p_{0},p_{\theta})>\epsilon_{n}+t_{\epsilon}\mid\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n},d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\right)\mathbb{I}[\Omega_{n}(C,\varepsilon_{n})]\\\
\leq\frac{\int_{\Theta}P\left[\left|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K(p_{0},p_{\theta})\right|>t_{\epsilon}\mid
d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\right]d\Pi(\theta)}{(1-C)\times\Pi(B_{2}(p_{0},\epsilon_{n}/{2})}$
(32)

Using Theorem 2.1 again, we obtain an upper bound for the display above with

$\frac{\rho_{n}(t_{\epsilon};C_{n};\delta_{n})}{(1-C)\times\Pi(B_{2}(p_{0},\epsilon_{n}/{2})(1+o(1))}.$

Using the prior Assumption 3, we can choose $t_{\epsilon}$ such that

$\Big{(}\frac{C_{n}\delta_{n}}{t_{\epsilon}}+\frac{1}{nt_{\epsilon}^{2}}\Big{)}/\epsilon_{n}^{\kappa}=1/M_{n}$

for some arbitrarily slowly increasing sequence $M_{n}>0$. We choose
$t_{\epsilon}=M_{n}C_{n}\delta_{n}/\epsilon_{n}^{\kappa}+\sqrt{M}_{n}n^{-1/2}/\epsilon_{n}^{\kappa/2}$.
Since $\delta_{n}\gtrsim n^{-1/2}$ and
$\epsilon_{n}^{-\kappa}\geq\epsilon_{n}^{-\kappa/2}$, the overall rate is then
driven by
$\epsilon_{n}+t_{\epsilon}=\epsilon_{n}+\widetilde{M}_{n}\delta_{n}\epsilon_{n}^{-\kappa}$,
where $\widetilde{M}_{n}=M_{n}C_{n}$.

###  D.3 Proof of Theorem 3.2

Recall from Theorem 3.1 that with the accept-reject strategy, the true KL
divergence $K(p_{0},p_{\theta})$ is contracting at the rate
$\lambda_{n}=\epsilon_{n}+\widetilde{M}_{n}\delta_{n}\epsilon_{n}^{-\kappa}$,
where $\widetilde{M}_{n}$ is a slowly increasing sequence that diverges faster
than $C_{n}$. Consider the case when
$\epsilon_{n}\gg\widetilde{M}_{n}\delta_{n}\epsilon_{n}^{-\kappa}$ or,
equivalently, $\epsilon_{n}\gg\delta_{n}^{1/(\kappa+1)}$. Denote

$x(\theta)=\epsilon_{n}^{-1}K(p_{0},p_{\theta})\quad\text{ and}\quad
f_{n}(\theta-\theta_{0})=f\big{(}\epsilon_{n}^{-1/2}(\theta-\theta_{0})\big{)}.$

We express the ABC posterior expectation of $f_{n}(\theta-\theta_{0})$ for a
non-negative and bounded function $f_{n}(\cdot)$ by

$\displaystyle
P_{0}^{(n)}E_{\hat{\Pi}^{AR}_{\epsilon_{n}}}\Big{[}f_{n}(\theta-\theta_{0})\Big{]}=P_{0}^{(n)}\int
f_{n}(\theta-\theta_{0})\mathrm{d}\hat{\Pi}^{AR}_{\epsilon_{n}}(\theta\mid\bm{X})$
$\displaystyle\quad=\underbrace{P_{0}^{(n)}\int
f_{n}(\theta-\theta_{0})\mathbb{I}[K(p_{0},p_{\theta})\leq\lambda_{n},d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n})]\mathrm{d}\hat{\Pi}^{AR}_{\epsilon_{n}}(\theta\mid\bm{X})}_{(\text{I})}$
$\displaystyle\qquad+\underbrace{P_{0}^{(n)}\int
f_{n}(\theta-\theta_{0})\mathbb{I}[K(p_{0},p_{\theta})\leq\lambda_{n},d(\hat{D}^{\theta}_{n,m},D_{\theta})>C_{n}\delta_{n})]\mathrm{d}\hat{\Pi}^{AR}_{\epsilon_{n}}(\theta\mid\bm{X})}_{(\text{II})}$
$\displaystyle\qquad+\underbrace{P_{0}^{(n)}\int
f_{n}(\theta-\theta_{0})\mathbb{I}[K(p_{0},p_{\theta})>\lambda_{n}]\mathrm{d}\hat{\Pi}^{AR}_{\epsilon_{n}}(\theta\mid\bm{X})}_{(\text{III})}$

where the term $(\text{III})$ can be controlled using Fubini’s theorem and the
concentration result in Theorem 3.1 as follows

$\displaystyle(\text{III})$ $\displaystyle=\int
f_{n}(\theta-\theta_{0})P_{0}^{(n)}\mathbb{I}[K(p_{0},p_{\theta})>\lambda_{n}]\mathrm{d}\hat{\Pi}^{AR}_{\epsilon_{n}}(\theta\mid\bm{X})$
$\displaystyle\leq\left\|f\right\|_{\infty}P_{0}^{(n)}\Pi\Big{(}K(p_{0},p_{\theta})>\lambda_{n}\mid\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n}\Big{)}=o(1).$

The second term $(\text{II})$ can be bounded similarly as

$\displaystyle(\text{II})$
$\displaystyle=P_{0}^{(n)}\int_{K(p_{0},p_{\theta})\leq\lambda_{n}}f_{n}(\theta-\theta_{0})\mathbb{I}[d(\hat{D}^{\theta}_{n,m},D_{\theta})>C_{n}\delta_{n})]\mathrm{d}\hat{\Pi}^{AR}_{\epsilon_{n}}(\theta\mid\bm{X})$
$\displaystyle\leq\left\|f\right\|_{\infty}\int_{K(p_{0},p_{\theta})\leq\lambda_{n}}P(d(\hat{D}^{\theta}_{n,m},D_{\theta})>C_{n}\delta_{n})\mathrm{d}\Pi(\theta)$
$\displaystyle\leq\left\|f\right\|_{\infty}\sup_{\theta\in\Theta}P(d(\hat{D}^{\theta}_{n,m},D_{\theta})>C_{n}\delta_{n})=o(1)$

where we use the fact that
$\sup_{\theta\in\Theta}P(d(\hat{D}^{\theta}_{n,m},D_{\theta})>C_{n}\delta_{n})=o(1)$
from Lemma 2.1.

Thus, the asymptotic behavior is mainly determined by the term $(\text{I})$.
Combined with the continuity of $\pi(\theta)$ at $\theta_{0}$ we can re-write
$(\text{I})$ as

$\displaystyle(\text{I})$
$\displaystyle=P_{0}^{(n)}\frac{\int_{K(p_{0},p_{\theta})\leq\lambda_{n}}\pi(\theta)f_{n}(\theta-\theta_{0})\tilde{P}^{(m)}\big{[}\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n},d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\big{]}\mathrm{d}\theta}{\int_{K(p_{0},p_{\theta})\leq\lambda_{n}}\pi(\theta)\tilde{P}^{(m)}\big{[}\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n},d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\big{]}\mathrm{d}\theta}$
$\displaystyle=P_{0}^{(n)}\frac{\int_{K(p_{0},p_{\theta})\leq\lambda_{n}}f_{n}(\theta-\theta_{0})\tilde{P}^{(m)}\big{[}\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n},d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\big{]}\mathrm{d}\theta}{\int_{K(p_{0},p_{\theta})\leq\lambda_{n}}\tilde{P}^{(m)}\big{[}\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n},d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\big{]}\mathrm{d}\theta}(1+o(1))$

We have $x(\theta)\geq 0$ for all $\theta\in\Theta$ and since

$\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})=\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K(p_{0},p_{\theta})+\epsilon_{n}x(\theta),$

we can write

$\tilde{P}^{(m)}\big{[}\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n},d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\big{]}\\\
\geq\tilde{P}^{(m)}\left[\left|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K(p_{0},p_{\theta})\right|\leq\epsilon_{n}(1-x(\theta)),d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\right]\\\
=1-\tilde{P}^{(m)}\left[\left|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K(p_{0},p_{\theta})\right|>\epsilon_{n}(1-x(\theta)),d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\right].$ (33)

Denoting with

$\tilde{Z}(\bm{X})=\int_{x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}\tilde{P}^{(m)}\left[\left|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K(p_{0},p_{\theta})\right|>\epsilon_{n}\big{(}1-x(\theta)\big{)},d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\right]\mathrm{d}\theta.$

Using Markov’s inequality and Fubini’s theorem we have, for every
$\tilde{C}>0$,

$\displaystyle P_{0}^{(n)}(\tilde{Z}(\bm{X})>\tilde{C})$
$\displaystyle\quad\leq\frac{1}{\tilde{C}}\int\limits_{x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}P\left[\left|\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})-K(p_{0},p_{\theta})\right|>\epsilon_{n}\big{(}1-x(\theta)\big{)},d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\right]\mathrm{d}\theta$
$\displaystyle\quad\leq\frac{1}{\tilde{C}}\int_{x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}\rho_{n}\big{(}\epsilon_{n}(1-x(\theta));C_{n};\delta_{n}\big{)}\mathrm{d}\theta$
$\displaystyle\quad\leq\frac{1}{\tilde{C}}\times\int_{x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}\mathrm{d}\theta\times\sup_{x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}\rho_{n}\big{(}\epsilon_{n}(1-x(\theta));C_{n};\delta_{n}\big{)}$
$\displaystyle\quad\leq\frac{\int_{x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}\mathrm{d}\theta}{\tilde{C}}\times\rho_{n}(M_{n}C_{n}\delta_{n};C_{n};\delta_{n})$
$\displaystyle\quad=\frac{\int_{x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}\mathrm{d}\theta}{\tilde{C}}\times
O\left(\frac{1}{M_{n}}+\frac{1}{nM_{n}^{2}C_{n}^{2}\delta_{n}^{2}}\right).$

where $\rho_{n}(\cdot;C_{n};\delta_{n})$ is defined in Theorem 2.1. Thus, we
can define a set $\tilde{\Omega}_{n}(\tilde{C}_{n})$, for some for some
arbitrarily slowly increasing sequence $\tilde{C}_{n}>0$, and
$\tilde{C}_{n}=O(1/M_{n})$, as

$\tilde{\Omega}_{n}(\tilde{C}_{n})=\left\\{\bm{X}:\tilde{Z}(\bm{X})\leq\tilde{C}_{n}\times\int_{x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}\mathrm{d}\theta\right\\}.$

Then we have that $P_{0}^{(n)}(\tilde{\Omega}_{n}(\tilde{C}_{n})^{c})=o(1)$.
Recall the inequality in (33), on $\tilde{\Omega}_{n}(\tilde{C}_{n})$, we have
that

$\int_{x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}\tilde{P}^{(m)}\big{[}\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n},d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\big{]}\mathrm{d}\theta\\\ \geq\int_{x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}\mathrm{d}\theta-\tilde{Z}(\bm{X})=\int_{x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}\mathrm{d}\theta\times(1-o(1))$

And this quantity is upper-bounded by $\int_{x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}\mathrm{d}\theta$. Therefore, we
can conclude that

$\int_{x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}\tilde{P}^{(m)}\big{[}\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n},d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\big{]}\mathrm{d}\theta=\int_{x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}\mathrm{d}\theta\big{(}1+o(1)\big{)}.$

Note that $x(\theta)\leq 1-M_{n}C_{n}\delta_{n}/\epsilon_{n}$ implies that
$K(p_{0},p_{\theta})\leq\epsilon_{n}-M_{n}C_{n}\delta_{n}$. On this set
$\tilde{\Omega}_{n}(\tilde{C}_{n})$, we can further lower bound the
denominator as

$\displaystyle\int_{K(p_{0},p_{\theta})\leq\lambda_{n}}\tilde{P}^{(m)}\big{[}\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n},d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\big{]}\mathrm{d}\theta$
$\displaystyle\qquad=\underbrace{\int_{K(p_{0},p_{\theta})\leq\epsilon_{n}-M_{n}C_{n}\delta_{n}}\mathrm{d}\theta}_{D_{1}}(1+o(1))$
$\displaystyle\qquad\quad+\underbrace{\int_{\epsilon_{n}-M_{n}C_{n}\delta_{n}<K(p_{0},p_{\theta})\leq\lambda_{n}}\tilde{P}^{(m)}\big{[}\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n},d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\big{]}\mathrm{d}\theta}_{D_{2}}.$

Next, we show that the second term $D_{2}$ is $o(D_{1})$. Let
$u=\epsilon_{n}^{-1/2}(\theta-\theta_{0})$. Under Assumption 5, we have
$K(p_{0},p_{\theta})=\frac{1}{2}(\theta-\theta_{0})^{\prime}I(\theta_{0})(\theta-\theta_{0})\\{1+o(1)\\}$,
which is
$x(\theta_{0}+\sqrt{\epsilon}_{n}u)=\frac{1}{2}u^{\prime}I(\theta_{0})u$.
Since $I(\theta_{0})$ is positive definite we can write

$\displaystyle\frac{D_{2}}{D_{1}}$
$\displaystyle\leq\frac{\int_{1-M_{n}C_{n}\delta_{n}/\epsilon_{n}<x(\theta)\leq
1+\widetilde{M}_{n}\delta_{n}/\epsilon_{n}^{\kappa+1}}\mathrm{d}\theta}{\int_{x(\theta)\leq
1-M_{n}C_{n}\delta_{n}/\epsilon_{n}}\mathrm{d}\theta}$ (34)
$\displaystyle\leq\frac{\int_{2(1-M_{n}C_{n}\delta_{n}/\epsilon_{n})\leq
u^{\prime}I(\theta_{0})u\leq
2(1+\widetilde{M}_{n}\delta_{n}/\epsilon_{n}^{\kappa+1})}\mathrm{d}\theta}{\int_{u^{\prime}I(\theta_{0})u\leq
2(1-M_{n}C_{n}\delta_{n}/\epsilon_{n})}\mathrm{d}\theta}\lesssim\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}+\frac{\widetilde{M}_{n}\delta_{n}}{\epsilon_{n}^{\kappa+1}}=o(1).$
(35)

where the approximation follows from the fact that
$\epsilon_{n}\gg\delta_{n}^{1/(\kappa+1)}$ and because the denominator is
approximating the integral $\int_{u^{\prime}I(\theta_{0})u\leq
2}\mathrm{d}\theta$ and the numerator is the length of shrinking intervals.
Combining the above results, we find that, on the set
$\tilde{\Omega}_{n}(\tilde{C}_{n})$, the denominator can be lower-bounded by
$\int_{x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}\mathrm{d}\theta\\{1+o(1)\\}.$

This implies

$\displaystyle(\text{I})$
$\displaystyle=\frac{\int_{K(p_{0},p_{\theta})\leq\lambda_{n}}f_{n}(\theta-\theta_{0})P\big{[}\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n},d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\big{]}\mathrm{d}\theta}{(1+o(1))\int_{x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}\mathrm{d}\theta}(1+o(1))+o(1)$
$\displaystyle=(N_{1}+N_{2})\\{1+o(1)\\}+o(1).$ (36)

with

$N_{1}\equiv{\frac{\int\limits_{x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}f\Big{(}\epsilon_{n}^{-1/2}(\theta-\theta_{0})\Big{)}\mathrm{d}\theta}{\int\limits_{x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}\mathrm{d}\theta}}$ (37)

and

$N_{2}\equiv\frac{\int\limits_{K(p_{0},p_{\theta})\leq\lambda_{n}}\mathbb{I}\left[{x(\theta)>1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}\right]f\Big{(}\epsilon_{n}^{-1/2}(\theta-\theta_{0})\Big{)}P\Big{[}\hat{K}(\bm{X},\tilde{\bm{X}}^{\theta})\leq\epsilon_{n},d(\hat{D}^{\theta}_{n,m},D_{\theta})\leq
C_{n}\delta_{n}\Big{]}\mathrm{d}\theta}{\int\limits_{x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}\mathrm{d}\theta},$ (38)

where the second equality follows from the fact that $x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}$ leads to
$K(p_{0},p_{\theta})\leq\epsilon_{n}-M_{n}C_{n}\delta_{n}$ and then
$K(p_{0},p_{\theta})\leq\lambda_{n}$ is trivially satisfied and where the last
$o(1)$ comes from the set $\tilde{\Omega}_{n}(\tilde{C}_{n})^{c}$.

Since we have $\epsilon_{n}\gg\delta_{n}$, with
$u=\epsilon_{n}^{-1/2}(\theta-\theta_{0})$, the first term is approximately
equal to

$N_{1}=\frac{\int_{K(p_{0},p_{\theta_{0}+\sqrt{\epsilon_{n}}u})\leq\epsilon_{n}}f(u)\mathrm{d}u}{\int_{K(p_{0},p_{\theta_{0}+\sqrt{\epsilon_{n}}u})\leq\epsilon_{n}}\mathrm{d}u}\left(1+o(1)\right).$

Using Assumption 5 again, we have

$\int_{K(p_{0},p_{\theta_{0}+\sqrt{\epsilon_{n}}u})\leq\epsilon_{n}}\mathrm{d}u=\int_{\frac{1}{2}\sqrt{\epsilon_{n}}u^{\prime}I(\theta)\sqrt{\epsilon_{n}}u\leq\epsilon_{n}}\mathrm{d}u+o(1)=\int_{u^{\prime}I(\theta_{0})u\leq
2}\mathrm{d}u+o(1).$

This leads to

$N_{1}=\frac{\int_{u^{\prime}I(\theta_{0})u\leq
2}f(u)\mathrm{d}u}{\int_{u^{\prime}I(\theta_{0})u\leq
2}\mathrm{d}u}\left(1+o(1)\right).$

Next, we show that the fraction $N_{2}$ converges to 0. The numerator can be
simplified as an integral over
$1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}<x(\theta)\leq
1+\frac{\widetilde{M}_{n}\delta_{n}}{\epsilon_{n}^{\kappa+1}}$, which can be
bounded by (34) as

$\displaystyle N_{2}$
$\displaystyle\leq\frac{\left\|f\right\|_{\infty}\int_{1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}<x(\theta)\leq
1+\frac{\widetilde{M}_{n}\delta_{n}}{\epsilon_{n}^{\kappa+1}}}\mathrm{d}\theta}{\int_{x(\theta)\leq
1-\frac{M_{n}C_{n}\delta_{n}}{\epsilon_{n}}}\mathrm{d}\theta}$
$\displaystyle\leq\left\|f\right\|_{\infty}\frac{D_{2}}{D_{1}}=o(1).$

Since we have $\epsilon_{n}^{\kappa+1}\gg\delta_{n}$, putting all the terms
together, we obtain that the $P_{0}^{(n)}$-averaged ABC posterior distribution
of $\epsilon_{n}^{-1/2}(\theta-\theta_{0})$ is asymptotically uniform over the
ellipsoid $\\{u:u^{\prime}I(\theta_{0})u\leq 2\\}$.