Simulation-Based Inference of Bayesian Hierarchical Models While Checking for Model Misspecification ^†

Abstract

This paper presents recent methodological advances for performing simulation-based inference (SBI) of a general class of Bayesian hierarchical models (BHMs) while checking for model misspecification. Our approach is based on a two-step framework. First, the latent function that appears as a second layer of the BHM is inferred and used to diagnose possible model misspecification. Second, target parameters of the trusted model are inferred via SBI. Simulations used in the first step are recycled for score compression, which is necessary for the second step. As a proof of concept, we apply our framework to a prey–predator model built upon the Lotka–Volterra equations and involving complex observational processes.

1. Introduction

Model misspecification is a long-standing problem for Bayesian inference: when the model differs from the actual data-generating process, posteriors tend to be biased and/or overly concentrated. In this paper, we are interested the problem of model misspecification for a particular, but common, class of Bayesian hierarchical models (BHMs): those that involve a latent function, such as the primordial power spectrum in cosmology (e.g., [1]) or the population model in genetics (e.g., [2]).

Simulation-based inference (SBI) only provides the posterior of top-level target parameters and marginalizes over all other latent variables of the BHM. Alone, it is therefore unable to diagnose whether the model is misspecified. Key insights regarding the issue of model misspecification can usually be obtained from the posterior distribution of the latent function, as there often exists an independent theoretical understanding of its values. An approximate posterior for the latent function (a much higher-dimensional quantity than the target vector of parameters) can be obtained using selfi (simulator expansion for likelihood-free inference, [1]), an approach based on the likelihood of an alternative parametric model, constructed by linearizing model predictions around an expansion point.

This paper presents a framework that combines selfi and SBI while recycling the necessary simulations. The simulator is first linearized to obtain the selfi posterior of the latent function. Next, the same simulations are used for data compression to the score function (the gradient of the log-likelihood with respect to the parameters), and the final SBI posterior of target parameters is obtained.

2. Method

2.1. Bayesian Hierarchical Models with a Latent Function

In this paper, we assume a given BHM consisting of the following variables: $\mathsf{\omega }\in {\mathbb{R}}^N$ (vector of N target parameters), $\mathsf{\theta }\in {\mathbb{R}}^S$ (vector containing the values of the latent function $\theta$ at S support points), $\mathsf{\Phi }\in {\mathbb{R}}^P$ (data vector of P components), and $\tilde{\mathsf{\omega }}\in {\mathbb{R}}^N$ (compressed data vector of size N). We typically expect $N\sim \mathcal{O}(5-10)$ target parameters, $S\sim \mathcal{O}({10}^2-{10}^3)$ support points; P can be any number and as large as $\mathcal{O}\left({10}^7\right)$ for complex data models. We further assume that $\mathsf{\omega }$ and $\mathsf{\theta }$ are linked by a deterministic function T, usually theoretically well-understood and numerically cheap. Therefore, the expensive and potentially misspecified part of the BHM is the probabilistic simulator linking the latent function $\mathsf{\theta }$ to the data $\mathsf{\Phi }$, $P\left(\mathsf{\Phi }\right|\mathsf{\theta })$. The deterministic compression step C linking $\mathsf{\Phi }$ to $\tilde{\mathsf{\omega }}$ is discussed in Section 2.4.

2.2. Latent Function Inference with SELFI

The first part of the framework proposed in this paper is to infer the latent function $\mathsf{\theta }$ conditional on observed data ${\mathsf{\Phi }}_{\mathrm{O}}$. This is an inference problem in high dimension (S, the number of support points for the latent function $\theta$), which means that usual SBI frameworks, allowing a general exploration of parameter space, will fail and that stronger assumptions are required. selfi [1] relies upon the simplification of the inference problem around an expansion point ${\mathsf{\theta }}_0$.

