Misspecification-robust Sequential Neural Likelihood for Simulation-based Inference

The traditional statistical approach to conducting inference on $\bm{\theta}$ in this setting is to use ABC methods. Using the assumed DGP, these methods search for values of $\bm{\theta}$ that produce pseudo-data $\bm{x}$ which is “close enough” to $\bm{y}$, and then retain these values to build an approximation to the posterior. The comparison is generally carried out using summaries of the data to ensure the problem is computationally feasible. Let $S:\mathbb{R}^{n}\rightarrow\mathbb{R}^{d}$, denote the vector summary statistic mapping used in the analysis, where $d\geq d_{\bm{\theta}}$, and $n\geq d$.

Two prominent statistical approaches for SBI are ABC and BSL. ABC approximates the likelihood for the summaries via the following:

where $\rho\{S(\bm{y}),S(\bm{x})\}$ measures the discrepancy between observed and simulated summaries and $K_{\epsilon}(\cdot)$ is a kernel that allocates higher weight to smaller $\rho$. The bandwidth of the kernel, $\epsilon$, is often referred to as the tolerance in the ABC literature. The above integral is intractable, but can be estimated unbiasedly by drawing $m\geq 1$ mock datasets $\bm{x_{1}},\ldots,\bm{x_{m}}\sim P^{(n)}_{\bm{\theta}}$ and computing

It is common to set $m=1$ and choose the indicator kernel function,

Using arguments from the exact-approximate literature (Andrieu & Roberts, 2009), unbiasedly estimating the ABC likelihood leads to a Bayesian algorithm that samples from the approximate posterior proportional to $g_{\epsilon}({S(\bm{y})\mid\bm{\theta}})\pi(\bm{\theta})$. As is evident from the above integral estimator, ABC non-parametrically estimates the summary statistic likelihood. Unfortunately, this non-parametric approximation causes ABC to exhibit the “curse of dimensionality”, meaning that the probability of a proposed parameter being accepted decreases dramatically as the dimension of the summary statistics increases (Barber et al., 2015; Csilléry et al., 2010).

In contrast, BSL uses a parametric estimator. The most common BSL approach approximates $g_{n}(\cdot\mid\bm{\theta})$ using a Gaussian:

where $\mu(\bm{\theta})=\mathbb{E}[S(\bm{x})\mid\bm{\theta}]$ and $\Sigma(\bm{\theta})=\mbox{Var}(S(\bm{x})\mid\bm{\theta})$ denote the mean and variance of the model summary statistic at $\bm{\theta}$. In almost all practical cases $\mu(\bm{\theta})$ and $\Sigma(\bm{\theta})$ are unknown, but we can replace these quantities with those estimated from $m$ independent model simulations, using for example the sample mean and variance:

and where each simulated data set $\bm{x^{i}}$, $i=1,\dots,m$, is generated i.i.d. from $P^{(n)}_{\bm{\theta}}$. The synthetic likelihood is then approximated as

Unlike ABC, $\hat{g}_{A}({S(\bm{y})\mid\bm{\theta}})$ is not an unbiased estimator of $g_{A}(S(\bm{y})\mid\bm{\theta})$. David T. Frazier & Kohn (2023) demonstrate that if the summary statistics are sub-Gaussian, then the choice of $m$ is immaterial so long as $m$ diverges as $n$ diverges. The insensitivity to $m$ is supported empirically in Price et al. (2018), provided that $m$ is chosen large enough so that the plug-in synthetic likelihood estimator has a small enough variance to ensure that MCMC mixing is not adversely affected. The Gaussian assumption can be limiting; however, neural likelihoods provide a more flexible alternative to BSL while retaining the advantages of a parametric approximation.

