# Detecting Model Misspecification in Amortized Bayesian Inference with Neural
Networks: An Extended Investigation

Marvin Schmitt<EMAIL_ADDRESS>
Cluster of Excellence SimTech, University of Stuttgart, Germany Paul-Christian
Bürkner<EMAIL_ADDRESS>
Department of Statistics, TU Dortmund University, Germany Ullrich Köthe
<EMAIL_ADDRESS>
Visual Learning Lab, Heidelberg University, Germany Stefan T. Radev
<EMAIL_ADDRESS>
Cognitive Science, Rensselaer Polytechnic Institute, United States

###### Abstract

Recent advances in probabilistic deep learning enable efficient amortized
Bayesian inference in settings where the likelihood function is only
implicitly defined by a simulation program (simulation-based inference; SBI).
But how faithful is such inference if the simulation represents reality
somewhat inaccurately, that is, if the true system behavior at test time
deviates from the one seen during training? We conceptualize the types of such
model misspecification arising in SBI and systematically investigate how the
performance of neural posterior approximators gradually deteriorates as a
consequence, making inference results less and less trustworthy. To notify
users about this problem, we propose a new misspecification measure that can
be trained in an unsupervised fashion (i.e., without training data from the
true distribution) and reliably detects model misspecification at test time.
Our experiments clearly demonstrate the utility of our new measure both on toy
examples with an analytical ground-truth and on representative scientific
tasks in cell biology, cognitive decision making, disease outbreak dynamics,
and computer vision. We show how the proposed misspecification test warns
users about suspicious outputs, raises an alarm when predictions are not
trustworthy, and guides model designers in their search for better simulators.

This article is an extended version of the conference paper Detecting Model
Misspecification in Amortized Bayesian Inference with Neural Networks (Schmitt
et al., 2024, Link).

##  1 Introduction

Computer simulations play a fundamental role in many fields of science and
engineering. However, the associated inverse problems of finding simulation
parameters that accurately reproduce or predict real-world behavior are
generally difficult and often analytically intractable. Here, we consider
_simulation-based inference_ (SBI; Cranmer et al., 2020) as a general approach
to address this problem within a Bayesian inference framework. That is, given
an assumed generative model $\mathcal{M}$ (as represented by the simulation
program, see Section 3.2 for details) and observations $\boldsymbol{x}$ (real
or simulated outcomes), we estimate the posterior distribution
$p(\boldsymbol{\theta}\,|\,\boldsymbol{x},\mathcal{M})$ of the simulation
parameters $\boldsymbol{\theta}$ that would reproduce the observed
$\boldsymbol{x}$. The recent introduction of efficient neural network
approximators for this task has inspired a rapidly growing literature on SBI
solutions for various application domains (e.g., Butter et al., 2022;
Lueckmann et al., 2021a; Gonçalves et al., 2020; Shiono, 2021; Bieringer et
al., 2021; von Krause et al., 2022; Ghaderi-Kangavari et al., 2023). These
empirical successes call for a systematic investigation of the trustworthiness
of SBI, see Figure 1.

We conduct an extensive analysis of neural posterior estimation (NPE) and
sequential neural posterior estimation (SNPE), two deep learning approaches to
approximate the posterior distribution
$p(\boldsymbol{\theta}\,|\,\boldsymbol{x},\mathcal{M})$. In particular, we
study their accuracy under model misspecification, where the generative model
$\mathcal{M}^{*}$ at test time (the “true data generating process”) deviates
from the one assumed during training (i.e., $\mathcal{M}^{*}\neq\mathcal{M}$),
a situation commonly known as simulation gap. As a consequence of a simulation
gap, the observed data of interest might lie outside of the simulated data
from the training phase of SBI. Paralleling the notion of “out-of-
distribution” in anomaly detection and representation learning (Yang et al.,
2021), simulation gaps may lead to “out-of-simulation” samples, and ultimately
to silent errors in posterior estimates.

Our experiments illustrate the difference between the two situations (see
Figure 2 for a quick preview): When the simulation model is well-specified,
the posterior estimates of a neural posterior estimator and a gold-standard
MCMC sampler (HMC-MCMC; Stan Development Team, 2022) are essentially equal,
whereas both approaches disagree considerably under model misspecification. In
order to avoid drawing incorrect conclusions from misspecified models and
incorrect posteriors, it is of crucial importance to detect whether
$\mathcal{M}^{*}\neq\mathcal{M}$ and to quantify the severity of the mismatch.
However, this is difficult in practice because the true data generating
process $\mathcal{M}^{*}$ is generally unknown (except in controlled
validation settings).

Training Set$\boldsymbol{x}\sim p(\boldsymbol{x}\,|\,\mathcal{M})$Generative
Model $\mathcal{M}$$\accentset{\text{o}}{\boldsymbol{x}}\sim
p^{*}(\boldsymbol{x})$$\huge{h_{\boldsymbol{\psi}}}$SummaryNetworkSummary
Space$h_{\boldsymbol{\psi}}(\boldsymbol{x})$$h_{\boldsymbol{\psi}}(\accentset{\text{o}}{\boldsymbol{x}})$$\huge{f_{\boldsymbol{\phi}}}$InferenceNetworkCorrectPosteriorIncorrectPosteriorMisspecificationMisspecification
detected

Figure 1: Conceptual overview of our neural approach. The summary network
$h_{\psi}$ maps observations $\boldsymbol{x}$ to summary statistics
$h_{\psi}(\boldsymbol{x})$, and the inference network $f_{\phi}$ estimates the
posterior $p(\boldsymbol{\theta}\,|\,\boldsymbol{x},\mathcal{M})$ from the
summary statistics. The generative model $\mathcal{M}$ creates training data
$\boldsymbol{x}$ in the green region, and the networks learn to map these data
to well-defined summary statistics and posteriors (green regions/dot/box). If
the generative model $\mathcal{M}$ is misspecificed, real observations
$\accentset{\text{o}}{\boldsymbol{x}}$ fall outside the training region and
are therefore mapped to outlying summary statistics and potentially incorrect
posteriors (red dots/box). Since our learning approach enforces a known inlier
summary distribution (e.g., Gaussian), misspecification can be detected by a
distribution mismatch in summary space, as signaled by a high maximum mean
discrepancy score (Gretton et al., 2012).

In this work, we propose a new misspecification measure that can be trained in
an unsupervised fashion (i.e., without knowledge of $\mathcal{M}^{*}$ or
training data from the true data distribution) and reliably quantifies by how
much $\mathcal{M}^{*}$ deviates from $\mathcal{M}$ at test time. Our
experiments clearly demonstrate the power of our new measure both on toy
examples with an analytical ground-truth, and on representative scientific
tasks in cell biology, cognitive decision making, and disease outbreak
dynamics. We show how simulation-based posterior inference gradually
deteriorates as the simulation gap widens and how the proposed
misspecification test warns users about suspicious outputs, raises an alarm
when predictions are not trustworthy, and guides model designers in their
search for better simulators.

Our investigations complement existing work on deep amortized SBI, whose main
focus has been on network architectures and training algorithms for high
accuracy in the well-specified case $\mathcal{M}^{*}=\mathcal{M}$ (Ramesh et
al., 2022; Pacchiardi & Dutta, 2022a; Durkan et al., 2020; Radev et al.,
2020a; Greenberg et al., 2019; Lueckmann et al., 2017; Papamakarios & Murray,
2016), with a steadily growing interest in inference algorithms that are
robust to certain types of model misspecification (Huang et al., 2023; Ward et
al., 2022; Dellaporta et al., 2022; Leclercq, 2022; Kelly et al., 2023;
Gloeckler et al., 2023). In particular, our paper makes the following key
contributions:

  1. (i)

We systematically conceptualize different sources of model misspecification in
amortized Bayesian inference with neural networks and propose a new detection
criterion that is widely applicable to different model structures, inputs, and
outputs.

  2. (ii)

We incorporate this criterion into existing neural posterior estimation
methods, with hand-crafted and learned summary statistics, with sequential or
amortized inference regimes, and we extend the associated learning algorithms
in a largely non-intrusive manner.

  3. (iii)

We conduct a systematic empirical evaluation of our detection criterion, the
influence of the summary space dimension, and the relationship between summary
outliers and posterior distortion under various types and strengths of model
misspecification.

(a) Well-specified: Similar posteriors.

(b) Misspecified: Dissimilar posteriors.

Figure 2: Preview of Experiment 3 on reaction time modeling in psychological
experiments. Posteriors obtained via non-amortized MCMC (HMC in Stan) and
amortized simulation-based inference (NPE in BayesFlow) are very similar when
the model is well-specified (left). However, a simulation gap (here: not
accounting for occasional slow responses due to mind wandering) leads to
considerable disagreement between these methods (right).

##  2 Related Work

Model misspecification has been studied both in the context of standard
Bayesian inference and generalizations thereof (Knoblauch et al., 2019; Schmon
et al., 2021). To alleviate model misspecification in generalized Bayesian
inference, researchers have investigated probabilistic classifiers (Thomas &
Corander, 2019), second-order PAC-Bayes bounds (Masegosa, 2020), scoring rules
(Giummolè et al., 2019), priors over a class of predictive models (Loaiza-Maya
et al., 2021), or Stein discrepancy as a loss function (Matsubara et al.,
2022). Notably, these approaches deviate from the standard Bayesian
formulation and investigate alternative schemes for belief updating and
learning (e.g., replacing the likelihood function with a generic loss
function). In contrast, our method remains grounded in the standard Bayesian
framework embodying an implicit likelihood principle (Berger & Wolpert, 1988;
Cranmer et al., 2020). Differently, power scaling methods incorporate a
modified likelihood (raised to a power $0<\alpha<1)$ in order to prevent
potentially overconfident Bayesian updating (Grünwald et al., 2017; Holmes &
Walker, 2017). However, the SBI setting assumes that the likelihood function
is not available in closed-form, which makes an explicit modification of the
implicitly defined likelihood less obvious.

Neural approaches to amortized SBI can be categorized as either targeting the
posterior (Radev et al., 2020a; Greenberg et al., 2019), the likelihood
(Papamakarios et al., 2019; Hermans et al., 2020), or both (Wiqvist et al.,
2021; Radev et al., 2023a). These methods employ simulations for training
neural approximators which can either generate samples from the posterior
directly (Radev et al., 2020a; Greenberg et al., 2019; Wiqvist et al., 2021)
or in tandem with Markov chain Monte Carlo (MCMC) sampling algorithms
(Papamakarios et al., 2019; Hermans et al., 2020). Since the behavior of these
methods depends on the fidelity of the simulations used as training data, we
hypothesize that their estimation quality will be, in general, unpredictable,
when faced with atypical real-world data. Indeed, the critical impact of model
misspecification in neural SBI has been commonly acknowledged in the
scientific research community (Cannon et al., 2022; Alquier & Ridgway, 2019;
Zhang & Gao, 2020; Frazier et al., 2020; Frazier & Drovandi, 2021; Pacchiardi
& Dutta, 2022b; Ward et al., 2022; Kelly et al., 2023; Gloeckler et al., 2023;
Huang et al., 2024).

Recent approaches to detect model misspecification in simulation-based
inference are usually based on the obtained approximate posterior distribution
(Dellaporta et al., 2022; Leclercq, 2022). However, we explicitly show in
Experiment 1 and Experiment 3 that the approximate posteriors in simulation-
based inference tend to show pathological behavior under misspecified models.
Posteriors from misspecified models may erroneously look legitimate, rendering
diagnostic methods on their basis unreliable. Moreover, the same applies for
approaches based on the posterior predictive distribution (Bürkner et al.,
2020; Gabry et al., 2019; Vehtari & Ojanen, 2012) since these also rely on the
fidelity of the posterior distribution and can therefore only serve as an
indirect measure of misspecification.

A few novel techniques aim to _mitigate_ model misspecification in simulation-
based inference to achieve robust inference. Ward et al. (2022) explore a way
to alleviate model misspecification with two neural approximators and
subsequent MCMC. While both approaches are appealing in theory, the
computational burden of MCMC sampling contradicts the idea of amortized
inference and prohibits their use in complex applications with learned summary
statistics and large amounts of data. In fact, von Krause et al. (2022) used
amortized neural SBI on more than a million data sets of multiple observations
each and demonstrated that an alternative inference method involving non-
amortized MCMC would have taken years of sampling. Huang et al. (2024) propose
a method to learn summary statistics that are robust to unmodeled data
contamination for sequential (i.e., non-amortized) neural posterior
estimation. Finally, Gloeckler et al. (2023) propose to penalize the Fisher
information of neural posterior estimators to increase their robustness to
adversarial attacks (which can be considered a form of subtle model
misspecification).

For robust non-amortized ABC samplers, the possibility of utilizing hand-
crafted summary statistics as an important element of misspecification
analysis has already been explored (Frazier et al., 2020; Frazier & Drovandi,
2021). Our work parallels these ideas and extends them to the case of
learnable summary statistics in amortized SBI on potentially massive data
sets, where ABC becomes infeasible. However, we show in Experiment 2 that our
method also works with hand-crafted summary statistics.

Finally, from the perspective of deep anomaly detection, our approach for
learning informative summary statistics can be viewed as a special case of
_generic normality feature learning_ (Pang et al., 2022). Standard learned
summary statistics are optimized with a generic feature learning objective
which is not primarily designed for anomaly detection (Radev et al., 2020a;
Chen et al., 2020). However, since learned summary statistics are also
optimized to be maximally informative for posterior inference, they will
likely capture underlying data regularities (Pang et al., 2022) and serve as a
natural target for detecting model misspecification.

##  3 Method

For simulation-based Bayesian inference, we define a generative model as a
triple
$\mathcal{M}=\big{(}g(\boldsymbol{\theta},\boldsymbol{\xi}),p(\boldsymbol{\xi}\,|\,\boldsymbol{\theta}),p(\boldsymbol{\theta})\big{)}$.
A generative model $\mathcal{M}$ generates data $\boldsymbol{x}\in\mathcal{X}$
according to the system

$\boldsymbol{x}=g(\boldsymbol{\theta},\boldsymbol{\xi})\quad\textrm{with}\quad\boldsymbol{\xi}\sim
p(\boldsymbol{\xi}\,|\,\boldsymbol{\theta}),\;\boldsymbol{\theta}\sim
p(\boldsymbol{\theta}),$ (1)