The first assumption is a Taylor expansion (linearization) of the mean data model around ${\mathsf{\theta }}_0$. Namely, if ${\widehat{\mathsf{\Phi }}}_{\mathsf{\theta }}\equiv \mathrm{E}\left[{\mathsf{\Phi }}_{\mathsf{\theta }}\right]$ is the expectation value of ${\mathsf{\Phi }}_{\mathsf{\theta }}$, where ${\mathsf{\Phi }}_{\mathsf{\theta }}$ are simulations of $\mathsf{\Phi }$ given $\mathsf{\theta }$ (i.e., ${\mathsf{\Phi }}_{\mathsf{\theta }}⤺P\left(\mathsf{\Phi }\right|\mathsf{\theta })$), we assume that

$${\widehat{\mathsf{\Phi }}}_{\mathsf{\theta }}\approx {\mathbf{f}}_0+\nabla {\mathbf{f}}_0·(\mathsf{\theta }-{\mathsf{\theta }}_0)\equiv \mathbf{f}\left(\mathsf{\theta }\right),$$

(1)

where ${\mathbf{f}}_0\equiv {\widehat{\mathsf{\Phi }}}_{{\mathsf{\theta }}_0}$ is the mean data model at the expansion point ${\mathsf{\theta }}_0$, and $\nabla {\mathbf{f}}_0$ is the gradient of ${\mathbf{f}}_0$ at the expansion point (for simplification, we note $\nabla {\mathbf{f}}_0={\nabla }_{\mathsf{\theta }}{\mathbf{f}}_0$, where the gradient is taken with respect to $\mathsf{\theta }$). The second assumption is that the (true) implicit likelihood of the problem is replaced by a Gaussian effective likelihood: $P\left({\mathsf{\Phi }}_{\mathrm{O}}\right|\mathsf{\theta })\equiv exp\left[{\widehat{\ell }}_{\mathsf{\theta }}\left(\mathsf{\theta }\right)\right]$ with

$$-2{\widehat{\ell }}_{\mathsf{\theta }}\left(\mathsf{\theta }\right)\approx log\left|2\pi {\mathbf{C}}_0\right|+{\left[{\mathsf{\Phi }}_{\mathrm{O}}-\mathbf{f}\left(\mathsf{\theta }\right)\right]}^{⊺}{\mathbf{C}}_0^{-1}\left[{\mathsf{\Phi }}_{\mathrm{O}}-\mathbf{f}\left(\mathsf{\theta }\right)\right],$$

(2)

where ${\mathbf{C}}_0$ is the data covariance matrix at the expansion point ${\mathsf{\theta }}_0$.

The selfi framework is fully characterized by ${\mathbf{f}}_0$, ${\mathbf{C}}_0$, and $\nabla {\mathbf{f}}_0$, which, if unknown, can be evaluated through forward simulations only. The numerical computation requires $N_0$ simulations at the expansion point (to evaluate the empirical mean ${\mathbf{f}}_0$ and empirical covariance matrix ${\mathbf{C}}_0$), and $N_s$ simulations in each direction of parameter space (to evaluate the empirical gradient $\nabla {\mathbf{f}}_0$ via first-order forward finite differences). The total is $N_0+N_s\times S$ simulations; $N_0$ and $N_s$ should be of the order of the dimensionality of the data space P, giving a total cost of $\mathcal{O}(\gtrsim P(S+1\left)\right)$ model evaluations.

To fully characterize the Bayesian problem, one requires a prior on $\mathsf{\theta }$, $P\left(\mathsf{\theta }\right)$. Any prior can be used if one is ready to use numerical techniques to explore the posterior (such as standard Markov Chain Monte Carlo), using the linearized data model and Gaussian effective likelihood. However, a remarkable analytic result with selfi is that, if the prior is Gaussian with a mean equal to the expansion point ${\mathsf{\theta }}_0$, i.e.,

$$-2logP\left(\mathsf{\theta }\right)\equiv log\left|2\pi \mathbf{S}\right|+{(\mathsf{\theta }-{\mathsf{\theta }}_0)}^{⊺}{\mathbf{S}}^{-1}(\mathsf{\theta }-{\mathsf{\theta }}_0),$$

(3)

