[description=True data generating
distribution]P$\operatorname{\operatorname{\mathbb{P}}^{\ast}}$
[description=Observations
$x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}$]I$x_{1:n}$
[description=Empirical measure
$\mathbb{P}_{n}=\frac{1}{n}\sum_{i=1}^{n}\delta_{x_{i}}$]F$\mathbb{P}_{n}$
[description=Centering measure of the Dirichlet Process prior
$\operatorname{\operatorname{DP}}(\alpha,\operatorname{\mathbb{F}})$]G$\operatorname{\mathbb{F}}$
[description=Centering measure of the Dirichlet Process posterior
$\operatorname{\operatorname{DP}}(\alpha^{\prime},\operatorname{\mathbb{F}}^{\prime})$
defined as
$\operatorname{\mathbb{F}}^{\prime}=\frac{\alpha}{\alpha+n}\operatorname{\mathbb{F}}+\frac{n}{\alpha+n}\mathbb{P}_{n}$]H$\operatorname{\mathbb{F}}^{\prime}$
[description=Probability measure on $\mathcal{P}$ defined by the sampling
procedure of samples from the Dirichlet approximation of the DP posterior
$\operatorname{\operatorname{DP}}(\alpha^{\prime},\operatorname{\mathbb{F}}^{\prime})$
in (2.3)]A$\operatorname{\hat{\nu}}$ [description=Probability measure on
$\mathcal{P}$ defined by the stick-break process representing samples from the
exact DP posterior $\operatorname{\mathbb{P}}\given
x_{1:n}\sim\operatorname{\operatorname{DP}}(\alpha^{\prime},\operatorname{\mathbb{F}}^{\prime})$]B$\operatorname{\nu}$
[description=Data space]D$\mathcal{X}$ [description=Space of Borel
distributions on $\mathcal{X}$]E$\operatorname{\mathcal{P}}$
[description=Probability measure
$\operatorname{\mathbb{P}_{\theta}}\in\mathcal{P}$ indexed by parameter
$\theta$ ]C$\operatorname{\mathbb{P}_{\theta}}$ [description=Family of
parametric models
$\operatorname{\mathcal{P}_{\Theta}}=\\{\operatorname{\mathbb{P}}_{\theta}:\theta\in\Theta\\}\subset\operatorname{\mathcal{P}}$]X$\operatorname{\mathcal{P}_{\Theta}}$
[description=Map $\theta^{*}_{l}:\mathcal{P}\rightarrow\Theta$ indexed by loss
$l$ which takes as input a probability measure $\operatorname{\mathbb{P}}$ and
returns the value of $\theta\in\Theta$ which minimises the expected loss; i.e.
$\theta^{*}_{l}(\operatorname{\mathbb{P}})=\operatorname{arginf}_{\theta\in\Theta}\operatorname{\mathbb{E}}_{X\sim\operatorname{\mathbb{P}}}[l(X,\theta)]$
]K$\theta^{*}_{l}$ [description=Map
$\operatorname{\theta^{*}}:\mathcal{P}\rightarrow\Theta$ for the MMD-based
loss which takes as input a probability measure $\operatorname{\mathbb{P}}$
and returns the value of $\theta\in\Theta$ which minimises the MMD between
$\operatorname{\mathbb{P}_{\theta}}$ and $\operatorname{\mathbb{P}}$; i.e.
$\operatorname{\theta^{*}}(\operatorname{\mathbb{P}})=\operatorname{arginf}_{\theta\in\Theta}\operatorname{\operatorname{MMD}}^{2}(\operatorname{\mathbb{P}},\operatorname{\mathbb{P}_{\theta}})$
]L$\operatorname{\theta^{*}}$

Robust Bayesian Inference for Simulator-based Models
via the MMD Posterior Bootstrap

Charita Dellaporta Jeremias Knoblauch Theodoros Damoulas François-Xavier Briol

University of Warwick University College London University of Warwick
University College London

###### Abstract

Simulator-based models are models for which the likelihood is intractable but
simulation of synthetic data is possible. They are often used to describe
complex real-world phenomena, and as such can often be misspecified in
practice. Unfortunately, existing Bayesian approaches for simulators are known
to perform poorly in those cases. In this paper, we propose a novel algorithm
based on the posterior bootstrap and maximum mean discrepancy estimators. This
leads to a highly-parallelisable Bayesian inference algorithm with strong
robustness properties. This is demonstrated through an in-depth theoretical
study which includes generalisation bounds and proofs of frequentist
consistency and robustness of our posterior. The approach is then assessed on
a range of examples including a g-and-k distribution and a toggle-switch
model.

##  1 INTRODUCTION

A key assumption of standard Bayesian inference is that the data-generating
mechanism lies within the family of models defined by our choice of likelihood
[9]. In reality, this assumption is often violated, and this can lead to
inconsistent and misleading inference outcomes; see for example [34, 55] for
the cases of linear or nonparametric regression.

The need for robust inference methods is even more apparent when dealing with
simulator-based models. In these models, independent sampling is possible, but
evaluating the likelihood is intractable. Simulators are used in a range of
fields including population genetics [67], ecology [5] and telecommunication
engineering [12]; see [22] for a recent review. In each of these examples, the
model is at best a rough approximation of a complex physical or biological
phenomenon, and will most likely not capture all of the key characteristics of
the underlying data generating process. Unfortunately, the main approaches to
Bayesian inference for simulators, including approximate Bayesian computation
(ABC) [7], tend to approximate the classical Bayesian posterior, which itself
lacks robustness to misspecification [30]. This leads ABC posteriors to
exhibit poor coverage; see [37] for an empirical study.

One promising approach for addressing model misspecification is generalised
Bayesian inference (GBI; e.g. [[]]hooker2014bayesian, ghosh2016robust,
bissiri2016general,Jewson2018,Knoblauch2019) and the closely-related
approaches of quasi-Bayesian inference [[, e.g.]]Chernozhukov2003, PAC-
Bayesian inference [[, e.g.]]ShaweTaylor,catoniPACBound, Safe-Bayes [35],
Split-BOLFI approach for higher-dimensional inference [[]]thomas2020split and
the Focused Bayesian Prediction approach [[]]loaiza2021focused. These
approaches typically work directly with the original model, but attempt to
correct for possible misspecification by scoring observations robustly. Recent
work has focused on GBI for intractable likelihood models [69, 56, 75, 51].
While some of these techniques address robustness concerns in simulator-based
methods, they have two main drawbacks: Firstly, they require setting a crucial
hyperparameter that determines the relative importance of the prior for the
posterior inferences. Secondly, some make use of Monte Carlo methods, which
imposes a substantial computational burden.

Another recent approach for inference under misspecification is Bayesian
nonparametric learning (NPL) [[, see]]lyddon2018nonparametric,
lyddon2019general, fong2019scalable. Unlike GBI, NPL does not address
misspecification by robustly scoring the statistical model. Instead,
robustness is achieved by obtaining a non-parametric posterior directly on the
data-generating process. This posterior then implies a posterior on the
parameter of interest through the use of a robust loss function.

Our paper’s contribution is the first NPL-based algorithm for simulators.
Specifically, we leverage the NPL framework to obtain a posterior belief
distribution about the parameters that minimise the maximum mean discrepancy
(MMD) [33] between our model and the data-generating mechanism. The MMD is a
probability metric that is not only robust, but also easy to approximate
through simulation—and therefore suitable for simulators. Further, the MMD has
numerous desirable theoretical robustness and generalisation properties [[,
see]]briol2019statistical, cherief2019finite, cherief2020mmd. One of the main
achievements of this paper is to show that these hold for our method whenever
we use a bounded kernel. Additionally, unlike ABC or GBI, our approach is
computationally efficient: it is trivially parallelisable, and never requires
discarding parameter samples or using inherently sequential Monte Carlo
methods.

The paper is structured as follows. In Section 2, we recall the details of NPL
and introduce MMD estimators. In Section 3, we propose our novel algorithm
which combines these concepts to create a scalable approach for robust
Bayesian inference with simulators. Then, Section 4 provides theoretical
results including consistency and robustness. Finally, Section 5 studies the
algorithm on a range of benchmark problems for simulators including the
g-and-k distribution, and a toggle-switch model describing the interaction of
genes through time.

##  2 BACKGROUND

Let $\mathcal{X}$ denote our data space, $\operatorname{\mathcal{P}}$ the
space of Borel distributions on $\mathcal{X}$, and
$\operatorname{\operatorname{\mathbb{P}}^{\ast}}\in\operatorname{\mathcal{P}}$
the true data-generating mechanism of our data. In Bayesian statistics, given
observations
$x_{1:n}=x_{1},\ldots,x_{n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}$,
one chooses a parametric model
$\operatorname{\mathcal{P}_{\Theta}}=\\{\operatorname{\mathbb{P}}_{\theta}:\theta\in\Theta\\}\subset\operatorname{\mathcal{P}}$.
Given $\operatorname{\mathcal{P}_{\Theta}}$, the Bayesian now places a prior
on the parameter $\theta$, and then conditions on $x_{1:n}$ to obtain a
posterior on $\theta$. For standard Bayesian inference to be well-behaved, we
have to assume that $\operatorname{\mathcal{P}_{\Theta}}$ is well-specified;
i.e $\exists{\theta_{0}}\in\Theta$ such that
$\operatorname{\mathbb{P}}_{\theta_{0}}=\operatorname{\operatorname{\mathbb{P}}^{\ast}}$.
When this assumption is violated, we call the model misspecified.

Misspecification in the Bayesian context has recently seen increasing interest
since the Bayesian posterior does not provide robust parameter inferences [[,
see]]safeBayes. This has led to GBI and NPL approaches aimed at rectifying the
issue. Another approach, albeit not specifically for simulators, is the
BayesBag algorithm [[]]huggins2019robust. In this paper, we combine the
strengths of GBI and NPL approaches for robust inference with simulators.

###  2.1 Simulator-Based Inference

The problem of simulator-based models is a significant challenge for
Bayesians. Consider some
$\operatorname{\mathbb{P}}_{\theta}\in\operatorname{\mathcal{P}_{\Theta}}$
with fixed $\theta\in\Theta$, whose density is the likelihood
$p(\cdot\given\theta)$ and suppose we have a prior $\pi(\theta)$ on the
parameter. The corresponding posterior density is given by

$\textstyle\pi(\theta\given
x_{1:n})=\frac{\prod_{i=1}^{n}p(x_{i}\given\theta)\pi(\theta)}{\int_{\Theta}\prod_{i=1}^{n}p(x_{i}\given\theta)\pi(\theta)d\theta}\propto\prod_{i=1}^{n}p(x_{i}\given\theta)\pi(\theta).$

The latter is intractable whenever the likelihood $p(\cdot\given\theta)$
cannot be evaluated pointwise. This has led to the development of simulator-
based inference methods (sometimes also called likelihood-free inference
methods).

In this paper, we focus on simulator-based models (also called generative
models), which are parametric families. This means that
$\operatorname{\mathbb{P}}_{\theta}$ can be represented using a distribution
$\operatorname{\mathbb{U}}$ on a space $\operatorname{\mathcal{U}}$ and a
simulator $G_{\theta}:\operatorname{\mathcal{U}}\rightarrow\mathcal{X}$ so
that a sample $y\sim\operatorname{\mathbb{P}}_{\theta}$ from the model can be
obtained by first sampling $u\sim\operatorname{\mathbb{U}}$, and then applying
the simulator $y:=G_{\theta}(u)\in\mathcal{X}$.

##### Approximate Bayesian Computation

ABC algorithms are arguably the most popular family of techniques for tackling
Bayesian posteriors of simulator-based models [[, see
]]beaumont2002approximate, sisson2018abc, beaumont2019approximate. Most ABC
algorithms are a variation of the following steps:

  1. (i)

For $b=1,2,...B$ and prior $\pi$, sample
$\theta_{b}\operatorname{\overset{iid}{\sim}}\pi$;

  2. (ii)

For each $\theta_{b}$, sample $m$ realisations from
$\operatorname{\mathbb{P}}_{\theta_{b}}$ (i.e. simulate
$u_{1:m}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{U}}$ and set
$y^{(b)}_{i}=G_{\theta_{b}}(u_{i})$ $\forall i$);

  3. (iii)

Compare $y^{(b)}_{1:m}$ with the true data $x_{1:n}$ using a discrepancy
$D:\operatorname{\mathcal{P}}\times\operatorname{\mathcal{P}}\rightarrow\operatorname{\mathbb{R}}$:
$D(\operatorname{\mathbb{P}}_{n},\widehat{\operatorname{\mathbb{P}}}_{\theta_{b}})$
where $\operatorname{\mathbb{P}}_{n}=\frac{1}{n}\sum_{i=1}^{n}\delta_{x_{i}}$
and
$\widehat{\operatorname{\mathbb{P}}}_{\theta_{b}}=\frac{1}{m}\sum_{j=1}^{m}\delta_{y^{(b)}_{j}}$.

  4. (iv)

Weight $\theta_{b}$ as an approximate sample from the posterior according to
$D(\operatorname{\mathbb{P}}_{n},\widehat{\operatorname{\mathbb{P}}}_{\theta_{b}})$.

The function $D$ is often chosen to be a discrepancy comparing summary
statistics of the two datasets, but could also be a probability metric. For
example, [11] used the Wasserstein distance, whilst [60] used the MMD. Step
(iv) is usually implemented by verifying whether the discrepancy is smaller
than some threshold $\varepsilon$, accepting the sample $\theta_{b}$ if so,
and rejecting it otherwise. This leads to the following approximation of the
standard Bayesian posterior density:

$\textstyle\pi_{\text{ABC}}(\theta\mid x_{1:n})$
$\textstyle\propto\int_{\mathcal{X}}\ldots\int_{\mathcal{X}}\prod_{j=1}^{m}1_{\\{D(\operatorname{\mathbb{P}}^{\ast}_{n},(\operatorname{\mathbb{P}_{\theta}})_{m})<\varepsilon\\}}(\theta)$
$\textstyle\qquad\qquad\times p(y_{j}|\theta)\pi(\theta)dy_{1}\ldots dy_{m}.$

Smaller values of $\varepsilon$ imply an increased computational cost since
more parameter samples will be rejected. However, as $\varepsilon\rightarrow
0$, the ABC posterior approaches the standard Bayesian posterior
$\pi(\theta|x_{1:n})$. The latter is usually seen as a desirable property of
ABC, but it can lead to poor performance because the Bayesian posterior itself
lacks robustness to misspecification.

Although the prototype ABC algorithm above is parallelisable, it will tend to
be inefficient because many samples will be rejected whenever $\epsilon$ is
small. As a result, it is common to use inherently sequential algorithms such
as population Monte Carlo [4] or sequential Monte Carlo with adaptive
resampling [24].

##### Alternative Approaches

There are a number of perhaps less prominent alternative approaches for
Bayesian inference for simulators, including Bayesian synthetic likelihoods
[65] and techniques relying on neural density estimation [[, see
e.g.]]papamakarios2016fast, papamakarios19a, Cranmer2020frontier. Just like
ABC, all these approaches are non-robust since they approximate the standard
Bayesian posterior.

###  2.2 Generalised Bayesian Inference (GBI)

To address the robustness concern for standard Bayesian inference, GBI
approaches have recently been proposed. In GBI,
$l_{n}:\mathcal{X}^{n}\times\Theta\to\operatorname{\mathbb{R}}$ denotes any
(empirical) loss function and $\beta>0$ a learning rate. Then, the associated
GBI posterior’s density is

$\textstyle\pi_{\text{GBI}}(\theta\given x_{1:n})=\frac{\exp\\{-\beta
l_{n}(x_{1:n},\theta)\\}\pi(\theta)}{\int_{\Theta}\exp\\{-\beta
l_{n}(x_{1:n},\theta)\\}\pi(\theta)d\theta}.$

Here, the loss is typically chosen in relation to the likelihood model
$p(x_{1:n}\given\theta)=\prod_{i=1}^{n}p(x_{i}\given\theta)$. Indeed, it is
easy to see that for the standard Bayesian posterior, one chooses
$l_{n}(x_{1:n},\theta)=-\log p(x_{1:n}|\theta)$. From this, it also becomes
clear that the Bayes posterior’s lack of robustness is intimately related to
choosing the loss function $l_{n}(x_{1:n},\theta)=-\log p(x_{1:n}|\theta)$. To
address this, a host of generalised posteriors have been derived by choosing a
discrepancy $D$ with desirable robustness properties, and then seeking to find
a loss so that $l_{n}(x_{1:n},\theta)\approx
D(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}}_{\theta})$
[[, see]]Jewson2018.