where $g$ denotes a (randomized) simulator, $\boldsymbol{\xi}\in\Xi$ is a
source of randomness (i.e., noise) with density function
$p(\boldsymbol{\xi}\,|\,\boldsymbol{\theta})$, and $p(\boldsymbol{\theta})$
encodes prior knowledge about plausible simulation parameters
$\boldsymbol{\theta}\in\Theta$. Throughout the paper, we use the decorated
symbol $\accentset{\text{o}}{\boldsymbol{x}}$ to mark data that was in fact
_observed_ in the real world and not merely simulated by the assumed model
$\mathcal{M}$. The parameters $\boldsymbol{\theta}$ consist of hidden
properties whose role in $g$ we explicitly understand and model, and
$\boldsymbol{\xi}$ takes care of nuisance effects that we only treat
statistically. The abstract spaces $\mathcal{X},\Xi$, and $\Theta$ denote the
domain of possible output data (possible worlds), the scope of noise, and the
set of admissible model parameters, respectively. The distinction between
hidden properties $\boldsymbol{\theta}$ and noise $\boldsymbol{\xi}$ is not
entirely clear-cut, but depends on our modeling goals and may vary across
applications.

Our generative model formulation is equivalent to the standard factorization
of the Bayesian joint distribution into likelihood and prior,
$p(\boldsymbol{\theta},\boldsymbol{x}\,|\,\mathcal{M})=p(\boldsymbol{x}\,|\,\boldsymbol{\theta},\mathcal{M})\,p(\boldsymbol{\theta}\,|\,\mathcal{M})$,
where $\mathcal{M}$ expresses the prior knowledge and assumptions embodied in
the model. The likelihood is obtained by marginalizing the joint distribution
$p(\boldsymbol{\xi},\boldsymbol{x}\,|\,\boldsymbol{\theta},\mathcal{M})$ over
all possible values of the nuisance parameters $\boldsymbol{\xi}$, that is,
over all possible execution paths of the simulation program, for fixed
$\boldsymbol{\theta}$:

$p(\boldsymbol{x}\,|\,\boldsymbol{\theta},\mathcal{M})=\int_{\Xi}p(\boldsymbol{\xi},\boldsymbol{x}\,|\,\boldsymbol{\theta},\mathcal{M})\,\mathrm{d}\boldsymbol{\xi}.$
(2)

This integral is typically intractable (Cranmer et al., 2020), but we assume
that it exists and is non-degenerate, that is, it defines a proper density
over the constrained manifold
$\left(g(\boldsymbol{\theta},\boldsymbol{\xi}),\boldsymbol{\xi}\right)$, and
this density can be learned. A major challenge in Bayesian inference is
approximating the posterior distribution
$p(\boldsymbol{\theta}\,|\,\boldsymbol{x},\mathcal{M})\propto
p(\boldsymbol{x}\,|\,\boldsymbol{\theta},\mathcal{M})\,p(\boldsymbol{\theta}\,|\,\mathcal{M})$.
Below, we focus on _amortized_ posterior approximation with neural networks,
which aims to achieve zero-shot posterior sampling for any input data
$\boldsymbol{x}$ compatible with the reference model $\mathcal{M}$.111We
demonstrate in Experiment 1 that model misspecification also affects the
performance of non-amortized sequential neural posterior estimation, and our
method can also detect such misspecification when deployed with SNPE.

###  3.1 Neural Posterior Estimation

Neural Posterior Estimation (NPE) with learned summary statistics
$h_{\boldsymbol{\psi}}(\boldsymbol{x})$ involves a posterior network and a
summary (aka embedding) network. The summary network transforms input data
$\boldsymbol{x}$ of variable size and structure (e.g., multiple IID
measurements) to a fixed-length representation
$\boldsymbol{z}=h_{\boldsymbol{\psi}}(\boldsymbol{x})$. The inference network
$f_{\boldsymbol{\phi}}$ generates random draws from an approximate posterior
$q_{\boldsymbol{\phi}}$ via a conditional generative neural network. In the
following, we focus on normalizing flows realized as conditional invertible
neural networks (Ardizzone et al., 2019), but most of our discussion applies
to any generative backbone architecture that relies on (explicit or implicit)
data compression for flexible inference. Ideally, the two networks jointly
minimize the expected forward Kullback-Leibler (KL) divergence between
analytic and approximate posterior

$\displaystyle(\boldsymbol{\psi}^{*},\boldsymbol{\phi}^{*})$
$\displaystyle=\operatorname*{argmin}_{\boldsymbol{\psi},\boldsymbol{\phi}}\mathbb{E}_{p^{*}(\boldsymbol{x})}\Big{[}\mathbb{KL}\left[p(\boldsymbol{\theta}\,|\,\boldsymbol{x},\mathcal{M})\,||\,q_{\boldsymbol{\phi}}\big{(}\boldsymbol{\theta}\,|\,h_{\boldsymbol{\psi}}(\boldsymbol{x}),\mathcal{M}\big{)}\right]\Big{]}$
(3)
$\displaystyle=\operatorname*{argmin}_{\boldsymbol{\psi},\boldsymbol{\phi}}\mathbb{E}_{p^{*}(\boldsymbol{x})}\Big{[}\mathbb{E}_{p(\boldsymbol{\theta}\,|\,\boldsymbol{x},\mathcal{M})}\left[-\log
q_{\boldsymbol{\phi}}\big{(}\boldsymbol{\theta}\,|\,h_{\boldsymbol{\psi}}(\boldsymbol{x}),\mathcal{M}\big{)}\right]\Big{]},$
(4)

where the expectation runs over the unknown data-generating distribution
$p^{*}(\boldsymbol{x})$ and Eq. 4 follows from Eq. 3 since the analytic
posterior $p(\boldsymbol{\theta}\,|\,\boldsymbol{x},\mathcal{M})$ and its
entropy do not depend on the trainable neural network parameters
($\boldsymbol{\psi}$, $\boldsymbol{\phi}$). In theory, minimizing the above
criterion ensures that the networks learn to perform correct amortized
posterior estimation over any observable configuration $\boldsymbol{x}\sim
p^{*}(\boldsymbol{x})$ (Greenberg et al., 2019; Radev et al., 2020a). In
practice, however, Eq. 4 is hardly a feasible target, since we may not have
enough real data to approximate the outer expectation and the inner
expectation runs over the analytic posterior which we assume is intractable in
the first place.

To render Eq. 4 tractable for NPE, we replace the unknown
$p^{*}(\boldsymbol{x})$ in the outer expectation with the model-implied
distribution (aka evidence, marginal likelihood, or prior predictive
distribution)

$p(\boldsymbol{x}\,|\,\mathcal{M})=\int
p(\boldsymbol{x}\,|\,\boldsymbol{\theta},\mathcal{M})\,p(\boldsymbol{\theta}\,|\,\mathcal{M})\,\mathrm{d}\boldsymbol{x}.$
(5)

The marginal likelihood is a key quantity in Bayesian inference (MacKay, 2003;
Lotfi et al., 2022) and, informally speaking, determines the generative scope
of model $\mathcal{M}$ as an average over the joint model
$p(\boldsymbol{\theta},\boldsymbol{x}\,|\,\mathcal{M})=p(\boldsymbol{x}\,|\,\boldsymbol{\theta},\mathcal{M})\,p(\boldsymbol{\theta}\,|\,\mathcal{M})$.
Thus, we can now rephrase Eq. 4 solely in terms of model-implied quantities:

$\displaystyle(\boldsymbol{\psi}^{*},\boldsymbol{\phi}^{*})$
$\displaystyle=\operatorname*{argmin}_{\boldsymbol{\psi},\boldsymbol{\phi}}\mathbb{E}_{p(\boldsymbol{x}\,|\,\mathcal{M})}\Big{[}\mathbb{E}_{p(\boldsymbol{\theta}\,|\,\boldsymbol{x},\mathcal{M})}\left[-\log
q_{\boldsymbol{\phi}}\big{(}\boldsymbol{\theta}\,|\,h_{\boldsymbol{\psi}}(\boldsymbol{x}),\mathcal{M}\big{)}\right]\Big{]}$
(6)
$\displaystyle=\operatorname*{argmin}_{\boldsymbol{\psi},\boldsymbol{\phi}}\mathbb{E}_{p(\boldsymbol{\theta},\boldsymbol{x}\,|\,\mathcal{M})}\Big{[}-\log
q_{\boldsymbol{\phi}}\big{(}\boldsymbol{\theta}\,|\,h_{\boldsymbol{\psi}}(\boldsymbol{x}),\mathcal{M}\big{)}\Big{]}.$
(7)

This “trick” allows us to approximate the expectation in Eq. 7 via an “online”
stream or an “offline” data set of simulations from the generative model
$\mathcal{M}$ and repeat the process until convergence to an amortized
posterior approximator that is “competent” over
$p(\boldsymbol{x}\,|\,\mathcal{M})$. However, in doing so, we have tacitly
assumed a closed-world setting in which sampling from the prior predictive
distribution $p(\boldsymbol{x}\,|\,\mathcal{M})$ is equivalent to sampling
from the true data distribution $p^{*}(\boldsymbol{x})$. This critical
assumption makes NPE susceptible to posterior errors in the open world due to
model misspecification or “out-of-distribtuion” (OOD) data (Yang et al.,
2021), as we discuss next.

###  3.2 Model Misspecification in Simulation-Based Inference

When modeling a complex system or process, we typically assume an unknown
(true) generator $\mathcal{M}^{*}$, which yields an unknown (true)
distribution $\accentset{\text{o}}{\boldsymbol{x}}\sim p^{*}(\boldsymbol{x})$
and is available to the data analyst only via a finite realization (i.e.,
actually observed data $\accentset{\text{o}}{\boldsymbol{x}}$). According to a
common definition (Frazier et al., 2020; White, 1982; Masegosa, 2020; Lotfi et
al., 2022), the generative model $\mathcal{M}$ is well-specified if a “true”
parameter $\boldsymbol{\theta}^{*}\in\Theta$ exists, such that the likelihood
function matches the data-generating distribution,

$\exists\boldsymbol{\theta}^{*}:p(\boldsymbol{x}\,|\,\boldsymbol{\theta}^{*},\mathcal{M})=p^{*}(\boldsymbol{x}),$
(8)

and misspecified otherwise. This likelihood-centered definition is well-
established and sensible in many domains of traditional Bayesian inference.

In _simulation-based_ inference, however, there is an additional difficulty
regarding model specification: Simulation-based training (see Eq. 7) takes the
expectation with respect to the model-implied distribution
$p(\boldsymbol{x}\,|\,\mathcal{M})$, not necessarily the “true” distribution
$p^{*}(\boldsymbol{x})$ that generates the observed data. Thus, optimal
convergence does not imply correct amortized inference or faithful prediction
in the real world when there is a simulation gap, that is, when the assumed
training model $\mathcal{M}$ deviates critically from the unknown true
generative model $\mathcal{M}^{*}$. Crucially, even if the generative model
$\mathcal{M}$ is well-specified according to the likelihood-centered
definition in Eq. 8, finite training with respect to a “wrong” prior
(predictive) distribution will likely result in insufficient learning of
relevant parameter (and data) regions. This scenario could also be framed as
“out-of-simulation” (OOSim) by analogy with the out-of-distribution (OOD)
problem in machine learning applications (Yang et al., 2021). In fact, we
observe in Experiment 1 that a misspecified prior distribution worsens
posterior inference just like a misspecified likelihood function. What is
more, the same issues also apply when the observed data lies within the
support of the assumed model but only populates a narrow subspace. While this
is not strictly an out-of-simulation setting, the simulation-based neural
network training will possibly be inefficient because the generative forward
model rarely yields relevant training examples. Hence, posterior performance
may be impeded given small simulation budgets—a common scenario in real-life
applications where data is scarce and training times limited.

Thus, our adjusted definition of model misspecification _in the context of
simulation-based inference_ considers the entire prior predictive distribution
$p(\boldsymbol{x}\,|\,\mathcal{M})$ relative to the data-generating process
$p^{*}(\boldsymbol{x})$: A generative model $\mathcal{M}$ is well-specified if
the information loss through modeling $p^{*}(\boldsymbol{x})$ with
$p(\boldsymbol{x}\,|\,\mathcal{M})$ is below an acceptance threshold
$\vartheta$,

$\mathbb{D}\big{[}p(\boldsymbol{x}\,|\,\mathcal{M})\,||\,p^{*}(\boldsymbol{x})\big{]}<\vartheta,$
(9)

and misspecified otherwise. The symbol $\mathbb{D}$ denotes a divergence
metric between the data distributions implied by reality and by the model
(i.e., the marginal likelihood). A natural choice for $\mathbb{D}$ would stem
from the family of $\mathcal{F}$-divergences, such as the KL divergence. We
choose the Maximum Mean Discrepancy (MMD),

$\mathbb{MMD}^{2}\big{[}p^{*}(\boldsymbol{x})\,||\,p(\boldsymbol{x}\,|\,\mathcal{M})\big{]}=\mathbb{E}_{x,\boldsymbol{x}^{\prime}\sim
p^{*}(\boldsymbol{x})}\big{[}\kappa(\boldsymbol{x},\boldsymbol{x}^{\prime})\big{]}+\mathbb{E}_{\boldsymbol{x},\boldsymbol{x}^{\prime}\sim
p(\boldsymbol{x}\,|\,\mathcal{M})}\big{[}\kappa(\boldsymbol{x},\boldsymbol{x}^{\prime})\big{]}-2\mathbb{E}_{{\scriptstyle\begin{array}[]{RCL}\boldsymbol{x}}{&\sim&p^{*}(\boldsymbol{x})\\\
\boldsymbol{x}^{\prime}&\sim&p(\boldsymbol{x}\,|\,\mathcal{M})\end{array}}}\big{[}\kappa(\boldsymbol{x},\boldsymbol{x}^{\prime})\big{]},$
(10)

where $\kappa(\cdot,\cdot)$ is a positive definite kernel. We can tractably
estimate the MMD based on finite samples from
$p(\boldsymbol{x}\,|\,\mathcal{M})$ and $p^{*}(\boldsymbol{x})$ and its
analytic value equals zero if and only if the two densities are equal (Gretton
et al., 2012). In our empirical evaluation, we use sums of Gaussian kernels
with different widths $\sigma_{i}$ as an established and flexible universal
kernel (Muandet et al., 2017). However, Ardizzone et al. (2018) argue that
kernels with heavier tails may improve performance by yielding superior
gradients for outliers. Thus, we repeated all experiments with a sum of
inverse multiquadratic kernels (Tolstikhin et al., 2017), and find that the
results are essentially equal.

Our adjusted definition of model misspecification no longer assumes the
existence of a true parameter vector $\boldsymbol{\theta}^{*}$ (cf. Eq. 8).
Instead, we focus on the marginal likelihood
$p(\boldsymbol{x}\,|\,\mathcal{M})$ which represents the entire prior
predictive distribution of a model and does not commit to a single most
representative parameter vector. In this way, multiple models whose marginal
distributions are representative of $p^{*}(\boldsymbol{x})$ can be considered
well-specified without any reference to some hypothetical ground-truth
$\boldsymbol{\theta}^{*}$, which may not even exist for opaque systems with
unknown properties.

###  3.3 Structured Summary Statistics