then the effective posterior is also Gaussian:

$$-2logP\left(\mathsf{\theta }\right|{\mathsf{\Phi }}_{\mathrm{O}})\approx log\left|2\pi \mathsf{\Gamma }\right|+{(\mathsf{\theta }-\mathsf{\gamma })}^{⊺}{\mathsf{\Gamma }}^{-1}(\mathsf{\theta }-\mathsf{\gamma }).$$

(4)

The posterior mean and covariance matrix are given by

$$\begin{array}{ccc}\hfill \mathsf{\gamma }& \equiv & {\mathsf{\theta }}_0+\mathsf{\Gamma }{(\nabla {\mathbf{f}}_0)}^{⊺}{\mathbf{C}}_0^{-1}({\mathsf{\Phi }}_{\mathrm{O}}-{\mathbf{f}}_0),\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill \mathsf{\Gamma }& \equiv & {\left[{(\nabla {\mathbf{f}}_0)}^{⊺}{\mathbf{C}}_0^{-1}\nabla {\mathbf{f}}_0+{\mathbf{S}}^{-1}\right]}^{-1}\hfill \end{array}$$

(6)

(see [1] Appendix B, for a derivation). They are fully characterized by the expansion variables ${\mathsf{\theta }}_0$, ${\mathbf{f}}_0$, ${\mathbf{C}}_0$, and $\nabla {\mathbf{f}}_0$, as well as the prior covariance matrix $\mathbf{S}$.

2.3. Check for Model Misspecification

The selfi posterior can be used as a check for model misspecification. Visually checking the reconstructed $\mathsf{\gamma }$ and $\mathsf{\Gamma }$ can yield interesting insights, especially if the latent function has some properties (such as an expected shape, periodicity, etc.) to which the data model may be sensitive if misspecified (see Section 4.2).

If a quantitative check for model misspecification is desired, we propose using the Mahalanobis distance between the reconstruction $\mathsf{\gamma }$ and the prior distribution $P\left(\mathsf{\theta }\right)$, defined formally by

$$d_{\mathrm{M}}(\mathsf{\theta },{\mathsf{\theta }}_0|\mathbf{S})\equiv \sqrt{{\left(\mathsf{\theta }-{\mathsf{\theta }}_0\right)}^{⊺}{\mathbf{S}}^{-1}(\mathsf{\theta }-{\mathsf{\theta }}_0)}.$$

(7)

The value of $d_{\mathrm{M}}(\mathsf{\gamma },{\mathsf{\theta }}_0|\mathbf{S})$ for the selfi posterior mean $\mathsf{\gamma }$ can be compared to an ensemble of values of $d_{\mathrm{M}}({\mathsf{\theta }}_{\mathsf{\omega }},{\mathsf{\theta }}_0|\mathbf{S})$ for simulated latent functions ${\mathsf{\theta }}_{\mathsf{\omega }}=T\left(\mathsf{\omega }\right)$, where samples $\mathsf{\omega }$ are drawn from the prior $P\left(\mathsf{\omega }\right)$.

2.4. Score Compression and Simulation-Based Inference

Having checked the BHM for model misspecification, we now address the second part of the framework, aimed at inferring top-level parameters $\mathsf{\omega }$ given observations. SBI is known to be difficult when the dimensionality of the data space P is high. For this reason, data compression is usually necessary. Data compression can be thought of as an additional layer at the bottom of the BHM, made of a deterministic function C acting on $\mathsf{\Phi }$. In practical scenarios, data compression shall preserve as much information about $\mathsf{\omega }$ as possible, meaning that compressed summaries $C\left(\mathsf{\Phi }\right)$ shall be as close as possible to sufficient summary statistics of $\mathsf{\Phi }$, i.e., $P\left(\mathsf{\omega }\right|C\left(\mathsf{\Phi }\right))=P(\mathsf{\omega }\left|\mathsf{\Phi }\right)$.