GBI for robustness in simulator-based inference has been studied in [57] and
[69]. The main drawback of the proposals in both papers are two-fold: firstly,
the uncertainty quantification properties of $\pi_{\text{GBI}}(\theta\given
x_{1:n})$ depend on $\beta$, which can only be chosen based on rough
heuristics [[, see]]wu2020comparison. Furthermore, the methods require
techniques that are computationally more cumbersome: [57] use an inherently
sequential sampling method based on pseudo-marginal MCMC [2], and [69] relies
on standard ABC methods.

In this paper, we perform Bayesian inference under model misspecification
through a robust discrepancy like in GBI methods—but without relying on
computationally burdensome methods like in ABC.

###  2.3 Bayesian Nonparametric Learning (NPL)

While the Bayesian nonparametric learning (NPL) framework of [49] was
introduced to deal with misspecification, it also possesses attractive
computational properties. NPL defines a nonparametric prior directly on the
true data-generating process
$\operatorname{\operatorname{\mathbb{P}}^{\ast}}$, which in turn leads to a
nonparametric posterior on $\operatorname{\operatorname{\mathbb{P}}^{\ast}}$.
Any posterior on the parameter space $\Theta$ of
$\operatorname{\mathcal{P}_{\Theta}}$ induced by NPL thus derives from this
nonparametric posterior on $\operatorname{\mathcal{P}}$.

Figure 1: _Sketch of the posterior bootstrap._ Samples from the NPL posterior
on $\Theta$ are obtained by mapping realisations of the posterior
$\operatorname{\hat{\nu}}$ through the map $\operatorname{\theta^{*}}$.

Note that this is different, but closely related to, standard Bayesian
inference where conditioning occurs at the level of parameters as opposed to
the level of the data-generating process.

Following [49] and [28], we use a Dirichlet Process (DP) prior on the data-
generating process:
$\operatorname{\mathbb{P}}\sim\operatorname{\operatorname{DP}}(\alpha,\operatorname{\mathbb{F}})$.
Here, $\alpha>0$ is a concentration parameter and
$\operatorname{\mathbb{F}}\in\operatorname{\mathcal{P}}$ a centering measure.
Given
$x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}$,
it follows by conjugacy that

$\textstyle\begin{split}\textstyle\operatorname{\mathbb{P}}\given
x_{1:n}&\textstyle\sim\operatorname{\operatorname{DP}}\left(\alpha^{\prime},\operatorname{\mathbb{F}}^{\prime}\right),\\\
\textstyle\alpha^{\prime}=\alpha+n,\qquad&\textstyle\operatorname{\mathbb{F}}^{\prime}=\frac{\alpha}{\alpha+n}\operatorname{\mathbb{F}}+\frac{n}{\alpha+n}\mathbb{P}_{n}\end{split}$
(1)

where $\mathbb{P}_{n}:=\frac{1}{n}\sum_{i=1}^{n}\delta_{{x}_{i}}$ and
$\delta_{x}$ is a Dirac measure at $x\in\mathcal{X}$. The size of $\alpha$
quantifies our confidence in the quality of the prior centering measure
$\operatorname{\mathbb{F}}$ and regulates the influence of the prior on the
posterior. Accordingly, the limiting case of $\alpha=0$ corresponds to a non-
informative prior and results in the posterior
$\operatorname{\operatorname{DP}}(n,\mathbb{P}_{n})$. Hence, unlike the
hyperparameter $\beta$ in GBI, the heuristic $\alpha=0$ has a clear
interpretation and yields reliable uncertainty quantification [[,
see]]fong2019scalable, knoblauch2020robust, galvani2021bayesian.

The posterior on the data-generating process is readily translated into a
posterior on a parameter space $\Theta$. Suppose we know the true data-
generating process $\operatorname{\operatorname{\mathbb{P}}^{\ast}}$ and a
loss $l:\mathcal{X}\times\Theta\rightarrow\operatorname{\mathbb{R}}$. Clearly,
having access to $\operatorname{\operatorname{\mathbb{P}}^{\ast}}$ is
equivalent to having access to infinitely many data points. Consequently, no
Bayesian uncertainty is needed, and one can simply compute

$\textstyle\theta^{*}_{l}(\operatorname{\operatorname{\mathbb{P}}^{\ast}}):=\arg\inf_{\theta\in\Theta}\operatorname{\mathbb{E}}_{X\sim\operatorname{\operatorname{\mathbb{P}}^{\ast}}}[l(X,\theta)].$
(2)

In practice, we do not know the true data-generating process. However, we have
a posterior belief over it, so that the simple push-forward measure
$(\theta^{*}_{l})_{\\#}(\operatorname{\operatorname{DP}}(\alpha^{\prime},\operatorname{\mathbb{F}}^{\prime}))$
gives a posterior on $\Theta$ denoted by $\operatorname{\Pi_{\text{NPL}}}$.
This is illustrated in Figure 1.

For the reader who is not familiar with pushforward distributions,
$\operatorname{\Pi_{\text{NPL}}}$ is best understood through the following
sampling mechanism to obtain independent realisations from
$\operatorname{\Pi_{\text{NPL}}}$. At iteration $j$,

  1. 1.

Sample $\operatorname{\mathbb{P}}^{(j)}$ from the posterior DP in (1);

  2. 2.

Compute $\theta^{(j)}=\theta^{*}_{l}(\operatorname{\mathbb{P}}^{(j)})$ using
(2).

This procedure is trivially parallelisable and does not discard any samples.
It therefore overcomes the computational inefficiencies of ABC-based methods,
whose rejection rate for samples is typically quite high—especially as
$\varepsilon$ is moved closer to zero.

Exact sampling from a DP as in step 1 above is usually infeasible. The most
common approximation is the truncated stick-breaking procedure [70]. This
procedure in turn can be approximated by the Dirichlet approximation of the
stick breaking process [52, 41]. In our case, this leads to

$\textstyle\tilde{x}^{(j)}_{1:T}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}$
$\textstyle,\quad(w^{(j)}_{1:n},\tilde{w}^{(j)}_{1:T})\sim\operatorname{\operatorname{Dir}}(1,\dots,1,\frac{\alpha}{T},\dots,\frac{\alpha}{T}).$
$\textstyle\operatorname{\mathbb{P}}^{(j)}$
$\textstyle\approx\sum_{i=1}^{n}w_{i}^{(j)}\delta_{x_{i}}+\sum_{k=1}^{T}\tilde{w}_{k}^{(j)}\delta_{\tilde{x}_{k}^{(j)}}.$
(3)

$\operatorname{\hat{\nu}}$ denotes the probability measure on
$\operatorname{\mathcal{P}}$ defined by (2.3), so that
$\operatorname{\mathbb{P}}=\sum_{i=1}^{n}w^{(j)}\delta_{x_{i}}+\sum_{k=1}^{T}\tilde{w}_{k}^{(j)}\delta_{\tilde{x}_{k}^{(j)}}\sim\operatorname{\hat{\nu}}$.

It is generally not possible to obtain the minimiser in (2) in closed form,
and this objective may not even be convex. This necessitates the use of
numerical optimisers like stochastic gradient descent, so that step 2 above is
typically only performed approximately.

In the current paper, we use the NPL framework with a loss $l$ that
corresponds to the Maximum Mean Discrepancy (MMD)—a robust discrepancy popular
in GBI methods [20, 57]. This implies that computationally, the second step in
our NPL algorithm amounts to minimum distance estimation as introduced in
[17].

###  2.4 Minimum Distance Estimation with Robust Discrepancies

Since NPL depends on a minimisation step, we will revisit a branch of
frequentist statistics whose theory and methodology we extensively draw on for
our algorithm: Minimum distance estimators (MDEs) [61]. MDEs are a frequentist
approach to parameter estimation. Given a discrepancy
$D:\operatorname{\mathcal{P}}\times\operatorname{\mathcal{P}}\rightarrow\operatorname{\mathbb{R}}_{+}$,
the MDE is

$\textstyle\hat{\theta}_{n}:=\arg\inf_{\theta\in\Theta}D(\operatorname{\mathbb{P}}^{\ast}_{n},\operatorname{\mathbb{P}}_{\theta}).$
(4)

MDEs can be robust to model misspecification when the underlying discrepancy
is chosen with this property in mind. A common choice of discrepancy $D$ is
integral pseudo-probability metrics [53]:

IPM
$\textstyle(\operatorname{\mathbb{P}},\operatorname{\mathbb{Q}})=\sup_{f\in\mathcal{F}}\left|\int_{\mathcal{X}}f(x)\operatorname{\mathbb{P}}(dx)-\int_{\mathcal{X}}f(y)\operatorname{\mathbb{Q}}(dy)\right|.$

IPMs can be thought of as comparing a family of summary statistics indexed by
the class $\mathcal{F}$. There are two common IPMs in the context of
simulators:

(i) Wasserstein Distance Let
$c:\mathcal{X}\times\mathcal{X}\rightarrow[0,\infty)$ be a metric and let
$p\geq 1$. Furthermore, let
$\mathcal{P}_{c,p}=\\{\operatorname{\mathbb{P}}\in\operatorname{\mathcal{P}}\given\int_{\mathcal{X}}c^{p}(x,y)\operatorname{\mathbb{P}}(dx)<\infty\;\forall
y\in\mathcal{X}\\}$. Then the Wasserstein distance is a map
$W:\mathcal{P}_{c,p}\times\mathcal{P}_{c,p}\rightarrow\operatorname{\mathbb{R}}_{+}$
obtained by considering an IPM with
$\mathcal{F}_{W}:=\\{f:\mathcal{X}\rightarrow\operatorname{\mathbb{R}}\given\forall
x,y\in\mathcal{X},|f(x)-f(y)|\leq c(x,y)\\}$. The Wasserstein distance was
used for MDE by [3, 10] and, as previously mentioned, was used for ABC by
[11].

(ii) Maximum Mean Discrepancy Let $\mathcal{H}_{k}$ be a reproducing kernel
Hilbert space (RKHS) with kernel
$k:\mathcal{X}\times\mathcal{X}\rightarrow\operatorname{\mathbb{R}}$ and norm
$\|\cdot\|_{\mathcal{H}_{k}}$. Let
$\mathcal{P}_{k}=\\{\operatorname{\mathbb{P}}\in\operatorname{\mathcal{P}}\given\int_{\mathcal{X}}\sqrt{k(x,x)}\operatorname{\mathbb{P}}(dx)<\infty\\}$.
Then the MMD is a map
$\operatorname{\operatorname{MMD}}:\mathcal{P}_{k}\times\mathcal{P}_{k}\rightarrow\operatorname{\mathbb{R}}_{+}$
obtained by considering an IPM with
$\mathcal{F}_{\text{MMD}}:=\\{f:\mathcal{X}\rightarrow\operatorname{\mathbb{R}}\given\|f\|_{\mathcal{H}_{k}}\leq
1\\}$. When the kernel $k$ is characteristic [73], the MMD is a probability
metric. A common characteristic kernel is the Gaussian kernel
$k(x,x^{\prime})=\exp(-\|x-x^{\prime}\|_{2}^{2}/(2l^{2}))$ where $l>0$.

The MMD was first considered for MDE by [17] and was further explored in [19,
1, 54]. It is closely related to the use of MMD in generative adversarial
networks [25, 47], and was used for goodness-of-fit testing with composite
hypotheses in [43]. As previously mentioned, it has also been used for ABC
[60] and GBI [20, 57].

The Wasserstein distance and MMD are both popular in the context of simulators
because they can both be computed or approximated for empirical measures. One
can therefore use the simulator to sample $m$ realisations from
$\operatorname{\mathbb{P}_{\theta}}$, then use
$D(\widehat{\operatorname{\mathbb{P}}}_{\theta_{b}},\operatorname{\mathbb{P}}^{\ast}_{n})$
as an approximation of
$D(\operatorname{\mathbb{P}_{\theta}},\operatorname{\mathbb{P}}^{\ast}_{n})$
in (4). In this context, the MMD has the advantage over the Wasserstein
distance in that it is a robust distance [17, 19], and can be computed in
quadratic, rather than cubic, time in $n$ and $m$.

##  3 METHODOLOGY

Our paper proposes _Bayesian nonparametric learning with the MMD_ for robust
and scalable inference in simulator models. Compared to ABC, our method has
two main computational advantages: it does not discard any samples, and it is
trivially parallelisable. Furthermore, we will show in Section 4 that the
approach inherits the robustness properties of both the NPL framework and the
MMD; and therefore satisfies a number of desirable properties—including
finite-sample generalisation bounds, robustness guarantees, and frequentist
consistency. Notably, all of these guarantees hold under model
misspecification, highlighting the approach’s usefulness for simulator models
of complex data generating mechanisms.

Assume we have observed data
$x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}$
and are interested in inference with a parametric family
$\operatorname{\mathcal{P}_{\Theta}}$ of simulator-based models. Our approach
is to use the NPL framework with the loss given by the kernel scoring rule
$l_{k}$ [26, 23]:

$\textstyle l_{k}(x,\theta)$
$\textstyle=k(x,x)-2\int_{\mathcal{X}}k(x,y)\operatorname{\mathbb{P}}_{\theta}(dy)$
$\textstyle\qquad+\int_{\mathcal{X}\times\mathcal{X}}k(y,z)\operatorname{\mathbb{P}}_{\theta}(dy)\operatorname{\mathbb{P}}_{\theta}(dz).$

For this loss, the minimiser in (2) becomes the MMD estimator of [17] given
by:

$\textstyle\operatorname{\theta^{*}}(\operatorname{\mathbb{P}})$
$\textstyle:=\arg\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}^{2}(\operatorname{\mathbb{P}},\operatorname{\mathbb{P}}_{\theta})$
$\textstyle\begin{split}\textstyle&\textstyle=\arg\inf_{\theta\in\Theta}\int_{\mathcal{X}\times\mathcal{X}}k(x,y)\operatorname{\mathbb{P}}_{\theta}(dx)\operatorname{\mathbb{P}}_{\theta}(dy)\\\
\textstyle&\textstyle\qquad-\frac{2}{n}\sum_{i=1}^{n}\int_{\mathcal{X}}k(x_{i},x)\operatorname{\mathbb{P}}_{\theta}(dx).\end{split}$
(5)

This objective can easily be approximated by sampling from both measures and
using a U-statistic. The same holds for its gradient, which naturally leads to
a stochastic gradient descent algorithm. Full details are provided in Appendix
C.

input : $x_{1:n}$, $T$, $B$, $\alpha$, $\operatorname{\mathbb{F}}$,
$\operatorname{\mathbb{U}}$, $G_{\theta}$.

1 for  _$j\leftarrow 1$ to $B$_ do

2 Sample
$\tilde{x}_{1:T}^{(j)}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}$
and

3
$(w_{1:n}^{(j)},\tilde{w}_{1:T}^{(j)})\sim\operatorname{\operatorname{Dir}}\left(1,\dots,1,\frac{\alpha}{T},\dots,\frac{\alpha}{T}\right)$.

4 Set
$\operatorname{\mathbb{P}}^{(j)}=\sum_{i=1}^{n}w_{i}^{(j)}\delta_{x_{i}}+\sum_{k=1}^{T}\tilde{w}_{k}^{(j)}\delta_{\tilde{x}_{k}^{(j)}}$.

5 Obtain $\theta^{(j)}=\theta^{\ast}(\operatorname{\mathbb{P}}^{(j)})$ using
numerical optimisation.

6

7 end for

return Posterior bootstrap sample $\theta^{(1:B)}$

Algorithm 1 MMD Posterior Bootstrap

Pseudo code of our approach is given in Algorithm 1, and we call the resulting
procedure MMD posterior bootstrap to pay homage to the fact that in the
limiting case of $\alpha\rightarrow 0$, it is computationally similar to the
Generalised Bayesian Bootstrap of [50]. Incorporating prior information in the
Posterior Bootstrap has also recently been discussed in [63]. From a practical
standpoint, this limiting case is particularly interesting because it
eliminates dependence on the centering measure $\operatorname{\mathbb{F}}$.
This corresponds to a non-informative prior [[, see]]fong2019scalable,
knoblauch2020robust, and reflects our uncertainty about the model in a
misspecified setting.

##  4 THEORY

Before presenting our experiments, we provide a theoretical study for which we
impose:

######  Standing Assumption 1.

$|k(x,y)|\leq 1\quad\forall x,y\in\mathcal{X}$.

