Robust and Efficient Approximate Bayesian Computation: A Minimum Distance Approach

Robust and Efficient Approximate Bayesian Computation: A Minimum Distance Approach

Introduction

The use of complex models has led researchers to construct statistical methods that are capable of conducting reliable inference in situations where the likelihood function of such models is intractable. In Bayesian inference, the class of approximate Bayesian computation (ABC) methods, which replaces evaluation of the likelihood with model simulation, are often employed; for a review on ABC methods see the recent handbook sisson2018handbook.

ABC is predicated on the belief that the observed data is drawn from a given class of models from which it is feasible to generate synthetic data. This fact is then used to construct an approximate posterior for the model unknowns by comparing, in some given distance measure, quantities calculated using observed data and data simulated from the model. While such a procedure allows us to conduct inference in very complex models, the broad applicability of ABC methods has been hindered by at least two factors. First, the choice of summary measures, and as a by product the resulting loss of information due to such a choice. Two, the reliability and robustness of such procedures to deviations from the modeling assumptions.

In many cases, the quantities used to compare observed and synthetic data correspond to moments or sample quantiles. Currently, there is no singularly accepted choice for such “summary statistics”, as different problems require different approaches, and researchers are therefore left to fish for useful summaries. Regardless of the choice of summaries, with few exceptions, the use of such summary measures will inevitably result in a loss of information relative to a procedure that uses the entire sample, such as (infeasible) exact Bayes approaches. More formally, in a frequentest sense, point estimators obtained from ABC procedures will generally be inefficient (see, e.g., LF2016 and FMRR2016).

In addition, it is known since frazier2020model that different versions of ABC can deliver very different point estimators under even mild levels of model misspecification. Moreover, the results of these authors suggest that the more complicated the ABC approach, such as, e.g., nonlinear regression adjustments ABC (blumF2010non), the poorer the performance. In this way, the authors argue Occam's razor should be applied to ABC inference when the idea of model misspecification is entertained.

Motivated by these practical concerns of efficiency and robustness of ABC inference, in Section (ref) we propose a novel ABC approach that can deliver estimators that are as efficient as the infeasible exact Bayes approach, when the model generating synthetic data is correctly specified, but simultaneously displays robustness to deviations from the modeling assumptions. To avoid the use of summarization via summary statistics, we follow the idea of Bernton2017 and view the observed and synthetic data through their empirical probability measures, and propose the use of a norm based on these measures to measure the discrepancy between the observed and synthetic datasets. It is well known that such norms give rise to so-called minimum distance functionals, donoho1988automatic, and so we term this new approach to ABC, minimum distance ABC (MD-ABC). The MD-ABC posterior is shown to be asymptotically normal under regularity conditions that are known to be satisfied for specific choices of the norm. Moreover, if one takes as the norm the Hellinger distance, the resulting ABC approach delivers point estimators that are asymptotically efficient. The later result extends to the Bayesian framework the results on minimum Hellinger distance estimation (MHDE) pioneered by beran1977minimum.

Furthermore, in Section (ref) we demonstrate that, in the class of ABC point estimators, MD-ABC is optimally robust to certain classes of deviations from the assumed model. Using a local model misspecification framework (see, e.g., rieder2012robust), we demonstrate that MD-ABC can display the smallest possible bias and mean squared error in the class of ABC point estimators. This result therefore extends the known robustness of frequentest based MD estimators, see, e.g., donoho1988automatic, in two novel directions: 1) to situations where inference is carried out in the Bayesian framework; 2) to situations where MD estimation must be implemented via simulation (i.e., where the model is intractable). Simulated examples considered in Section (ref) demonstrate that in locally misspecified models, MD-ABC can produce superior inference to exact Bayes inference and competing ABC approaches, such as ABC based on the Wasserstein distance.

Implicit models, and ABC Inference

Implicit Models

Let $(\Omega,\mathcal{F},\mathbb{P})$ denote the intrinsic probability space, with associated expectation operator $\mathbb{E}$, on which all the random variables are defined. Denote by $\mathcal{P}(\mathcal{X})$ the set of probability measures on a space $\mathcal{X}$. We observe $y_{1:n} = (y_1,\ldots,y_n)'$ data points, taking values in $\mathcal{Y}^n$, and which is distributed according to some unknown measure $\mu_\star^{(n)}$.