Here, we propose to use score compression [3]. We make the assumption (for compression only, not for later inference) that $P\left(\mathsf{\Phi }\right|\mathsf{\omega })$ is Gaussian distributed: $P\left({\mathsf{\Phi }}_{\mathrm{O}}\right|\mathsf{\omega })\equiv exp\left[{\widehat{\ell }}_{\mathsf{\omega }}\left(\mathsf{\omega }\right)\right]$ where ${\widehat{\ell }}_{\mathsf{\omega }}\left(\mathsf{\omega }\right)={\widehat{\ell }}_{\mathsf{\theta }}\left(T\left(\mathsf{\omega }\right)\right)$ (see Equation (2)). The score function ${\nabla }_{\mathsf{\omega }}{\widehat{\ell }}_{\mathsf{\omega }0}$ is the gradient of this log-likelihood with respect to the parameters $\mathsf{\omega }$ at a fiducial point ${\mathsf{\omega }}_0$ in parameter space. Using as fiducial point the values that generate the selfi expansion point (i.e., ${\mathsf{\omega }}_0$ such that ${\mathsf{\theta }}_0=T\left({\mathsf{\omega }}_0\right)$), a quasi maximum-likelihood estimator for the parameters is ${\tilde{\mathsf{\omega }}}_{\mathrm{O}}\equiv {\mathsf{\omega }}_0+{\mathbf{F}}_0^{-1}{\nabla }_{\mathsf{\omega }}{\widehat{\ell }}_{\mathsf{\omega }0}$, where the Fisher matrix ${\mathbf{F}}_0$ and the gradient of the log-likelihood are evaluated at ${\mathsf{\omega }}_0$. Compression of ${\mathsf{\Phi }}_{\mathrm{O}}$ to ${\tilde{\mathsf{\omega }}}_{\mathrm{O}}$ yields N compressed statistics that are optimal in the sense that they preserve the Fisher information content of the data [3].

In our case, the covariance matrix ${\mathbf{C}}_0$ is assumed not to depend on parameters (${\nabla }_{\mathsf{\omega }}{\mathbf{C}}_0=0$), and the expression for $C\left(\mathsf{\Phi }\right)$ is therefore

$$C\left(\mathsf{\Phi }\right)=\tilde{\mathsf{\omega }}\equiv {\mathsf{\omega }}_0+{\mathbf{F}}_0^{-1}\left[{\left({\nabla }_{\mathsf{\omega }}{\mathbf{f}}_0\right)}^{⊺}{\mathbf{C}}_0^{-1}(\mathsf{\Phi }-{\mathbf{f}}_0)\right].$$

(8)

The Fisher matrix of the problem further takes a simple form:

$${\mathbf{F}}_0\equiv -\mathrm{E}\left[{\nabla }_{\mathsf{\omega }}{\nabla }_{\mathsf{\omega }}{\widehat{\ell }}_{\mathsf{\omega }0}\left(\mathsf{\omega }\right)\right]={\left({\nabla }_{\mathsf{\omega }}{\mathbf{f}}_0\right)}^{⊺}{\mathbf{C}}_0^{-1}{\nabla }_{\mathsf{\omega }}{\mathbf{f}}_0.$$

(9)

We therefore need to evaluate

$${\nabla }_{\mathsf{\omega }}{\mathbf{f}}_0=\nabla {\mathbf{f}}_0·{\left.\frac{\partial T\left(\mathsf{\omega }\right)}{\partial \mathsf{\omega }}\right|}_{\mathsf{\omega }={\omega }_0}.$$

(10)

Importantly, in Equations (8)–(10), ${\mathbf{C}}_0$ and $\nabla {\mathbf{f}}_0$ have already been computed for latent function inference with selfi. The only missing quantity is the second matrix in the right-hand side of Equation (10), that is, ${\nabla }_{\mathsf{\omega }}{T}_0$, the gradient of T evaluated at ${\mathsf{\omega }}_0$. If unknown, its computation (e.g., via finite differences) does not require any more simulation of $\mathsf{\Phi }$. It is usually easy, as there are only N directions in parameter space and T is the numerically cheap part of the BHM. We note that, because we have to calculate ${\mathbf{F}}_0$, we can easily get the Fisher–Rao distance between any simulated summaries $\tilde{\mathsf{\omega }}$ and the observed summaries ${\tilde{\mathsf{\omega }}}_{\mathrm{O}}$,