In simulation-based inference, summary statistics have a dual purpose because
(i) they are fixed-length vectors, even if the input data $\boldsymbol{x}$
have variable length; and (ii) they usually contain crucial features of the
data, which simplifies neural posterior inference. However, in complex real-
world scenarios such as COVID-19 modeling (see Experiment 4), it is not
feasible to rely on hand-crafted summary statistics. Thus, combining neural
posterior estimation with _learned summary statistics_ leverages the benefits
of summary statistics (i.e., compression to fixed-length vectors) while
avoiding the virtually impossible task of manually designing summariesfor
complex models.

In simulation-based inference, the summary network $h_{\boldsymbol{\psi}}$
acts as an interface between the data $\boldsymbol{x}$ and the inference
network $f_{\boldsymbol{\phi}}$. Its role is to learn maximally informative
summary vectors of fixed size $S$ from complex and structured observations
(e.g., sets of $i.\\-i.\\-d.\\-$ measurements or multivariate time series).
Since the learned summary statistics are optimized to be maximally informative
for posterior inference, they are forced to capture underlying data
regularities (see Section 2). Therefore, we deem the summary network’s
representation $\boldsymbol{z}=h_{\boldsymbol{\psi}}(\boldsymbol{x})$ as an
adequate target to detect simulation gaps.

Specifically, we prescribe an $S$-dimensional multivariate unit (aka.
standard) Gaussian distribution to the summary space,
$p\big{(}\boldsymbol{z}=h_{\boldsymbol{\psi}}(\boldsymbol{x})\,|\,\mathcal{M}\big{)}\approx\mathcal{N}(\boldsymbol{z}\,|\,\boldsymbol{0},\mathbb{I})$,
by minimizing the MMD between summary network outputs and random draws from a
unit Gaussian distribution. To ensure that the summary vectors comply with the
support of the Gaussian density, we use a linear (bottleneck) output layer
with $S$ units in the summary network. A random vector in summary space takes
the form
$h_{\boldsymbol{\psi}}(\boldsymbol{x})\equiv\boldsymbol{z}\equiv(z_{1},\ldots,z_{S})\in\mathbb{R}^{S}$.
The extended optimization objective follows as

$\boldsymbol{\psi}^{*},\boldsymbol{\phi}^{*}=\operatorname*{argmin}_{\boldsymbol{\psi},\boldsymbol{\phi}}\mathbb{E}_{p(\boldsymbol{\theta},\boldsymbol{x}\,|\,\mathcal{M})}\Big{[}-\log
q_{\boldsymbol{\phi}}\big{(}\boldsymbol{\theta}\,\big{|}\,h_{\boldsymbol{\psi}}(\boldsymbol{x}),\mathcal{M}\big{)}\Big{]}+\gamma\,\mathbb{MMD}^{2}\Big{[}p\big{(}h_{\boldsymbol{\psi}}(\boldsymbol{x})\,\big{|}\,\mathcal{M}\big{)}\,\Big{|}\Big{|}\,\mathcal{N}\big{(}\boldsymbol{0},\mathbb{I}\big{)}\Big{]}$
(11)

with a hyperparameter $\gamma$ to control the relative weight of the MMD term.
Intuitively, this objective encourages the approximate posterior
$q_{\boldsymbol{\phi}}\big{(}\boldsymbol{\theta}\,|\,h_{\boldsymbol{\psi}}(\boldsymbol{x}),\mathcal{M}\big{)}$
to match the correct posterior and the summary distribution
$\smash{p\big{(}h_{\boldsymbol{\psi}}(\boldsymbol{x})\,|\,\mathcal{M}\big{)}}$
to match a unit Gaussian. The extended objective does not constitute a
theoretical trade-off between the two terms, since the MMD merely reshapes the
summary distribution in an information-preserving manner. In practice, the
extended objective may render optimization of the summary network more
difficult, but our experiments suggest that it does not restrict the quality
of the amortized posteriors.

This method is also directly applicable to hand-crafted summary statistics,
which already have fixed length and usually contain rich information for
posterior inference. Thus, the task of the summary network $h_{\psi}$
simplifies to transforming the hand-crafted summary statistics to a unit
Gaussian (Eq. 11) to enable model misspecification via distribution matching
during test time (see below). We apply our method to a problem with hand-
crafted summary statistics in Experiment 2.

###  3.4 Detecting Model Misspecification with Finite Data

Once the simulation-based training phase is completed, we can generate $M$
validation samples
$\smash{\\{\boldsymbol{\theta}^{(m)},\boldsymbol{x}^{(m)}\\}}$ from our
generative model $\mathcal{M}$ and pass them through the summary network to
obtain a sample of latent summary vectors
$\smash{\\{\boldsymbol{z}^{(m)}\\}}$, where
$\smash{\boldsymbol{z}=h_{\boldsymbol{\psi}}(\boldsymbol{x})}$ denotes the
summary network output. This sample contains important convergence
information: If $\smash{\boldsymbol{z}}$ is approximately unit Gaussian, we
can assume a structured summary space given the training model $\mathcal{M}$.
This enables misspecification diagnostics via distribution checking during
inference on real data (algorithm in Appendix).

Let $\smash{\\{\accentset{\text{o}}{\boldsymbol{x}}^{(n)}\\}}$ be an
_observed_ sample, either simulated from a different generative model, or
arising from real-world observations with an unknown generator. Before
invoking the inference network, we pass this sample through the summary
network to obtain the summary statistics for the sample:
$\smash{\\{\accentset{\text{o}}{\boldsymbol{z}}^{(n)}\\}}$. We then compare
the validation summary distribution $\\{\boldsymbol{z}^{(m)}\\}$ with the
summary statistics of the observed data
$\smash{\\{\accentset{\text{o}}{\boldsymbol{z}}^{(n)}\\}}$ according to the
sample-based MMD estimate
$\smash{\widehat{\text{MMD}}^{2}\left(p(\boldsymbol{z})\,||\,p(\accentset{\text{o}}{\boldsymbol{z}})\right)}$
(Gretton et al., 2012). Importantly, we are not limited to pre-determined
sizes of simulated or real-world data sets, as the MMD estimator is defined
for arbitrary $M$ and $N$. To allow MMD estimation for data sets with single
instances ($N=1$ or $M=1$), we do not use the unbiased MMD version from
Gretton et al. (2012). Singleton data sets are an important use case for our
method in practice, and potential advantages of unbiased estimators do not
justify exclusion of such data. To enhance visibility, the figures in the
experimental section will plot the square root of the squared MMD estimate.

Whenever we estimate the MMD from finite data, its estimates vary according to
a sampling distribution and we can resort to a frequentist hypothesis test to
determine the probability of observed MMD values under well-specified models.
Although this sampling distribution under the null hypothesis is unknown, we
can estimate it from multiple sets of simulations from the generative model,
$\smash{\\{\boldsymbol{z}^{(m)}\\}}$ and $\smash{\\{\boldsymbol{z}^{(n)}\\}}$,
with $M$ large and $N$ equal to the number of real data sets. Based on the
estimated sampling distribution, we can obtain a critical MMD value for a
fixed Type I error probability $(\alpha)$ and compare it to the one estimated
with the observed data. In general, a larger $\alpha$ level corresponds to a
more conservative modeling approach: A larger type I error implies that more
tests reject the null hypothesis, which corresponds to more frequent model
misspecification alarms and a higher chance that incorrect models will be
recognised. The Type II error probability ($\beta$) of this test will
generally be high (i.e., the _power_ of the test will be low) whenever the
number of real data sets $N$ is very small. However, we show in Experiment 4
that even as few as $5$ real data sets suffice to achieve $\beta\approx 0$ for
a complex model on COVID-19 time series.

##  4 Experiments

###  4.1 Experiment 1: 2D Gaussian Means

Figure 3: Experiment 1. Summary space samples for the minimal sufficient summary network ($S=2$) from a well-specified model $\mathcal{M}$ (blue) and several misspecified configurations. Left: Prior misspecification can be detected. Right: Noise misspecification can be detected, while simulator scale misspecification is indistinguishable from the validation summary statistics. Model (MMS) | Prior | Likelihood
---|---|---
$\mathcal{M}$ (No MMS) | $\boldsymbol{\mu}\sim\mathcal{N}(\boldsymbol{0},\mathbb{I})$ | $\boldsymbol{x}_{k}\sim\mathcal{N}(\boldsymbol{\mu},\mathbb{I})$
$\mathcal{M}_{P}$ (Prior) | $\boldsymbol{\mu}\sim\mathcal{N}(\boldsymbol{\mu}_{0},\tau_{0}\mathbb{I}),\tau_{0}\in\mathbb{R}^{+}$ | $\boldsymbol{x}_{k}\sim\mathcal{N}(\boldsymbol{\mu},\mathbb{I})$
$\mathcal{M}_{S}$ (Simulator) | $\boldsymbol{\mu}\sim\mathcal{N}(\boldsymbol{0},\mathbb{I})$ | $\boldsymbol{x}_{k}\sim\mathcal{N}(\boldsymbol{\mu},\tau\mathbb{I}),\tau\in\mathbb{R}^{+}$
$\mathcal{M}_{N}$ (Noise) | $\boldsymbol{\mu}\sim\mathcal{N}(\boldsymbol{0},\mathbb{I})$ | $\boldsymbol{x}_{k}\sim\lambda\cdot\mathrm{Beta}(2,5)+(1-\lambda)\cdot\mathcal{N}(\boldsymbol{\mu},\mathbb{I})$
Table 1: Experiment 1. Investigated model misspecifications (MMS).
$\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ denotes a (multivariate)
Gaussian distribution with location $\boldsymbol{\mu}$ and covariance
$\boldsymbol{\Sigma}$.

We set the stage by estimating the means of a $2$-dimensional conjugate
Gaussian model with $K=100$ observations per data set and a known analytic
posterior in order to illustrate our method. This experiment contains the
Gaussian examples from Frazier et al. (2020) and Ward et al. (2022), and
extends them by (i) studying misspecifications beyond the likelihood variance
(see below); and (ii) implementing continuously widening simulation gaps, as
opposed to a single discrete misspecification. The data generating process is
defined as

$\boldsymbol{x}_{k}\sim\mathcal{N}(\boldsymbol{x}\,|\,\boldsymbol{\mu},\boldsymbol{\Sigma})\quad\text{for
}k=1,...,K\qquad\text{with
}\;\;\boldsymbol{\mu}\sim\mathcal{N}(\boldsymbol{\mu}\,|\,\boldsymbol{\mu}_{0},\boldsymbol{\Sigma}_{0}).$
(12)

As a summary network, we use a permutation invariant network (Bloem-Reddy &
Teh, 2020) with $S=2$ output dimensions, which equal the number of minimal
sufficient statistics implied by the analytic posterior. The terms “minimal”,
“sufficient”, and “overcomplete” refer to the inference task and _not_ to the
data. Thus, $S=2$ summary statistics are _sufficient_ to solve the inference
task, namely recover two means. For training the posterior approximator, we
set the prior of the generative model $\mathcal{M}$ to a unit Gaussian and the
likelihood covariance $\boldsymbol{\Sigma}$ to an identity matrix.

We induce prior misspecification by altering the prior location
$\boldsymbol{\mu}_{0}$ and covariance
$\boldsymbol{\Sigma}_{0}=\tau_{0}\mathbb{I}$ (diagonal covariance, controlled
through the factor $\tau_{0}$). Further, we achieve misspecified likelihoods
by manipulating the likelihood covariance $\boldsymbol{\Sigma}=\tau\mathbb{I}$
(diagonal covariance, controlled through $\tau$). We induce noise
misspecification by replacing a fraction $\lambda\in[0,1]$ of data
$\boldsymbol{x}$ with samples from a scaled $\text{Beta}(2,5)$ distribution.

Results. The neural posterior estimator trained to minimize the augmented
objective (Eq. 11) exhibits excellent recovery and calibration (Talts et al.,
2020; Säilynoja et al., 2021) in the well-specified case, as shown in the
Appendix. All prior misspecifications manifest themselves in anomalies in the
summary space which are directly detectable through visual inspection of the
$2$-dimensional summary space in Figure 3 (left). Note that the combined prior
misspecification (location and scale) exhibits a summary space pattern that
combines the location and scale of the respective location and scale
misspecifications. However, based on the $2$-dimensional summary space,
misspecifications in the fixed parameters of the likelihood ($\tau$) and
mixture noise are not detectable via an increased MMD (see 5(a), top right).

We further investigate the effect of an _overcomplete_ summary space, namely
$S=4$ learned summary statistics with an otherwise equal architecture. In
addition to prior misspecifications, the overcomplete summary network also
captures misspecifications in the noise and simulator via the MMD criterion
(see 5(a), bottom row). Furthermore, the misspecifications in the noise and
simulator are visually detectable in the summary space (see Appendix). Recall
that the $2$-dimensional summary space fails to capture these
misspecifications (see 5(a), top right).

|  | Model Misspecification |
---|---|---|---
|  | Prior |  Simulator & noise |
Summary Network |  minimal |  |  |
overcomplete |  |  |

(a) Summary MMD by misspecification severity.

|  | Model Misspecification |
---|---|---|---
|  | Prior | Simulator & noise |
Summary Network |  minimal |  |  |
overcomplete |  |  |

(b) Posterior error by misspecification severity.

Figure 5: Experiment 1. Summary space discrepancy (MMD to training
distribution) and posterior error (RMSE of correct vs. analytic posterior
means) as a function of misspecification severity. White stars indicate the
well-specified model configuration (i.e., equal to the training model
$\mathcal{M}$), where both MMD and posterior error are low.

Posterior Errors. While detecting model misspecification might be valuable by
itself for applied modelers, one fundamental question remains: Are the
inference results of NPE even impaired by model misspecification? To answer
this question, we compute the error in posterior recovery as a function of the
misspecification severity. To ease visualization, we use the RMSE of the
approximated posterior mean from $L$ samples
$\hat{\overline{\boldsymbol{\theta}}}^{(n)}=\frac{1}{L}\sum_{l=1}^{L}\boldsymbol{\theta}^{(l,n)}$
against the analytic posterior mean $\overline{\boldsymbol{\theta}}^{(n)}$
across $N$ data sets.222Since the approximate posterior in the Gaussian model
is likely to be symmetric—and the analytic posterior is symmetric by
definition—we deem the posterior mean as an appropriate evaluation target for
the RMSE across data sets. In fact, error metrics over several posterior
quantiles (i.e., $Q_{25}$, $Q_{50}$ and $Q_{75}$) in the place of posterior
means yield identical results. Other common metrics (i.e., MSE and MAE) yield
identical results.

5(b) illustrates that more severe model misspecifications generally coincide
with a larger error in posterior estimation across all model misspecifications
for both $S=2$ and $S=4$ learned summary statistics, albeit with fundamental
differences, as explained in the following. When processing data emerging from
models with misspecified noise and simulator (see 5(b), right column), minimal
and overcomplete summary networks exhibit a drastically different behavior:
While the minimal summary network cannot detect noise or simulator simulation
gaps, its posterior estimation performance is not heavily impaired either (see
5(b), top right). On the other hand, the overcomplete summary network is able
to capture noise and simulator misspecifications, but also incurs larger
posterior inference error (see 5(b), bottom right). This might suggest a
trade-off between model misspecification detection and posterior inference
error, depending on the number of learnable summary statistics.