In most cases, the space $\mathcal{P}(\mathcal{X})$ is specified to be a parametric class of models $$\mathcal{M}^{(n)}=\{\mu^{(n)}_\theta : \theta \in \Theta\} \subset \mathcal{P}(\mathcal{Y}^n),\text{ where }\Theta\subset \mathbb{R}^{d_\theta},$$ and where $\Pi\in\mathcal{P}(\Theta)$ quantifies our prior beliefs about the unknown parameters $\theta$. Regardless of the model complexity, it is often feasible to generate synthetic observations $z_{1:n}$ according to $\mu^{(n)}_\theta$, for any $\theta\in\Theta$, but there are many cases where it is not feasible to evaluate the likelihood associated to $\mu^{(n)}_\theta$. Even in cases where the likelihood is infeasible to calculate, useful information about the model can still be obtained by comparing the observed data, $y_{1:n}$, and the simulated data, $z_{1:n}$.

For $A\in\mathcal{F}$, the empirical measure of $A$ is $\hat{\mu}_n(A) = n^{-1}\sum_{i=1}^n \delta_{y_i}(A)$, where $\delta_y$ denotes the Dirac measure on $y\in\mathcal{Y}$. Likewise, for any $\theta\in\Theta$, let $\hat{\mu}_{\theta,n}(A)=n^{-1}\sum_{i=1}^n \delta_{z_i}(A)$ with $z_{1:n}\sim \mu_\theta^{(n)}$, denote the simulated empirical measure. Under weak conditions on the observed data and model class $\mathcal{M}^{(n)}$, the sequence $\{\hat{\mu}_{\theta,n}\}_{n\geq 1}$ of random probability measures, on $\mathcal{Y}^n$, will converge (weakly) toward $\mu^{(n)}_\theta\in\mathcal{P}(\mathcal{Y}^n)$, and, likewise, the sequence $\{\hat{\mu}_n\}_{n\geq1}$ will converge toward $\mu^{(n)}_\star\in \mathcal{P}(\mathcal{Y}^n)$. Even if the likelihood is intractable, $\hat\mu_n$ and $\hat\mu_{\theta,n}$ can be easily constructed.

Standard Approximate Bayesian Computation

The following algorithm describes the general approach to approximate Bayesian computation. Let $\epsilon>0$ denote some pre-defined threshold, and let $\mathcal{D}$ denotes an arbitrary distance measure between two data sets $y_{1:n}$ and $z_{1:n}$. The most

Accepted draws from the above can be shown to from the ABC posterior distribution: for $A\subset\Theta$,

where $\mathds{1}$ is the indicator function.

In most cases, using a distance calculated on the data directly is infeasible. The most common approach is to first degraded the data down to a low-dimensional vector of summary statistics and replace $\mathcal{D}(y_{1:n},z_{1:n})$ with a distance based on the summaries. For $\eta:\mathcal{Y}^n\rightarrow\mathcal{B}\subset\mathbb{R}^{d_\eta}$, with $d_\eta$ much smaller than $n$, ABC is most commonly implemented using $\mathcal{D}(y_{1:n},z_{1:n}) = \|\eta(y_{1:n})-\eta(z_{1:n})\|$, and where $\|\cdot\|$ denotes the Euclidean norm.

While practically useful, replacing the data with summaries creates several issues. Firstly, by the nature of the models to which ABC is generally applied, no set of sufficient statistics exist and the use of summaries inevitably leads to a loss of information. Consequently, as $\epsilon\rightarrow0$ the ABC posterior converges to the partial posterior, i.e., the posterior of $\theta$ given $\eta(y_{1:n})$, and, as shown in LF2016 and FMRR2016, point estimators for $\theta$ obtained using ABC based on summary statistics will always have a larger variance than those obtained by exact Bayes. Thirdly, it is well-known that in certain situations the use of different summaries in ABC can lead to significantly different results.

To circumvent the loss of information incurred by using summary statistics in ABC, Bernton2017 propose Wasserstein-based ABC (WABC): replace $\mathcal{D}(y_{1:n},z_{1:n})$ with the $p$-Wasserstein distance