$$d_{\mathrm{FR}}(\tilde{\mathsf{\omega }},{\tilde{\mathsf{\omega }}}_{\mathrm{O}})\equiv \sqrt{{\left(\tilde{\mathsf{\omega }}-{\tilde{\mathsf{\omega }}}_{\mathrm{O}}\right)}^{⊺}{\mathbf{F}}_0(\tilde{\mathsf{\omega }}-{\tilde{\mathsf{\omega }}}_{\mathrm{O}})},$$

(11)

which can be used by any non-parametric SBI method.

We specify a prior $P\left(\mathsf{\omega }\right)$ (typically peaking at or centered on ${\mathsf{\omega }}_0$, for consistency with the assumptions made for data compression). Having defined C, we now have a full BHM that maps $\mathsf{\omega }$ (of dimension N) to compressed summaries $\tilde{\mathsf{\omega }}$ (of size N) and has been checked for model misspecification for the part linking $\mathsf{\theta }$ to $\mathsf{\Phi }$. We can then proceed with SBI via usual techniques. These can include likelihood-free rejection sampling, but also more sophisticated techniques such as delfi (e.g., [4,5]) or bolfi (e.g., [6,7,8]).

3. Lotka–Volterra BHM

3.1. Lotka–Volterra Solver

The Lotka–Volterra equations describe the dynamics of an ecological system in which two species interact, as a pair of first-order non-linear differential equations:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}& =& \alpha x-\beta xy,\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\mathrm{d}y}{\mathrm{d}t}}& =& \delta xy-\gamma y.\hfill \end{array}$$

(13)

where $x\left(t\right)$ is the number of prey at time t, and $y\left(t\right)$ is the number of predators at time t. The model is characterized by $\mathsf{\omega }=(\alpha ,\beta ,\gamma ,\delta )$, a vector of four real parameters describing the interaction of the two species.

The initial conditions of the problem $\left\{x\left(0\right),y\left(0\right)\right\}=\left\{x_0,y_0\right\}$ are assumed to be exactly known. Throughout the paper, timestepping and number of timesteps are fixed: $t_i=i\mathsf{\Delta }t$ for $i\in 〚0,S/2〛$.

The expression T is an algorithm that numerically solves the ordinary differential equations. For simplicity, we choose an explicit Euler method: for all $i\in 〚0,S/2-1〛$,

$$\begin{array}{ccc}\hfill x\left(t_{i+1}\right)& =& x\left(t_i\right)\times \left[1+\alpha -\beta y\left(t_i\right)\right]\times {\mathsf{\Delta }}_t,\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill y\left(t_{i+1}\right)& =& y\left(t_i\right)\times \left[1+\delta x\left(t_i\right)-\gamma \right]\times {\mathsf{\Delta }}_t.\hfill \end{array}$$

(15)

The latent function $\theta \left(t\right)$ is a concatenation of $x\left(t\right)$ and $y\left(t\right)$ evaluated at the timesteps of the problem. The corresponding vector is $\mathsf{\theta }\equiv \left\{{\left\{x\left(t_i\right)\right\}}_{0\le i<S/2},{\left\{y\left(t_i\right)\right\}}_{0\le i<S/2}\right\}$ of size S.

3.2. Lotka–Volterra Observer

3.2.1. Full Data Model

To go from $\mathsf{\theta }$ to $\mathsf{\Phi }$, we assume a complex, probabilistic observational process of prey and predator populations, later referred to as “model A” and defined as follows.

Signal. The (unobserved) signal $s_z$ is a delayed and non-linearly perturbed observation of the true population function for species $z\in \left\{x,y\right\}$, modulated by some seasonal efficiency $e_z\left(t\right)$. Formally, $s_x\left(0\right)=x_0$, $s_y\left(0\right)=y_0$, and for $i\in 〚0,S/2-1〛$,