From a practical modeling perspective, researchers might wonder how to choose
the number of learnable summary statistics. While an intuitive heuristic might
suggest “the more, the merrier”, the observed results in this experiment beg
to differ depending on the modeling goals. If the focus in a critical
application lies in detecting potential simulation gaps, it might be
advantageous to utilize a large (overcomplete) summary vector. However,
modelers might also desire a network which is as robust as possible during
test time, opting for a smaller number of summary statistics. Experiments 4
addresses this question for a complex non-linear time series model where the
number of sufficient summary statistics is unknown.

Replication with the SNPE-C method. Our method successfully detects model
misspecification when we use non-amortized sequential neural posterior
estimation (SNPE-C; Greenberg et al., 2019) with a structured summary space
(see Appendix for details). The results are largely equivalent to those
obtained with amortized NPE, as implemented in the BayesFlow framework (Radev
et al., 2023c).

###  4.2 Experiment 2: Cancer and Stromal Cell Model

This experiment illustrates model misspecification detection in a marked point
process model of cancer and stromal cells (Jones-Todd et al., 2019). We use
the implementation of Ward et al. (2022) with hand-crafted summary statistics
and showcase the applicability of our method in scenarios where good summary
statistics are already known. The inference parameters are three Poisson rates
$\lambda_{c},\lambda_{p},\lambda_{d}$, and the setup in Ward et al. (2022)
extracts four hand-crafted summary statistics from the 2D plane data: (1–2)
number of cancer and stromal cells; (3–4) mean and maximum distance from
stromal cells to the nearest cancer cell.

We achieve misspecification during inference by mimicking necrosis, which
often occurs in core regions of tumors. A Bernoulli distribution with
parameter $\pi$ controls whether a cell is affected by necrosis or not.
Consequently, $\pi=0$ implies no necrosis (and thus no simulation gap), and
$\pi=1$ entails that all cells are affected. The experiments by Ward et al.
(2022) study a single misspecification scenario, namely the case $\pi=0.75$ in
our implementation. In order to employ our proposed method for model
misspecification detection, we add a small summary network
$h_{\boldsymbol{\psi}}:\mathbb{R}^{4}\to\mathbb{R}^{4}$ consisting of three
hidden fully connected layers with $64$ units each. This network
$h_{\boldsymbol{\psi}}$ merely transforms the hand-crafted summary statistics
into a $4$-D unit Gaussian, followed by NPE for posterior inference.

(a) MMD by misspecification severity.

(b) Power for $\pi=0.75$ is essentially 1.

Figure 6: Experiment 2. MMD increases with misspecification severity (6(a);
mean, SD of 20 repetitions). Our test easily detects the setting from Ward et
al. (2022) (6(b)).

Results. Our MMD misspecification score increases with increasingly severe
model misspecification (i.e., increasing necrosis rate $\pi$), see 6(a). What
is more, for the single misspecification $\pi=0.75$ studied by Ward et al.
(2022), we illustrate (i) the power of our proposed hypothesis test; and (ii)
the summary space distribution for misspecified data. The power ($1-\beta$)
essentially equals $1$, as shown in 6(b): The MMD sampling distributions under
the training model ($H_{0}$) and under the observed data generating process
($\mathcal{M}^{*}$) are completely separated.

###  4.3 Experiment 3: Drift-Diffusion Model of Decision Making

Figure 7: Experiment 3. In each trial, the drift-diffusion model (DDM) models
the information uptake in experimental condition $j\in\\{1,2\\}$ as a random
walk with drift $v_{j}$ and initial non-decision time $t_{0}$. Once the
evidence reaches a decision threshold $a_{j}$, the associated reaction time is
recorded. After $N$ trials, the $i.\,i.\,d.\,$ set of reaction times
$\boldsymbol{x}$ is used to infer the mechanistic parameters
$\boldsymbol{\theta}=(v_{1},v_{2},a_{1},a_{2},t_{0})$.

This experiment aims to (i) apply the new optimization objective to a complex
model of human decision making; (ii) tackle strategies to determine the
required number of learned summary statistics in more complex applications;
and (iii) compare the posterior estimation of amortized NPE under a
misspecified model with a reference posterior from the gold-standard non-
amortized Markov Chain Monte Carlo algorithm (HMC-MCMC).

We focus on the drift diffusion model (DDM)—a cognitive model describing
reaction times (RTs) in binary decision tasks (Ratcliff & McKoon, 2008) which
is well amenable to amortized inference (Radev et al., 2020b). The DDM assumes
that perceptual information for a choice alternative accumulates continuously
according to a Wiener diffusion process. Thus, the change in information in
experimental condition $j$ follows a random walk with drift and Gaussian
noise: $v\mathrm{d}t+\xi\sqrt{\mathrm{d}t}$ with $\xi\sim\mathcal{N}(0,1)$.
Our model implementation features five free parameters
$\boldsymbol{\theta}=(v_{1},v_{2},a_{1},a_{2},t_{0})$ which produce data
stemming from two simulated conditions. The summary network is a permutation-
invariant network which reduces $i.\,i.\,d.\ $RT data sets to $S=10$ summary
statistics each.

We realize a simulation gap by simulating typically observed contaminants:
fast guesses (e.g., due to inattention), very slow responses (e.g., due to
mind wandering), or a combination of the two. The parameter $\lambda$ controls
the fraction of the observed data which is contaminated (see Section E.3 in
the Appendix for implementation details). For the comparison with the MCMC
reference posterior from Stan, we simulate $100$ uncontaminated DDM data sets
from the true model $\mathcal{M}$, as well as three scenarios (fast guesses,
slow responses, fast and slow combined) with a fraction of $\lambda=10\%$ of
contaminated observations.

(a) MMD as a function of misspecification severity.

(b) Correlation between parameters $\boldsymbol{\theta}$ and principal
components (PCs) of learned summary statistics.

Figure 8: Experiment 3. Contamination of reaction times is detectable with our
method (8(a)), and the principal components from the learned summary
statistics coincide with the true parameters $\boldsymbol{\theta}$ (8(b)).
White star indicates the configuration for the well-specified model (i.e.,
training model $\mathcal{M}$ without contamination).

Results. During inference, our criterion reliably detects the induced
misspecifications: Increasing fractions $\lambda$ of contaminants (fast, slow,
and combined) manifest themselves in increasing MMD values (see 8(a)). The
results of applying PCA to the summary network outputs
$\\{\accentset{\text{o}}{\boldsymbol{z}}^{(n)}\\}$ for the well-specified
model (no contamination) are illustrated in 8(b). We observe that the first
five principal components exhibit a large overlap with the true model
parameters $\boldsymbol{\theta}$ and jointly account for 85% of the variance
in the summary output. Furthermore, the drift rates and decision thresholds
within conditions are entangled (i.e., $v_{1},a_{1}$ and $v_{2},a_{2}$). This
entanglement mimics the strong posterior correlations observed between these
two parameters. In practical applications, dimensionality reduction might act
as a guideline for determining the number of minimally sufficient summary
statistics or parameter redundancies for a given problem.

Model (Contamination) | Posterior error | Summary space MMD
---|---|---
$\mathcal{M}_{\ }$: True model | $0.25\,[0.13,0.56]$ | $0.45\,[0.42,0.52]$
$\mathcal{M}_{1}$: Fast contaminants | $2.66\,[1.44,3.40]$ | $2.68\,[2.61,2.74]$
$\mathcal{M}_{2}$: Slow contaminants | $0.55\,[0.23,1.01]$ | $1.18\,[1.13,1.26]$
$\mathcal{M}_{3}$: Fast & slow contaminants | $1.90\,[0.83,3.18]$ | $2.33\,[2.19,2.43]$
Table 2:  Experiment 3. Posterior error of the approximate neural posterior
(MMD to reference MCMC posterior; median and 95% CI). The bootstrapped MMD
values for the summary statistics of the $100$ investigated data sets and
$1\,000$ samples from the true uncontaminated model $\mathcal{M}$ illustrate
that posterior errors are mirrored by anomalies in the summary space and thus
detectable.

For the comparison with the MCM reference posterior from Stan, we juxtapose
$4\,000$ samples from the approximate neural posterior with $4\,000$ samples
obtained from the Stan sampler after ensuring MCMC convergence and sufficient
sampling efficiency for each data set (see Figure 2 for an illustration).
Because Stan is currently considered state-of-the-art for likelihood-based
Bayesian inference, we assume the Stan samples are representative of the
_correct posterior under the potentially misspecified model_ (see Section B).
In order to quantify the difference between the posterior samples from NPE
(BayesFlow) and the MCMC reference (Stan), we use the MMD criterion as well.
When no model misspecification is present, the posterior samples from NPE
(BayesFlow) and MCMC (Stan) match almost perfectly (see 2(a)). This means that
our augmented optimization objective still enables correct posterior
approximation under well-specified models. In contrast, the results in 2(b)
and Table 2 clearly indicate that the amortized posteriors deteriorate as a
result of the induced misspecification. Moreover, these results closely mirror
the overall detectability of misspecification obtained by matching the summary
representations of $1000$ data sets from the true uncontaminated process with
the representations of the $100$ data sets for each of the above scenarios via
MMD (see Table 2). This confirms the hypothesis that our proposed MMD
criterion for quantifying model misspecification acts as a proxy for the
posterior estimation error, which is of hallmark importance but typically
unknown in practice due to a lack of ground-truths.

###  4.4 Experiment 4: Epidemiological Model for COVID-19

(a) Representation of Germany’s COVID-19 time series w.r.t. the MMD
distribution under the null hypothesis
$H_{0}:p^{*}(\boldsymbol{x})=p(\boldsymbol{x}\,|\,\mathcal{M})$.

(b) Cumulative explained variance ratio for the COVID-19 summary statistics as
a function of the number of principal components (PCs).

Figure 9: Experiment 4. Our model misspecification hypothesis test indicates
that the assumed model is well-specified for the real observed COVID data in
Germany (a). Further, linear dimensionality reduction reveals that 40
principal components can explain 95% of the variance in learned summary
statistics.

As a challenging real-world example, we treat a high-dimensional compartmental
model representing the early months of the COVID-19 pandemic in Germany (Radev
et al., 2021). We investigate the utility of our method to detect simulation
gaps in a much more realistic and non-trivial extension of the SIR settings in
Lueckmann et al. (2021b) and Ward et al. (2022) with substantially increased
complexity. Moreover, we perform inference on real COVID-19 data from Germany
and use our new method to test whether the model used in Radev et al. (2021)
is misspecified, possibly undermining the trustworthiness of political
conclusions that are based on the inferred posteriors. To achieve this, we
train an NPE setup identical to Radev et al. (2021) but using our new
optimization objective (Eq. 11) to encourage a structured summary space. We
then simulate $1000$ time series from the training model $\mathcal{M}$ and
$1000$ time series from three misspecified models: (i) a model
$\mathcal{M}_{1}$ without an intervention sub-model; (ii) a model
$\mathcal{M}_{2}$ without an observation sub-model; (iii) a model
$\mathcal{M}_{3}$ without a “carrier” compartment (Dehning et al., 2020).

| Bootstrap MMD | Power ($1-\beta$)
---|---|---
Model | $N=1$ | $N=2$ | $N=5$ | $N=1$ | $N=2$ | $N=5$
$\mathcal{M}_{\ }$ | $3.70\,[3.65,3.79]$ | $2.61\,[2.54,2.91]$ | $1.66\,[1.59,1.84]$ | — | — | —
$\mathcal{M}_{1}$ | $3.76\,[3.72,3.80]$ | $2.86\,[2.62,3.16]$ | $2.11\,[1.82,2.50]$ | $.998$ | $.958$ | $\approx 1.0$
$\mathcal{M}_{2}$ | $3.80\,[3.73,3.83]$ | $2.81\,[2.65,3.00]$ | $2.01\,[1.82,2.19]$ | $.789$ | $.804$ | $\approx 1.0$
$\mathcal{M}_{3}$ | $3.78\,[3.74,3.83]$ | $2.81\,[2.68,3.11]$ | $2.07\,[1.92,2.41]$ | $.631$ | $.690$ | $\approx 1.0$
Table 3:  Experiment 4. Results for different variations of the COVID-19
compartmental model (median and 95% CI of 100 bootstrap samples). Our proposed
detection criterion reliably flags data from misspecified models, with
essentially no false-negatives for as few as $N=5$ observations.

Results. Table 3 shows the MMD between the summary representation of $N=1,2,5$
bootstrapped time series from each model and the summary representation of the
$1000$ time series from model $\mathcal{M}$ (see the Appendix for
bootstrapping details). We also calculate the power ($1-\beta$) of our
hypothesis test for each misspecified model under the sampling distribution
estimated from $1\,000$ samples of the $1\,000$ time series from $\mathcal{M}$
at a type I error probability of $\alpha=.05$. We observe that the power of
the test rapidly increases with more data sets and the Type II error
probability ($\beta$) is essentially zero for as few as $N=5$ time series.

As a next step, we pass the reported COVID-19 data between 1 March and 21
April 2020 (Dong et al., 2020) through the summary network and compute the
critical MMD value for a sampling-based hypothesis test with an $\alpha$ level
of $.05$ (see 9(a)). The MMD of the Germany data is well below the critical
MMD value—it essentially lies in the bulk of the distribution. Hence, we
conclude that the assumed training model $\mathcal{M}$ is well-specified for
this time period. Finally, we perform linear dimensionality reduction (PCA) on
the summary space and find that the first 40 principal components jointly
explain $95\%$ of the variance in the $192$-dimensional summary space outputs
(see 9(b)). Thus, a $40$-dimensional learned summary vector might provide a
good approximation of the true (unknown) minimally sufficient summary
statistics and render inference less fragile under potential misspecifications
(cf. Experiment 1).

###  4.5 Experiment 5: Style Inference of EMNIST Image Simulator

Figure 10: Experiment 5. Overview of the training objective. The latent style
vector $\boldsymbol{\theta}\in\mathbb{R}^{20}$ and the digit label $C$
parameterize the GIN decoder (Sorrenson et al., 2020) which simulates a
realistic hand-drawn digit $\boldsymbol{x}$. The inverse problem aims to
recover the posterior distribution $p(\boldsymbol{\theta}\,|\,\boldsymbol{x})$
of plausible styles given an observed image $\boldsymbol{x}$. Our method can
detect anomalies in the structured learned summary statistics that originate
from (i) atypical styles due to a wider prior $p(\boldsymbol{\theta})$; and
(ii) domain shifts to new digits. In Experiment 2, we additionally pass the
embedded digit label to the inference network to mimic knowledge about a
domain shift in practical applications.