While the kernel is required to be bounded, the choice of upper bound $1$ is
without loss of generality, since we only consider the minimiser of the MMD,
which does not depend on the bound’s magnitude. Note that the condition also
ensures that $\operatorname{\mathbb{F}}$,
$\operatorname{\operatorname{\mathbb{P}}^{\ast}}$, and any element of
$\operatorname{\mathcal{P}_{\Theta}}$ are in $\mathcal{P}_{k}$. While no other
assumption on the kernel is needed for our theory, it is desirable for the
kernel to be characteristic—as this guarantees that the MMD is a metric on
$\mathcal{P}_{k}$ and hence that
$\operatorname{\operatorname{\mathbb{P}}^{\ast}}$ can be recovered in well-
specified models.

Since the conjugacy of the DP posterior in the NPL setting relies on the
assumption that the data are i.i.d, we also inherit this requirement:

######  Standing Assumption 2.

$x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}$.

Our results provide the first generalisation, robustness, and consistency
guarantees for NPL posteriors. Our first result is a generalisation bound in
terms of the MMD expected under $\operatorname{\hat{\nu}}$ (Footnote 1).
Beyond that, we also prove that consistency in the frequentist sense (Theorem
4), and robustness to outliers (Corollary 5). A particular characteristic is
that all of these results hold for the misspecified setting where
$\operatorname{\operatorname{\mathbb{P}}^{\ast}}\not\in\operatorname{\mathcal{P}_{\Theta}}$.
The well-specified counterparts are Corollaries 12, 13, and 15 in Appendix B,
and are obtained by noting that in this case
$\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}_{\theta}},\operatorname{\operatorname{\mathbb{P}}^{\ast}})=0$.

###  4.1 Generalisation Error

First, we bound the generalisation error of our procedure, i.e. the error we
expect based on unseen data from the true data-generating mechanism. Here, the
notion of error considered is the MMD to be expected under the (approximated)
posterior, which is given by
$\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})})\right]$.
We therefore bound the expected value of this quantity under unseen data from
the data-generating mechanism:

######  Theorem 3.

111An earlier version of the paper contained a small error in this bound. The
overall rate of the bound and discussions hold the same.

$\textstyle\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})})\right]\right]$
$\textstyle\leq\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}}_{\theta})+\frac{2}{\sqrt{n}}+2\sqrt{\frac{2(n-1)+\alpha(\alpha+1)}{(\alpha+n)(n+\alpha+1)}}$
$\textstyle\quad\quad+2\sqrt{\frac{\alpha(1+\alpha)}{(\alpha+n)(\alpha+n+1)}}.$

Since the expectation over
$\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})})\right]$
is lower-bounded by
$\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})$,
this result tells us that, in expectation, the error is at most
${2}/{\sqrt{n}}+{2\sqrt{2(n-1)+\alpha(\alpha+1)}/{\sqrt{(\alpha+n)(\alpha+n+1)}}}+{2\sqrt{\alpha(1+\alpha)}}/{\sqrt{(\alpha+n)(\alpha+n+1)}}$.
All terms vanish as $n\to\infty$, which implies a $\sqrt{n}$-rate as the first
two terms dominate the third. This is the same rate as for the frequentist MMD
estimators of [17] (Theorem 1), as well as the MMD-based GBI of [19] (Theorem
3.1).

Figure 2: Posterior marginal distributions for the Gaussian location model in
$d=4$. The true parameter $\theta_{0}$ is indicated by dotted lines. (Left &
middle) Comparison against WABC and NPL-WLL methods for the well-specified
case and $\epsilon=0.1$. (Right) Comparison against NPL with the Wasserstein
distance for $\epsilon=0.1$. We note that in the low middle panel, the NPL-WLL
method is not visible as the samples lie significantly away from $\theta_{0}$.
Figure 3: Marginal densities of the posterior obtained by NPL-MMD and ABC with
the MMD for the Gaussian location model in $d=4$ for $\epsilon=0.1$. The true
parameter $\theta_{0}$ is indicated by dotted lines.

###  4.2 Posterior Consistency

To deepen our analysis, we consider consistency in the frequentist sense. This
guarantees that the posterior contracts around the optimal value
$\theta^{\ast}(\operatorname{\operatorname{\mathbb{P}}^{\ast}})$ as
$n\to\infty$. For standard Bayesian inference with posterior measure $\Pi_{n}$
(whose density was previously denoted $\pi(\cdot\given x_{1:n})$) defined on
$\Theta$ directly, and where our model is misspecified so that $\exists C>0$
such that
$\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})=C$,
this would amount to

$\textstyle\Pi_{n}\left(\theta\in\Theta:\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})>C+\frac{M_{n}}{n^{1/2}}\right)\overset{n\to\infty}{\longrightarrow}0$
(6)

for a sequence $M_{n}\rightarrow+\infty$ so that
$M_{n}n^{-\frac{1}{2}}\rightarrow 0$ as $n\rightarrow\infty$. In other words,
the posterior measure over regions of $\Theta$ that induces large values for
$\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})$
goes to $C$ as we obtain more data so that it must ultimately concentrate
around increasingly small neighbourhoods of optimal MMD-minimising values for
$\theta$.

In our case, the posterior $\operatorname{\Pi_{\text{NPL}}}$ is defined
implicitly: given the function $\theta^{\ast}(\operatorname{\mathbb{P}})$ and
the approximate posterior $\operatorname{\hat{\nu}}$ on
$\operatorname{\mathbb{P}}$ constructed via equation (2.3), our posterior on
$\theta$ is obtained by the push-forward operation. Thus, the equivalent
statement for our case concerns $\operatorname{\hat{\nu}}$:

######  Theorem 4.

Suppose our model is misspecified so that for some $C>0$ we have
$\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})=C$.
Then, we have that for any $M_{n}\rightarrow+\infty$:

$\textstyle\operatorname{\hat{\nu}}\left(\operatorname{\mathbb{P}}\in\operatorname{\mathcal{P}}:\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})},\operatorname{\operatorname{\mathbb{P}}^{\ast}})>C+\frac{M_{n}}{n^{1/2}}\right)\overset{n\to\infty}{\longrightarrow}0.$

###  4.3 Robustness to Outliers

To assess robustness against the presence of outliers in the dataset, we
consider Huber’s contamination model [39]. In this setting, a proportion
$1-\varepsilon$ of the observed data is generated from the distribution of
interest $\tilde{\operatorname{\mathbb{P}}}$, and the rest follows a
contaminating noise distribution; i.e.
$\operatorname{\operatorname{\mathbb{P}}^{\ast}}=(1-\epsilon)\tilde{\operatorname{\mathbb{P}}}+\epsilon\operatorname{\mathbb{Q}}$
for
$\tilde{\operatorname{\mathbb{P}}},\operatorname{\mathbb{Q}}\in\operatorname{\mathcal{P}}$
and $\epsilon\in[0,1]$. Here, $\operatorname{\mathbb{Q}}$ is the contaminant,
and so the goal is to place most posterior mass on values of $\theta$ for
which
$\operatorname{\mathbb{P}}_{\theta}\approx\tilde{\operatorname{\mathbb{P}}}$,
where closeness is measured via the MMD.

######  Corollary 5.

Suppose
$\operatorname{\operatorname{\mathbb{P}}^{\ast}}=(1-\epsilon)\tilde{\operatorname{\mathbb{P}}}+\epsilon\operatorname{\mathbb{Q}}$.
Then

$\textstyle\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\big{[}\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\big{[}\operatorname{\operatorname{MMD}}\big{(}\tilde{\operatorname{\mathbb{P}}},\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})}\big{)}\big{]}\big{]}$
$\textstyle\leq\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\tilde{\operatorname{\mathbb{P}}},\operatorname{\mathbb{P}}_{\theta})+4\epsilon+\frac{2}{\sqrt{n}}+2\sqrt{\frac{2(n-1)+\alpha(\alpha+1)}{(\alpha+n)(n+\alpha+1)}}$
$\textstyle\quad\quad+2\sqrt{\frac{\alpha(1+\alpha)}{(\alpha+n)(\alpha+n+1)}}.$

Similarly to the generalisation bound discussed above, the rate at which this
bound goes to zero is $\max\\{n^{-\frac{1}{2}},\epsilon\\}$. Since there are
at most $\lceil\epsilon n\rceil$ contaminated data points in a dataset, the
maximum number of outliers the dataset can have while maintaining the same
rate is of order $\sqrt{n}$, which agrees with [19], Corollary 3.4, who
studied the frequentist minimum MMD estimator.

##  5 EXPERIMENTS

We now study the performance of our method using three examples. Throughout,
we use the Gaussian kernel, which satisfies Standing Assumption 1. We also use
a non-informative prior by setting $\alpha=0$. Further experimental details
and results are reported in Appendix C. Appendix D provides additional
experiments which examine sensitivity to hyperparameters and provide
comparison of our method with the MMD-Bayes method in [57]. We further
consider an example of misspecification which is not based in a contamination
model by wrongly fitting a Gaussian location model to Cauchy distributed data.
The code for all experiments can be found at
https://github.com/haritadell/npl_mmd_project.git.

###  5.1 Gaussian Location Model

We start by considering a toy example, the Gaussian location model. While the
likelihood of this model is available, we treat it as a simulator to study the
properties of our proposed method. We take
$\operatorname{\mathbb{P}}_{\theta}=\mathcal{N}(\theta,I_{d\times d})$ and use
a true data generating process
$\operatorname{\operatorname{\mathbb{P}}^{\ast}}=(1-\epsilon)\operatorname{\mathbb{P}}_{\theta_{0}}+\epsilon\operatorname{\mathbb{P}}_{\theta^{\prime}}$,
with $\theta_{0}=(1,\dots,1)\in\operatorname{\mathbb{R}}^{d}$ and
$\theta^{\prime}=(20,\dots,20)\in\operatorname{\mathbb{R}}^{d}$ in $d=4$. We
assess robustness by considering both the well-specified case $\epsilon=0$ and
the case $\epsilon=0.1$ and $n=200$ realisations from
$\operatorname{\operatorname{\mathbb{P}}^{\ast}}$.

Our simulation study will illustrate how robustness is inherited from both the
NPL framework, which is more robust to model misspecification than standard
Bayes or ABC, and the MMD being more robust than alternative losses such as
the Wasserstein distance or negative log-likelihood.

##### Wasserstein ABC

We first compare our method against the Wasserstein ABC (WABC) [11] algorithm,
which uses Sequential Monte Carlo (SMC). We compare against WABC since it is a
popular algorithm which, unlike standard ABC, does not require hand-crafting
of summary statistics. As observed in Figure 2, NPL-MMD outperforms WABC in
both the well-specified and misspecified settings. The former is explained by
the fact that the Wasserstein distance exhibits poor sample complexity for
$d>1$ [29], whereas the MMD can be estimated at a $\sqrt{n}-$rate. The latter
is explained by the fact that the WABC is an approximation of the exact
Bayesian posterior, which is not robust.

##### NPL with Wasserstein distance

Secondly, we consider the NPL framework using the Wasserstein distance (NPL-
WAS) such that for
$W:\operatorname{\mathbb{P}}^{2}\to\operatorname{\mathbb{R}}_{\geq 0}$ the
2-Wasserstein distance,
$\theta^{*}_{W}(\operatorname{\mathbb{P}}):=\arg\inf_{\theta\in\Theta}W(\operatorname{\mathbb{P}},\operatorname{\mathbb{P}}_{\theta}).$
We use the POT package [27] for approximating the Wasserstein distance and the
Powell optimiser from Scipy [77]. We focus specifically on the misspecified
setting, and notice that NPL-MMD outperforms the NPL-WAS in that case. This
clearly demonstrates the advantage of using a robust loss function, even when
using a robust inference framework such as NPL.

Figure 4: Posterior marginal distributions for $\theta$ in the univariate
g-and-k model for an increasing percentage of outliers present in the dataset.
The true parameter $\theta_{0}$ is indicated by the dotted lines.

##### NPL with log-likelihood

Thirdly, the availability of the likelihood in this toy problem allows for
comparison with the original NPL posterior bootstrap using
$l(x,\theta)=\frac{1}{\sqrt{2\pi}}\exp\\{-\frac{1}{2}(x-\theta)^{2}\\}$ as in
[49], which we call the weighted log-likelihood NPL (NPL-WLL). Figure 2 shows
that NPL-WLL and NPL-MMD perform similarly in the well-specified case, but
NPL-MMD significantly outperforms NPL-WLL in the misspecified case. This is
once again due to the fact that the MMD is a robust loss, whereas the negative
log-likelihood is not.

##### ABC with the MMD

Finally, we further explore ABC method with the SMC samplers considered in
[11], this time using the MMD instead of the Wasserstein distance. The ABC
with the MMD has also been previously explored in [[]]park2016k2. The
posterior marginals are visualised in Figure 3 for $\epsilon=0.1$. This figure
clearly demonstrates the advantage of the NPL framework over ABC.

###  5.2 Simulator-based Models

We now consider two more complex numerical examples, for which likelihood-
based inference is impossible. We compare solely against the Wasserstein ABC
for simplicity by directly using the experimental setup in [11].

Figure 5: Posterior marginal distributions for $\theta$ in the toggle switch
model with $10\%$ of noisy data in the dataset. Bottom right figure presents
2000 samples from the well-specified toggle switch model.

##### G-and-k Distribution Model

First, consider $G_{\theta}:[0,1]^{2}\rightarrow\operatorname{\mathbb{R}}$
with $\theta=(a,b,g,k)$ where

$\textstyle G_{\theta}(u)$
$\textstyle=a+b\left(1+0.8\frac{1-\exp(-gz(u))}{1+\exp(-gz(u))}\right)(1+z(u)^{2})^{k}z(u),$

$z(u)=\sqrt{-2\log(u_{1})}\cos(2\pi u_{2})$ and
$\operatorname{\mathbb{U}}=\text{Unif}([0,1]^{2})$. Parameters $a,b,g$ and $k$
control the location, scale, skewness and kurtosis respectively. For
computational convenience, we reparametrise the last parameter by setting
$k^{\prime}=\exp(k)$. Although $\operatorname{\mathbb{P}}_{\theta}$ is one-
dimensional, it is a popular baseline for simulator-based models [64] because
of the challenge of inferring the four parameters simultaneously. It has also
been used extensively in applications; for example to model the price of
AirBnB rentals [68], the air pollution [66] or for non-life insurance
modelling [62].

Our data consists of $n=2^{11}$ realisations from
$\operatorname{\operatorname{\mathbb{P}}^{\ast}}=(1-\epsilon)\operatorname{\mathbb{P}}_{\theta_{0}}+\epsilon\operatorname{\mathbb{Q}}$
where $\operatorname{\mathbb{P}}_{\theta_{0}}$ denotes the g-and-k with
$\theta_{0}=(3,1,1,-\log(2))$, and $\operatorname{\mathbb{Q}}$ is the shifted
distribution
$\operatorname{\mathbb{Q}}=\operatorname{\mathbb{P}}_{\theta_{0}}\pm 50$ with
an equal number of points shifted to either direction. The resulting
posteriors are shown in Figure 4 displayed as bivariate plots. The WABC method
appears more sensitive to contamination in the dataset—in contrast to the NPL-
MMD, which concentrates significantly closer to the true parameter values for
an increasing proportion of outliers. This is particularly the case for the
last two parameters $g$ and $\log(k)$ which are well-known to be more
challenging to estimate.

##### Toggle Switch Model

Finally, we consider the toggle-switch model arising in Systems Biology [15,
14]. This is a dynamic model used to study cellular networks; more precisely,
the network describes the interaction of two genes $u$ and $v$ over time. The
simulator is too complex to include in the main text, but is given in Appendix
C. The data is one-dimensional, but the model has $7$ parameters and the
latent space is $601$-dimensional.

We consider inference on $\theta$ for $n=2000$ data points simulated from the
toggle-switch model with true parameter
$\theta_{0}=(22,12,4,4.5,325,0.25,0.15)$ in which $10\%$ of the data have some
added Cauchy noise of location parameter 0 and scale parameter 10. Such noise
can be interpreted as measurement error in the collection of data. The
posterior marginal distributions for $\theta$ are shown in Figure 5, and
indicate that the NPL-MMD method is successful in concentrating around
$\theta_{0}$ despite the Cauchy noise.

###  5.3 Computational Time

We provide the computational time recorded for each experiment over multiple
independent runs. The Gaussian location and G-and-k experiments were repeated
10 times and the Toggle-switch was repeated 5 times due to the higher
computational cost. For each run, a new dataset was generated and $B$
posterior samples were obtained for $\theta$. The average time for the
generation of $B$ samples is recorded for each method. Our method is
compatible with the use of GPU so in Table 1 we provide the average time
recorded when run in Google Colaboratory Pro.