where $\mathcal{S}_n$ is the set of permutations of $\{1, \dots,n\}$. In the case where the data is scalar valued, and if $\gamma(y_i,z_j) = |y_i - z_j|$, the infimum is achieved by sorting $y_{1:n}$ and $z_{1:n}$ in increasing order and matching the order statistics. The use of order statistics within ABC is well-founded, see e.g. FP2012, and in the multivariate case Wassertein ABC can be interpreted as generalizing the use of order statistics within ABC to arbitrary dimension in a principled manner.

While undoubtedly novel, practical application of WABC can be hindered by the specialized tools required to manipulate the Wasserstein distance in cases other than $p=1$ and in the multivariate case. In addition, the nature of the Wasserstein distance makes analyzing the theoretical behavior of WABC difficult. In particular, the asymptotic shape of the WABC posterior is currently unknown, as is the behavior, and in particular the asymptotic efficiency of, point estimators obtained by from WABC. Lastly, it is unclear precisely how the deviations from the model affect the resulting inference in WABC.

Regardless of the above, the ability of WABC to use the entire data has been shown to lead to strong efficiency gains over summary based ABC. Herein, we propose an alternative to WABC that obviates the need to calculate summaries of the data. However, unlike WABC, or ABC based on summary statistics, this new approach can deliver efficient inference in the case of correct model specification, and is simultaneously robust to deviations from the underlying modeling assumptions.

Minimum Distance Based ABC Inference

Replacing the dataset with a low-dimensional summary statistic inevitable leads to an information loss, and a significant reduction in the efficiency of the resulting inference. Moreover, as discussed in frazier2020model, different choices of summary statistics can lead to substantially different behavior when the assumed model is not correct, i.e., inference based on summary statistics are not “stable” under model misspecification. With these two features in mind, efficiency and stability of the resulting inferences, we propose a new approach to Bayesian inference that does not resort to a low-dimensional summary statistics but is instead based on distances between observed and simulated empirical measures. Let $\|\cdot\|_\mathcal{P}:\mathcal{P}(\mathcal{Y})\times \mathcal{P}(\mathcal{Y})\rightarrow\mathbb{R}_{+}$ denote a pseudo-norm on $\mathcal{P}(\mathcal{Y})$. Replacing $\mathcal{D}$ by $\|\hat{\mu}_{\theta,n}-\hat{\mu}_{n}\|_\mathcal{P}$ in equation (ref) yields the ABC posterior

It is important to point out that MD-ABC is not limited to synthetic datasets that are of the same dimension as the observed data. Indeed, there is no reason we could not generate synthetic data $z_{1:m}$, with $m\geq n$, and with corresponding empirical measure $\hat\mu_{\theta,m}$. Such a mechanism is potentially useful to alleviate simulation noise associated with $\hat\mu_{\theta,n}$. When such a distinction between the observed and simulate size is necessary, we will write the simulated empirical measure as $\hat\mu_{\theta,m}$.

While there are many possible norms to choose from, and thus many different possible posteriors on which we can base our inferences, for reasons or efficiency and robustness we restrict our attention to classes of norms that yield “minimum distance functionals”. Optimization of such norms is known to yield point estimators that are inherently robust and potentially efficient (see, e.g., donoho1988automatic). Arguably, the two most commonly used “robust” distances for parameter inference are the Hellinger and Cramer-von Mises (hereafter, CvM) distances.

Example: Mixture Model

We now consider a simplistic example that allows us to easily compare MD-ABC, WABC and exact Bayes (hereafter, EB) inference. We observe a sequence $y_{1:n}=(y_1,\dots,y_n)$ generated independent and identically distributed (iid) from the following two-component mixture model: for $\omega\in[0,1]$, and $\varphi(\cdot;\mu,\sigma^2)$ denoting a normal density with mean $\mu$ and variance $\sigma^2$, $$ f_\theta(\cdot):=(1-\omega)\varphi(\cdot;\mu,\sigma^2_1)+\omega\varphi(\cdot;-\mu,\sigma^2_2),\;\;\theta:=(\mu,\omega,\sigma_1,\sigma_2)'. $$ From this model, we generate $n=100$ observations where the unknown parameters $\theta$ are fixed at $\theta=(-2,0.5,1,1)'$. Our prior beliefs over $\theta$ are $$ \mu\sim\mathcal{N}(0,1),\;\omega\sim\mathcal{U}[0,1],\;\sigma_1\sim\mathcal{U}[0,10],\;\sigma_2\sim\mathcal{U}[0,10]. $$