$$\begin{array}{ccc}\hfill s_x\left(t_{i+1}\right)& =& e_x\left(t_i\right)\left[x\left(t_i\right)-px\left(t_i\right)y\left(t_i\right)+qx{\left(t_i\right)}^2\right],\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill s_y\left(t_{i+1}\right)& =& e_y\left(t_i\right)\left[y\left(t_i\right)+px\left(t_i\right)y\left(t_i\right)-qy{\left(t_i\right)}^2\right].\hfill \end{array}$$

(17)

These equations involve two parameters: p accounts for hunts between $t_i$ and $t_{i+1}$ (temporarily making prey more likely to hide and predators more likely to be visible), and q accounts for the gregariousness of prey and independence of predators. The free functions $e_x\left(t\right)$ and $e_y\left(t\right)$, valued in $[0,1]$, describe how prey and predators are likely to be detectable at any time, accounting, for example, for seasonal variation (hibernation, migration).

Noise. The signal $s_z$ is subject to additive noise, giving a noisy signal $u_z\left(t\right)=s_z\left(t\right)+n_z^{\mathrm{D}}\left(t\right)+n_z^{\mathrm{O}}\left(t\right)$, where the noise has two components:

• Demographic Gaussian noise with zero mean and variance proportional to the true underlying population, i.e., $n_x^{\mathrm{D}}\left(t\right)⤺G\left[0,rx\left(t\right)\right]$ and $n_y^{\mathrm{D}}\left(t\right)⤺G\left[0,ry\left(t\right)\right]$. The parameter r gives the strength of demographic noise.
• Observational Gaussian noise that accounts for observer efficiency, coupling prey and predators such that

  $$\left(\begin{array}{c}{n}_x^{\mathrm{O}}\left(t\right)\\ n_y^{\mathrm{O}}\left(t\right)\end{array}\right)⤺G\left[\left(\begin{array}{c}0\\ 0\end{array}\right),s\left(\begin{array}{cc}y\left(t\right)& t\sqrt{x\left(t\right)y\left(t\right)}\\ t\sqrt{x\left(t\right)y\left(t\right)}& x\left(t\right)\end{array}\right)\right].$$

  (18)
  The parameter s gives the overall amplitude of observational noise, and the parameter t controls the strength of the non-diagonal component (it should be chosen such that the covariance matrix appearing in Equation (18) is positive semi-definite).

Censoring. Finally, observed data are a censored and thresholded version of the noisy signal: for each timestep $t_i$, ${\mathsf{\Phi }}_z\left(t_i\right)=m_z\left(t_i\right)\times min\left[u_z\left(t_i\right),M_z\right]$, where $M_z$ is the maximum number of prey or predators that can be detected by the observer, and $m_z$ is a mask (taking either the value 0 or 1). Masked data points are discarded. The data vector is $\mathsf{\Phi }=\left\{\left\{{\mathsf{\Phi }}_x\left(t_i\right)\right\},\left\{{\mathsf{\Phi }}_y\left(t_i\right)\right\}\right\}$. It contains $P\le S$ elements depending on the number of masked timesteps for each species z (formally, $P={\sum }_{i=0}^{S/2-1}\left({\mathsf{\delta }}_{\mathrm{K}}^{m_x\left(t_i\right),1}+{\mathsf{\delta }}_{\mathrm{K}}^{m_y\left(t_i\right),1}\right)$, where ${\mathsf{\delta }}_{\mathrm{K}}$ is a Kronecker delta symbol).

All of the free parameters (p, q, r, s, t, $M_x$, $M_y$) and free functions ($e_x\left(t\right)$, $e_y\left(t\right)$, $m_x\left(t\right)$, $m_y\left(t\right)$) appearing in the Lotka–Volterra observer data model described in this section are assumed known and fixed throughout the paper. Parameters used are $x_0=10$, $y_0=5$, $p=0.05$, $q=0.01$, $r=0.15$, $s=0.05$, $t=0.2$.