Unfortunately, both ABC and BSL are inefficient in terms of the number of model simulations required to produce posterior approximations. In particular, most algorithms for ABC and BSL are wasteful in the sense that they use a relatively large number of model simulations associated with rejected parameter proposals. In contrast, methods have been developed in machine learning that utilise all model simulations to learn either the likelihood (e.g. Papamakarios et al. (2019)), posterior (e.g. Greenberg et al. (2019); Papamakarios & Murray (2016)) or likelihood ratio (e.g. Durkan et al. (2020); Hermans et al. (2020); Thomas et al. (2022)). We consider the SNL method of Papamakarios et al. (2019) in more detail in Section 3.

Neural SBI methods have seen rapid advancements, with various approaches approximating the likelihood (Boelts et al., 2022; Wiqvist et al., 2021). These methods include diffusion models for approximating the score of the likelihood (Simons et al., 2023), energy-based models for surrogating the likelihood (Glaser et al., 2023) and a “neural conditional exponential family” trained via score matching (Pacchiardi & Dutta, 2022). Bon et al. (2022) present a method for refining approximate posterior samples to minimise bias and enhance uncertainty quantification by optimising a transform of the approximate posterior, which maximises a scoring rule.

One way to categorise neural SBI methods is by differentiating between amortised and sequential sampling schemes. These methods differ in their proposal distribution for $\bm{\theta}$. Amortised methods estimate the neural density for any $\bm{x}$ within the support of the prior predictive distribution. This allows the trained flow to approximate the posterior for any observed statistic, making it efficient when analysing multiple datasets. However, this requires using the prior as the proposal distribution. When the prior and posterior differ significantly, there will be few training samples of $\bm{x}$ close to $\bm{y}$, resulting in the trained flow potentially being less accurate in the vicinity of the observed statistic.

When we apply a summary statistic mapping, we need to redefine the usual notion of model misspecification, i.e., no value of $\bm{\theta}\in\Theta$ such that $P^{(n)}_{\bm{\theta}}=P_{0}^{(n)}$, as it is still possible for $P^{(n)}_{\bm{\theta}}$ to generate summary statistics that match the observed statistic even if the model is incorrect (Frazier et al., 2020b). We define $\bm{b}(\bm{\theta})=\mathbb{E}[S(\bm{x})\mid\bm{\theta}]$ and $\bm{b_{0}}=\mathbb{E}[S(\bm{y})]$ as the expected values of the summary statistic with respect to the probability measures $P^{(n)}_{\bm{\theta}}$ and $P^{(n)}_{0}$, respectively. The meaningful notion of misspecification in SBI is when no $\bm{\theta}\in\Theta$ satisfies $\bm{b}(\bm{\theta})=\bm{b_{0}}$, implying there is no parameter value for which the expected simulated and observed summaries match. This is the definition of incompatibility proposed in Marin et al. (2014).

Wilkinson (2013) reframes the approximate ABC posterior as an exact result for a different distribution that incorporates an assumption of model error (i.e. model misspecification). Ratmann et al. (2009) proposes an ABC method that augments the likelihood with unknown error terms, allowing model criticism by detecting summary statistics that require high tolerances. Frazier et al. (2020a) incorporate adjustment parameters, inspired by RBSL, in an ABC context.

The behaviour of ABC and BSL under incompatibility is now well understood. In the context of ABC, we consider the model to be misspecified if

for some metric $\rho$, and the corresponding pseudo-true parameter is defined as

Frazier et al. (2020b) show, under various conditions, the ABC posterior concentrates onto $\bm{\theta_{0}}$ for large sample sizes, providing an inherent robustness to model misspecification. However, they also demonstrate that the asymptotic shape of the ABC posterior is non-Gaussian and credible intervals lack valid frequentist coverage; i.e., confidence sets do not have the correct level under $P^{(n)}_{0}$.

In the context of BSL, Frazier et al. (2021) show that when the model is incompatible, i.e. $b(\bm{\theta})\neq\bm{b_{0}}$ $\forall\bm{\theta}\in\Theta$, the Kullback-Leibler divergence between the true data generating distribution and the Gaussian distribution associated with the synthetic likelihood diverges as $n$ diverges. In BSL, we say that the model is incompatible if