In this experiment, it is simple to construct the likelihood and sample from the exact posterior. To this end, EB inference is implemented using the NUTS sampler via Rstan (carpenter2017stan). {The NUTS sample was run for 10,000 iterations with an initial warmup period of 1,000 iterations.}

We compare EB with ABC based on the Wasserstein distance (WABC), and MD-ABC using the CvM distance (CvM-ABC) and the Hellinger distance (H-ABC). We set the number of simulated sample points for CvM-ABC and WABC to be $n$, i.e., we take $m=n$, while for H-ABC we use $m=2\cdot n$. All ABC posteriors are sampled using the same approach described in the previous examples.

The resulting posteriors for EB, WABC, and MD-ABC are given in Figure (ref). All posteriors are centered over, or close to, the true values used to generated the data, and there is remarkable agreement across the posteriors. Figure (ref) compares the posterior means of the four different approaches across fifty replicated data sets (generated using the same specification), and demonstrates that, in repeated samples, there is substantial agreement between the point estimators from these different methods.

Lastly, Table (ref) compares the posterior mean, (Means), standard deviation (Std), Monte Carlo coverage (Cov), and root mean squared error (RMSE) for the posterior mean across the replications.\footnote{Monte Carlo coverage is calculated as the percentage of the replications where the true value $\theta_\star$ is contained in the marginal 95% credible set.} The results illustrate that all methods generally give reliable point estimates of the unknown parameters, with H-ABC displaying the smallest bias and RMSE across all parameters. In terms of posterior uncertainty, WABC displays the largest posterior uncertainty, as measured by the posterior standard deviation, while H-ABC displays the smallest amount of posterior uncertainty, followed closely by EB. In terms of Monte Carlo coverage, we note that H-ABC is arguably closest to the nominal 95% level across the different parameters, while WABC yields 100% coverage for all parameters except $\pi$.

The repeated sampling results demonstrate that MD-ABC can produce inferences that are competitive with those obtained using EB, and also yields point estimators that outperform WABC according to the chosen loss measures.

Theoretical properties

In this section, we analyze the asymptotic behavior of MD-ABC: we prove posterior concentration, a Bernstein von-Mises result, and document the behavior of the MD-ABC posterior mean. In addition, we rigorously analyze the robustness of MD-ABC using a local perturbation scheme. The results demonstrate that MD-ABC can deliver point estimators that are as efficient as those obtained by exact Bayes, for certain choices of $\|\cdot\|_{\mathcal{P}}$, and is simultaneously robust to deviations from the underlying modeling assumptions.

Before rigorously analyzing the behavior of this new ABC approach, we first define the following additional notations that will be used throughout the remainder of the paper. For sequences $\{a_{n}\}$ and $\{b_{n}\}$, real valued, $a_{n}\lesssim b_{n}$ denotes $a_{n}\leq Cb_{n}$ for some $C>0$, $a_{n}\asymp b_{n}$ denotes equivalent order of magnitude, $a_{n}{\gg }b_{n}$ indicates a larger order of magnitude and the symbols $o_{p}(a_{n}),\;O_{p}(a_{n})$ have their usual meaning. We reserve $\Rightarrow $ to denote weak convergence under some stated measure. Whenever $\mu_\star$ and $\mu_\theta$ is used, it is assumed that these objects exist. Unless otherwise noted, all limits are taken as $n\rightarrow\infty$.

Posterior Concentration

For a sequence of tolerance values $\{\epsilon_n\}_{n\ge1}$, we first prove that the MD-ABC posterior concentrates onto $ \epsilon_\star:=\inf_{\theta\in\Theta}\|\mu_{\theta}-\mu_\star\|_\mathcal{P}\geq0. $ Since $\epsilon_\star\geq0$, this posterior concentration result will not require that the model be correctly specified. We say the model is correctly specified if there exists some $\theta\in\Theta$, possibly not unique, such that $\mu_{\theta_\star}=\mu_\star$.

We maintain the following assumptions, which are similar to those considered in FMRR2016:

The proof of the following follows the same lines as Theorem 1 of FMRR2016, and Proposition 3.2 in Bernton2017 and hence is omitted for brevity.

The result of Theorem (ref) does not require the parameters to be identifiable, in which case posterior concentrates is onto the set $\{\theta\in\Theta:\epsilon_\star=\|\mu_\theta-\mu_\star\|_\mathcal{P}\}$. However, in many cases the model $\mathcal{M}^{(n)}$ will be identifiable and it would be convenient to state conditions under which the MD-ABC posterior concentrates onto a Dirac mass.

Asymptotic Behavior of the Posterior and the Posterior Mean

Under the case of correct model specification, beran1977minimum demonstrated that the minimum Hellinger distance estimator, the minimizer of $\mathcal{H}\left(\mu_\theta,\hat{\mu}_{n}\right)$, was asymptotically equivalent to the MLE under correct model specification. In addition, ozturk1997generalised have shown that the minimum CvM distance estimator is asymptotically normal and can be nearly as efficient as MLE. In this section, we extend such results to the Bayesian paradigm in the case of intractable likelihoods. In particular, we demonstrate that under correct model specification: 1) the MD-ABC posterior is asymptotically Gaussian ; 2) the MD-ABC posterior mean is asymptotically equivalent to the minimizer of the map $\theta\mapsto \|\mu_{\theta}-\hat{\mu}_{n}\|_{\mathcal{P}}$. The latter result directly implies that, in the class of regular models, MD-ABC based on the Hellinger distance delivers point estimators that are asymptotically as efficient as the MLE .

In order to compare the efficiency of MD-ABC with other approaches, we explicitly consider the case of correctly specified models. In addition, we consider that MD-ABC is implemented with $m=n$ simulations. In what follows consider that $\|\cdot\|_\mathcal{P}$ can be associated to the inner product $\langle\cdot,\cdot\rangle$.

Throughout the remainder, we employ the following explicit assumptions.

Define the following point estimator of $\theta_\star$

where we highlight that $\mu_\theta$ in (ref) is not estimated via simulation but is known. As such, in many models, there is no reason to suspect that $\hat\theta_n$ is feasible to obtain in practice.

Robustness

Recall that the dual objective of MD-ABC was to deduce a method that was more efficient than summary-statistic based ABC, but which simultaneously displays meaningful robustness to deviations from the underlying modeling assumption. Corollary (ref) demonstrates that MD-ABC satisfies the first of these tenants, and in this section we formally demonstrate that MD-ABC point estimators display robustness to perturbations from the assumed model. Throughout the remainder, we refer to such a setting as local model misspecification.\footnote{We refer the reader to frazier2020model for an analysis of ABC methods under global model misspecification.}

Let $T=\{T_n\}_{n\ge1}$, $T_n:\mathcal{M}^{(n)}\rightarrow\Theta$ denote a functional of interest. For all $Q\in\mathcal{M}^{(n)}$, define the map $Q\mapsto \bar{T}(Q)$ by $$ \bar{T}(Q):=\operatornamewithlimits{argmin\,}_{t\in\Theta} \int \frac{(\sqrt{n}\{t-\theta\})\mathrm{d}\Pi(\theta)\int_{\mathcal{Y}^n}\mathds{1}\left[\|\hat{\mu}_{\theta,n}-Q\|_\mathcal{P}\leq\epsilon\right]\mathrm{d}\mu_\theta^{(n)}(z_{1:n})}{\int_{{\Theta}}\mathrm{d}\Pi(\theta)\int_{\mathcal{Y}^n}\mathds{1}\left[\|\hat{\mu}_{\theta,n}-Q\|_\mathcal{P}\leq\epsilon\right]\mathrm{d}\mu_\theta^{(n)}(z_{1:n})}. $$ The MD-ABC posterior mean corresponds to $\bar{\theta}_n=\bar{T}_{}(\hat{\mu}_n)$.