As a final experiment, we demonstrate the sensitivity of our detection method
on a pre-trained generative network which we treat as an image simulator
(i.e., an implicit likelihood model), see Figure 10 for an overview. We use a
general incompressible-flow network (GIN; Sorrenson et al., 2020) which has
been trained on the digits of the EMNIST data set (Cohen et al., 2017). The
GIN learns a class-conditional bijective mapping from the image domain to a
disentangled latent space, which can be interpreted as a nonlinear independent
component analysis (ICA). We treat the decoder (latent $\mapsto$ image) of the
pre-trained GIN as a simulation program and only use the 20 dimensions with
the highest variance of the latent space, which encode global (across all
digits) and local (digit-specific) style features (see Sorrenson et al., 2020,
for a detailed analysis). This introduces uncertainty, because we omit
information from the GIN’s learned bijective mapping. We use NPE to compensate
for the induced uncertainty by approximating the inverse problem, namely,
performing probabilistic inference on the 20 latent style dimensions
(inference targets $\boldsymbol{\theta}$) based on simulated images
(observables $\boldsymbol{x}$). The summary network is a simple CNN that
learns a 64-dimensional representation of the observable image data
$\boldsymbol{x}$. The inference network is tasked with learning the Bayesian
posterior distribution $p(\boldsymbol{\theta}\,|\,\boldsymbol{x})$ of the
latent styles $\boldsymbol{\theta}\in\mathbb{R}^{20}$ conditional on the
observable images $\boldsymbol{x}\in\mathbb{R}^{784}$. We realize the
inference network as an affine coupling flow with 10 affine coupling layers,
512 units per layer, $10^{-4}$ kernel regularization, and dropout layers with
a $0.15$ dropout rate.

In this experiment, we induce model misspecification via two strategies:
First, we train the neural approximator on a narrow subset of possible styles
by halving the dispersion of the style prior $p(\boldsymbol{\theta})$. During
inference, out-of-distribution styles can easily be generated by sampling from
a wider style prior. Second, we only train the neural approximator on a subset
of all “domains”, namely the digits 0–4. During inference, we achieve out-of-
domain samples for the digits 5–8.

We further investigate two variations of this experiment. In Variation 1
(class-agnostic), the neural posterior approximator is trained and evaluated
without any explicit label information (digits 0–9). In Variation 2 (class-
informed), a simple 1D embedding of the digit index is passed as an additional
condition to the invertible inference network.

Results. NPE exhibits excellent recovery of the latent style vectors of in-
domain samples with our modified misspecification-aware optimization objective
in Eq. 11 (see Appendix J for detailed recovery plots of all inferred style
dimensions). Moreover, NPE shows a heavily impaired performance (RMSE and bias
to ground-truth) under extreme styles and domain shifts (see Table 4).
Crucially, our proposed detection criterion successfully detects extreme style
shifts (wider prior $p(\boldsymbol{\theta})$) as well as domain shifts to
unseen digits (5–9) via a significantly increased MMD score compared to the
in-domain MMD score. These effects persist regardless of whether the neural
approximators are class-agnostic (variation 1) or class-informed (variation
2). These results suggest that our proposed detection method is a reliable
proxy to gauge the trustworthiness of NPE under a potentially misspecified
probabilistic model and alert the user about unfaithful inference in case of
unknown (variation 1) or known (variation 2) domain shifts.

Table 4: Experiment 5. Our method consistently detects atypically extreme styles, as well as domain shifts to the unseen digits 5–9, as indexed by a large MMD score compared to the in-distribution baseline (first row). The neural posterior approximator shows impaired performance in these unknown domains, as quantified by higher RMSE and bias relative to the ground-truth styles.  | Variation 1: Class-agnostic | Variation 2: Class-informed
---|---|---
| MMD | RMSE | Bias | MMD | RMSE | Bias
In distribution, digits 0–4 | 0.020 | 0.109 | 0.004 | 0.026 | 0.124 | 0.006
Extreme styles, digits 0–4 | 1.592 | 0.418 | 0.073 | 0.165 | 0.896 | 0.086
New domain: Digit 5 | 4.031 | 0.860 | 0.565 | 2.576 | 1.050 | 0.774
New domain: Digit 6 | 3.485 | 0.973 | 0.859 | 2.248 | 1.223 | 1.372
New domain: Digit 7 | 3.786 | 0.993 | 0.915 | 2.033 | 1.021 | 0.968
New domain: Digit 8 | 3.873 | 0.844 | 0.549 | 2.393 | 0.917 | 0.573
New domain: Digit 9 | 5.702 | 0.745 | 0.506 | 2.470 | 0.693 | 0.440

##  5 Conclusions

This paper approached a fundamental problem in amortized simulation-based
Bayesian inference, namely, flagging potential posterior errors due to model
misspecification. We argued that misspecified models might cause so-called
_simulation gaps_ , resulting in deviations between simulations during
training time and actual observed data at test time. We further showed that
simulation gaps can be detrimental for the performance and faithfulness of
simulation-based inference relying on neural networks. We proposed to increase
the networks’ awareness of posterior errors by compressing simulations into a
structured latent summary space induced by a modified optimization objective
in an unsupervised fashion. We then applied the maximum mean discrepancy (MMD)
estimator, equipped with a sampling-based hypothesis test, as a criterion to
spotlight discrepancies between model-implied and actually observed
distributions in summary space. While we focused on the application to NPE
(BayesFlow implementation; Radev et al., 2020a; 2023b) and SNPE (sbi
implementation; Tejero-Cantero et al., 2020), our proposed method can be
easily integrated into other inference algorithms and frameworks as well. Most
notably, while we focus on detection, our reliable misspecification detection
method can be easily combined with techniques that aim to mitigate model
misspecification. Our method can pinpoint individual data sets where inference
is potentially unfaithful, and local non-amortized techniques may be used as a
targeted post-hoc correction for these problematic data sets. Such a two-step
approach might consist of (i) detection with our method; and (ii) targeted
mitigation with local methods, and constitute a promising alternative to
current global regularization techniques. Our methods are implemented in the
open-source BayesFlow library (www.bayesflow.org) and can be seamlessly
integrated into an end-to-end workflow for amortized simulation-based
inference.

### Acknowledgments

We acknowledge support by Cyber Valley CyVy-RF-2021-16 (MS, PCB), Deutsche
Forschungsgemeinschaft EXC-2075 - 390740016 (MS, PCB), the Informatics for
Life initiative (UK), and Deutsche Forschungsgemeinschaft EXC-2181 - 390900948
(STR).

### Data Availability

Our methods are implemented in the free open-source library BayesFlow
(www.bayesflow.org) and can be seamlessly integrated into amortized Bayesian
workflows.