3.2.2. Simplified Data Model

In this section, we introduce “model B”, a simplified (misspecified) data model linking $\mathsf{\theta }$ to $\mathsf{\Phi }$. Model B assumes that underlying functions are directly observed, i.e., $s_z\left(t\right)=z\left(t\right)$. It omits observational noise, such that $u_z\left(t\right)=s_z\left(t\right)+n_z^{\mathrm{D}}\left(t\right)$. In model B, parameters p, q, s, and t are not involved, and the value of r (strength of demographic noise) can be incorrect (we used $r=0.105$). Finally, model B fails to account for the thresholds: ${\mathsf{\Phi }}_z\left(t\right)=m_z\left(t\right)u_z\left(t\right)$.

4. Results

In this section, we apply the two-step inference method described in Section 2 to the Lotka–Volterra BHM introduced in Section 3. We generate mock data ${\mathsf{\Phi }}_{\mathrm{O}}$ from model A, using ground truth parameters ${\mathsf{\omega }}_{\mathrm{gt}}=({\alpha }_{\mathrm{gt}},{\beta }_{\mathrm{gt}},{\gamma }_{\mathrm{gt}},{\delta }_{\mathrm{gt}})=(0.55,0.2,0.2,0.05)$. We assume that ground truth parameters are known a priori with a precision of approximately $3\%$. Consistently, we choose a Gaussian prior $P\left(\mathsf{\omega }\right)$ with mean ${\mathsf{\omega }}_0=(0.5768,0.1963,0.1968,0.0484)$ and diagonal covariance matrix $\mathrm{diag}(0.{0173}^2,0.{0059}^2,0.{0059}^2,0.{0015}^2)$.

4.1. Inference of Population Functions with SELFI

We first seek to reconstruct the latent population functions $x\left(t\right)$ and $y\left(t\right)$, conditional on the data ${\mathsf{\Phi }}_{\mathrm{O}}$, using selfi. We choose as an expansion point the population functions simulated from the mean of the prior on $\mathsf{\omega }$, i.e., ${\mathsf{\theta }}_0=T\left({\mathsf{\omega }}_0\right)$. We use $N_0=150$ and $N_s=100$; the computational workload is therefore a fixed number of $10,150$ simulations for each model. It is known a priori and perfectly parallel.

We adopt a Gaussian prior $P\left(\mathsf{\theta }\right)$ and combine it with the effective likelihood to obtain the selfi effective posterior $P\left(\mathsf{\theta }\right|{\mathsf{\Phi }}_{\mathrm{O}})$. Figure 1 (left panels) shows the inferred population functions $\mathsf{\gamma }$ in comparison with the prior mean and expansion point ${\mathsf{\theta }}_0$ and the ground truth ${\mathsf{\theta }}_{\mathrm{gt}}$. The figure shows $2\sigma$ credible regions for the prior and the posterior (i.e., $2\sqrt{\mathrm{diag}\left(\mathbf{S}\right)}$ and $2\sqrt{\mathrm{diag}\left(\mathsf{\Gamma }\right)}$, respectively). The full posterior covariance matrix $\mathsf{\Gamma }$ for each model is shown in the rightmost column of Figure 1.

4.2. Check for Model Misspecification

The inferred population functions allow us to check for model misspecification. From Figure 1, it is clear that model B fails to produce a plausible reconstruction of population functions: model B breaks the (pseudo-)periodicity of the predator population function $y\left(t\right)$, which is a property required by the model. In the bottom left-hand panels, the red lines differ in shape from fiducial functions $T\left(\mathsf{\omega }\right)$ (grey lines), and the credible intervals exclude the expansion point. On the contrary, with model A, the reconstructed population functions are consistent with the expansion point. The inference is unbiased, as the ground truth typically lies within the $2\sigma$ credible region of the reconstruction.