We define

The behaviour of BSL under misspecification depends on the number of roots of $M_{n}(\bm{\theta})=\bm{0}$. If there is a single solution, and under various assumptions, the BSL posterior will concentrate onto the pseudo-true parameter $\bm{\theta_{0}}$ with an asymptotic Gaussian shape, and the BSL posterior mean satisfies a Bernstein von-Mises result. However, if there are multiple solutions to $M_{n}(\bm{\theta})=\bm{0}$, then the BSL posterior will asymptotically exhibit multiple modes that do not concentrate on $\bm{\theta_{0}}$. The number of solutions to $M_{n}(\bm{\theta})=\bm{0}$ for a given problem is not known a priori and is very difficult to explore.

In addition to the theoretical issues faced by BSL under misspecification, there are also computational challenges. Frazier & Drovandi (2021) point out that under incompatibility, since the observed summary lies in the tail of the estimated synthetic likelihood for any value of $\bm{\theta}$, the Monte Carlo estimate of the likelihood suffers from high variance. Consequently, a significantly large value of $m$ is needed to enable the MCMC chain to mix and avoid getting stuck, which is computationally demanding.

Due to the undesirable properties of BSL under misspecification, Frazier & Drovandi (2021) propose RBSL as a way to identify incompatible statistics and make inferences more robust simultaneously. RBSL is a model expansion that introduces auxiliary variables, represented by the vector $\bm{\Gamma}=(\gamma_{1},\ldots,\gamma_{d})^{\top}$. These variables shift the means (RBSL-M) or inflate the variances (RBSL-V) in the Gaussian approximation, ensuring that the extended model is compatible by absorbing any misspecification. This approach guarantees that the observed summary does not fall far into the tails of the expanded model. However, the expanded model is now overparameterised since the dimension of the combined vector $(\bm{\theta},\bm{\Gamma})^{\top}$ is larger than the dimension of the summary statistics.

To regularise the model, Frazier & Drovandi (2021) impose a prior distribution on $\bm{\Gamma}$ that favours compatibility. However, each component of $\bm{\Gamma}$’s prior has a heavy tail, allowing it to absorb the misspecification for a subset of the summary statistics. This method identifies the statistics the model is incompatible with while mitigating their influence on the inference. Frazier & Drovandi (2021) demonstrate that under compatibility, the posterior for $\bm{\Gamma}$ is the same as its prior, so that incompatibility can be detected by departures from the prior.

The mean adjusted synthetic likelihood is denoted

where $\sigma(\bm{\theta})=\mathrm{diag}\{\Sigma(\bm{\theta})\}$ is the vector of estimated standard deviations of the model summary statistics, and $\circ$ denotes the Hadamard (element-by-element) product. The role of $\sigma(\bm{\theta})$ is to ensure that we can treat each component of $\bm{\Gamma}$ as the number of standard deviations that we are shifting the corresponding model summary statistic.

Frazier & Drovandi (2021) suggest using a prior where $\bm{\theta}$ and $\bm{\Gamma}$ are independent, with the prior density for $\bm{\Gamma}$ being a Laplace prior with scale $\lambda$ for each $\gamma_{j}$. The prior is chosen because it is peaked at zero but has a moderately heavy tail. Sampling the joint posterior can be done using a component-wise MCMC algorithm that iteratively updates using the conditionals $\bm{\theta}\mid S,\bm{\Gamma}$ and $\bm{\Gamma}\mid S,\bm{\theta}$. The update for $\bm{\Gamma}$ holds the $m$ model simulations fixed and uses a slice sampler, resulting in an acceptance rate of one without requiring tuning a proposal distribution. Frazier & Drovandi (2021) find empirically that sampling over the joint space $(\bm{\theta},\bm{\Gamma})^{\top}$ does not slow down mixing on the $\bm{\theta}$-marginal space. In fact, in cases of misspecification, the mixing is substantially improved as the observed value of the summaries no longer falls in the tail of the Gaussian distribution.