We now analyze the behavior of the map $Q\mapsto \bar{T}_{}(Q)$ as we vary the true process that has generated the data in a neighborhood of the assumed model. To this end, following, among others, rieder2012robust, we assume that the observed data is generated from a member in the class of simply perturbation from the assumed model $\mu^{(n)}_\theta$ along the tangent $\zeta_n$: for and $t\in\Theta$, $\theta_n:=\theta_\star+t/\sqrt{n}$, observed data is generated according to\footnote{Similar to the analysis in Section (ref), in this section we only consider the specific case where $v_n=v_m=\sqrt{n}$ to simplify the analysis.} $$ \mathrm{d} \mu^{(n)}_{\theta_n,\zeta_n}=\left(1+\frac{t^{\intercal}\zeta_n}{\sqrt{n}}\right)\mathrm{d}\mu^{(n)}_{\theta_\star},\quad\lim _{n \rightarrow\infty} \int(\zeta_n-\zeta)^2\mathrm{d}\mu_\theta=0, \quad \left|\zeta_{n}\right|=o(\sqrt{n}). $$It is important to note that the form of model misspecification induced by the simple perturbations $\mu_{\theta_n,\zeta_n}$ excludes the case of global model misspecification analyzed in frazier2020model. Unlike the clase of global model misspecification, the “modelling bias” associated with the simple perturbations eventually converges to zero.

Define the neighborhood of radius $r/\sqrt{n}$ around $\mu_\star$, based on the distance $d(\cdot,\cdot)$, as follows: $$\mathcal{B}_{d}\left(\mu_\star, r / \sqrt{n}\right)=\left\{\mu \in \mathcal{M}^{}: d\left(\mu^{}, \mu_\star^{}\right) \leq r / v_n\right\}.$$ Implicitly, the neighborhood $\mathcal{B}_{d}\left(\mu_\star, r / \sqrt{n}\right) $ sets the classes of models under which we will analyze the behavior of functionals $Q\mapsto T(Q)$ and statistics obtained from $T(\cdot)$. Optimality of these functionals is defined relative to a given neighborhood system in which we measure the deviations from the model.

Herein, we consider two common types of neighborhoods: Hellinger and Cramer-von Mises (CvM) neighborhoods. The Hellinger neighborhood replaces the generic distance $d$ in $\mathcal{B}_d$ with the Hellinger distance, see Definition (ref), and we denote the neighborhood by $\mathcal{B}_{\mathcal{H}}$. To define the CvM neighborhood, let $\lambda$ be a $\sigma$-finite measure, and let $Q(y)$ and $P(y)$ denote the CDF of $Q$ and $P$ at the point $y\in\mathcal{Y}$. Then, the CvM neighbourhood is defined as $$ \mathcal{B}_{\mathcal{C}}\left(\mu_\star, r / \sqrt{n}\right),\text{ where } \mathcal{C}(Q,P):=\int_{\mathcal{Y}} \left|Q(y)-P(y)\right|\lambda(\mathrm{d} y). $$

We now analyze the robustness of MD-ABC to deviations from the model that lie in $\mathcal{B}_{d}$, for $d\in\{\mathcal{H},\mathcal{C}\}$, and compare this behavior to any other functional $\mu\mapsto T_\eta(\mu)$ computed from summary statistics $\eta$. We carry out this analysis in the case where the functionals are regular and Fisher consistent. Let $T_a\in\{\bar{T},T_\eta\}$, $T_a:\mathcal{M}\mapsto\Theta$ be an estimator of $\theta_\star$.

Let $\tau:\Theta\mapsto\mathbb{R}$ be a (possibly) nonlinear transformation of the parameter with a bounded first derivative.\footnote{The use of $\tau$ simplifies the calculation of bias and variance terms in the following calculations.} The following result describes the robustness of MD-ABC functionals to deviations from the underlying modelling assumptions.

Additional Examples: Robustness

The results in Section (ref) suggest that the MD-ABC approaches have useful robustness properties that may not be replicated by alternative inference approaches. In this section, we investigate this claim by comparing the robustness of MD-ABC, WABC and exact Bayes to deviations from the underlying modeling assumptions.

Example: Mixture Model Revisited