As a quantitative check, we compute the Mahalanobis distance between $\mathsf{\gamma }$ and $P\left(\mathsf{\theta }\right)$ (Equation (7)) for each model. We find that $d_{\mathrm{M}}(\mathsf{\gamma },{\mathsf{\theta }}_0|\mathbf{S})$ is much smaller for model A than for model B ($5.35$ versus $12.54$). The numbers can be compared to the empirical mean among our set of fiducial populations functions, $〈d_{\mathrm{M}}(T\left({\mathsf{\omega }}_n\right),{\mathsf{\theta }}_0|\mathbf{S})〉=9.43$.

At this stage, we therefore consider that model B is excluded, and we proceed further with model A.

4.3. Score Compression

As T is numerically cheap, we get ${\nabla }_{\mathsf{\omega }}{T}_0$ via sixth-order central finite differences around ${\mathsf{\omega }}_0$, then obtain ${\nabla }_{\mathsf{\omega }}{\mathbf{f}}_0$ using Equation (10). This does not require any further evaluation of the data model $P\left(\mathsf{\Phi }\right|\mathsf{\theta })$, as $\nabla {\mathbf{f}}_0$ has already been computed.

Using Equations (8) and (9), we compress ${\mathsf{\Phi }}_{\mathrm{O}}$ and obtain ${\tilde{\mathsf{\omega }}}_{\mathrm{O}}=(0.7050,0.2287,0.1471,0.0415)$.

4.4. Inference of Parameters Using Likelihood-Free Rejection Sampling

As a last step, we infer top-level parameters $\mathsf{\omega }$ given compressed summaries ${\tilde{\mathsf{\omega }}}_{\mathrm{O}}$. As the problem studied in this paper is sufficiently simple, we rely on the simplest solution for SBI, namely likelihood-free rejection sampling (sometimes also known as approximate Bayesian computation, e.g., [9]). To do so, we use the Fisher–Rao distance between simulated $\tilde{\mathsf{\omega }}$ and observed ${\tilde{\mathsf{\omega }}}_{\mathrm{O}}$, which comes naturally from score compression (see Equation (11)), and we set a threshold $\epsilon =2$. We draw samples from the prior $P\left(\mathsf{\omega }\right)$, simulate $\tilde{\mathsf{\omega }}$, then accept $\mathsf{\omega }$ as a sample of $P\left(\mathsf{\omega }\right|{\tilde{\mathsf{\omega }}}_{\mathrm{O}})$ if $d_{\mathrm{FR}}(\tilde{\mathsf{\omega }},{\tilde{\mathsf{\omega }}}_{\mathrm{O}})<\epsilon$, and reject it otherwise.

In Figure 2, we find that the inference of top-level parameters is unbiased, with the ground truth ${\mathsf{\omega }}_{\mathrm{gt}}$ (dashed lines) lying within the $2\sigma$ credible region of the posterior. We observe that the data correctly drive some features that are not built into the prior, for instance, the degeneracy between $\alpha$ and $\gamma$, respectively, the reproduction rate of prey and the mortality rate of predators.

5. Conclusions

One of the biggest challenges in statistical data analysis is checking data models for misspecification, so as to obtain meaningful parameter inferences. In this work, we described a novel two-step simulation-based Bayesian approach, combining selfi and SBI, which can be used to tackle this issue for a large class of models. BHMs to which the approach can be applied involve a latent function depending on parameters and observed through a complex probabilistic process. They are ubiquitous, e.g., in astrophysics and ecology.

In this paper, we introduced a prey–predator model, consisting of a numerical solver of the Lotka–Volterra system of equations and of a complex observational process of population functions. As a proof of concept, we applied our technique to this model and to a simplified (misspecified) version of it. We demonstrated successful identification of the misspecified model and unbiased inference of the parameters of the correct model.

In conclusion, the method developed constitutes a computationally efficient and easily applicable framework to perform SBI of BHMs while checking for model misspecification. It allows one to infer the latent function as an intermediate product, then to perform score compression at no additional simulation cost. This study opens up a new avenue to increase the robustness and reliability of Bayesian data analysis using fully non-linear, simulator-based models.