Table 1: Average Clock Time In Seconds. Method | Gaussian | G-and-k | Toggle-S
---|---|---|---
NPL-MMD | $1.09\times 10^{2}$ | $3.02\times 10$ | $7.42\times 10^{3}$
WABC | $3.34\times 10^{3}$ | $1.05\times 10^{2}$ | $3.08\times 10^{4}$

##  6 CONCLUSION

Our paper proposes the very first posterior bootstrap method in the simulator
setting. Unlike ABC, our method does not discard any samples and can run in
parallel—leading to a substantially decreased computational burden. Further,
the approach is based on MMD estimators, and therefore inherits both the
robustness of the NPL framework and that of the MMD. We support this claim
with three theoretical results that hold even in the presence of
misspecification. The results include a generalisation bound, a robustness
guarantee, and a consistency result. We verified the practical utility of our
theory through deploying it on three models, which highlighted both robustness
and computational advantages of our method.

A particular strength of the method is that no assumption is required on
$\mathcal{X}$. Hence, the results apply directly to any set $\mathcal{X}$ on
which one can define kernels (e.g. graphs, strings, or discrete spaces).

In future work, we would like to tackle the challenges posed by the
optimisation step. Specifically, since the MMD-objective is usually non-
convex, a particularly interesting question would be if the kernel can be used
to enforce a more well-behaved objective. Another direction for improvement is
the cost of computing the MMD minimiser. While naively, this scales
quadratically in the number of observations, it could be reduced to linear
time using approaches like in Lemma 14 of [33].

#### Acknowledgements

CD was funded by EPSRC grant [EP/T51794X/1] as part of the Warwick CDT in
Mathematics and Statistics. JK was funded through the Biometrika Fellowship
courtesy of the Biometrika Trust. TD acknowledges support from a UKRI Turing
AI Fellowship [EP/V02678X/1]. TD and FXB were supported by the Lloyd’s
Register Foundation Programme on Data-Centric Engineering and The Alan Turing
Institute under the EPSRC grant [EP/N510129/1]. The authors thank all
reviewers for their useful comments.

#### References

## References

  * [1] P. Alquier and M. Gerber  “Universal robust regression via maximum mean discrepancy”  In _arXiv:2006.00840_ , 2020
  * [2] Christophe Andrieu and Gareth O. Roberts  “The pseudo-marginal approach for efficient Monte Carlo computations”  In _Annals of Statistics_ 37.2, 2009, pp. 697–725
  * [3] F. Bassetti, A. Bodini and E. Regazzini  “On minimum Kantorovich distance estimators”  In _Statistics and Probability Letters_ 76, 2006, pp. 1298–1302
  * [4] M.. Beaumont, J-M. Cornuet, J-M. Marin and C.. Robert  “Adaptive approximate Bayesian computation”  In _Biometrika_ 96.4, 2009, pp. 983–990
  * [5] Mark A Beaumont  “Approximate Bayesian computation in evolution and ecology”  In _Annual Review of Ecology, Evolution, and Systematics_ 41, 2010, pp. 379–406
  * [6] Mark A Beaumont  “Approximate Bayesian computation”  In _Annual Review of Statistics and Its Application_ 6, 2019, pp. 379–403
  * [7] Mark A Beaumont, Wenyang Zhang and David J Balding  “Approximate Bayesian computation in population genetics”  In _Genetics_ 162.4 Oxford University Press, 2002, pp. 2025–2035
  * [8] Alain Berlinet and Christine Thomas-Agnan  “Reproducing kernel Hilbert spaces in probability and statistics”  Springer Science & Business Media, 2011
  * [9] José M Bernardo and Adrian FM Smith  “Bayesian theory”  John Wiley & Sons, 2009
  * [10] E. Bernton, P.. Jacob, M. Gerber and C.. Robert  “Inference in generative models using the Wasserstein distance”  In _Information and Inference: A Journal of the IMA_ 8.4, 2017, pp. 657–676
  * [11] Espen Bernton, Pierre E Jacob, Mathieu Gerber and Christian P Robert  “Approximate Bayesian computation with the Wasserstein distance”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 81.2 Wiley Online Library, 2019, pp. 235–269
  * [12] Ayush Bharti, Francois-Xavier Briol and Troels Pedersen  “A general method for calibrating stochastic radio channel models with kernels”  In _IEEE Transactions on Antennas and Propagation_ IEEE, 2021
  * [13] Pier Giovanni Bissiri, Chris C Holmes and Stephen G Walker  “A general framework for updating belief distributions”  In _Journal of the Royal Statistical Society. Series B, Statistical methodology_ 78.5 Wiley-Blackwell, 2016, pp. 1103
  * [14] Fernando V Bonassi and Mike West  “Sequential Monte Carlo with adaptive weights for approximate Bayesian computation”  In _Bayesian Analysis_ 10.1 International Society for Bayesian Analysis, 2015, pp. 171–187
  * [15] Fernando V Bonassi, Lingchong You and Mike West  “Bayesian learning from marginal data in bionetwork models”  In _Statistical applications in genetics and molecular biology_ 10.1 De Gruyter, 2011
  * [16] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin and Skye Wanderman-Milne  “JAX: composable transformations of Python+ NumPy programs, 2018”  In _URL http://github. com/google/jax_ 4, 2020, pp. 16
  * [17] Francois-Xavier Briol, Alessandro Barp, Andrew B Duncan and Mark Girolami  “Statistical inference for generative models with maximum mean discrepancy”  In _arXiv:1906.05944_ , 2019
  * [18] O Catoni  “Pac-Bayesian supervised classification: the thermodynamics of statistical learning”  In _Institute of Mathematical Statistics Lecture Notes—Monograph Series 56_ , 2007
  * [19] Badr-Eddine Chérief-Abdellatif and Pierre Alquier  “Finite sample properties of parametric MMD estimation: robustness to misspecification and dependence”  In _arXiv preprint arXiv:1912.05737_ , 2019
  * [20] Badr-Eddine Chérief-Abdellatif and Pierre Alquier  “MMD-Bayes: robust Bayesian estimation via maximum mean discrepancy”  In _Symposium on Advances in Approximate Bayesian Inference_ , 2020, pp. 1–21  PMLR
  * [21] V. Chernozhukov and H. Hong  “An MCMC approach to classical estimation”  In _Journal of Econometrics_ 115.2, 2003, pp. 293–346
  * [22] K. Cranmer, J. Brehmer and G. Louppe  “The frontier of simulation-based inference”  In _Proceedings of the National Academy of Sciences of the United States of America_ 117.48, 2020
  * [23] A.. Dawid  “The geometry of proper scoring rules”  In _Annals of the Institute of Statistical Mathematics_ 59.1, 2007, pp. 77–93
  * [24] Pierre Del Moral, Arnaud Doucet and Ajay Jasra  “On adaptive resampling strategies for sequential Monte Carlo methods”  In _Bernoulli_ 18.1, 2012, pp. 252–278
  * [25] Gintare Karolina Dziugaite, Daniel M Roy and Zoubin Ghahramani  “Training generative neural networks via maximum mean discrepancy optimization”  In _Proceedings of the Thirty-First Conference on Uncertainty in Artificial Intelligence_ , 2015, pp. 258–267
  * [26] M.. Eaton  “A method for evaluating improper prior distributions”  In _Statistical Decision Theory and Related Topics III_ , 1982, pp. 329–352
  * [27] Rémi Flamary et al.  “Pot: Python optimal transport”  In _Journal of Machine Learning Research_ 22.78, 2021, pp. 1–8
  * [28] Edwin Fong, Simon Lyddon and Chris Holmes  “Scalable nonparametric sampling from multimodal posteriors with the posterior bootstrap”  In _International Conference on Machine Learning_ , 2019, pp. 1952–1962  PMLR
  * [29] N. Fournier and A. Guillin  “On the rate of convergence in Wasserstein distance of the empirical measure”  In _Probability Theory and Related Fields_ 162.3-4, 2015, pp. 707–738  DOI: 10.1007/s00440-014-0583-7
  * [30] D.. Frazier, C.. Robert and J. Rousseau  “Model misspecification in ABC: consequences and diagnostics”  In _Journal of the Royal Statistical Society B: Statistical Methodology_ 82.2, 2020, pp. 421–444
  * [31] Marta Galvani, Chiara Bardelli, Silvia Figini and Pietro Muliere  “A Bayesian nonparametric learning approach to ensemble models using the proper Bayesian bootstrap”  In _Algorithms_ 14.1 Multidisciplinary Digital Publishing Institute, 2021, pp. 11
  * [32] Abhik Ghosh and Ayanendranath Basu  “Robust Bayes estimation using the density power divergence”  In _Annals of the Institute of Statistical Mathematics_ 68.2 Springer, 2016, pp. 413–437
  * [33] Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard Schölkopf and Alexander Smola  “A kernel two-sample test”  In _The Journal of Machine Learning Research_ 13.1 JMLR. org, 2012, pp. 723–773
  * [34] P Grünwald and Thijs Van Ommen  “Inconsistency of Bayesian inference for misspecified linear models, and a proposal for repairing it”  In _Bayesian Analysis_ 12.4 International Society for Bayesian Analysis, 2017, pp. 1069–1103
  * [35] Peter Grünwald  “Safe learning: bridging the gap between Bayes, MDL and statistical learning theory via empirical convexity”  In _Journal of Machine Learning Research_ 19, 2011, pp. 397–419
  * [36] Peter Grünwald  “The safe Bayesian”  In _International Conference on Algorithmic Learning Theory_ , 2012, pp. 169–183  Springer
  * [37] J. Hermans, A. Delaunoy, F. Rozet, A. Wehenkel and G. Louppe  “Averting a crisis in simulation-based inference”  In _arXiv:2110.06581_ , 2021
  * [38] Giles Hooker and Anand N Vidyashankar  “Bayesian model robustness via disparities”  In _Test_ 23.3 Springer, 2014, pp. 556–584
  * [39] Peter J Huber  “Robust estimation of a location parameter”  In _Breakthroughs in statistics_ Springer, 1992, pp. 492–518
  * [40] Jonathan H Huggins and Jeffrey W Miller  “Robust inference and model criticism using bagged posteriors”  In _arXiv preprint arXiv:1912.07104_ , 2019
  * [41] Hemant Ishwaran and Mahmoud Zarepour  “Exact and approximate sum representations for the Dirichlet process”  In _Canadian Journal of Statistics_ 30.2 Wiley Online Library, 2002, pp. 269–283
  * [42] J. Jewson, J.. Smith and C. Holmes  “Principled Bayesian minimum divergence inference”  In _Entropy_ 20.6, 2018, pp. 442
  * [43] O. Key, T. Fernandez, A. Gretton and F-X. Briol  “Composite goodness-of-fit tests with kernels”  In _arXiv:2111.10275_ , 2021
  * [44] Diederik P Kingma and Jimmy Ba  “Adam: a method for stochastic optimization”  In _arXiv preprint arXiv:1412.6980_ , 2014
  * [45] J. Knoblauch, J. Jewson and T. Damoulas  “Generalized variational inference: three arguments for deriving new posteriors”  In _arXiv:1904.02063_ , 2019
  * [46] Jeremias Knoblauch and Lara Vomfell  “Robust Bayesian inference for discrete outcomes with the total variation distance”  In _arXiv:2010.13456_ , 2020
  * [47] Yujia Li, Kevin Swersky and Rich Zemel  “Generative moment matching networks”  In _International Conference on Machine Learning_ , 2015, pp. 1718–1727  PMLR
  * [48] Ruben Loaiza-Maya, Gael M Martin and David T Frazier  “Focused Bayesian prediction”  In _Journal of Applied Econometrics_ 36.5 John Wiley & Sons, Ltd., 2021, pp. 517–543
  * [49] Simon Lyddon, Stephen Walker and Chris Holmes  “Nonparametric learning from Bayesian models with randomized objective functions”  In _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , 2018, pp. 2075–2085
  * [50] SP Lyddon, CC Holmes and SG Walker  “General Bayesian updating and the loss-likelihood bootstrap”  In _Biometrika_ 106.2 Oxford University Press, 2019, pp. 465–478
  * [51] T. Matsubara, J. Knoblauch, F.-X. Briol and C.. Oates  “Robust generalised Bayesian inference for intractable likelihoods”  In _arXiv:2104.07359_ , 2021
  * [52] Pietro Muliere and Piercesare Secchi  “Bayesian nonparametric predictive inference and bootstrap techniques”  In _Annals of the Institute of Statistical Mathematics_ 48.4 Springer, 1996, pp. 663–673
  * [53] Alfred Müller  “Integral probability metrics and their generating classes of functions”  In _Advances in Applied Probability_ 29.2 Cambridge University Press, 1997, pp. 429–443
  * [54] Z. Niu, J. Meier and F-X Briol  “Discrepancy-based inference for intractable generative models using quasi-Monte Carlo”  In _arXiv:2106.11561_ , 2021
  * [55] Houman Owhadi, Clint Scovel and Tim Sullivan  “Brittleness of Bayesian inference under finite information in a continuous world”  In _Electronic Journal of Statistics_ 9.1 Institute of Mathematical StatisticsBernoulli Society, 2015, pp. 1–79
  * [56] Lorenzo Pacchiardi and Ritabrata Dutta  “Score matched conditional exponential families for likelihood-free inference”  In _arXiv:2012.10903_ , 2020
  * [57] Lorenzo Pacchiardi and Ritabrata Dutta  “Generalized Bayesian likelihood-free inference using scoring rules estimators”  In _arXiv:2104.03889_ , 2021
  * [58] George Papamakarios and Iain Murray  “Fast $\varepsilon$-free inference of simulation models with Bayesian conditional density estimation”  In _Advances in Neural Information Processing Systems_ , 2016, pp. 1028–1036
  * [59] George Papamakarios, David Sterratt and Iain Murray  “Sequential neural likelihood: fast likelihood-free inference with autoregressive flows” 89, International Conference on Artificial Intelligence and Statistics, 2019, pp. 837–848
  * [60] Mijung Park, Wittawat Jitkrittum and Dino Sejdinovic  “K2-ABC: approximate Bayesian computation with kernel embeddings”  In _Artificial Intelligence and Statistics_ , 2016, pp. 398–407  PMLR
  * [61] W C Parr and W R Schucany  “Minimum distance and robust estimation”  In _Journal of the American Statistical Association_ 75.371 Taylor & Francis, 1980, pp. 616–624
  * [62] G. Peters, W. Chen and R. Gerlach  “Estimating quantile families of loss distributions for non-life insurance modelling via L-moments”  In _Risks_ 4.2, 2016, pp. 14
  * [63] Emilia Pompe  “Introducing prior information in Weighted Likelihood Bootstrap with applications to model misspecification”  In _arXiv preprint arXiv:2103.14445_ , 2021
  * [64] Dennis Prangle  “gk: an R Package for the g-and-k and Generalised g-and-h distributions”  In _The R Journal_ 12.1, 2020, pp. 7–20
  * [65] L.. Price, C.. Drovandi, A. Lee and D.. Nott  “Bayesian synthetic likelihood”  In _Journal of Computational and Graphical Statistics_ 8600, 2017
  * [66] G.. Rayner and H.. MacGillivray  “Numerical maximum likelihood estimation for the g-and-k and generalized g-and-h distributions”  In _Statistics and Computing_ 12.1, 2002, pp. 57–75
  * [67] Adam J Riesselman, John B Ingraham and Debora S Marks  “Deep generative models of genetic variation capture the effects of mutations”  In _Nature methods_ 15.10 Nature Publishing Group, 2018, pp. 816–822
  * [68] G.. Rodrigues, David J. Nott and S.. Sisson  “Likelihood-free approximate Gibbs sampling”  In _Statistics and Computing_ 30.4, 2020, pp. 1057–1073
  * [69] Sebastian M Schmon, Patrick W Cannon and Jeremias Knoblauch  “Generalized Posteriors in approximate Bayesian computation”  In _Third Symposium on Advances in Approximate Bayesian Inference_ , 2020
  * [70] Jayaram Sethuraman  “A constructive definition of Dirichlet priors”  In _Statistica sinica_ JSTOR, 1994, pp. 639–650
  * [71] John Shawe-Taylor and Robert C Williamson  “A PAC analysis of a Bayesian estimator”  In _Proceedings of the tenth annual conference on Computational learning theory_ 6.09, 1997, pp. 2–9
  * [72] SA Sisson and Y Fan  “ABC samplers”  In _Handbook of Approximate Bayesian Computation_ ChapmanHall/CRC, 2018, pp. 87–123
  * [73] Bharath K Sriperumbudur, Arthur Gretton, Kenji Fukumizu, Bernhard Schölkopf and Gert RG Lanckriet  “Hilbert space embeddings and metrics on probability measures”  In _The Journal of Machine Learning Research_ 11 JMLR. org, 2010, pp. 1517–1561
  * [74] Danica J. Sutherland, Hsiao-Yu Tung, Heiko Strathmann, Soumyajit De, Aaditya Ramdas, Alexander J. Smola and Arthur Gretton  “Generative models and model criticism via optimized maximum mean discrepancy”  In _International Conference on Learning Representations_ , 2017
  * [75] Owen Thomas, Henri Pesonen and Jukka Corander  “Generalised Bayes Updates with $f$-divergences through Probabilistic Classifiers”  In _arXiv preprint arXiv:2007.04358_ , 2020
  * [76] Owen Thomas, Henri Pesonen, Raquel Sá-Leão, Hermı́nia Lencastre, Samuel Kaski and Jukka Corander  “Split-BOLFI for for misspecification-robust likelihood free inference in high dimensions”  In _arXiv preprint arXiv:2002.09377_ , 2020
  * [77] Pauli Virtanen et al.  “SciPy 1.0: fundamental algorithms for scientific computing in Python”  In _Nature methods_ 17.3 Nature Publishing Group, 2020, pp. 261–272
  * [78] Pei-Shien Wu and Ryan Martin  “A comparison of learning rate selection methods in generalized Bayesian inference”  In _arXiv:2012.11349_ , 2020