Recent developments in SBI have focused on detecting and addressing model misspecification for both neural posterior estimation (Ward et al., 2022) and for amortised and sequential inference (Schmitt et al., 2023a). Schmitt et al. (2023a) employs a maximum mean discrepancy (MMD) estimator to detect a “simulation gap” between observed and simulated data, while Ward et al. (2022) detects and corrects for model misspecification by introducing an error model $\pi(S(\bm{y})\mid S(\bm{x}))$. Noise is added directly to the summary statistics during training in Bernaerts et al. (2023). Huang et al. (2023) noted that as incompatibility is based on the choice of summary statistics, if the summary statistics are learnt via a NN, training this network with a regularised loss function that penalises statistics with a mismatch between the observed and simulated values will lead to robust inference. Glöckler et al. (2023), again focused on the case of learnt (via NN) summaries, propose a scheme for robustness against adversarial attacks (i.e. small worst-case perturbations) on the observed data. Schmitt et al. (2023c) introduce a meta-uncertainty framework that blends real and simulated data to quantify uncertainty in posterior model probabilities, applicable to SBI with potential model misspecification

Generalised Bayesian inference (GBI) is an alternative class of methods suggested to handle model misspecification better than standard Bayesian methods (Knoblauch et al., 2022). Instead of targeting the standard Bayesian posterior, the GBI framework targets a generalised posterior,

where $\ell(\bm{y},\bm{\theta})$ is some loss function and $w$ is a tuning parameter that needs to be calibrated appropriately (Bissiri et al., 2016). Various approaches have applied GBI to misspecified models with intractable likelihoods (Chérief-Abdellatif & Alquier, 2020; Matsubara et al., 2022; Pacchiardi & Dutta, 2021; Schmon et al., 2021). Gao et al. (2023) extends GBI to the amortised SBI setting, using a regression neural network to approximate the loss function, achieving favourable results for misspecified examples. Dellaporta et al. (2022) employ similar ideas to GBI for an MMD posterior bootstrap.

Here we consider the contaminated normal example from Frazier & Drovandi (2021) to assess how SNL and RSNL perform under model misspecification. In this example, the DGP is assumed to follow:

where $i=1,\ldots,n$. However, the actual DGP follows:

The sufficient statistic for $\theta$ under the assumed DGP is the sample mean, $S_{1}(\bm{y})=\frac{1}{n}\sum_{i=1}^{n}y_{i}$. For demonstration purposes, let us also include the sample variance, $S_{2}(\bm{y})=\frac{1}{n-1}\sum_{i=1}^{n}(y_{i}-S_{1}(\bm{y}))^{2}$. When $\sigma_{\epsilon}\neq 1$, we are unable to replicate the sample variance under the assumed model. We use the prior, $\theta\sim\mathcal{N}(0,10^{2})$ and $n=100$. The actual DGP is set to $\omega=0.8$ and $\sigma_{\epsilon}=2.5$, and hence the sample variance is incompatible. Since $S_{1}(\bm{y})$ is sufficient, so is $S(\bm{y})$, and one might still be optimistic that useful inference will result. We thus want our robust algorithm to concentrate the posterior around the sample mean. Under the assumed DGP we have that $b(\theta)=(\theta,1)^{\top}$, for all $\theta\in\mathbb{R}$. We thus have $\inf_{\theta\in\mathbb{R}}||b(\theta)-b_{0}||>0$ meaning our model is misspecified. We include additional results for the contaminated normal, included in Appendix D, due to the ease of obtaining an analytical true posterior. We also consider this example with no summarisation (i.e. using 100 draws directly) to see how RSNL scales as we increase the number of summaries, with results found in Appendix F.