## References

  * Alquier & Ridgway (2019) Pierre Alquier and James Ridgway.  Concentration of tempered posteriors and of their variational approximations.  _arXiv:1706.09293 [cs, math, stat]_ , April 2019.  arXiv: 1706.09293.
  * Ardizzone et al. (2018) Lynton Ardizzone, Jakob Kruse, Sebastian J. Wirkert, Daniel Rahner, Eric W. Pellegrini, Ralf S. Klessen, Lena Maier-Hein, Carsten Rother, and Ullrich Köthe.  Analyzing Inverse Problems with Invertible Neural Networks.  _CoRR_ , abs/1808.04730, 2018.
  * Ardizzone et al. (2019) Lynton Ardizzone, Carsten Lüth, Jakob Kruse, Carsten Rother, and Ullrich Köthe.  Guided Image Generation with Conditional Invertible Neural Networks, 2019.
  * Barnard et al. (2000) John Barnard, Robert McCulloch, and Xiao-Li Meng.  Modelling Covariance Matrices in Terms of Standard Deviations and Correlations, with Application To Shrinkage.  _Statistica Sinica_ , 10:1281–1311, 10 2000.
  * Berger & Wolpert (1988) James O. Berger and Robert L. Wolpert.  _The likelihood principle_.  Number v. 6 in Lecture notes-monograph series. Institute of Mathematical Statistics, Hayward, Calif, 2nd ed edition, 1988.  ISBN 978-0-940600-13-3.
  * Betancourt (2017) Michael Betancourt.  A conceptual introduction to Hamiltonian Monte Carlo.  _arXiv preprint_ , 2017.
  * Bieringer et al. (2021) Sebastian Bieringer, Anja Butter, Theo Heimel, Stefan Höche, Ullrich Köthe, Tilman Plehn, and Stefan T Radev.  Measuring QCD Splittings with Invertible Networks.  _SciPost Physics Proceedings_ , 10(6), 2021.
  * Bloem-Reddy & Teh (2020) Benjamin Bloem-Reddy and Yee Whye Teh.  Probabilistic Symmetries and Invariant Neural Networks.  _J. Mach. Learn. Res._ , 21:90–1, 2020.
  * Butter et al. (2022) Anja Butter, Tilman Plehn, Steffen Schumann, Simon Badger, Sascha Caron, Kyle Cranmer, Francesco Armando Di Bello, Etienne Dreyer, Stefano Forte, Sanmay Ganguly, et al.  Machine learning and lhc event generation.  _arXiv preprint arXiv:2203.07460_ , 2022.
  * Bürkner et al. (2020) Paul-Christian Bürkner, Jonah Gabry, and Aki Vehtari.  Approximate leave-future-out cross-validation for Bayesian time series models.  _Journal of Statistical Computation and Simulation_ , 90(14):2499–2523, September 2020.  doi: 10.1080/00949655.2020.1783262.  arXiv:1902.06281 [stat].
  * Cannon et al. (2022) Patrick Cannon, Daniel Ward, and Sebastian M. Schmon.  Investigating the Impact of Model Misspecification in Neural Simulation-based Inference, September 2022.  arXiv:2209.01845 [cs, stat].
  * Chen et al. (2020) Yanzhi Chen, Dinghuai Zhang, Michael Gutmann, Aaron Courville, and Zhanxing Zhu.  Neural approximate sufficient statistics for implicit models.  _arXiv preprint arXiv:2010.10079_ , 2020.
  * Cohen et al. (2017) Gregory Cohen, Saeed Afshar, Jonathan Tapson, and Andre Van Schaik.  Emnist: Extending mnist to handwritten letters.  In _2017 international joint conference on neural networks (IJCNN)_ , pp. 2921–2926. IEEE, 2017.
  * Cranmer et al. (2020) Kyle Cranmer, Johann Brehmer, and Gilles Louppe.  The frontier of simulation-based inference.  _Proceedings of the National Academy of Sciences_ , 117(48):30055–30062, 2020.
  * Dehning et al. (2020) Jonas Dehning, Johannes Zierenberg, F Paul Spitzner, Michael Wibral, Joao Pinheiro Neto, Michael Wilczek, and Viola Priesemann.  Inferring change points in the spread of COVID-19 reveals the effectiveness of interventions.  _Science_ , 369(6500), 2020.
  * Dellaporta et al. (2022) Charita Dellaporta, Jeremias Knoblauch, Theodoros Damoulas, and François-Xavier Briol.  Robust Bayesian Inference for Simulator-based Models via the MMD Posterior Bootstrap, 2022.
  * Dong et al. (2020) Ensheng Dong, Hongru Du, and Lauren Gardner.  An interactive web-based dashboard to track COVID-19 in real time.  _The Lancet Infectious Diseases_ , 20(5):533–534, May 2020.  doi: 10.1016/S1473-3099(20)30120-1.
  * Durkan et al. (2020) Conor Durkan, Iain Murray, and George Papamakarios.  On contrastive learning for likelihood-free inference.  In _International Conference on Machine Learning_ , pp. 2771–2781. PMLR, 2020.
  * Frazier & Drovandi (2021) David T. Frazier and Christopher Drovandi.  Robust Approximate Bayesian Inference With Synthetic Likelihood.  _Journal of Computational and Graphical Statistics_ , 30(4):958–976, October 2021.  doi: 10.1080/10618600.2021.1875839.
  * Frazier et al. (2020) David T. Frazier, Christian P. Robert, and Judith Rousseau.  Model misspecification in approximate Bayesian computation: consequences and diagnostics.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 82(2):421–444, April 2020.  doi: 10.1111/rssb.12356.
  * Gabry et al. (2019) Jonah Gabry, Daniel Simpson, Aki Vehtari, Michael Betancourt, and Andrew Gelman.  Visualization in Bayesian workflow.  _Journal of the Royal Statistical Society: Series A (Statistics in Society)_ , 182(2):389–402, 2019.  Publisher: Wiley Online Library.
  * Ghaderi-Kangavari et al. (2023) Amin Ghaderi-Kangavari, Jamal Amani Rad, and Michael D. Nunez.  A general integrative neurocognitive modeling framework to jointly describe EEG and decision-making on single trials.  _Computational Brain and Behavior_ , 2023.  doi: 10.1007/s42113-023-00167-4.
  * Giummolè et al. (2019) Federica Giummolè, Valentina Mameli, Erlis Ruli, and Laura Ventura.  Objective bayesian inference with proper scoring rules.  _Test_ , 28(3):728–755, 2019.
  * Gloeckler et al. (2023) Manuel Gloeckler, Michael Deistler, and Jakob H Macke.  Adversarial robustness of amortized bayesian inference.  _arXiv preprint arXiv:2305.14984_ , 2023.
  * Gonçalves et al. (2020) Pedro J Gonçalves, Jan-Matthis Lueckmann, Michael Deistler, Marcel Nonnenmacher, Kaan Öcal, Giacomo Bassetto, Chaitanya Chintaluri, William F Podlaski, Sara A Haddad, Tim P Vogels, et al.  Training deep neural density estimators to identify mechanistic models of neural dynamics.  _Elife_ , 9:e56261, 2020.
  * Greenberg et al. (2019) David Greenberg, Marcel Nonnenmacher, and Jakob Macke.  Automatic posterior transformation for likelihood-free inference.  In _International Conference on Machine Learning_ , pp. 2404–2414. PMLR, 2019.
  * Gretton et al. (2012) A Gretton, K. Borgwardt, Malte Rasch, Bernhard Schölkopf, and AJ Smola.  A Kernel Two-Sample Test.  _The Journal of Machine Learning Research_ , 13:723–773, 03 2012.
  * Grünwald et al. (2017) Peter Grünwald, Thijs Van Ommen, et al.  Inconsistency of Bayesian inference for misspecified linear models, and a proposal for repairing it.  _Bayesian Analysis_ , 12(4):1069–1103, 2017.
  * Hermans et al. (2020) Joeri Hermans, Volodimir Begy, and Gilles Louppe.  Likelihood-free mcmc with amortized approximate ratio estimators.  In _International Conference on Machine Learning_ , pp. 4239–4248. PMLR, 2020.
  * Holmes & Walker (2017) Chris C Holmes and Stephen G Walker.  Assigning a value to a power likelihood in a general bayesian model.  _Biometrika_ , 104(2):497–503, 2017.
  * Huang et al. (2023) Daolang Huang, Ayush Bharti, Amauri H Souza, Luigi Acerbi, and Samuel Kaski.  Learning robust statistics for simulation-based inference under model misspecification.  In _Thirty-seventh Conference on Neural Information Processing Systems_ , 2023.  URL https://openreview.net/forum?id=STrXsSIEiq.
  * Huang et al. (2024) Daolang Huang, Ayush Bharti, Amauri Souza, Luigi Acerbi, and Samuel Kaski.  Learning robust statistics for simulation-based inference under model misspecification.  _Advances in Neural Information Processing Systems_ , 36, 2024.
  * Jones-Todd et al. (2019) Charlotte M. Jones-Todd, Peter Caie, Janine B. Illian, Ben C. Stevenson, Anne Savage, David J. Harrison, and James L. Bown.  Identifying prognostic structural features in tissue sections of colon cancer patients using point pattern analysis: Point pattern analysis of colon cancer tissue sections.  _Statistics in Medicine_ , 38(8):1421–1441, April 2019.  doi: 10.1002/sim.8046.
  * Kelly et al. (2023) Ryan P Kelly, David J Nott, David T Frazier, David J Warne, and Chris Drovandi.  Misspecification-robust sequential neural likelihood.  _arXiv preprint arXiv:2301.13368_ , 2023.
  * Knoblauch et al. (2019) Jeremias Knoblauch, Jack Jewson, and Theodoros Damoulas.  Generalized variational inference: Three arguments for deriving new posteriors.  _arXiv preprint arXiv:1904.02063_ , 2019.
  * Leclercq (2022) Florent Leclercq.  Simulation-based inference of Bayesian hierarchical models while checking for model misspecification, September 2022.  arXiv:2209.11057 [astro-ph, q-bio, stat].
  * Loaiza-Maya et al. (2021) Ruben Loaiza-Maya, Gael M Martin, and David T Frazier.  Focused bayesian prediction.  _Journal of Applied Econometrics_ , 36(5):517–543, 2021.
  * Lotfi et al. (2022) Sanae Lotfi, Pavel Izmailov, Gregory Benton, Micah Goldblum, and Andrew Gordon Wilson.  Bayesian model selection, the marginal likelihood, and generalization.  _arXiv preprint arXiv:2202.11678_ , 2022.
  * Lueckmann et al. (2017) Jan-Matthis Lueckmann, Pedro J Goncalves, Giacomo Bassetto, Kaan Öcal, Marcel Nonnenmacher, and Jakob H Macke.  Flexible statistical inference for mechanistic models of neural dynamics.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Lueckmann et al. (2021a) Jan-Matthis Lueckmann, Jan Boelts, David Greenberg, Pedro Goncalves, and Jakob Macke.  Benchmarking simulation-based inference.  In _International Conference on Artificial Intelligence and Statistics_ , pp. 343–351. PMLR, 2021a.
  * Lueckmann et al. (2021b) Jan-Matthis Lueckmann, Jan Boelts, David Greenberg, Pedro Goncalves, and Jakob Macke.  Benchmarking Simulation-Based Inference.  In Arindam Banerjee and Kenji Fukumizu (eds.), _Proceedings of The 24th International Conference on Artificial Intelligence and Statistics_ , volume 130 of _Proceedings of Machine Learning Research_ , pp. 343–351. PMLR, April 2021b.
  * MacKay (2003) David JC MacKay.  _Information theory, inference and learning algorithms_.  Cambridge university press, 2003.
  * Mardia et al. (1979) Kantilal Varichand Mardia, John T. Kent, and John M. Bibby.  _Multivariate analysis_.  Probability and mathematical statistics. Acad. Press, London, 1979.
  * Masegosa (2020) Andres Masegosa.  Learning under model misspecification: Applications to variational and ensemble methods.  _Advances in Neural Information Processing Systems_ , 33:5479–5491, 2020.
  * Matsubara et al. (2022) Takuo Matsubara, Jeremias Knoblauch, François-Xavier Briol, and Chris J. Oates.  Robust Generalised Bayesian Inference for Intractable Likelihoods, January 2022.  arXiv:2104.07359 [math, stat].
  * Muandet et al. (2017) Krikamol Muandet, Kenji Fukumizu, Bharath Sriperumbudur, and Bernhard Schölkopf.  Kernel Mean Embedding of Distributions: A Review and Beyond.  _Foundations and Trends® in Machine Learning_ , 10(1-2):1–141, 2017.  doi: 10.1561/2200000060.
  * Murphy (2007) Kevin Murphy.  Conjugate Bayesian analysis of the Gaussian distribution, 11 2007.
  * Pacchiardi & Dutta (2022a) Lorenzo Pacchiardi and Ritabrata Dutta.  Likelihood-free inference with generative neural networks via scoring rule minimization.  _arXiv preprint arXiv:2205.15784_ , 2022a.
  * Pacchiardi & Dutta (2022b) Lorenzo Pacchiardi and Ritabrata Dutta.  Score Matched Neural Exponential Families for Likelihood-Free Inference, January 2022b.  arXiv:2012.10903 [stat].
  * Pang et al. (2022) Guansong Pang, Chunhua Shen, Longbing Cao, and Anton van den Hengel.  Deep Learning for Anomaly Detection: A Review.  _ACM Computing Surveys_ , 54(2):1–38, March 2022.  doi: 10.1145/3439950.  arXiv:2007.02500 [cs, stat].
  * Papamakarios & Murray (2016) George Papamakarios and Iain Murray.  Fast $\varepsilon$-free inference of simulation models with bayesian conditional density estimation.  _Advances in neural information processing systems_ , 29, 2016.
  * Papamakarios et al. (2019) George Papamakarios, David Sterratt, and Iain Murray.  Sequential neural likelihood: Fast likelihood-free inference with autoregressive flows.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pp. 837–848. PMLR, 2019.
  * Radev et al. (2020a) Stefan T Radev, Ulf K Mertens, Andreas Voss, Lynton Ardizzone, and Ullrich Köthe.  BayesFlow: Learning complex stochastic models with invertible neural networks.  _IEEE Transactions on Neural Networks and Learning Systems_ , 2020a.
  * Radev et al. (2020b) Stefan T. Radev, Andreas Voss, Eva Marie Wieschen, and Paul-Christian Buerkner.  Amortized Bayesian Inference for Models of Cognition, 2020b.
  * Radev et al. (2021) Stefan T Radev, Frederik Graw, Simiao Chen, Nico T Mutters, Vanessa M Eichel, Till Bärnighausen, and Ullrich Köthe.  OutbreakFlow: Model-based Bayesian inference of disease outbreak dynamics with invertible neural networks and its application to the COVID-19 pandemics in Germany.  _PLOS Computational Biology_ , 17(10):e1009472, 2021.
  * Radev et al. (2023a) Stefan T. Radev, Marvin Schmitt, Valentin Pratz, Umberto Picchini, Ullrich Köthe, and Paul-Christian Bürkner.  JANA: Jointly Amortized Neural Approximation of Complex Bayesian Models.  In Robin J. Evans and Ilya Shpitser (eds.), _Proceedings of the 39th Conference on Uncertainty in Artificial Intelligence_ , volume 216 of _Proceedings of Machine Learning Research_ , pp. 1695–1706. PMLR, 2023a.
  * Radev et al. (2023b) Stefan T. Radev, Marvin Schmitt, Lukas Schumacher, Lasse Elsemüller, Valentin Pratz, Yannik Schälte, Ullrich Köthe, and Paul-Christian Bürkner.  Bayesflow: Amortized bayesian workflows with neural networks.  _Journal of Open Source Software_ , 8(89):5702, 2023b.  ISSN 2475-9066.  doi: 10.21105/joss.05702.  URL http://dx.doi.org/10.21105/joss.05702.
  * Radev et al. (2023c) Stefan T. Radev, Marvin Schmitt, Lukas Schumacher, Lasse Elsemüller, Valentin Pratz, Yannik Schälte, Ullrich Köthe, and Paul-Christian Bürkner.  Bayesflow: Amortized bayesian workflows with neural networks.  _Journal of Open Source Software_ , 8(89):5702, 2023c.  doi: 10.21105/joss.05702.  URL https://doi.org/10.21105/joss.05702.
  * Ramesh et al. (2022) Poornima Ramesh, Jan-Matthis Lueckmann, Jan Boelts, Álvaro Tejero-Cantero, David S Greenberg, Pedro J Gonçalves, and Jakob H Macke.  Gatsbi: Generative adversarial training for simulation-based inference.  _arXiv preprint arXiv:2203.06481_ , 2022.
  * Ratcliff & McKoon (2008) Roger Ratcliff and Gail McKoon.  The Diffusion Decision Model: Theory and Data for Two-Choice Decision Tasks.  _Neural Computation_ , 20(4):873–922, April 2008.  doi: 10.1162/neco.2008.12-06-420.
  * Schmitt et al. (2024) Marvin Schmitt, Paul-Christian Bürkner, Ullrich Köthe, and Stefan T. Radev.  Detecting model misspecification in amortized bayesian inference with neural networks.  In Ullrich Köthe and Carsten Rother (eds.), _Pattern Recognition_ , pp. 541–557. Springer Nature, 2024.
  * Schmon et al. (2021) Sebastian M. Schmon, Patrick W. Cannon, and Jeremias Knoblauch.  Generalized Posteriors in Approximate Bayesian Computation, February 2021.  arXiv:2011.08644 [stat].
  * Shiono (2021) Takashi Shiono.  Estimation of agent-based models using Bayesian deep learning approach of BayesFlow.  _Journal of Economic Dynamics and Control_ , 125:104082, 2021.
  * Sorrenson et al. (2020) Peter Sorrenson, Carsten Rother, and Ullrich Köthe.  Disentanglement by nonlinear ica with general incompressible-flow networks (gin).  In _International Conference on Learning Representations_ , 2020.
  * Stan Development Team (2022) Stan Development Team.  The Stan Core Library, 2022.  Version 2.30.0.
  * Stine (1989) Robert Stine.  An Introduction to Bootstrap Methods.  _Sociological Methods & Research_, 18(2-3):243–291, November 1989.  doi: 10.1177/0049124189018002003.
  * Säilynoja et al. (2021) Teemu Säilynoja, Paul-Christian Bürkner, and Aki Vehtari.  Graphical Test for Discrete Uniformity and its Applications in Goodness of Fit Evaluation and Multiple Sample Comparison, November 2021.  arXiv:2103.10522 [stat].
  * Talts et al. (2020) Sean Talts, Michael Betancourt, Daniel Simpson, Aki Vehtari, and Andrew Gelman.  Validating Bayesian Inference Algorithms with Simulation-Based Calibration, October 2020.  arXiv:1804.06788 [stat].
  * Tejero-Cantero et al. (2020) Alvaro Tejero-Cantero, Jan Boelts, Michael Deistler, Jan-Matthis Lueckmann, Conor Durkan, Pedro J Gonçalves, David S Greenberg, and Jakob H Macke.  Sbi–a toolkit for simulation-based inference.  _arXiv preprint arXiv:2007.09114_ , 2020.
  * Thomas & Corander (2019) Owen Thomas and Jukka Corander.  Diagnosing model misspecification and performing generalized Bayes’ updates via probabilistic classifiers.  _arXiv preprint arXiv:1912.05810_ , 2019.
  * Tolstikhin et al. (2017) Ilya Tolstikhin, Olivier Bousquet, Sylvain Gelly, and Bernhard Schoelkopf.  Wasserstein auto-encoders, 2017.
  * Vehtari & Ojanen (2012) Aki Vehtari and Janne Ojanen.  A survey of Bayesian predictive methods for model assessment, selection and comparison.  _Statistics Surveys_ , 6(none), January 2012.  doi: 10.1214/12-SS102.
  * von Krause et al. (2022) Mischa von Krause, Stefan T. Radev, and Andreas Voss.  Mental speed is high until age 60 as revealed by analysis of over a million participants.  _Nature Human Behaviour_ , 6(5):700–708, May 2022.  doi: 10.1038/s41562-021-01282-7.
  * Ward et al. (2022) Daniel Ward, Patrick Cannon, Mark Beaumont, Matteo Fasiolo, and Sebastian M. Schmon.  Robust Neural Posterior Estimation and Statistical Model Criticism, October 2022.  arXiv:2210.06564 [cs, stat].
  * White (1982) Halbert White.  Maximum likelihood estimation of misspecified models.  _Econometrica_ , 50(1):1–25, 1982.
  * Wiqvist et al. (2021) Samuel Wiqvist, Jes Frellsen, and Umberto Picchini.  Sequential Neural Posterior and Likelihood Approximation.  _arXiv preprint arXiv:2102.06522_ , 2021.
  * Yang et al. (2021) Jingkang Yang, Kaiyang Zhou, Yixuan Li, and Ziwei Liu.  Generalized out-of-distribution detection: A survey.  _arXiv:2110.11334_ , 2021.
  * Zhang & Gao (2020) Fengshuo Zhang and Chao Gao.  Convergence rates of variational posterior distributions.  _The Annals of Statistics_ , 48(4), August 2020.  doi: 10.1214/19-AOS1883.

Appendix

##  Appendix A Theoretical Implications

Attaining the global minimum of the neural posterior estimation loss with an
arbitrarily expressive neural architecture
$\\{h_{\boldsymbol{\psi}^{*}},f_{\boldsymbol{\phi}^{*}},\mathcal{M}\\}$
implies that (i) the inference and summary network jointly amortize the
analytic posterior $p(\boldsymbol{\theta}\,|\,\boldsymbol{x},\mathcal{M})$;
and (ii) the summary network transforms the data
$p(\boldsymbol{x}\,|\,\mathcal{M})$ into a unit Gaussian in summary space:
$p\big{(}\boldsymbol{z}=h_{\boldsymbol{\psi}^{*}}(\boldsymbol{x})\big{)}=\mathcal{N}(\boldsymbol{z}\,|\,\mathbf{0},\mathbb{I})$.
According to (i), the set of inference network parameters
$\boldsymbol{\phi}^{*}$ is a minimizer of the (expected) negative log
posterior learned by the inference network,

$\boldsymbol{\phi}^{*}=\operatorname*{argmin}_{\boldsymbol{\phi}}\mathbb{E}_{p(\boldsymbol{z})}\mathbb{E}_{p(\boldsymbol{\theta}\,|\,\boldsymbol{z})}\Big{[}-\log
q_{\boldsymbol{\phi}}\big{(}\boldsymbol{\theta}\,|\,\boldsymbol{z}\big{)}\Big{]}.$
(13)

At the same time, (ii) ensures that deviances in the summary space (according
to MMD) imply differences in the data generating processes,

$\underbrace{\mathbb{MMD}^{2}\big{[}p^{*}(\boldsymbol{z})\,||\,p(\boldsymbol{z}\,|\,\mathcal{M})\big{]}>0}_{\text{summary
space
difference}}\implies\underbrace{\mathbb{MMD}^{2}\big{[}p^{*}(\boldsymbol{x})\,||\,p(\boldsymbol{x}\,|\,\mathcal{M})\big{]}>0}_{\text{data
space difference}},$ (14)