Supplementary Material:
Robust Bayesian Inference for Simulator-based Models
via the MMD Posterior Bootstrap

In Appendix A we summarise the notation used throughout the paper. In Appendix
B we prove all theoretical results. In Appendix C we provide details on the
experiments introduced in the main text and in Appendix D we provide some
additional experiments.

##  Appendix A NOTATION

In the following section, we recall the notation used in the paper:

[type=symbols, title=]

##  Appendix B PROOFS

We now provide proofs for all theoretical results presented in the main text
as well as some additional theory. We start with the necessary background
material on the MMD that we use in our proofs. We continue by proving our
generalisation error bound both in the case where the DP posterior sample is
represented exactly by an infinite sum of the stick-breaking process (Theorem
9) and using the approximation of algorithm 1 (Theorem 1). We also provide the
corresponding corollaries for the well-specified case (Corollaries 12 and 10)
for reference. We then proceed by proving the general posterior consistency
result of Theorem 4; along with the special case of a well-specified model
(Corollary 13). Finally, we prove the results in the case of a contaminated
data generating process in both the misspecified and well-specified case
(Corollaries 5 and 15). Recall that throughout all our theoretical results we
impose Standing Assumptions 1 and 2.

###  B.1 MMD Through Kernel Mean Embeddings

In the following proofs we will use an equivalent definition of the MMD
through kernel mean embeddings. Consider an RKHS $\mathcal{H}_{k}$ indexed by
a reproducing kernel $k$. We say $k$ is a reproducing kernel if (i)
$k(\cdot,x)\in\mathcal{H}_{k}$, $\forall x\in\mathcal{X}$, (ii)
$\left<f,k(\cdot,x)\right>_{\mathcal{H}_{k}}=f(x)$ $\forall x\in\mathcal{X}$
and $\forall f\in\mathcal{H}_{k}$ (reproducing property); and inner product
$\left<\cdot,\cdot\right>_{\mathcal{H}_{k}}$; see [8]. For a function
$f\in\mathcal{H}_{k}$ the mean embedding
$\mu_{\operatorname{\mathbb{P}}}\in\mathcal{H}_{k}$ with respect to
probability measure $\operatorname{\mathbb{P}}\in\mathcal{P}$ is defined as

$\textstyle\mu_{\operatorname{\mathbb{P}}}(\cdot):=\mathbb{E}_{X\sim\operatorname{\mathbb{P}}}[k(X,\cdot)]\in\mathcal{H}_{k}$

and satisfies

$\textstyle\mathbb{E}_{X\sim\operatorname{\mathbb{P}}}[f(X)]=\left<f,\mu_{\operatorname{\mathbb{P}}}\right>_{\mathcal{H}_{k}}.$

The MMD between probability measures $\operatorname{\mathbb{P}}$ and
$\operatorname{\mathbb{Q}}$ takes the form [[, see]Lemma 4]gretton2012kernel

$\textstyle\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}},\operatorname{\mathbb{Q}})=\|\mu_{\operatorname{\mathbb{P}}}-\mu_{\operatorname{\mathbb{Q}}}\|_{\mathcal{H}_{k}}.$

###  B.2 Generalisation Error for the Stick-Breaking Process

We start with a bound of the generalisation error for the case where the DP
sample is represented exactly—i.e., by the infinite sum obtained from the
stick-breaking process. We denote by $\nu$ the probability measure on
$\operatorname{\mathcal{P}}$ corresponding to the stick-breaking process given
as

$\textstyle w_{1:\infty}$
$\textstyle\sim\operatorname{\operatorname{GEM}}(\alpha^{\prime}),\quad\alpha^{\prime}:=\alpha+n,$
(7) $\textstyle z_{1:\infty}$
$\textstyle\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}:=\frac{\alpha}{\alpha+n}\operatorname{\mathbb{F}}+\frac{n}{n+\alpha}\mathbb{P}_{n},$
(8) $\textstyle\operatorname{\mathbb{P}}$
$\textstyle=\sum_{i=1}^{\infty}w_{i}\delta_{z_{i}}\sim\operatorname{\nu}.$

Note that instead of writing the expectation directly over
$\operatorname{\nu}$, we instead often use the separate expectations
$\operatorname{\mathbb{E}}_{w_{1:\infty}\sim\operatorname{\operatorname{GEM}}(\alpha^{\prime})}$
and
$\operatorname{\mathbb{E}}_{z_{1:\infty}\sim\operatorname{\mathbb{F}}^{\prime}}$
induced by $\operatorname{\nu}$ in the proofs below.

Before stating and proving the main result of this section, we provide three
lemmas that bound the expected MMD between the true data generating mechanism
$\operatorname{\operatorname{\mathbb{P}}^{\ast}}$ and $\mathbb{P}_{n}$, the
centering measure of the DP posterior $\operatorname{\mathbb{F}}^{\prime}$ and
$\operatorname{\mathbb{P}}$ and lastly $\operatorname{\mathbb{F}}^{\prime}$
and $\mathbb{P}_{n}$. The reason for this is that the main proof will use a
decomposition of the MMD using the triangle inequality; and the following
technical Lemmas bound three of the trickier terms arising from said triangle
inequality.

######  Lemma 6.

For
$x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}$
we have

$\textstyle\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\operatorname{MMD}}(\mathbb{P}_{n},\operatorname{\operatorname{\mathbb{P}}^{\ast}})\right]\leq\frac{1}{\sqrt{n}}.$

where $\mathbb{P}_{n}$ denotes the empirical measure of the sample data, i.e.
$\mathbb{P}_{n}=\frac{1}{n}\sum_{i=1}^{n}\delta_{x_{i}}$.

###### Proof.

Proof follows directly from the proof of [19], Lemma 7.1 and using the
Jensen’s inequality to obtain

$\textstyle\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\operatorname{MMD}}(\mathbb{P}_{n},\operatorname{\operatorname{\mathbb{P}}^{\ast}})\right]\leq\sqrt{\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\operatorname{MMD}}^{2}(\mathbb{P}_{n},\operatorname{\operatorname{\mathbb{P}}^{\ast}})\right]}\leq\frac{1}{\sqrt{n}}.$

∎

######  Lemma 7.

Let
$\operatorname{\mathbb{P}}=\sum_{i=1}^{\infty}w_{i}\delta_{z_{i}}\sim\operatorname{\nu}$
denote a sample from the DP posterior, where
$z_{1:\infty}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}$
and $w_{1:\infty}\sim\operatorname{\operatorname{GEM}}(\alpha^{\prime})$, then

$\textstyle\operatorname{\mathbb{E}}_{w_{1:\infty}\sim\operatorname{\operatorname{GEM}}(\alpha^{\prime})}\left[\operatorname{\mathbb{E}}_{z_{1:\infty}\sim\operatorname{\mathbb{F}}^{\prime}}\left[\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}},\operatorname{\mathbb{F}}^{\prime})\right]\right]\leq\frac{1}{\sqrt{\alpha+n+2}}.$

###### Proof.

First note that

$\textstyle\operatorname{\operatorname{MMD}}^{2}(\operatorname{\mathbb{P}},\operatorname{\mathbb{F}}^{\prime})$
$\textstyle=\|\mu_{\operatorname{\mathbb{P}}}-\mu_{\operatorname{\mathbb{F}}^{\prime}}\|_{\mathcal{H}_{k}}^{2}$
$\textstyle=\|\sum_{i=1}^{\infty}w_{i}k(z_{i},\cdot)-\mu_{\operatorname{\mathbb{F}}^{\prime}}\|_{\mathcal{H}_{k}}^{2}$
$\textstyle=\|\sum_{i=1}^{\infty}w_{i}\left[k(z_{i},\cdot)-\mu_{\operatorname{\mathbb{F}}^{\prime}}\right]\|_{\mathcal{H}_{k}}^{2}$
$\textstyle=\sum_{i=1}^{\infty}w_{i}^{2}\|k(z_{i},\cdot)-\mu_{\operatorname{\mathbb{F}}^{\prime}}\|_{\mathcal{H}_{k}}^{2}+2\sum_{i\neq
j}w_{i}w_{j}\left<k(z_{i},\cdot)-\mu_{\operatorname{\mathbb{F}}^{\prime}},k(z_{j},\cdot)-\mu_{\operatorname{\mathbb{F}}^{\prime}}\right>$

Note that since
$z_{1:\infty}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}$
we have that for any $i\neq j$:

$\textstyle\operatorname{\mathbb{E}}_{z_{i},z_{j}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}\left[\left<k(z_{i},\cdot)-\mu_{\operatorname{\mathbb{F}}^{\prime}},k(z_{j},\cdot)-\mu_{\operatorname{\mathbb{F}}^{\prime}}\right>\right]$
$\textstyle=\operatorname{\mathbb{E}}_{z_{i},z_{j}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}\left[\left<k(z_{i},\cdot)-\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}[k(z_{i},\cdot)],k(z_{j},\cdot)-\operatorname{\mathbb{E}}_{z_{j}\sim\operatorname{\mathbb{F}}^{\prime}}[(k(z_{j},\cdot)]\right>\right]$
$\textstyle=\operatorname{\mathbb{E}}_{z_{i},z_{j}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}[k(z_{i},z_{j})]-\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}[k(z_{i},\cdot)]\operatorname{\mathbb{E}}_{z_{j}\sim\operatorname{\mathbb{F}}^{\prime}}[k(z_{j},\cdot)]$
$\textstyle\qquad-\operatorname{\mathbb{E}}_{z_{j}\sim\operatorname{\mathbb{F}}^{\prime}}[k(z_{j},\cdot)]\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}[k(z_{i},\cdot)]+\operatorname{\mathbb{E}}_{z_{j}\sim\operatorname{\mathbb{F}}^{\prime}}[k(z_{j},\cdot)]\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}[k(z_{i},\cdot)]$
$\textstyle=\operatorname{\mathbb{E}}_{z_{i},z_{j}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}[k(z_{i},z_{j})]-\operatorname{\mathbb{E}}_{z_{i},z_{j}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}[k(z_{i},z_{j})]-\operatorname{\mathbb{E}}_{z_{i},z_{j}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}[k(z_{i},z_{j})]+\operatorname{\mathbb{E}}_{z_{i},z_{j}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}[k(z_{i},z_{j})]$
$\textstyle=0.$ (9)

Moreover, for any $i=1,2,\dots$

$\textstyle\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}\left[\|k(z_{i},\cdot)-\mu_{\operatorname{\mathbb{F}}^{\prime}}\|_{\mathcal{H}_{k}}^{2}\right]$
$\textstyle=\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}\left[\|k(z_{i},\cdot)-\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}[k(z_{i},\cdot)]\|_{\mathcal{H}_{k}}^{2}\right]$
$\textstyle=\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}\left[\|k(z_{i},\cdot)\|_{\mathcal{H}_{k}}^{2}-2\left<k(z_{i},\cdot),\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}[k(z_{i},\cdot)]\right>+\|\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}[k(z_{i},\cdot)]\|_{\mathcal{H}_{k}}^{2}\right]$
$\textstyle=\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}\left[\|k(z_{i},\cdot)\|_{\mathcal{H}_{k}}^{2}\right]-2\left<\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}[k(z_{i},\cdot)],\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}[k(z_{i},\cdot)]\right>+\|\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}[k(z_{i},\cdot)]\|_{\mathcal{H}_{k}}^{2}$
$\textstyle=\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}\left[\|k(z_{i},\cdot)\|_{\mathcal{H}_{k}}^{2}\right]-2\|\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}[k(z_{i},\cdot)]\|_{\mathcal{H}_{k}}^{2}+\|\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}[k(z_{i},\cdot)]\|_{\mathcal{H}_{k}}^{2}$
$\textstyle=\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}[\|k(z_{i},\cdot)\|_{\mathcal{H}_{k}}^{2}]-\|\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}[k(z_{i},\cdot)]\|_{\mathcal{H}_{k}}^{2}$
$\textstyle\leq\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}[\|k(z_{i},\cdot)\|_{\mathcal{H}_{k}}^{2}]$
$\textstyle=\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}[|k(z_{i},z_{i})|]$
(10)

where the last equality follows from the fact that

$\|k(z_{i},\cdot)\|_{\mathcal{H}_{k}}^{2}=|\left<k(z_{i},\cdot),k(\cdot,z_{i})\right>|=|k(z_{i},z_{i})|$

using the reproducing property of the RKHS which says that $\forall
f\in\mathcal{H}_{k}\left<f,k(x,\cdot)\right>_{\mathcal{H}_{k}}=f(x)$. We then
obtain:

$\textstyle\operatorname{\mathbb{E}}_{w_{1:\infty}\sim\operatorname{\operatorname{GEM}}(\alpha^{\prime})}\left[\operatorname{\mathbb{E}}_{z_{1:\infty}\sim\operatorname{\mathbb{F}}^{\prime}}\left[\operatorname{\operatorname{MMD}}^{2}(\operatorname{\mathbb{P}},\operatorname{\mathbb{F}}^{\prime})\right]\right]$
$\textstyle\leq\sum_{i=1}^{\infty}\operatorname{\mathbb{E}}_{w_{1:\infty}\sim\operatorname{\operatorname{GEM}}(\alpha^{\prime})}[w_{i}^{2}]\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}[\|k(z_{i},\cdot)-\mu_{\operatorname{\mathbb{F}}^{\prime}}\|_{\mathcal{H}_{k}}^{2}]+2\sum_{i\neq
j}0$
$\textstyle\leq\sum_{i=1}^{\infty}\operatorname{\mathbb{E}}_{w_{1:\infty}\sim\operatorname{\operatorname{GEM}}(\alpha^{\prime})}[w_{i}^{2}]\operatorname{\mathbb{E}}_{z_{i}\sim\operatorname{\mathbb{F}}^{\prime}}[|k(z_{i},z_{i})|]$
$\textstyle\leq\sum_{i=1}^{\infty}\operatorname{\mathbb{E}}_{w_{1:\infty}\sim\operatorname{\operatorname{GEM}}(\alpha^{\prime})}[w_{i}^{2}].$

The first inequality above follows from equation (9), the second inequality
follows from equation (10), and the last inequality follows from the
boundedness of the kernel in standing assumption 1.

From the properties of the GEM and Beta distributions we have that since
$w_{k}=\beta_{k}\prod_{i=1}^{k-1}(1-\beta_{i})$ where
$\beta_{1:\infty}\operatorname{\overset{iid}{\sim}}\text{Beta}(1,\alpha+n)$
then

$\textstyle\operatorname{\mathbb{E}}_{w_{1:\infty}\sim\operatorname{\operatorname{GEM}}(\alpha^{\prime})}[w_{i}]=\frac{(\alpha+n)^{i-1}}{(1+\alpha+n)^{i}}\qquad\text{
and
}\qquad\operatorname{\mathbb{E}}_{w_{1:\infty}\sim\operatorname{\operatorname{GEM}}(\alpha^{\prime})}[w_{i}^{2}]=\frac{2(\alpha+n)^{i-1}}{(\alpha+n+2)^{i}(\alpha+n+1)}.$