For the contaminated normal example, Figure 1 demonstrates that RSNL yields reliable inference, with high posterior density around the true parameter value, $\theta=1$. In contrast, SNL provides unreliable inference, offering minimal support near the true value. The posteriors for the components of $\bm{\Gamma}$ are also shown. The prior and posterior for $\gamma_{1}$ (linked to the compatible summary statistic) are nearly identical, aligning with RBSL behaviour. For $\gamma_{2}$ (associated with the incompatible statistic), misspecification is detected as the posterior exhibits high density away from 0. A visual inspection suffices for modellers to identify the misspecified summary and offers guidance for adjusting the model. Further, we observed that in the well-specified case, the effect of the adjustment parameters on the coverage is minimal (see Appendix D).

We follow the misspecified moving average (MA) of order 1 example in Frazier & Drovandi (2021), where the assumed DGP is an MA(1) model, $y_{t}=w_{t}+\theta w_{t-1}$, $-1\leq\theta\leq 1$ and $w_{t}\overset{i.i.d.}{\sim}\mathcal{N}(0,1)$. However, the true DGP is a stochastic volatility model of the form:

where $0<\kappa,\sigma_{v}<1$, and $u_{t},v_{t}\overset{i.i.d.}{\sim}\mathcal{N}(0,1)$. We generate the observed data using the parameters $\omega=-0.76$, $\kappa=0.90$ and $\sigma_{v}=0.36$. The data is summarised using the autocovariance function, $\zeta_{j}(\bm{x})=\frac{1}{T}\sum_{i=j}^{T}x_{i}x_{i-j-1}$, where $T$ is the number of observations and $j\in\{0,1\}$ is the lag. We use the prior $\theta\sim\mathcal{U}(-1,1)$ and set $T=100$.

It can be shown that for the assumed DGP that $b(\theta)=(1+\theta^{2},\theta)^{\top}$. Under the true DGP, $\bm{b_{0}}=(\exp(\frac{\omega}{1-\kappa}+\frac{\sigma_{v}^{2}}{2(1-\kappa^{2})}),0)^{\top}\approx(0.0007,0)^{\top}$. As evidently $\inf_{\theta\in\Theta}||b(\theta)-\bm{b_{0}}||>0$, the model is misspecified. We also have a unique pseudo-true value with $||b(\theta)-\bm{b_{0}}||$ minimised at $\theta=0$. The desired behaviour for our robust algorithm is to detect incompatibility in the first summary statistic and centre the posterior around this pseudo-true value. As the first element of $b_{0}$ goes from $1\rightarrow 0$, $||b(\theta)-\bm{b_{0}}||$ increases, and the impact of model misspecification becomes more pronounced. We consider a specific run where we observe $S(\bm{y})=(0.013,0.006)^{\top}$.

Figure 3 shows that RSNL both detects the incompatible sample variance statistic and ensures that the approximate posterior concentrates onto the parameter value that favours matching of the compatible statistic, i.e. $\theta=0$. RNPE has a similar concentration around the pseudo-true value. SNL, however, is biased and has less support for the pseudo-true value.

As expected, $\gamma_{1}$ (corresponding to the incompatible statistic) has significant posterior density away from 0 as seen in Figure 3. Also, the posterior for $\gamma_{2}$ (corresponding to the compatible statistic) closely resembles the prior. The computational price of making inferences robust for the misspecified MA(1) model is minimal, with RSNL taking around 20 minutes to run and SNL taking around 10 minutes.

We may be concerned with the under-confidence of RSNL (and RNPE) on the misspecified MA(1) example as illustrated in the coverage plot in Figure 2. But from the posterior plots in Figure 3, we see that, for the misspecified MA(1) benchmark example, SNL is not only over-confident but it is over-confident for a point in the parameter space where the actual summaries and observed summaries are very different. In contrast, RSNL is under-confident, in that its posteriors are inflated relative to the standard level of frequentist coverage, however, it is under-confident for the right point: the values over which the RSNL posterior are centred deliver simulated summaries that are as close as possible to the observed summaries in the Euclidean norm. We also note that, in general when models are misspecified, Bayesian inference does not deliver valid frequentist coverage (Kleijn & van der Vaart, 2012).