We now revisit the mixture model to analyze the robustness of the various inference approaches in this context. The assumed data generating process (DGP) for $z_1,\dots,z_n$ is independent and identically distributed (iid) from a two-component mixture: for $\omega\in[0,1]$, recalling that $\varphi(\cdot;\mu,\sigma^2)$ denotes a normal density with mean $\mu$ and variance $\sigma^2$, $$ f_\theta(\cdot):=(1-\omega)\varphi(\cdot;\mu,\sigma^2_1)+\omega\varphi(\cdot;-\mu,\sigma^2_2),\;\;\theta:=(\mu,\omega,\sigma^2_1,\sigma^2_2)'. $$Following cutler1996minimum, we examine the situation where the observed data $y_1,\dots,y_n$ is generated iid from the contaminated normal mixture: for $\alpha\in[0,1]$, $$ f_\star(\cdot):=(1-\alpha)f_\theta(\cdot) + \alpha \varphi(\cdot;\zeta,\nu). $$The value $\alpha$ determines the level of contamination in the observed data and for $\alpha>0$ the assumed DGP is misspecified. The impact of this assumption on the resulting inferences can be controlled by choosing $\alpha,\zeta,\nu$.

In this example we compare the robustness of exact Bayesian (EB) inference, ABC based on the Wasserstein distance (WABC), and MD-ABC based on the CvM and Hellinger distances. We fix $\alpha=0.05$, $\nu=0.01$, and we consider a grid of values for $\zeta\in[-10,10]$. The parameters $\theta$ are fixed at $\theta_\star=(-2,0.5,1,1)'$. The observed sample size is again $n=100$, and we set the number of simulated sample points for CvM-ABC to be $n$, i.e., we take $m=n$, and for H-ABC we set $m=2\cdot n$.

Across all experiments, EB inference is implemented using the NUTS sampler: {NUTS was run for 10,000 iterations with an initial warm-up period of 1,000 iterations.} For all ABC approaches, we again target the posterior using the sequential Monte Carlo (SMC) approach of del2012adaptive with $N=1,024$ particles, and with all procedures using $2\cdot10^5$ simulations from the model.

First, we analyze the robustness properties of these approaches by comparing the posteriors when the observed data is simulated using a particularly large value of $\zeta$, namely $\zeta=9$. Figure (ref) displays the posteriors for the four procedures in this example. In this case, we see that the ABC posteriors are similar to those obtained under correct specification (see Figure (ref)) and signal that these methods are somewhat insensitive to the contamination in the data. In contrast, the contamination leads to significantly different outcomes for EB: for all parameters except $\sigma_2$, the EB posteriors have fatter tails and are farther away from the true values used to generate the data than for the ABC posteriors.

In Figure (ref) we graphically compare the evolution of the posteriors means for these four approaches across a grid of values $\zeta\in[-10,10]$, with incitements of $0.50$. The results suggest that EB can be highly sensitive to the data contamination, while MD-ABC and WABC display significantly less sensitivity to model misspecification. Across the different values of $\zeta$, and across the different methods, MD-ABC based on the CvM distance displays the least amount of sensitivity to model misspecification across the parameters.

Figure (ref) displays the marginal posterior means of the four inference approaches across fifty data sets generated from the above design with $\zeta=9$. Collectively, the results in Figures (ref)-(ref) demonstrate that EB is the most sensitive of the methods to this particular form of model misspecification, while MD-ABC based on the CvM distance displays the least amount of sensitive to misspecification, followed closely by MD-ABC based on the Hellinger distance.

Table (ref) compares the posterior means (Means), standard deviation (Std), root mean squared error (RMSE) of the posterior mean and the Monte Carlo coverage (COV) from the replications associated with the results in Figure (ref). As suggested by Figures (ref) and (ref), MD-ABC based on the CvM and Hellinger distances displays the least amount of sensitivity to model misspecification, while EB is the most sensitive. Indeed, the posterior means of EB are generally far away from $\theta_\star$, and its standard deviation is generally quite large. Moreover, due to the large bias, EB displays very poor coverage for $\sigma_1$, with the 95% credible set covering the true value in only 2% of the repeated samples. In contrast to EB, the marginal MD-ABC posteriors are significantly less sensitive to the model misspecification, with MD-ABC displaying better behavior than EB or WABC across virtually all the accuracy measures.

The repeated sampling results demonstrate that while EB is the gold standard under correct model specification, robust Bayesian approaches should be entertained when model misspecification is a possibility. The robustness of MD-ABC in this example, relative to EB and WABC, is a direct consequence of the methods focus on distance metrics which are known to yield inferences that are less sensitive to perturbations from the underlying modeling assumptions.