Therefore

$\textstyle\operatorname{\mathbb{E}}_{w_{1:\infty}\sim\operatorname{\operatorname{GEM}}(\alpha^{\prime})}\left[\operatorname{\mathbb{E}}_{z_{1:\infty}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}\left[\operatorname{\operatorname{MMD}}^{2}(\operatorname{\mathbb{P}},\operatorname{\mathbb{F}}^{\prime})\right]\right]$
$\textstyle\leq\sum_{i=1}^{\infty}\operatorname{\mathbb{E}}_{w_{1:\infty}\sim\operatorname{\operatorname{GEM}}(\alpha^{\prime})}[w_{i}^{2}]$
$\textstyle=\frac{2}{(\alpha+n+2)(\alpha+n+1)}\sum_{i=1}^{\infty}\left(\frac{\alpha+n}{\alpha+n+2}\right)^{i-1}$
$\textstyle=\frac{2}{(\alpha+n+2)(\alpha+n+1)}\frac{1}{1-\frac{\alpha+n}{\alpha+n+2}}$
$\textstyle=\frac{1}{\alpha+n+2}$

Finally, by Jensen’s inequality,

$\textstyle\operatorname{\mathbb{E}}_{w_{1:\infty}\sim\operatorname{\operatorname{GEM}}(\alpha^{\prime})}\left[\operatorname{\mathbb{E}}_{z_{1:\infty}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}\left[\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}},\operatorname{\mathbb{F}}^{\prime})\right]\right]\leq\sqrt{\operatorname{\mathbb{E}}_{w_{1:\infty}\sim\operatorname{\operatorname{GEM}}(\alpha^{\prime})}\left[\operatorname{\mathbb{E}}_{z_{1:\infty}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}\left[\operatorname{\operatorname{MMD}}^{2}(\operatorname{\mathbb{P}},\operatorname{\mathbb{F}}^{\prime})\right]\right]}\leq\frac{1}{\sqrt{\alpha+n+2}}.$

∎

######  Lemma 8.

Let $\mathbb{P}_{n}$ denote the empirical measure
$\mathbb{P}_{n}=\frac{1}{n}\sum_{i=1}^{n}\delta_{x_{i}}$ then

$\textstyle\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{z_{1:\infty}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}\left[\operatorname{\operatorname{MMD}}(\mathbb{P}_{n},\operatorname{\mathbb{F}}^{\prime})\right]\right]\leq\frac{2\alpha}{\alpha+n}.$

###### Proof.

By definition of $\operatorname{\mathbb{F}}^{\prime}$ and linearity of
expectations, it follows that
$\operatorname{\mathbb{E}}_{\operatorname{\mathbb{F}}^{\prime}}[\cdot]=\frac{\alpha}{\alpha+n}\operatorname{\mathbb{E}}_{\operatorname{\mathbb{F}}}[\cdot]+\frac{n}{n+\alpha}\operatorname{\mathbb{E}}_{\mathbb{P}_{n}}[\cdot]$
and hence

$\textstyle\mu_{\operatorname{\mathbb{F}}^{\prime}}=\operatorname{\mathbb{E}}_{\operatorname{\mathbb{F}}^{\prime}}[k(z,\cdot)]=\frac{\alpha}{\alpha+n}\operatorname{\mathbb{E}}_{\operatorname{\mathbb{F}}}[k(z,\cdot)]+\frac{n}{n+\alpha}\operatorname{\mathbb{E}}_{\mathbb{P}_{n}}[k(z,\cdot)]=\frac{\alpha}{\alpha+n}\mu_{\operatorname{\mathbb{F}}}+\frac{n}{n+\alpha}\mu_{\mathbb{P}_{n}}.$
(11)

Therefore,

$\textstyle\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{z_{1:\infty}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}\left[\operatorname{\operatorname{MMD}}^{2}(\mathbb{P}_{n},\operatorname{\mathbb{F}}^{\prime})\right]\right]$
$\textstyle=\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{z_{1:\infty}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}\left[\|\mu_{\mathbb{P}_{n}}-\mu_{\operatorname{\mathbb{F}}^{\prime}}\|^{2}_{\mathcal{H}_{k}}\right]\right]$
$\textstyle=\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{z_{1:\infty}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}\left[\left\|\mu_{\mathbb{P}_{n}}-\frac{\alpha}{\alpha+n}\mu_{\operatorname{\mathbb{F}}}-\frac{n}{n+\alpha}\mu_{\mathbb{P}_{n}}\right\|^{2}_{\mathcal{H}_{k}}\right]\right]$
$\textstyle=\frac{\alpha^{2}}{(\alpha+n)^{2}}\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{z_{1:\infty}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}\left[\|\mu_{\mathbb{P}_{n}}-\mu_{\operatorname{\mathbb{F}}}\|^{2}_{\mathcal{H}_{k}}\right]\right]$
$\textstyle=\frac{\alpha^{2}}{(\alpha+n)^{2}}\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{z_{1:\infty}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}\left[\left\|\frac{1}{n}\sum_{i=1}^{n}k(x_{i},\cdot)\right\|^{2}_{\mathcal{H}_{k}}-\frac{2}{n}\sum_{i=1}^{n}\left<k(x_{i},\cdot),\operatorname{\mathbb{E}}_{z\sim\operatorname{\mathbb{F}}}[k(z,\cdot)]\right>+\left\|\mu_{\operatorname{\mathbb{F}}}\right\|^{2}_{\mathcal{H}_{k}}\right]\right]$
$\textstyle=\frac{\alpha^{2}}{(\alpha+n)^{2}}\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{z_{1:\infty}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}\left[\frac{1}{n^{2}}\sum_{i=1}^{n}k(x_{i},x_{i})+\frac{2}{n^{2}}\sum_{i\neq
j}k(x_{i},x_{j})-\frac{2}{n}\sum_{i=1}^{n}\operatorname{\mathbb{E}}_{z\sim\operatorname{\mathbb{F}}}[k(x_{i},z)]\right.\right.$
$\textstyle\quad\quad\left.\left.+\|\operatorname{\mathbb{E}}_{z\sim\operatorname{\mathbb{F}}}[k(z,\cdot)]\|^{2}_{\mathcal{H}_{k}}\right]\right]$
$\textstyle=\frac{\alpha^{2}}{(\alpha+n)^{2}}\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{z_{1:\infty}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}\left[\frac{1}{n^{2}}\sum_{i=1}^{n}k(x_{i},x_{i})+\frac{2}{n^{2}}\sum_{i\neq
j}k(x_{i},x_{j})-\frac{2}{n}\operatorname{\mathbb{E}}_{z\sim\operatorname{\mathbb{F}}}\left[\sum_{i=1}^{n}k(x_{i},z)\right]\right.\right.$
$\textstyle\quad\quad\left.\left.+\operatorname{\mathbb{E}}_{z,z^{\prime}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}}[k(z,z^{\prime})]\right]\right]$
$\textstyle\leq\frac{\alpha^{2}}{(\alpha+n)^{2}}\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{z_{1:\infty}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}\left[\frac{1}{n}+\frac{n(n-1)}{n^{2}}+2+1\right]\right]$
$\textstyle=\frac{4\alpha^{2}}{(\alpha+n)^{2}}.$

The first three equalities above follow from the definition of the MMD in
terms of the kernel mean embeddings and equation (11). The following three
equalities use the definition of the kernel mean embedding as well as the
characteristic property which ensures that

$\textstyle\|\operatorname{\mathbb{E}}_{z\sim\operatorname{\mathbb{F}}}[k(z,\cdot)]\|^{2}_{\mathcal{H}_{k}}$
$\textstyle=\left<\operatorname{\mathbb{E}}_{z\sim\operatorname{\mathbb{F}}}k(z,\cdot),\operatorname{\mathbb{E}}_{z\sim\operatorname{\mathbb{F}}}k(z,\cdot)\right>_{\mathcal{H}_{k}}$
$\textstyle=\left<\operatorname{\mathbb{E}}_{z\sim\operatorname{\mathbb{F}}}k(z,\cdot),\operatorname{\mathbb{E}}_{z^{\prime}\sim\operatorname{\mathbb{F}}}k(\cdot,z^{\prime})\right>_{\mathcal{H}_{k}}$
$\textstyle=\operatorname{\mathbb{E}}_{z,z^{\prime}\sim\operatorname{\mathbb{F}}}\left<k(z,\cdot),k(\cdot,z^{\prime})\right>_{\mathcal{H}_{k}}$
$\textstyle=\operatorname{\mathbb{E}}_{z,z^{\prime}\sim\operatorname{\mathbb{F}}}[k(z,z^{\prime})].$

For the inequality above we used again the standing assumption 1 about the
boundedness of the kernel. To conclude, by Jensen’s inequality we have

$\textstyle\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{z_{1:\infty}\sim\operatorname{\mathbb{F}}^{\prime}}\left[\operatorname{\operatorname{MMD}}(\mathbb{P}_{n},\operatorname{\mathbb{F}}^{\prime})\right]\right]\leq\sqrt{\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{z_{1:\infty}\sim\operatorname{\mathbb{F}}^{\prime}}\left[\operatorname{\operatorname{MMD}}^{2}(\mathbb{P}_{n},\operatorname{\mathbb{F}}^{\prime})\right]\right]}\leq\frac{2\alpha}{\alpha+n}.$

∎

We can now formulate and prove the generalisation error bound below.

######  Theorem 9.

Assume
$x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}$
and let $\operatorname{\mathbb{P}}$ be a sample from the exact DP posterior
with law $\operatorname{\nu}$. Then

$\textstyle\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\nu}}\left[\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})})\right]\right]\leq\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})+\frac{2}{\sqrt{n}}+\frac{2}{\sqrt{\alpha+n+2}}+\frac{4\alpha}{\alpha+n}.$

###### Proof.

Using the triangle inequality we have that for any $\theta\in\Theta$:

$\textstyle\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})})$
$\textstyle\leq\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})},\operatorname{\mathbb{P}})+\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}})$
$\textstyle\leq\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\mathbb{P}})+\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}})$
$\textstyle\leq\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})+2\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}})$
$\textstyle\leq\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})+2\operatorname{\operatorname{MMD}}(\mathbb{P}_{n},\operatorname{\operatorname{\mathbb{P}}^{\ast}})+2\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}},\mathbb{P}_{n})$
$\textstyle\leq\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})+2\operatorname{\operatorname{MMD}}(\mathbb{P}_{n},\operatorname{\operatorname{\mathbb{P}}^{\ast}})+2\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}},\operatorname{\mathbb{F}}^{\prime})+2\operatorname{\operatorname{MMD}}(\mathbb{P}_{n},\operatorname{\mathbb{F}}^{\prime}).$

For the first step above we used the triangle inequality and for the second we
used the fact that since
$\theta^{\ast}=\arg\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}}_{\theta})$
it follows that
$\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})},\operatorname{\mathbb{P}})\leq\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}_{\theta}},\operatorname{\mathbb{P}})$
for any $\theta\in\Theta$. The remaining inequalities are obtained by
reapplying the triangle inequality on
$\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}_{\theta}},\operatorname{\mathbb{P}})$,
$\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}})$
and
$\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}},\mathbb{P}_{n})$
respectively. Since the above is true for any $\theta\in\Theta$ it follows
that

$\textstyle\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})})\leq\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})+2\operatorname{\operatorname{MMD}}(\mathbb{P}_{n},\operatorname{\operatorname{\mathbb{P}}^{\ast}})+2\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}},\operatorname{\mathbb{F}}^{\prime})+2\operatorname{\operatorname{MMD}}(\mathbb{P}_{n},\operatorname{\mathbb{F}}^{\prime}).$

Taking expectations on both sides we obtain

$\textstyle\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{w_{1:\infty}\sim\operatorname{\operatorname{GEM}}(\alpha^{\prime})}\left[\operatorname{\mathbb{E}}_{z_{1:\infty}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}\left[\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}},\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})})\right]\right]\right]\leq$
$\textstyle\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})+2\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\operatorname{MMD}}(\mathbb{P}_{n},\operatorname{\operatorname{\mathbb{P}}^{\ast}})\right]$
$\textstyle+2\operatorname{\mathbb{E}}_{w_{1:\infty}\sim\operatorname{\operatorname{GEM}}(\alpha^{\prime})}\left[\operatorname{\mathbb{E}}_{z_{1:\infty}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}\left[\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}},\operatorname{\mathbb{F}}^{\prime})\right]\right]$
$\textstyle+2\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{z_{1:\infty}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}^{\prime}}\left[\operatorname{\operatorname{MMD}}(\mathbb{P}_{n},\operatorname{\mathbb{F}}^{\prime})\right]\right].$

The result now follows by Lemmas 6, 7 and 8.

∎

The well-specified case is an immediate consequence of the Theorem:

######  Corollary 10.

Suppose that the model is well-specified, i.e.
$\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})=0$
then

$\textstyle\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\nu}}\left[\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})})\right]\right]\leq\frac{2}{\sqrt{n}}+\frac{2}{\sqrt{\alpha+n+2}}+\frac{4\alpha}{\alpha+n}.$

###### Proof.

The corollary follows directly from Theorem 9 by setting
$\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})=0$.
∎

###  B.3 Theorem 1

We start with a lemma bounding the expected MMD between the empirical measure
and the approximated DP posterior sample. Recall that the probability measure
$\operatorname{\hat{\nu}}$ refers to the law on $\operatorname{\mathcal{P}}$
induced by the sampling process:

$\textstyle(w_{1:n},\tilde{w}_{n+T})$
$\textstyle\sim\operatorname{\operatorname{Dir}}\left(1,\dots,1,\frac{\alpha}{T},\dots,\frac{\alpha}{T}\right)$
(12) $\textstyle\tilde{x}_{1:T}$
$\textstyle\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}$ (13)
$\textstyle\operatorname{\mathbb{P}}:=\sum_{i=1}^{n}w_{i}\delta_{x_{i}}+\sum_{k=1}^{T}\tilde{w}_{k}\delta_{\tilde{x}_{k}}$
$\textstyle\sim\operatorname{\hat{\nu}}.$

For clarity we use the individual expectations
$\operatorname{\mathbb{E}}_{(w_{1:n},\tilde{w}_{n+T})\sim\operatorname{\operatorname{Dir}}\left(1,\dots,1,\frac{\alpha}{T},\dots,\frac{\alpha}{T}\right)}$,
which we denote by
$\operatorname{\mathbb{E}}_{w\sim\operatorname{\operatorname{Dir}}}$, and
$\operatorname{\mathbb{E}}_{\tilde{x}_{1:T}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{F}}}$
arising from Equations (12)-(13) in the proofs below—rather than a single
expectation over $\operatorname{\hat{\nu}}$.

######  Lemma 11.

Let $\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}$ denote a single
draw from the approximated DP posterior, i.e.

$\textstyle\operatorname{\mathbb{P}}=\sum_{i=1}^{n}w_{i}\delta_{x_{i}}+\sum_{k=1}^{T}\tilde{w}_{k}\delta_{\tilde{x}_{k}}.$

Then

$\textstyle\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}},\mathbb{P}_{n})\right]\right]\leq
2\sqrt{\frac{\alpha(1+\alpha)}{(\alpha+n)(\alpha+n+1)}}$

###### Proof.

The mean embedding for $\operatorname{\mathbb{P}}$ is
$\mu_{\operatorname{\mathbb{P}}}=\sum_{i=1}^{n}w_{i}k(x_{i},\cdot)+\sum_{j=1}^{T}\tilde{w}_{j}k(\tilde{x}_{j},\cdot)$.
Hence using the triangle inequality:

$\textstyle\operatorname{\operatorname{MMD}}\left(\operatorname{\mathbb{P}},\mathbb{P}_{n}\right)$
$\textstyle=\|\mu_{\operatorname{\mathbb{P}}}-\mu_{\mathbb{P}_{n}}\|_{\mathcal{H}_{k}}$
$\textstyle=\left\|\sum_{i=1}^{n}w_{i}k(x_{i},\cdot)+\sum_{j=1}^{T}\tilde{w}_{j}k(\tilde{x}_{j},\cdot)-\mu_{\mathbb{P}_{n}}\right\|_{\mathcal{H}_{k}}$
$\textstyle\leq\left\|\sum_{i=1}^{n}w_{i}k(x_{i},\cdot)-\mu_{\mathbb{P}_{n}}\right\|_{\mathcal{H}_{k}}+\left\|\sum_{j=1}^{T}\tilde{w}_{j}k(\tilde{x}_{j},\cdot)\right\|_{\mathcal{H}_{k}}.$
(14)