since a deviation of
$p^{*}(\boldsymbol{z}=h_{\boldsymbol{\psi}^{*}}(\boldsymbol{x}))$ from a unit
Gaussian means that the summary network is no longer transforming samples from
$p(\boldsymbol{x}\,|\,\mathcal{M})$.

Accordingly, the LHS of Eq. 14 no longer guarantees that the inference network
parameters $\boldsymbol{\phi}^{*}$ are maximally informative for posterior
inference. The preceding argumentation also motivates our augmented objective,
since a divergence of summary statistics for observed data
$\accentset{\text{o}}{\boldsymbol{z}}=h_{\boldsymbol{\psi}}(\accentset{\text{o}}{\boldsymbol{x}})$
from a unit Gaussian signalizes a deficiency in the assumed generative model
$\mathcal{M}$ and a need to revise the latter. We also hypothesize and show
empirically that we can successfully detect simulation gaps in practice even
when the summary network outputs have not exactly converged to a unit Gaussian
(e.g., in the presence of correlations in summary space, cf. Experiment 3 and
Experiment 4).

However, the converse of Eq. 14 is not true in general. In other words, a
discrepancy in data space (non-zero MMD on the RHS of Eq. 14) does _not_
generally imply a difference in summary space (non-zero MMD on the LHS of Eq.
14). To illustrate this via a counter-example, consider the assumed Gaussian
generative model $\mathcal{M}$ defined by

$\displaystyle\mu$ $\displaystyle\sim p(\mu),$ (15) $\displaystyle
x_{1},x_{2}$ $\displaystyle\sim\mathcal{N}(\mu,\sigma^{2}=2),$

for $N=2$ observations and a summary network with a single-output ($S=1$).
Since the variance is fixed, the only inference target is the mean $\mu$.

Then, an optimal summary network $\boldsymbol{\psi}^{*}$ outputs the minimal
sufficient summary statistic for recovering the mean, namely the empirical
average:
$h_{\boldsymbol{\psi}^{*}}(x_{1},x_{2})=\bar{x}\equiv(x_{1}+x_{2})/2$.
Consequently, the distribution in summary space is given as
$p(\bar{x})=\mathcal{N}(0,1)$.333This follows from the property
$\text{Var}(\bar{x})=\text{Var}((x_{1}+x_{2})/2)=(\text{Var}(x_{1})+\text{Var}(x_{2}))/2^{2}=1$
for independent $x_{1},x_{2}$. In terms of the MMD criterion, we see that

$\mathbb{MMD}^{2}\big{[}\underbrace{p(h_{\boldsymbol{\psi}^{*}}(x_{1},x_{2}))}_{=\mathcal{N}(0,1)}\,||\,\mathcal{N}(0,1)\big{]}=0.$
(16)

Now, suppose that the real data are actually generated by a different model
$\mathcal{M}^{*}$ with $\accentset{\text{o}}{x}\sim p^{*}(\boldsymbol{x})$
given by

$\displaystyle\mu$ $\displaystyle\sim p(\mu),$ (17)
$\displaystyle\accentset{\text{o}}{x}_{1}$
$\displaystyle\sim\mathcal{N}(\mu,\sigma^{2}=1),$
$\displaystyle\accentset{\text{o}}{x}_{2}$
$\displaystyle\sim\mathcal{N}(\mu,\sigma^{2}=3).$

Clearly, this process $p^{*}(\boldsymbol{x})$ differs from
$p(\boldsymbol{x}\,|\,\mathcal{M})$ (Eq. 15) on the data domain according to
the MMD metric:
$\mathbb{MMD}^{2}\big{[}p^{*}(\boldsymbol{x})\,||\,p(\boldsymbol{x}\,|\,\mathcal{M})\big{]}>0$.
However, using the same calculations as above, we find that the summary space
for the process $p^{*}(\boldsymbol{x})$ also follows a unit Normal
distribution: $p^{*}(\bar{x})=\mathcal{N}(0,1)$. Thus, the processes
$p^{*}(\boldsymbol{x})$ and $p(\boldsymbol{x}\,|\,\mathcal{M})$ are
indistinguishable in the summary space despite the fact that the first
generative model $\mathcal{M}$ is clearly misspecified.

The above example shows that learning _minimal sufficient summary statistics_
for solving the inference task (i.e., the mean in this example) might not be
optimal for detecting simulation gaps. In fact, increasing the output
dimensions $S$ of the summary network $h_{\boldsymbol{\psi}}$ would enable the
network to learn structurally richer (overcomplete) sufficient summary
statistics. The latter would be invariant to fewer misspecifications and thus
more useful for uncovering simulation gaps. In the above example, an
_overcomplete_ summary network with $S=2$ which simply copies and scales the
two variables by their corresponding variances is able to detect the
misspecification. Experiment 1 studies the influence of the number of summary
statistics in a controlled setting and provides empirical illustrations.
Experiments 3 and 4 further address the choice of the number of summary
statistics in more complex models of decision making and disease outbreak.
Next, we describe how to practically detect simulation gaps during inference
using only _finite realizations_ from $\mathcal{M}$ and $\mathcal{M}^{*}$.

##  Appendix B Posterior Inference Errors

Given a generative model $\mathcal{M}$, the _analytic posterior under the
potentially misspecified model_
$p(\boldsymbol{\theta}\,|\,\boldsymbol{x},\mathcal{M})$ always exists, even if
$\mathcal{M}$ is misspecified for the data $\boldsymbol{x}$. Obtaining a
trustworthy approximation of the analytic posterior is the fundamental basis
for any follow-up inference (e.g., parameter estimation or model comparison)
and must be at least an intermediate goal in real world applications. Assuming
optimal convergence under a misspecified model $\mathcal{M}$, the amortized
posterior
$q_{\boldsymbol{\phi}}\big{(}\boldsymbol{\theta}\,|\,\boldsymbol{z}=h_{\boldsymbol{\psi}}(\boldsymbol{x}),\mathcal{M}\big{)}$
still corresponds to the analytic posterior
$p(\boldsymbol{\theta}\,|\,\boldsymbol{x},\mathcal{M})$, as any transformed
$\accentset{\text{o}}{\boldsymbol{x}}$ arising from $p^{*}(\boldsymbol{x})$
has non-zero density in the latent Gaussian summary space.444We assume that we
have no hard limits in the prior or simulator in $\mathcal{M}$. Thus, the
posterior approximator should still be able to obtain the correct pushforward
density under $\mathcal{M}$ for any query
$\accentset{\text{o}}{\boldsymbol{x}}$. However, optimal convergence can never
be achieved after finite training time, so we need to address its implications
for the validity of amortized simulation-based posterior inference in
practice.

Given finite training data, the summary and inference networks will mostly see
simulations from the _typical set_
$\mathcal{T}(\mathcal{M})\subset\mathcal{X}$ of the generative model
$\mathcal{M}$, that is, training instances whose self-information $-\log
p(\boldsymbol{x}\,|\,\mathcal{M})$ is close to the entropy
$\mathbb{E}\big{[}-\log p(\boldsymbol{x}\,|\,\mathcal{M})\big{]}$. In high
dimensional problems, the typical set will comprise a rather small subset of
the possible outcome space, determined by a complex interaction between the
components of $\mathcal{M}$ (Betancourt, 2017). Accordingly, good convergence
in practice may mean that i) only observations from $\mathcal{T}(\mathcal{M})$
actually follow the approximate Gaussian in latent summary space and ii) the
inference network has only seen enough training examples in
$\mathcal{T}(\mathcal{M})$ to learn accurate posteriors for observables
$\boldsymbol{x}\in\mathcal{T}(\mathcal{M})$, but remains inaccurate for well-
specified, but rare $\boldsymbol{x}\not\in\mathcal{T}(\mathcal{M})$.

Since atypical or improbable outcomes occur rarely during simulation-based
training, they have negligible influence on the neural networks’ loss.
Consequently, posterior approximation errors for observations outside of
$\mathcal{T}(\mathcal{M})$ can be large, simply because the networks have not
yet converged in these unusual regions, and the highly non-linear mapping of
the inference network still deviates considerably from the true solution.
Better training methods might resolve this problem in the future, but for now
our proposed MMD criterion reliably signals low fidelity posterior estimates
by quantifying the “distance from the typical generative set”
$\mathcal{T}(\mathcal{M})$ in the structured summary space.

Moreover, we hypothesize and demonstrate empirically in the following
experiments that the difference between the correct posterior
$p(\boldsymbol{\theta}\,|\,\boldsymbol{x},\mathcal{M})$ and the approximate
posterior
$q_{\boldsymbol{\phi}}(\boldsymbol{\theta}\,|\,h_{\boldsymbol{\psi}}(\boldsymbol{x}),\mathcal{M})$
for misspecified models increases as a function of MMD, and thus the latter
also measures the amount of misspecification. Therefore, our MMD criterion
serves a dual purpose in practice: It can uncover potential simulation gaps
and, simultaneously, signal errors in posterior estimation of rare (but valid)
events.

##  Appendix C Description of the Method as Algorithm

Algorithm 1 Misspecification-aware amortized Bayesian inference. The algorithm
illustrates online learning for simplicity. The training phase can utilize any
other learning paradigm as well (e.g., round-based). Thus, details like batch
size or round indices are omitted in the training phase of the algorithm.

1:Training phase:

2:repeat

3: Sample parameters and data $(\boldsymbol{\theta},\boldsymbol{x})$ from the
specified generative model $\mathcal{M}$.

4: Pass the data $\boldsymbol{x}$ through the summary network:
$\boldsymbol{z}=h_{\boldsymbol{\psi}}(\boldsymbol{x})$.

5: Compute $-\log
q_{\boldsymbol{\phi}}\big{(}\boldsymbol{\theta}\,|\,\boldsymbol{z},\mathcal{M}\big{)}$.

6: Compute loss function.

7: Update neural network parameters $\boldsymbol{\psi},\boldsymbol{\phi}$ via
backpropagation.

8:until convergence to
$\widehat{\boldsymbol{\psi}},\widehat{\boldsymbol{\phi}}$

9:

10:Inference phase (given $N$ observed or test data sets
$\\{\accentset{\text{o}}{\boldsymbol{x}}^{(n)}\\}$, query $L$ draws from the
amortized posterior, use $M$ draws from the validation summary distribution):

11:for $m=1,\ldots,M$ do

12: Re-use data set $\boldsymbol{x}^{(m)}$ from the training phase or simulate
a new one from $\mathcal{M}$.

13: Pass the data set $\boldsymbol{x}^{(m)}$ through the converged summary
network:
$\boldsymbol{z}^{(m)}=h_{\widehat{\boldsymbol{\psi}}}(\boldsymbol{x}^{(m)})$.

14:end for

15:for $n=1,\ldots,N$ do

16: Pass the observed data set $\accentset{\text{o}}{\boldsymbol{x}}^{(n)}$
through the summary network:
$\accentset{\text{o}}{\boldsymbol{z}}^{(n)}=h_{\widehat{\boldsymbol{\psi}}}(\accentset{\text{o}}{\boldsymbol{x}}^{(n)})$.

17: Pass $\accentset{\text{o}}{\boldsymbol{z}}^{(n)}$ through the posterior
approximator for $L$ draws $\\{\boldsymbol{\theta}^{(n,l)}\\}$ from
$q_{\widehat{\boldsymbol{\phi}}}(\boldsymbol{\theta}\,|\,\accentset{\text{o}}{\boldsymbol{z}}^{(n)})$

18:end for

19:Estimate the MMD distance of the inference data
$\\{\accentset{\text{o}}{\boldsymbol{z}}^{(n)}\\}$ from the validation summary
space under the generative model $\mathcal{M}$ from training:
$\widehat{\mathrm{MMD}}^{2}\left(\left\\{\accentset{\text{o}}{\boldsymbol{z}}^{(n)}\right\\}\,\middle|\middle|\,\left\\{\boldsymbol{z}^{(m)}\right\\}\right)$.

20:return Draws from the approximate posterior distribution for each queried
data set $\accentset{\text{o}}{\boldsymbol{x}}^{(n)}$
$\\{\boldsymbol{\theta}^{(n,l)}\\}$ and
$\widehat{\mathrm{MMD}}^{2}\left(\left\\{\accentset{\text{o}}{\boldsymbol{z}}^{(n)}\right\\}\,\middle|\middle|\,\left\\{\boldsymbol{z}^{(m)}\right\\}\right)$.

##  Appendix D Additional Experiments

###  D.1 Experiment 6: 5D Gaussian Means and Covariance

(a) Prior misspecification.

(b) Simulator and noise misspecification.

Figure D.1: Experiment 6. MMD as a function of model misspecification severity
in the prior (left) as well as simulator and noise (right). All induced model
misspecifications are detectable through an increased MMD. White stars
represent the configuration for the well-specified model (i.e., training model
$\mathcal{M}$). The deviating pattern in the top-most row of (1(b)) is caused
by the infinite variance of the Student $t$ likelihood with $\nu=2$ degrees of
freedom: This is an extreme simulation gap with respect to the unit Gaussian
model $\mathcal{M}$, and consequently detected as such.

This experiments extends the Gaussian conjugate model to higher dimensions and
a more difficult task, i.e., recovering the means and the full covariance
matrix of a $5$-dimensional Gaussian. There is a total of $20$ inference
parameters—5 means and 15 (co-)variances—meaning that $20$ summary statistics
would suffice to solve the inference task. The mean vector $\boldsymbol{\mu}$
and the covariance matrix $\boldsymbol{\Sigma}$ are drawn from a joint prior,
namely a normal-inverse-Wishart distribution ($\text{N-}\mathcal{W}^{-1}$;
Barnard et al., 2000). The normal-inverse-Wishart prior
$\text{N-}\mathcal{W}^{-1}(\boldsymbol{\mu},\boldsymbol{\Sigma}\,|\,\boldsymbol{\mu}_{0},\lambda_{0},\boldsymbol{\Psi}_{0},\nu_{0})$
implies a hierarchical process, and the implementation details are described
in Section E.1 in the Appendix. We set the number of summary statistics to
$S=40$ to balance the trade-off between posterior error and misspecification
detection (see Experiment 1). The model $\mathcal{M}$ used for training the
networks as well as the induced model misspecifications (prior, simulator, and
noise) are detailed in the Appendix.

Results. The converged posterior approximator can successfully recover the
analytic posterior for all inference parameters when no model misspecification
is present. Thus, our method does not impede posterior inference when the
models are well-specified. Since the summary space comprises $S=40$
dimensions, visual inspection is no longer feasible, and we resort to the
proposed MMD criterion. Both induced prior misspecifications—i.e., location
and variance—are detectable through an increased MMD (see 1(a)). Model
misspecifications via a heavy-tailed simulator—i.e., Student-$t$ with $\nu=2$
degrees of freedom—, as well as Beta noise, are also detectable with our MMD
criterion (see 1(b)).