The simple likelihood complex posterior (SLCP) model devised in Papamakarios et al. (2019) is a popular example in the SBI literature. The assumed DGP is a bivariate normal distribution with the mean vector, $\bm{\mu_{\theta}}=(\theta_{1},\theta_{2})^{\top}$, and covariance matrix:

where $s_{1}=\theta_{3}^{2}$, $s_{2}=\theta_{4}^{2}$ and $\rho=\tanh({\theta_{5}})$. This results in a nonlinear mapping from $\bm{\theta}=(\theta_{1},\theta_{2},\theta_{3},\theta_{4},\theta_{5})\in\mathbb{R}^{5}\rightarrow\bm{y_{j}}\in\mathbb{R}^{2}$, for $j=1,\ldots,5$. The posterior is “complex", having multiple modes due to squaring and vertical cutoffs from the uniform prior that we define in more detail later. Hence, the likelihood is expected to be easier to emulate than the posterior, making it suitable for an SNL approach. Four draws are generated from this bivariate distribution giving the likelihood, $g(\bm{y}\mid\bm{\theta})=\prod_{j=1}^{4}\mathcal{N}(\bm{y_{j}};\bm{\mu_{\theta}},\bm{\Sigma_{\theta}})$ for $\bm{y}=(\bm{y_{1}},\bm{y_{2}},\bm{y_{3}},\bm{y_{4}})$. No summarisation is done, and the observed data is used in place of the summary statistic. We generate the observed data at parameter values, $\bm{\theta}=(0.7,-2.9,-1.0,-0.9,0.6)^{\top}$, and place an independent $\mathcal{U}(-3,3)$ prior on each component of $\bm{\theta}$.

To impose misspecification on this illustrative example, we draw a contaminated 5-th observation, $\bm{y_{5}}$ and use the observed data $\bm{y}=(\bm{y_{1}},\bm{y_{2}},\bm{y_{3}},\bm{y_{4}},\bm{y_{5}})$. Contamination is introduced by applying the (stochastic) misspecification transform considered in Cannon et al. (2022), $\bm{y_{5}}=\bm{x_{5}}+100\bm{z_{5}}$, where $\bm{x_{5}}\sim\mathcal{N}(\mu_{\bm{\theta}},\Sigma_{\bm{\theta}})$, and $\bm{z_{5}}\sim\mathcal{N}((0,0)^{\top},100\mathbb{I}_{2})$. The assumed DGP is incompatible with this contaminated observation, and ideally, the approximate posterior would ignore the influence of this observation. Due to the stochastic transform, there is a small chance that the contaminated draw is compatible with the assumed DGP. However, the results presented here consider a specific example, where the observed contaminated draw is $\bm{y_{5}}=(23.41,-178.90)^{\top}$, which is very unlikely under the assumed DGP.

We thus want our inference to only use information from the four draws from the true DGP. The aim is to closely resemble the SNL posterior, where the observed data is the four non-contaminated draws. Figure 4 shows the estimated posterior densities for SNL (for both compatible and incompatible summaries) and RSNL for the contaminated SLCP example. When including the contaminated 5-th draw, SNL produces a nonsensical posterior with little useful information. Similarly, RNPE cannot correct the contaminated draw and does not give useful inference. Conversely, the RSNL posterior has reasonable density around the true parameters and has identified the separate modes.

The first eight compatible statistics are shown in Figure 6. The prior and posteriors reasonably match each other. In contrast, the observed data from the contaminated draw is recognised as being incompatible and has significant density away from 0, as evident in Figure 6. Again, there is no significant computational burden induced to estimate the adjustment parameters, with a total computational time of around 6 hours to run RSNL and around 4 hours for SNL.

Figure 2 shows that while RSNL is reasonably well-calibrated, RNPE does not give reliable inference. This is because RNPE could not correct for the high model misspecification of the contaminated draw. Additionally, approaches that emulate the likelihood rather than the posterior typically give better inference for this example (Lueckmann et al., 2021).