Example: Stochastic Volatility

In this example, we analyze the stability of MD-ABC inference in a time series context where data is generated according to an autoregressive log-normal stochastic volatility model. The observed series $y_{t}$ is generated according to $$

$$where $v_t,u_t$ are, for simplicity, $\mathcal{N}(0,1)$ iid random variables. The unknown parameters $\theta=(\pi,\beta,\sigma_u)'$ have the following prior specification $$ \pi\sim\mathcal{N}(0,1),\delta\sim\mathcal{U}[0,1],\sigma_u\sim\mathcal{U}[0,5]. $$

To analyze the stability of MD-ABC to local deviations from the SV model, we consider a sequence of models where the error process $v_t$ is generated from the normal mixture: $$ f_{v_t}(\cdot):=(1-\alpha)\varphi(\cdot;0,1)+\alpha\varphi(\cdot;\zeta,.001),\text{ for }\alpha\in[0,1]. $$This particular error process emulates the “jumps” that are observed in financial return series, where the parameter $\alpha$ proxies as “jump intensity/frequency” and the parameter $\zeta$ proxies for “jump magnitude”. We generate realizations from the above process with $\alpha=0.05$, and where $\zeta$ varies on a uniform grid from -5 to 0, with increments of 0.50. Following the simulation design of jacquier2002bayesian, we fix $\theta=(-0.736,0.900,0.363)'$. The sample size is set to $n$ = 500 across all experiments.

To assure convergence to the stationary distribution, the return sequence is generated using a sample size of 1000, where the first 500 observations are discarded. The same procedure is carried out when simulating samples for ABC.

In this example, we compare the posterior means and standard deviations of MD-ABC, based on the CvM and Hellinger distance, against those from the exact Bayes (EB), and WABC. The exact posterior is estimated using the NUTS sampler (with the same specification as in Example (ref)). For MD-ABC based on the Hellinger distance, we use $m=10\cdot n$ simulated observations, and for MD-ABC based on the CvM distance, we use $m=n$. All ABC approaches used the sequential Monte Carlo (SMC) approach of del2012adaptive with $N=1,024$ particles, and $2\cdot10^5$ simulations from the model.

Figures (ref) and (ref) present the results for the posterior means and standard deviations, respectively, across the difference values of $\zeta$. The results demonstrate that there is reasonable agreement across all methods when $\zeta$ is small, and indicates that MD-ABC can be reliably applied to time-series settings. However, for smaller values of $\zeta$, i.e., more negative values, the posterior means from EB are severely affected by the jump components and depart significantly from the true values used to generate the data. In addition, the posterior standard deviations for EB fluctuate significantly over the grid of values for $\zeta$. In contrast, MD-ABC point estimates of $\pi$ and $\beta$, as well as their posterior standard deviations, are virtually unaffected by the jumps in the data; inference on $\sigma$ is also less sensitive but does display some sensitivity. These results reinforce that the fact that even under relatively mild model misspecification, the reliability of EB-based inference can rapidly degraded; even when only five percent of the observed data does not agree with the assumed model, as in this example.

Discussion

We have proposed a new approach to Bayesian inference, MD-ABC, based on comparing differences between empirical probability measures calculated from observed and simulated data. This method is particularly useful in models for which the calculation of likelihood functions is intractable. For particular choices of the norm, such as when the norm is chosen to be the Hellinger distance, MD-ABC point estimators are asymptotically efficient and, thus, are equivalent to those obtained from exact Bayesian methods (under the assumption that the model is correctly specified).

While originally motivated by situations where the likelihood is intractable, we have demonstrated that MD-ABC methods are a useful alternative to exact Bayesian methods if one is concerned with deviations from the underlying modeling assumptions. Moreover, we have shown that if the model is correctly specified the differences between MD-ABC and exact inference will generally be minor (depending on the choice of distance used in MD-ABC). As such, MD-ABC represents an approach to inference that is as efficient, or nearly so, as exact methods, but is simultaneously robust to certain types of deviations from the underlying modeling assumptions.

\noindentAcknowledgments: David T. Frazier gratefully acknowledges support by the Australian Research Council through grant DE200101070. We would also like to thank Professor Pierre Jacob for making the code needed to implement Wasserstein ABC publicly available.