To bound
$\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\|\sum_{i=1}^{n}w_{i}k(x_{i},\cdot)-\mu_{\mathbb{P}_{n}}\|_{\mathcal{H}_{k}}\right]\right]$,
first we note that from the moments of a Dirichlet distribution we have that
for any $i,j\in\\{1,\dots,n\\}$:

$\textstyle\operatorname{\mathbb{E}}_{w\sim\operatorname{\operatorname{Dir}}}[w_{i}]=\frac{1}{n+\alpha},\qquad\operatorname{\mathbb{E}}_{w\sim\operatorname{\operatorname{Dir}}}[w_{i}^{2}]=\frac{2}{(\alpha+n+1)(\alpha+n)}\qquad\text{
and
}\qquad\operatorname{\mathbb{E}}_{w\sim\operatorname{\operatorname{Dir}}}[w_{i}w_{j}]$
$\textstyle=\frac{1}{(\alpha+n)(\alpha+n+1)}.$ (15)

Hence

$\textstyle\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\left\|\frac{1}{n}\sum_{i=1}^{n}k(x_{i},\cdot)-\sum_{i=1}w_{i}k(x_{i},\cdot)\right\|_{\mathcal{H}_{k}}^{2}\right]\right]$
$\textstyle=\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\left\|\sum_{i=1}^{n}\left(\frac{1}{n}-w_{i}\right)k(x_{i},\cdot)\right\|_{\mathcal{H}_{k}}^{2}\right]\right]$
$\textstyle=\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\sum_{i=1}^{n}\left(\frac{1}{n}-w_{i}\right)^{2}|k(x_{i},x_{i})|+2\sum_{1\leq
i<j\leq
n}\left(\frac{1}{n}-w_{i}\right)\left(\frac{1}{n}-w_{j}\right)k(x_{i},x_{j})\right]\right]$
$\textstyle\leq\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\sum_{i=1}^{n}\left(\frac{1}{n}-w_{i}\right)^{2}\right]+2\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\sum_{1\leq
i<j\leq
n}\left(\frac{1}{n}-w_{i}\right)\left(\frac{1}{n}-w_{j}\right)k(x_{i},x_{j})\right]\right]$
$\textstyle=\sum_{i=1}^{n}\left[\frac{1}{n^{2}}-\frac{2}{n}\operatorname{\mathbb{E}}_{w\sim\operatorname{\operatorname{Dir}}}(w_{i})+\operatorname{\mathbb{E}}_{w\sim\operatorname{\operatorname{Dir}}}(w_{i}^{2})\right]$
$\textstyle\quad\quad+2\sum_{1\leq i<j\leq
n}\left[\left(\frac{1}{n^{2}}-\frac{1}{n}\operatorname{\mathbb{E}}_{w\sim\operatorname{\operatorname{Dir}}}(w_{i})-\frac{1}{n}\operatorname{\mathbb{E}}_{w\sim\operatorname{\operatorname{Dir}}}(w_{j})+\operatorname{\mathbb{E}}_{w\sim\operatorname{\operatorname{Dir}}}(w_{i}w_{j})\right)\mathbb{E}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}(k(x_{i},x_{j}))\right]$
$\textstyle=\frac{1}{n}-\frac{2n}{n(n+\alpha)}+\frac{2n}{(n+\alpha)(n+\alpha+1)}+2\sum_{1\leq
i<j\leq
n}\left[\left(\frac{1}{n^{2}}-\frac{2}{n(n+\alpha)}+\frac{1}{(n+\alpha)(n+\alpha+1)}\right)\mathbb{E}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}(k(x_{i},x_{j}))\right]$
$\textstyle=\frac{n^{2}-n+\alpha^{2}+\alpha}{n(n+\alpha)(n+\alpha+1)}+2\sum_{1\leq
i<j\leq
n}\left[\frac{\alpha^{2}+\alpha-n}{n^{2}(n+\alpha)(n+\alpha+1)}\mathbb{E}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}(k(x_{i},x_{j}))\right]$
$\textstyle\leq\frac{n^{2}-n+\alpha^{2}+\alpha}{n(n+\alpha)(n+\alpha+1)}+2\sum_{1\leq
i<j\leq
n}\left[\frac{|\alpha^{2}+\alpha-n|}{n^{2}(n+\alpha)(n+\alpha+1)}\right]$
$\textstyle\leq\frac{n^{2}-n+\alpha^{2}+\alpha}{n(n+\alpha)(n+\alpha+1)}+2\sum_{1\leq
i<j\leq
n}\left[\frac{\alpha^{2}+\alpha+n}{n^{2}(n+\alpha)(n+\alpha+1)}\right]$
$\textstyle=\frac{2(n-1)+\alpha(\alpha+1)}{(n+\alpha)(n+\alpha+1)}$

Here the inequalities follow again from standing assumption 1 and the result
is obtained by using the expectations in (15). Hence, by Jensen’s inequality,

$\textstyle\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\left\|\sum_{i=1}^{n}w_{i}k(x_{i},\cdot)-\mu_{\mathbb{P}_{n}}\right\|_{\mathcal{H}_{k}}\right]\right]\leq\sqrt{\frac{2(n-1)+\alpha(\alpha+1)}{(\alpha+n)(n+\alpha+1)}}.$
(16)

Similarly due to standing assumption 1, to bound
$\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}[\|\sum_{j=1}^{T}\tilde{w}_{j}k(\tilde{x}_{j},\cdot)\|_{\mathcal{H}_{k}}]$,
we note that

$\textstyle\left\|\sum_{j=1}^{T}\tilde{w}_{j}k(\tilde{x}_{j},\cdot)\right\|_{\mathcal{H}_{k}}^{2}$
$\textstyle=\sum_{j=1}^{T}\tilde{w}_{j}^{2}k(\tilde{x}_{j},\tilde{x}_{j})+2\sum_{1\leq
j<k\leq T}\tilde{w}_{j}\tilde{w}_{j}k(\tilde{x}_{j},\tilde{x}_{k})$
$\textstyle\leq\sum_{j=1}^{T}\tilde{w}_{j}^{2}+2\sum_{1\leq j<k\leq
T}\tilde{w}_{j}\tilde{w}_{k}$

hence

$\textstyle\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\|\sum_{j=1}^{T}\tilde{w}_{j}k(\tilde{x}_{j},\cdot)\|^{2}_{\mathcal{H}_{k}}\right]$
$\textstyle\leq\operatorname{\mathbb{E}}_{w\sim\operatorname{\operatorname{Dir}}}\left[\sum_{j=1}^{T}\tilde{w}_{j}^{2}\right]+\operatorname{\mathbb{E}}_{w\sim\operatorname{\operatorname{Dir}}}\left[2\sum_{1\leq
j<k\leq T}\tilde{w}_{j}\tilde{w}_{k}\right].$ (17)

Now from the moments of a Dirichlet distribution we have that for any
$j,k\in\\{1,\dots,T\\}$:

$\textstyle\operatorname{\mathbb{E}}_{w\sim\operatorname{\operatorname{Dir}}}[\tilde{w}_{j}]=\frac{\alpha}{T(n+\alpha)},\qquad\operatorname{\mathbb{E}}_{w\sim\operatorname{\operatorname{Dir}}}[\tilde{w}_{j}^{2}]=\frac{\alpha^{2}+T\alpha}{T^{2}(\alpha+n+1)(\alpha+n)}\qquad\text{
and
}\qquad\operatorname{\mathbb{E}}_{w\sim\operatorname{\operatorname{Dir}}}[\tilde{w}_{j}\tilde{w}_{k}]$
$\textstyle=\frac{\alpha^{3}+\alpha^{2}n}{T^{2}(\alpha+n)^{2}(\alpha+n+1)}.$

Substituting these in equation (17) we obtain

$\textstyle\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\left\|\sum_{j=1}^{T}\tilde{w}_{j}k(\tilde{x}_{j},\cdot)\right\|^{2}_{\mathcal{H}_{k}}\right]$
$\textstyle\leq\frac{\alpha^{2}+T\alpha}{T(\alpha+n+1)(\alpha+n)}+(T-1)\frac{\alpha^{3}+\alpha^{2}n}{T(\alpha+n)^{2}(\alpha+n+1)}$
$\textstyle=\frac{T(\alpha+n)(\alpha+1)\alpha}{T(\alpha+n+1)(\alpha+n)}$
$\textstyle=\frac{\alpha(1+\alpha)}{(\alpha+n)(\alpha+n+1)}.$

Hence, by Jensen’s inequality

$\textstyle\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\left\|\sum_{j=1}^{T}\tilde{w}_{j}k(\tilde{x}_{j},\cdot)\right\|_{\mathcal{H}_{k}}\right]\leq\sqrt{\frac{\alpha(1+\alpha)}{(\alpha+n)(\alpha+n+1)}}$
(18)

for some fixed $\alpha$. Therefore, by substituting equations (16) and (18) in
equation (14) we obtain the required result:

$\textstyle\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}},\mathbb{P}_{n})\right]\right]$
$\textstyle\leq\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\left\|\sum_{i=1}^{n}w_{i}k(x_{i},\cdot)-\mu_{\mathbb{P}_{n}}\right\|_{\mathcal{H}_{k}}\right]+\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\left\|\sum_{j=1}^{T}\tilde{w}_{j}k(\tilde{x}_{j},\cdot)\right\|_{\mathcal{H}_{k}}\right]$
$\textstyle\leq\sqrt{\frac{2(n-1)+\alpha(\alpha+1)}{(\alpha+n)(n+\alpha+1)}}+\sqrt{\frac{\alpha(1+\alpha)}{(\alpha+n)(\alpha+n+1)}}.$

∎

We now proceed with the proof of the generalisation error in Theorem 1.

###### Proof.

We have that for any $\theta\in\Theta$:

$\textstyle\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})})$
$\textstyle\leq\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})},\operatorname{\mathbb{P}})+\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}})$
$\textstyle\leq\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\mathbb{P}})+\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}})$
$\textstyle\leq\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})+2\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}})$
$\textstyle\leq\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})+2\operatorname{\operatorname{MMD}}(\mathbb{P}_{n},\operatorname{\operatorname{\mathbb{P}}^{\ast}})+2\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}},\mathbb{P}_{n})$

For the first inequality above we used the triangle inequality and for the
second inequality we used the fact that
$\theta^{\ast}(\operatorname{\mathbb{P}})=\arg\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}},\operatorname{\mathbb{P}}_{\theta})$.
The last two inequalities follow again by repeatedly applying the triangle
inequality to
$\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}_{\theta}},\operatorname{\mathbb{P}})$
and
$\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}})$.
The above statements hold for any $\theta\in\Theta$ hence we have

$\textstyle\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})})\leq\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})+2\operatorname{\operatorname{MMD}}(\mathbb{P}_{n},\operatorname{\operatorname{\mathbb{P}}^{\ast}})+2\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\mathbb{P}_{n}).$

Taking double expectation on both sides we obtain

$\textstyle\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})})\right]\right]\leq$
$\textstyle\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})+2\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\operatorname{MMD}}(\mathbb{P}_{n},\operatorname{\operatorname{\mathbb{P}}^{\ast}})\right]$
$\textstyle+2\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}},\mathbb{P}_{n})\right]\right].$

The result then follows from Lemmas 6 and 11. ∎

The well-specified case is an immediate consequence of the Theorem:

######  Corollary 12.

Suppose that the model is well-specified, i.e.
$\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})=0$
then

$\textstyle\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})})\right]\right]\leq\frac{2}{\sqrt{n}}+2\sqrt{\frac{2(n-1)+\alpha(\alpha+1)}{(\alpha+n)(n+\alpha+1)}}+2\sqrt{\frac{\alpha(1+\alpha)}{(\alpha+n)(\alpha+n+1)}}.$

###### Proof.

The corollary follows directly from Theorem 1 by setting
$\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})=0$.
∎

###  B.4 Theorem 4

We can prove that the posterior consistency result indeed holds by adapting
the arguments in [20] (Theorem 2) from equation (6) to the required statement
of Theorem 4.

###### Proof.

Using Proposition 1 and Markov’s inequality we have:

$\textstyle\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\hat{\nu}}\left(\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})},\operatorname{\operatorname{\mathbb{P}}^{\ast}})-\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})>M_{n}\cdot
n^{-1/2}\right)\right]$
$\textstyle\leq\frac{\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})},\operatorname{\operatorname{\mathbb{P}}^{\ast}})-C\right]\right]}{M_{n}\cdot
n^{-1/2}}=\frac{\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})},\operatorname{\operatorname{\mathbb{P}}^{\ast}})\right]\right]-C}{M_{n}\cdot
n^{-1/2}}$ $\textstyle\leq\frac{2n^{-1/2}}{M_{n}\cdot
n^{-1/2}}+\frac{2\sqrt{2n(n-1)+\alpha
n(\alpha+1)}}{M_{n}\sqrt{(\alpha+n)(\alpha+n+1)}}+\frac{2\sqrt{n\alpha(1+\alpha)}}{M_{n}\sqrt{(\alpha+n)(\alpha+n+1)}}\overset{n\to\infty}{\longrightarrow}0.$

∎

Again, the well-specified case is a special case of Theorem 4 for $C=0$.

######  Corollary 13.

Suppose we have a well-specified model, i.e.
$\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})=0$
and let
$\operatorname{\mathbb{P}}=\sum_{i=1}^{n}w_{i}\delta_{x_{i}}+\sum_{k=1}^{T}\tilde{w}_{k}\delta_{\tilde{x}_{k}}\sim\operatorname{\hat{\nu}}$
be a sample from the approximate posterior. Then for any
$M_{n}\rightarrow+\infty$ such that $M_{n}\cdot n^{-1/2}\rightarrow 0$ as
$n\rightarrow\infty$

$\textstyle\operatorname{\hat{\nu}}\biggl{(}\operatorname{\mathbb{P}}\in\operatorname{\mathcal{P}}:$
$\textstyle\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})})>\frac{M_{n}}{n^{1/2}}\biggr{)}\overset{n\to\infty}{\longrightarrow}0.$

###### Proof.

Using Theorem 1 in the well specified case and Markov’s inequality we have
that

$\textstyle\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\hat{\nu}}\left(\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})})>M_{n}\cdot
n^{-1/2}\right)\right]$
$\textstyle\leq\frac{\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})})\right]\right]}{M_{n}\cdot
n^{-1/2}}$ $\textstyle\leq\frac{2n^{-1/2}}{M_{n}\cdot
n^{-1/2}}+\frac{2\sqrt{2n(n-1)+\alpha
n(\alpha+1)}}{M_{n}\sqrt{(\alpha+n)(\alpha+n+1)}}+\frac{2\sqrt{n\alpha(1+\alpha)}}{M_{n}\sqrt{(\alpha+n)(\alpha+n+1)}}\overset{n\to\infty}{\longrightarrow}0.$

∎

###  B.5 Corollary 5

Before we present the generalisation bound we will use the following lemma
from [19] which bounds the absolute difference of the MMD between the
parametric model and each of the contaminated and non-contaminated probability
measures.

######  Lemma 14 ([19] Lemma 3.3).

For any $\theta\in\Theta$,

$\textstyle|\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})-\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\tilde{\operatorname{\mathbb{P}}})|\leq
2\epsilon.$

Using this lemma and the previously obtained generalisation bounds we can show
the required statement of Corollary 5.

###### Proof.

From Lemma 14 we have that for any $\theta\in\Theta$:

$\textstyle\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\tilde{\operatorname{\mathbb{P}}})\leq
2\epsilon+\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\mathbb{P}}^{\ast}_{0})$
(19)
$\textstyle\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\operatorname{\mathbb{P}}^{\ast}})\leq
2\epsilon+\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\tilde{\operatorname{\mathbb{P}}})$
(20)

Hence, using equations 19, 20 and Proposition 1 we have that

$\textstyle\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\operatorname{\operatorname{MMD}}(\tilde{\operatorname{\mathbb{P}}},\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})})\right]\right]$
$\textstyle\leq
2\epsilon+\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})})\right]\right]$
$\textstyle\leq
2\epsilon+\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\operatorname{\operatorname{\mathbb{P}}^{\ast}},\operatorname{\mathbb{P}}_{\theta})+\frac{2}{\sqrt{n}}+2\sqrt{\frac{2(n-1)+\alpha(\alpha+1)}{(\alpha+n)(n+\alpha+1)}}+2\sqrt{\frac{\alpha(1+\alpha)}{(\alpha+n)(\alpha+n+1)}}$
$\textstyle\leq
4\epsilon+\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\tilde{\operatorname{\mathbb{P}}},\operatorname{\mathbb{P}}_{\theta})+\frac{2}{\sqrt{n}}+2\sqrt{\frac{2(n-1)+\alpha(\alpha+1)}{(\alpha+n)(n+\alpha+1)}}+2\sqrt{\frac{\alpha(1+\alpha)}{(\alpha+n)(\alpha+n+1)}}.$