We follow the misspecified susceptible-infected-recovered (SIR) model described in Ward et al. (2022). SIR models are a simple compartmental model for the spread of infectious diseases. The standard SIR model has an infection rate, $\beta$, and recovery rate, $\eta$. The population consists of three states: susceptible ($S$), infected ($I$) and recovered ($R$). The standard SIR model is defined by:

The assumed DGP is a stochastic extension of the deterministic SIR model with a time-varying infection rate, $\tilde{\beta}_{t}$. This allows the SIR model to better capture heterogeneous disease outcomes, virus mutations and mitigation strategies (Spannaus et al., 2022). In addition to the equations in 7, the assumed SIR model is also parameterised by a time-varying effective reproduction number, $R_{et}=\frac{\tilde{\beta}_{t}}{\eta}$,

where $\upsilon$ is the reversion to $R_{et}$, $\sigma$ is the volatility and $W_{t}$ is Brownian motion. We set $\upsilon=0.5$ and $\sigma=0.5$ as in Ward et al. (2022), and use priors, $\eta\sim\mathcal{U}(0,0.5)$ and $\beta\mid\eta\sim\mathcal{U}(\eta,0.5)$. We use $\beta=0.15$ and $\eta=0.1$ for all observed data. The assumed DGP is run for 365 days, and only the daily number of infected individuals is considered. The initial number of infected individuals is 100. The infected counts are scaled by 100,000 to represent a larger population. A visualisation of the observed DGP is shown in Appendix D.

The true DGP has a reporting lag in recorded infections. Weekend days have the number of recorded infections reduced by 5%, being recorded on Monday, which sees an increase of 10%. Six summary statistics are considered: mean, median, max, max day (day of max infections), half day (day when half of the total infections was reached), and the autocorrelation of lag 1.

The model is misspecified, as the assumed SIR model cannot replicate the observed autocorrelation summary. We want our robust algorithm to detect misspecification in the autocorrelation summary and deliver useful inferences. We consider a specific example where the observed autocorrelation is 0.9957. Under the assumed DGP, the simulated autocorrelation is tightly centred around 0.9997. Despite the seemingly minor difference between the observed and simulated summaries, the observed summary lies far in the tails post standardisation.

As shown in Figure 7, RSNL produces useful inference with high density around the true parameters. Inspecting the adjustment parameters in Figure 8, we can see that the misspecification has been detected for the autocorrelation summary statistic and adjusted accordingly. Consequently, the modeller can further refine the simulation model to better capture the observed autocorrelation. This refinement could be achieved by recognising that the assumed DGP fails to account for the reporting lag evident in the observed data. RNPE behaves similarly to RSNL and concentrates around the true data-generating parameters. SNL, however, focuses overconfidently in an inconsequential region of the parameter space.

As the main computational cost in this example is to run the SIR simulations, there is negligible impact from the addition of adjustment parameters. In the example considered in Figures 7 and 8, SNL and RSNL both took approximately 24 hours.

We consider here the animal movement model by Marchand et al. (2017) to simulate the dispersal of Fowler’s toads (Anaxyrus fowleri). This is an individual-based model that encapsulates two main behaviours that have been observed in amphibians: high site fidelity and a small possibility of long-distance movement. The assumed behaviour is for each toad to act independently, stay at a refuge site during the day and move to forage at night. After foraging, the toad either stays in its current location or returns to a previous refuge site. A toad returning to a previous refuge site occurs with a constant probability $p_{0}$. We concentrate on “model 2” in Marchand et al. (2017), which models the behaviour of returning to its nearest previous refuge. This specific model was chosen as there is evidence of model misspecification, allowing us to assess our method on a misspecified example with real data.