##  Appendix E Implementation Details

###  E.1 5D Gaussian with Covariance (Experiment 6)

The normal-inverse-Wishart prior
$\text{N-}\mathcal{W}^{-1}(\boldsymbol{\mu},\boldsymbol{\Sigma}\,|\,\boldsymbol{\mu}_{0},\lambda_{0},\boldsymbol{\Psi}_{0},\nu_{0})$
implies a hierarchical prior. Suppose the covariance matrix
$\boldsymbol{\Sigma}$ has an inverse Wishart distribution
$\mathcal{W}^{-1}(\boldsymbol{\Sigma}\,|\,\boldsymbol{\Psi}_{0},\nu_{0})$ and
the mean vector $\boldsymbol{\mu}$ has a multivariate normal distribution
$\mathcal{N}(\boldsymbol{\mu}\,|\,\boldsymbol{\mu}_{0},\frac{1}{\lambda_{0}}\boldsymbol{\Sigma})$,
then the tuple $(\boldsymbol{\mu},\boldsymbol{\Sigma})$ has a normal-inverse-
Wishart distribution
$\text{N-}\mathcal{W}^{-1}(\boldsymbol{\mu},\boldsymbol{\Sigma}\,|\,\boldsymbol{\mu}_{0},\lambda_{0},\boldsymbol{\Psi}_{0},\nu_{0})$.
Finally, the likelihood is
Gaussian:$\boldsymbol{x}_{k}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})\;\text{for}\;k=1,\ldots,K$.

For a multivariate Gaussian with unknown mean and unknown covariance matrix,
the analytic joint posterior
$p(\boldsymbol{\mu}_{p},\boldsymbol{\Sigma}_{p}\,|\,\boldsymbol{x})$ follows a
normal-inverse Wishart distribution again:

$\displaystyle(\boldsymbol{\mu}_{p},\boldsymbol{\Sigma}_{p}\,|\,\boldsymbol{x})$
$\displaystyle\sim\text{N-}\mathcal{W}^{-1}(\boldsymbol{\mu}_{p},\boldsymbol{\Sigma}_{p}\,|\,\boldsymbol{\mu}_{K},\lambda_{K},\boldsymbol{\Psi}_{K},\nu_{K})\quad\text{with}$
(18) $\displaystyle\boldsymbol{\mu}_{K}$
$\displaystyle=\dfrac{\lambda_{0}\boldsymbol{\mu}_{0}+K\bar{\boldsymbol{x}}}{\lambda_{0}+K},\quad\lambda_{K}=\lambda_{0}+K,\quad\nu_{K}=\nu_{0}+K,$
$\displaystyle\boldsymbol{\Psi}_{K}$
$\displaystyle=\boldsymbol{\Psi}_{0}+\sum\limits_{k=1}^{K}(\boldsymbol{x}_{k}-\bar{\boldsymbol{x}})(\boldsymbol{x}_{k}-\bar{\boldsymbol{x}})^{T}+\dfrac{\lambda_{0}K}{\lambda_{0}+K}(\bar{\boldsymbol{x}}-\boldsymbol{\mu}_{0})(\bar{\boldsymbol{x}}-\boldsymbol{\mu}_{0})^{T}$

The marginal posteriors for $\boldsymbol{\mu}_{p}$ and
$\boldsymbol{\Sigma}_{p}$ then follow as (Murphy, 2007):

$\displaystyle\boldsymbol{\mu}_{p}$ $\displaystyle\sim
t_{\nu_{K}-D-1}\Big{(}\boldsymbol{\mu}_{p}\,\Big{|}\,\boldsymbol{\mu}_{K},\frac{\boldsymbol{\Psi}_{K}^{-1}}{\lambda_{K}(\nu_{K}-D+1)}\Big{)}$
(19) $\displaystyle\boldsymbol{\Sigma}_{p}$
$\displaystyle\sim\mathcal{W}^{-1}(\boldsymbol{\Sigma}_{p}\,|\,\boldsymbol{\Psi}_{K},\nu_{K})$

The model $\mathcal{M}$ used for training the networks as well as the types of
induced model misspecifications are outlined in Table 5.

Model | Prior | Likelihood
---|---|---
$\mathcal{M}$ (No MMS) | $\big{(}\boldsymbol{\mu},\boldsymbol{\Sigma}\big{)}\sim\text{N-}\mathcal{W}^{-1}(\boldsymbol{0},5,\mathbb{I},10)$ | $\boldsymbol{x}_{k}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$
$\mathcal{M}_{P}$ (Prior) | $\big{(}\boldsymbol{\mu},\boldsymbol{\Sigma}\big{)}\sim\text{N-}\mathcal{W}^{-1}(\boldsymbol{\mu}_{0},5,\tau_{0}\mathbb{I},10)$ | $\boldsymbol{x}_{k}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$
$\mathcal{M}_{S}$ (Simulator) | $\big{(}\boldsymbol{\mu},\boldsymbol{\Sigma}\big{)}\sim\text{N-}\mathcal{W}^{-1}(\boldsymbol{0},5,\mathbb{I},10)$ | $\boldsymbol{x}_{k}\sim t_{\text{df}}(\boldsymbol{\mu},\boldsymbol{\Sigma}),\quad\text{df}\in\mathbb{N}_{>0}$
$\mathcal{M}_{N}$ (Noise) | $\big{(}\boldsymbol{\mu},\boldsymbol{\Sigma}\big{)}\sim\text{N-}\mathcal{W}^{-1}(\boldsymbol{0},5,\mathbb{I},10)$ | $\boldsymbol{x}_{k}\sim\lambda\cdot\mathrm{Beta}(2,5)+(1-\lambda)\cdot\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$
Table 5: Investigated model misspecifications (MMS) for the $5$-dimensional
Gaussian with fully estimated covariance matrix.
$\text{N-}\mathcal{W}^{-1}(\boldsymbol{\mu}_{0},\lambda_{0},\boldsymbol{\Psi},\nu)$
denotes a multivariate Normal-Inverse-Wishart distribution with location
$\boldsymbol{\mu}_{0}$, precision scaling factor $\lambda_{0}$, precision
$\boldsymbol{\Psi}$, and $\nu$ degrees of freedom.

In the evaluation, we compare the means of the approximate posterior samples
with the first moment of the respective marginal analytic posterior from Eq.
19. We evaluate correlation matrices with standard deviations on the diagonal.
For the $t$ distributed posterior mean and inverse-Wishart distributed
posterior covariance, we obtain (Mardia et al., 1979):

$\mathbb{E}(\boldsymbol{\mu}_{p})=\boldsymbol{\mu}_{K},\quad\mathbb{E}(\boldsymbol{\Sigma}_{p})=\frac{\boldsymbol{\Psi}_{K}}{\nu_{K}-D-1}$
(20)

###  E.2 Cancer and Stromal Cells (Experiment 2)

As implemented by Ward et al. (2022), the CS model simulates the development
of cancer and stromal cells. The respective cell counts (total $N_{c}$,
unobserved parents $N_{p}$, daughters for each parent $N_{d}^{(i)}$) are
stochastically defined through Poisson distributions

$N_{c}\sim\mathrm{Poisson}(\lambda_{c}),\quad
N_{p}\sim\mathrm{Poisson}(\lambda_{c}),\quad
N_{d}^{(i)}\sim\mathrm{Poisson}(\lambda_{c}),\quad i=1,\ldots,N_{p}.$ (21)

The distance based metrics for the hand-crafted summary statistics 3–4 are
empirically approximated from 50 stromal cells to avoid computing the full
distance matrix (Ward et al., 2022). The prior distributions are specified as

$\lambda_{c}\sim\mathrm{Gamma}(25,0.03),\quad\lambda_{d}\sim\mathrm{Gamma}(5,0.5),\quad\lambda_{p}\sim\mathrm{Gamma}(45,3)$
(22)

where $\mathrm{Gamma}(a,b)$ denotes the Gamma distribution with location $a$
and rate $b$. To induce model misspecification through necrosis, a Bernoulli
variable $w_{i}\sim\mathrm{Bernoulli}(\pi)$ is sampled, removing cancer cells
within a specified radius around the parent cell (Ward et al., 2022). Thus,
the Bernoulli parameter $\pi$ controls the degree of misspecification, ranging
from no misspecification ($\pi=0$; no necrosis) to maximal misspecification
with respect to that parameter ($\pi=1$; necrosis of all susceptible cells).

###  E.3 DDM (Experiment 3)

The starting point of the evidence accumulation process is unbiased,
$x_{t=0}=\frac{a}{2}$. During training, all parameters are drawn from Gamma
prior distributions:

$v_{1},v_{2},a_{1},a_{2},t_{0}\sim\Gamma(5,0.5).$ (23)

We first generate uncontaminated data $\boldsymbol{x}^{*}$ from the well-
specified generative model $\mathcal{M}$. Second, we randomly choose a
fraction $\lambda\in[0,1]$ of the data $\boldsymbol{x}^{*}$. Third, we replace
this fraction with data-dependent contaminants $\xi$,

$\displaystyle\text{Fast guesses:}\quad\xi$
$\displaystyle\sim\mathcal{U}\big{(}0.1,Q_{10}(\boldsymbol{x}^{*})\big{)}$
(24) $\displaystyle\text{Slow responses:}\quad\xi$
$\displaystyle\sim\mathcal{U}\big{(}Q_{75}(\boldsymbol{x}^{*}),10\big{)},$

where $Q_{k}(\boldsymbol{x}^{*})$ denotes the $k^{\text{th}}$ percentile of
$\boldsymbol{x}^{*}$. The asymmetry in percentiles between fast and slow
responses arises from the inherent positive skewness of reaction time
distributions. The fixed upper limit of slow response contamination is
motivated by the maximum number of iterations of the utilized diffusion model
simulator. The contamination procedure is executed separately for each
condition and response type. If an experiment features both fast and slow
contamination, the fraction $\lambda$ is equally split between fast and slow
contamination. The uncontaminated data set is generated once and acts as a
baseline for all analyses of an experiment, resulting in a baseline MMD of 0
since $\boldsymbol{x}^{*}$ is unaltered if $\lambda=0$.

##  Appendix F Bootstrapping Procedure

In Experiment 4, we estimate a sampling distribution of the MMD between
samples from the specified training model $\mathcal{M}$ with
$\boldsymbol{x}\sim p(\boldsymbol{x}\,|\,\mathcal{M}^{*})$ and samples from
the (opaque) observed model $\mathcal{M}^{*}$ with
$\accentset{\text{o}}{x}\sim p^{*}(\boldsymbol{x})$. Since simulating time
series from the compartmental models is time-consuming, we opt for
bootstrapping (Stine, 1989) on $M=1\,000$ pre-simulated time series
$\\{\boldsymbol{x}^{(m)}\\}$ from $\mathcal{M}$ and $N=1\,000$ pre-simulated
time series $\\{\accentset{\text{o}}{\boldsymbol{x}}^{(n)}\\}$ from
$\mathcal{M}^{*}$. In each bootstrapping iteration, we draw $M=1\,000$ samples
with replacement from $\\{\boldsymbol{x}^{(m)}\\}$ as well as
$N_{B}\in\\{1,2,5\\}$ samples (with replacement) from
$\\{\accentset{\text{o}}{\boldsymbol{x}}^{(n)}\\}$ and calculate the MMD
between the sets of bootstrap samples.

##  Appendix G 2D Gaussian Means: Further Results

###  G.1 Performance under well-specified model

(a) Excellent parameter recovery.

(b) Excellent calibration.

Figure G.1: Experiment 1. When the assumed model is well-specified for the
observed data ($\mathcal{M}=\mathcal{M}^{*}$), both posterior recovery (1(a))
and calibration (1(b)) remain excellent with our adjusted optimization
objective.

###  G.2 Overcomplete Summary Statistics

Figure G.2 shows the latent summary space when overcomplete sufficient summary
statistics ($S=4$) are used in Experiment 1 to recover the means of a
$2$-dimensional Gaussian. Model misspecification with respect to both
simulator and noise is detectable through anomalies in the latent summary
space. A network with $S=2$ summary statistics and otherwise equivalent
architecture could not capture these types of model misspecification.

Figure G.2: Pairplot of $10\,000$ latent summary space samples from the
overcomplete summary network. Both noise (orange) and simulator (pink)
misspecifications are distinguishable from the typical latent generative space
(blue).

##  Appendix H Replication of Experiment 1 with SNPE-C

Here, we show the results of repeating Experiment 1 with SNPE-C (Greenberg et
al., 2019) in the sbi implementation (Tejero-Cantero et al., 2020) for
posterior inference, instead of NPE in the BayesFlow implementation (Radev et
al., 2020a). The results are largely equivalent to those obtained with NPE in
BayesFlow.

Figure H.1: Experiment 1, SNPE-C. Summary space samples for the minimal
sufficient summary network ($S=2$) from a well-specified model $\mathcal{M}$
(blue) and misspecified configurations. Left: Prior misspecification can be
detected. Right: Simulator scale misspecification is indistinguishable from
the validation summary statistics.

|  | Model Misspecification |
---|---|---|---
|  |  Prior ($\mathcal{M}_{P}$) |  Simulator ($\mathcal{M}_{S}$) & noise ($\mathcal{M}_{N}$) |
Summary Network |  minimal |  |  |
overcomplete |  |  |

Figure H.2: Experiment 1, SNPE-C. Summary space MMD as a function of
misspecification severity. White stars indicate the well-specified model
configuration (i.e., equal to the training model $\mathcal{M}$).

##  Appendix I Epidemiological COVID-19 Model: Power Analysis

| $N=1$ | $N=2$ | $N=5$
---|---|---|---
$\mathcal{M}_{1}$ |  |  |
$\mathcal{M}_{2}$ |  |  |
$\mathcal{M}_{3}$ |  |  |

Figure I.1: Experiment 4. Detailed illustration of the power analysis. As few
as $N=5$ observed data sets suffice to achieve a negligible type II error
$\beta\approx 0$.

##  Appendix J Style-Inference of EMNIST Image Simulator: Additional Results

(a) In distribution (digits 0–4)

(b) Extreme styles (digits 0–4)

(c) Domain shift (digit 5)

(d) Domain shift (digit 6)

(e) Domain shift (digit 7)

(f) Domain shift (digit 8)

(g) Domain shift (digit 9)

Figure J.1: Experiment 5, variation 1 (class-agnostic), parameter recovery.

(a) In distribution (digits 0–4)

(b) Extreme styles (digits 0–4)

(c) Domain shift (digit 5)

(d) Domain shift (digit 6)

(e) Domain shift (digit 7)

(f) Domain shift (digit 8)

(g) Domain shift (digit 9)

Figure J.2: Experiment 5, variation 2 (class-informed), parameter recovery.