∎

The well-specified case is an immediate consequence:

######  Corollary 15 (Well-specified Case).

Suppose the model is well specified in terms of the target
$\tilde{\operatorname{\mathbb{P}}}$, i.e. $\exists\theta\in\Theta$ such that
$\operatorname{\mathbb{P}}_{\theta}=\tilde{\operatorname{\mathbb{P}}}$. Then

$\textstyle\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\left[\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\left[\operatorname{\operatorname{MMD}}(\tilde{\operatorname{\mathbb{P}}},\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})})\right]\right]\leq
4\epsilon+\frac{2}{\sqrt{n}}+2\sqrt{\frac{2(n-1)+\alpha(\alpha+1)}{(\alpha+n)(n+\alpha+1)}}+2\sqrt{\frac{\alpha(1+\alpha)}{(\alpha+n)(\alpha+n+1)}}.$

###### Proof.

The corollary follows from Corollary 5 noting that in the well-specified case
$\inf_{\theta\in\Theta}\operatorname{\operatorname{MMD}}(\operatorname{\mathbb{P}}_{\theta},\tilde{\operatorname{\mathbb{P}}})=0$.
∎

##  Appendix C EXPERIMENTAL DETAILS

We now provide further experimental details of our method. We first explain
how gradient-based numerical optimisation can be used to minimise the MMD
objective in algorithm 1. Next, we provide the experimental setup necessary
for the reproduction of the experiments presented in the main text. We then
give some additional details on the G-and-k and Toggle-Switch models. We
finally provide results of all experiments over a number of independent runs.

###  C.1 Numerical Optimisation for the MMD Objective

We first explain how the MMD objective in equation (5) can be approximated and
minimized through gradient-based methods. For two probability measures
$\operatorname{\mathbb{P}},\operatorname{\mathbb{Q}}\in\mathcal{P}$ recall
that the squared MMD is defined as

$\textstyle\operatorname{\operatorname{MMD}}^{2}(\operatorname{\mathbb{P}},\operatorname{\mathbb{Q}})=\int_{\mathcal{X}\times\mathcal{X}}$
$\textstyle
k(x,y)\operatorname{\mathbb{P}}(dx)\operatorname{\mathbb{P}}(dy)-2\int_{\mathcal{X}\times\mathcal{X}}k(x,y)\operatorname{\mathbb{P}}(dx)\operatorname{\mathbb{Q}}(dy)$
$\textstyle+\int_{\mathcal{X}\times\mathcal{X}}k(x,y)\operatorname{\mathbb{Q}}(dx)\operatorname{\mathbb{Q}}(dy).$

Since these integrals are usually intractable we can easily approximate this
by sampling
$x_{1:N}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{P}}$ and
$y_{1:m}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{Q}}$ and using
the U-statistic:

$\textstyle\widehat{\operatorname{\operatorname{MMD}}}^{2}(\operatorname{\mathbb{P}},\operatorname{\mathbb{Q}})=\frac{1}{N(N-1)}$
$\textstyle\sum_{i\neq
i^{\prime}}^{N}k(x_{i},x_{i^{\prime}})-\frac{2}{Nm}\sum_{i=1}^{N}\sum_{j=1}^{m}k(x_{i},y_{j})$
$\textstyle+\frac{1}{m(m-1)}\sum_{j\neq j^{\prime}}k(y_{j},y_{j^{\prime}}).$

In our case we are interested in minimising the MMD between the parametric
model
$\operatorname{\mathbb{P}_{\theta}}\in\operatorname{\mathcal{P}_{\Theta}}$
defined through the simulator
$G_{\theta}:\operatorname{\mathbb{U}}\rightarrow\mathcal{X}$ and an
approximated sample from the DP posterior $\operatorname{\mathbb{P}}$. Hence,
for samples
$y_{1:N}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{P}}$ and
$u_{1:M}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{U}}$ the
approximated squared MMD is

$\textstyle\widehat{\operatorname{\operatorname{MMD}}}^{2}(\operatorname{\mathbb{P}},\operatorname{\mathbb{P}_{\theta}})=\frac{1}{N(N-1)}$
$\textstyle\sum_{i\neq
i^{\prime}}^{N}k(y_{i},y_{i^{\prime}})-\frac{2}{NM}\sum_{j=1}^{M}\sum_{i=1}^{N}k(G_{\theta}(u_{j}),y_{i})$
$\textstyle+\frac{1}{M(M-1)}\sum_{j\neq
j^{\prime}}k(G_{\theta}(u_{j}),G_{\theta}(u_{j^{\prime}})).$ (21)

The gradient of equation (21) can also be approximated using a U-statistic as
follows: assuming that the generator is differentiable with respect to
$\theta$ with gradient $\nabla_{\theta}G_{\theta}$ then the U-statistic of the
gradient of the squared MMD is

$\textstyle\hat{J}_{\theta}(u_{1:M},y_{1:N}):=\nabla_{\theta}\widehat{\operatorname{\operatorname{MMD}}}^{2}(\operatorname{\mathbb{P}},\operatorname{\mathbb{P}}_{\theta})=\frac{2}{M(M-1)}$
$\textstyle\sum_{j\neq
j^{\prime}}^{M}\nabla_{\theta}G_{\theta}(u_{j})\nabla_{1}k(G_{\theta}(u_{j}),G_{\theta}(u_{j^{\prime}}))$
$\textstyle-\frac{2}{NM}\sum_{j=1}^{M}\sum_{i=1}^{N}\nabla_{\theta}G_{\theta}(u_{j})\nabla_{1}k(G_{\theta}(u_{j}),y_{i})$

where $\nabla_{1}$ denotes the partial derivative with respect to the first
argument and by noting that the first term in (21) does not depend on
$\theta$. This is an unbiased statistic in the sense that

$\textstyle\operatorname{\mathbb{E}}_{u_{1:M}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{U}}}\left[\operatorname{\mathbb{E}}_{y_{1:N\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{P}}}}[\nabla_{\theta}\widehat{\operatorname{\operatorname{MMD}}}^{2}(\operatorname{\mathbb{P}},\operatorname{\mathbb{P}_{\theta}})]\right]=\nabla_{\theta}\operatorname{\operatorname{MMD}}^{2}(\operatorname{\mathbb{P}},\operatorname{\mathbb{P}_{\theta}}).$

Using this approximation we can use gradient-based methods to minimise the
required objective. For example, one can use stochastic gradient descent as
follows; for learning rate $\eta$, at each iteration $k=1,2,\dots$:

  1. 1.

Sample $u_{1:M}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{U}}$
and $y_{1:N}\operatorname{\overset{iid}{\sim}}\operatorname{\mathbb{P}}$

  2. 2.

Set $\theta_{k}=\theta_{k-1}-\eta\hat{J}_{\theta_{k-1}}(u_{1:M},y_{1:N})$

###  C.2 Experimental Setup

We provide details on the experimental setup required to produce figures 2, 4
and 5. For the WABC method we make use of the winference R package and the
experimental setup provided in [11]. For all experiments with our method we
use JAX [16] to parallelize the bootstrap sampling and perform the
optimisation step.

For all three experiments we use the Adam optimiser [44] for the minimization
of the MMD at each bootstrap estimation. The learning rate is set to
$\eta=0.1$ for the Gaussian location and g-and-k models and to $\eta=0.04$ for
the toggle-switch model. To ensure convergence, we perform 2000 optimizations
steps in the toggle-switch experiment and 1000 steps for the rest.

The prior distributions set in the WABC method serve as an indication for the
initialisation point of the optimization. In the Gaussian location model, a
centered normal prior is imposed on the parameter, hence we initialize the
optimization step at $\theta=0$ at each bootstrap iteration. Similarly, the
g-and-k model uses a Uniform prior on $[0,10]$ for all parameters so we
initialise at the midpoint $(5,5,5,5)$. Finally, for the toggle-switch model
we employ a random restart method as follows; 500 random values are sampled
from the prior uniform distributions set in [14] and [11]. The approximated
MMD loss is computed for all of them and the three initial values with the
smallest loss are chosen as initial points. The optimization is performed for
all three starting points for each bootstrap iteration and the estimator with
the smallest loss is retained.

Finally, the length scale $l$ of the Gaussian kernel is set using the median
heuristic $l=\sqrt{\text{median}_{1\leq 1,j\leq n}(\|x_{i}-x_{j}\|^{2}_{2})}$
suggested in [33] and [25] in the Gaussian location example while for the
g-and-k model we set $l=0.15$. For the toggle-switch model, we suggest an
unweighted mixture of Gaussian kernels discussed in [74], i.e.
$k(x,x^{\prime})=\sum_{i=1}^{I}\exp(-\|x-x^{\prime}\|_{2}^{2}/(2l_{i}^{2}))$
for a range of values $l_{i}\in\\{1,10,20,40,80,100,130,200,400,800,1000\\}$.

###  C.3 The G-and-k Distribution

We visualise the observed data from the G-and-k distribution used in the
contaminated model example of section 5 for an increasing number of outliers
with $\theta_{0}=(3,1,1,-\log(2))$. We further provide the density obtained by
generating 2000 samples from the model using the posterior mean
$\theta=\frac{1}{B}\sum_{j=1}^{B}\theta^{(j)}$ for each method. The outliers
of the contaminated models are visualised on the histograms of figure 6. The
densities of the fitted models show the sensitivity of the two methods to the
two degrees of outliers; the densities fitted by the NPL-MMD remain similar
and close to $\operatorname{\mathbb{P}}_{\theta_{0}}$ whereas the densities
fitted by WABC are affected by the outliers in the misspecifed cases.

Figure 6: Histogram of 2000 samples from the g-and-k distribution with
$\theta_{0}=(3,1,1,-\log(2))$ (blue) and 2000 samples from the fitted models
using the posterior mean (orange) for the NPL-MMD method (above) and the WABC
method (below).

Figure 7 visualises the densities obtained at each bootstrap iteration of the
NPL-MMD method. Here, each density corresponds to 2000 samples from the
g-and-k distribution and parameter $\theta^{(j)}$ for $j=1,\dots,B$. It is
hence clear how each sample $\operatorname{\mathbb{P}}^{(j)}$ from the DP
posterior is mapped through $\theta^{*}(\cdot)$ to a value $\theta^{(j)}$
which gives rise to a different density. As we would expect, there is more
variability in the densities for higher contamination levels.

Figure 7: Realisations of the G-and-k distribution with each $\theta^{(j)}$
where $j=1,\dots,B$ obtained from the bootstrap iteration of the MMD posterior
bootstrap algorithm for an increasing degree of contamination in the dataset.

Finally, we illustrate the rate obtained in Corollary 5 in the absence of
outliers for the G-and-k distribution example. Since we take
$\epsilon=\alpha=0$, the Corollary implies that the expected MMD is bounded
above by $\frac{2}{\sqrt{n}}$. To obtain an approximation of
$\operatorname{\mathbb{E}}_{x_{1:n}\operatorname{\overset{iid}{\sim}}\operatorname{\operatorname{\mathbb{P}}^{\ast}}}\big{[}\operatorname{\mathbb{E}}_{\operatorname{\mathbb{P}}\sim\operatorname{\hat{\nu}}}\big{[}\operatorname{\operatorname{MMD}}^{2}\big{(}\tilde{\operatorname{\mathbb{P}}},\operatorname{\mathbb{P}}_{\theta^{\ast}(\operatorname{\mathbb{P}})}\big{)}\big{]}\big{]}$
we perform 10 runs of the algorithm and take the posterior mean
$\hat{\theta}=\frac{1}{B}\sum_{i=1}^{B}\theta^{(i)}$ as the estimate of
$\theta$. We then sample $15\times 10^{3}$ instances from
$\operatorname{\operatorname{\mathbb{P}}^{\ast}}=\operatorname{\mathbb{P}}_{\theta_{0}}$
and $\operatorname{\mathbb{P}}_{\hat{\theta}}$ and estimate the squared MMD
using the U-statistic in equation (21). We repeat this experiment for ten
values of sample size $n\in[250,4000]$ and plot
$\sqrt{\operatorname{\mathbb{E}}[\widehat{\operatorname{\operatorname{MMD}}}^{2}(\operatorname{\mathbb{P}}_{\hat{\theta}},\operatorname{\operatorname{\mathbb{P}}^{\ast}})]}$
against $\frac{2}{\sqrt{n}}$ in figure 8. We observe that the estimate of the
square root of the expected squared MMD is bounded above by the curve
$\frac{2}{\sqrt{n}}$ and by Jensen’s inequality it is implied that the
estimate of the approximate MMD will also be bounded above since:

$\operatorname{\mathbb{E}}[\widehat{\operatorname{\operatorname{MMD}}}(\operatorname{\mathbb{P}}_{\hat{\theta}},\operatorname{\operatorname{\mathbb{P}}^{\ast}})]\leq\sqrt{\operatorname{\mathbb{E}}[\widehat{\operatorname{\operatorname{MMD}}}^{2}(\operatorname{\mathbb{P}}_{\hat{\theta}},\operatorname{\operatorname{\mathbb{P}}^{\ast}})]}\leq\frac{2}{\sqrt{n}}.$

###  C.4 The Toggle-Switch Model

####  C.4.1 Description of the Simulator

For cell $i$ and unknown parameters
$\theta=(\alpha_{1},\alpha_{2},\beta_{1},\beta_{2},\mu,\sigma,\gamma)^{\top}$,
the simulator input is
$u_{i}=(u_{i,1,1},u_{i,1,2},\ldots,u_{i,T,1},u_{i,T,2},u_{i,T+1,1})^{\top}\sim\text{Unif}([0,1]^{2T+1})$.
The simulator $G_{\theta}$ is defined through:

$\textstyle G_{\theta}(u_{i})$
$\textstyle=\Phi^{-1}\Big{(}\Phi\Big{(}\frac{-(\mu+v_{i,T})v_{i,T}^{\gamma}}{\mu\sigma}\Big{)}+u_{i,T+1,1}\Big{(}1-\Phi\Big{(}\frac{-(\mu+v_{i,T})v_{i,T}^{\gamma}}{\mu\sigma}\Big{)}\Big{)}\Big{)}\frac{\mu\sigma}{v_{i,T}^{\gamma}}+(\mu+v_{i,T})$

where for $t=1,\dots,T-1$, we have

$\textstyle\tilde{v}_{i,t+1}$
$\textstyle=v_{i,t}+\frac{\alpha_{1}}{1+w_{i,t}^{\beta_{1}}}-(1+0.03v_{i,t})$
$\textstyle\tilde{w}_{i,t+1}$
$\textstyle=w_{i,t}+\frac{\alpha_{2}}{1+v_{i,t}^{\beta_{2}}}-(1+0.03w_{i,t})$
$\textstyle v_{i,t+1}$
$\textstyle=\tilde{v}_{i,t+1}+0.5\Phi^{-1}\big{(}\Phi(-2\tilde{v}_{i,t+1})+u_{i,t,1}(1-\Phi(-2\tilde{v}_{i,t+1}))\big{)}$
$\textstyle w_{i,t+1}$
$\textstyle=\tilde{w}_{i,t+1}+0.5\Phi^{-1}\big{(}\Phi(-2\tilde{w}_{i,t+1})+u_{i,t,2}(1-\Phi(-2\tilde{w}_{i,t+1}))\big{)}$

and $\Phi$ denotes the CDF of the standard Gaussian distribution. We use the
initial conditions $v_{i,0}=10$, $w_{i,0}=10$.

Figure 8: Illustration of the upper bound of the generalisation error in
Corollary 5 in the absence of outliers for the G-and-k distribution.

###  C.5 Results Over Multiple Independent Runs

We provide results for a number of independent runs of the experiments in
section 5. For each run, a new dataset was generated and $B$ posterior samples
were obtained for each parameter $\theta$. For each such sample, the mean
estimator is recorded and the normalized mean squared error between the
estimator and the true value of $\theta$ are presented with their standard
deviations in tables 2, 4 and 4 for the Gaussian location, G-and-k and Toggle-
Switch models respectively.

Table 2: Experiment results for the Gaussian model over 10 runs Method | NMSE (std)
---|---