The dispersal distance is modelled using a Lévy alpha-stable distribution, parameterised by a stability factor $\alpha_{\text{toad}}$ and scale factor $\delta$. This distribution was chosen for its heavy tails, which allow for occasional long-distance movement while still being symmetric around zero. Although the Lévy alpha-stable distribution lacks a closed form, it is straightforward to simulate. Thus the model is governed by three parameters: $\bm{\theta}=[\alpha_{\text{toad}},\delta,p_{0}]$. We assume the following uniform prior distributions for the model parameters: $\alpha_{\text{toad}}\sim\mathcal{U}(1,2),\delta\sim\mathcal{U}(20,70),$ and $p_{0}\sim\mathcal{U}(0.4,0.9)$.

The Marchand et al. (2017) GPS data was collected from 66 individual toads, with the daytime location (i.e. while resting refuge) being recorded. The number of recorded days varied across toads, with a maximum of 63 days. The two-dimensional GPS data is converted to a one-dimensional movement component, resulting in an observed $(63\times 66)$ matrix. The observed matrix was summarised using four sets of displacement vectors with time lags of 1, 2, 4 and 8 days. For each lag, the number of absolute displacements less than 10m, the median of the absolute displacements greater than 10m, and the log difference of the $0,0.1,\ldots,1$ of the absolute displacements greater than 10m are calculated, resulting in a total of 48 summary statistics (12 for each time lag).

In addition to SNL and RNPE, we evaluate our method against RBSL, as we are assessing performance on the real observed data, and the ground truth is unavailable for direct comparison. We examined plots for RSNL, SNL, RNPE and RBSL using real and simulated data from the toad example. In Figure 9, the estimated model parameter posteriors are conditioned on real observed summary statistics, with RBSL as the baseline for comparison with RSNL. The marginal plots reveal that RSNL closely resembles RBSL, while SNL differs, showing minimal density around the RSNL and RBSL maximum a posteriori estimates. RNPE also gives similar marginal posteriors to RSNL and RBSL. The 95% credible intervals for each parameter are displayed in Table 1. However, when considering simulated data (see Appendix D), SNL yields similar inferences to the robust methods. This suggests that SNL’s differing results in the real data scenario arise from incompatible summary statistics.

The most incompatible summary statistics were identified as the number of returns with a lag of 1 and the first quantile differences for lags 4 and 8. These were determined via MCMC output and visual inspection. We depict the posteriors for the adjustment parameters corresponding to these incompatible summaries in Figure 10, alongside the first three posteriors for compatible summaries, which are expected to closely match their priors. This example illustrates the advantage of carefully selecting summary statistics that hold intrinsic meaning for domain experts. For example, the insight that the current model cannot adequately capture the low number of observed toad returns—particularly while also fitting other movement behaviours—has direct meaning to the domain expert modeller, enabling them to refine the model accordingly.

In Figure 11, we present distributions for the log distance travelled by non-returning toads derived from the (pre-summarised) observation matrix. The simulated distances closely align with, and are tightly distributed around, the observed distances. Differences between the simulated and observed log distances appear most noticeable at lags 4 and 8 for shorter distances, aligning with the identified incompatible summaries. Figure 12 displays boxplots for the number of returns across the four lags. As confirmed by the first incompatible summary, the model has difficulties replicating the observed number of toad returns while accurately capturing other summary statistics. Overall, the posterior predictive distributions largely agree with observed data, with discrepancies coinciding with the incompatible summaries identified.

For the toad movement model, as the focus is on the performance of the methods on real data (with no ground-truth parameter), we instead consider the posterior predictive distribution. To compare predictive performance across methods, we computed the MMD between the observed summary statistic and samples from the posterior predictive, as presented in Table 1. Implementation details for the MMD can be found in Appendix D. We highlight that RSNL has a lower discrepancy than SNL and similar results to RNPE. While RBSL is the best-performing method in this example, it requires orders of magnitude more model simulations. We also note that the most incompatible summary statistics identified via the adjustment parameters (see results in Appendix D) agree with those found in Frazier & Drovandi (2021).