Likelihood-Free Parameter Estimation with Neural Bayes Estimators

Abstract

Neural Bayes estimators are neural networks that approximate Bayes estimators. They are fast, likelihood-free, and amenable to rapid bootstrap-based uncertainty quantification. In this article, we aim to increase the awareness of statisticians to this relatively new inferential tool, and to facilitate its adoption by providing user-friendly open-source software. We also give attention to the ubiquitous problem of estimating parameters from replicated data, which we address using permutation-invariant neural networks. Through extensive simulation studies we demonstrate that neural Bayes estimators can be used to quickly estimate parameters in weakly identified and highly parameterized models with relative ease. We illustrate their applicability through an analysis of extreme sea-surface temperature in the Red Sea where, after training, we obtain parameter estimates and bootstrap-based confidence intervals from hundreds of spatial fields in a fraction of a second.

Keywords:

• Amortized inference
• Deep learning
• Exchangeable data
• Extreme-value model
• Permutation invariant
• Point estimation
• Spatial statistics

1 Introduction

The most popular methods for estimating parameters in parametric statistical models are those based on the likelihood function. However, it is not always possible to formulate the likelihood function (Diggle and Gratton 1984; Lintusaari et al. 2017) and, even when it is available, it may be computationally intractable. For example, popular models for spatial extremes are max-stable processes, for which the number of terms involved in the likelihood function grows more than exponentially fast with the number of observations (Padoan 2010; Huser et al. 2019). One common workaround is to replace the full likelihood with a composite likelihood (e.g., Cox and Reid 2004; Varin and Vidoni 2005; Varin, Reid, and Firth 2011), but this usually results in a loss of statistical efficiency (Huser and Davison 2013; Castruccio, Huser, and Genton 2016), and it is not always clear how one should construct the composite likelihood. The related Vecchia approximation (Vecchia 1988) has been applied successfully both in Gaussian (e.g., Stein, Chi, and Welty 2004) and max-stable (Huser, Stein, and Zhong 2022) settings, but this approximation still trades statistical efficiency for computational efficiency.

To bypass these challenges, several model-fitting methods have been developed that preclude the need to evaluate the likelihood function. The most popular of these so-called likelihood-free methods is approximate Bayesian computation (ABC; Beaumont, Zhang, and Balding 2002; Sisson, Fan, and Beaumont 2018). In its simplest form, ABC involves sampling parameters from the prior, simulating from the model, and retaining parameters as an approximate sample from the posterior if the simulated data are “similar” to the observed data, with similarity typically assessed by comparing low-dimensional summary statistics. ABC methods are sensitive to the choice of summary statistics, and they are notoriously difficult to calibrate: for example, the number of summary statistics and the tolerance used when comparing statistics affect both the computational efficiency and the statistical efficiency of ABC.

Recently, neural networks have emerged as a promising approach to likelihood-free inference. In this work, we focus on neural networks that map data to parameter point estimates; we refer to such neural networks as neural point estimators. Neural point estimators date back to at least Chon and Cohen (1997), but they have only been adopted widely in recent years, for example in applications involving models for stock returns (Creel?); population-genetics (Flagel, Brandvain, and Schrider 2018); time series (Rudi, Julie, and Lenzi 2021); spatial fields (Gerber and Nychka 2021; Banesh et al. 2021; Lenzi et al. 2023); spatio-temporal fields (Zammit-Mangion and Wikle 2020); and agent-based models (Gaskin, Pavliotis, and Girolami 2023). The computational bottleneck in neural point estimation is the training procedure, whereby the neural network “learns” a useful mapping between the sample space and the parameter space. Importantly though, this training cost is amortized: for a given statistical model and prior distribution over the parameters, a trained neural point estimator with sufficiently flexible architecture can be re-used for new datasets at almost no computational cost, provided that each data replicate has the same format as those used to train the point estimator (e.g., images of a pre-specified width and height). Uncertainty quantification of the estimates proceeds naturally through the bootstrap distribution, which is essentially available “for free” with neural point estimators since they are so fast to evaluate. As we shall show, neural point estimators can be trained to approximate Bayes estimators and, in this case, we refer to them as neural Bayes estimators.

Parameter estimation from replicated data is commonly required in statistical applications, and we therefore give this topic particular attention. Neural point estimation with replicated data is not straightforward; for example, Gerber and Nychka (2021) considered two approaches to handling replicated data, both with drawbacks. In their first approach, they trained a neural point estimator for a single spatial field, applied it to each field independently, and averaged the resulting estimates; we call this the “one-at-a-time” approach. This approach does not reduce the bias commonly seen in small-sample estimators as the sample size increases, and it is futile when the parameters are unidentifiable from a single replicate. In their second approach, they adapted the size of the neural estimator’s input layer to match the number of independent replicates. The resulting estimator is not invariant to the ordering of the replicates; it results in an explosion of neural-network parameters with increasing number of replicates; and it requires a different architecture for every possible sample size, which reduces its applicability and generality.

In this article, we first clarify the connection between neural point estimators and classical Bayes estimators, which is sometimes misconstrued or ignored in the literature. Second, we propose a novel way to perform neural Bayes estimation from independent replicates by leveraging permutation-invariant neural networks, constructed using the DeepSets framework (Zaheer et al. 2017). To the best of our knowledge, this article is the first to explore its use and the related practical considerations in a point estimation context. We show that these architectures lead to a substantial improvement in estimation accuracy when compared to those that do not account for replication appropriately. Third, we discuss important practicalities for designing and training neural Bayes estimators and, in particular, we describe a way to construct an estimator that is approximately Bayes for any sample size. For illustration, we estimate parameters in Gaussian and max-stable processes, as well as the highly-parameterized spatial conditional extremes model (Wadsworth and Tawn 2022). We use the latter model in the analysis of sea-surface temperature extremes in the Red Sea where, using a neural Bayes estimator, we obtain estimates and bootstrap confidence intervals from hundreds of spatial fields in a fraction of a second. A primary motivation of this work is to facilitate the adoption of neural point estimation by statisticians and, to this end, we accompany the article with user-friendly open-source software in the Julia (Bezanson et al. 2017) and R (R Core Team 2023) programming languages: our software package, NeuralEstimators, is available at https://github.com/msainsburydale/NeuralEstimators, and can be used in a wide range of applied settings.

The remainder of this article is organized as follows. In Section 2, we outline the theory underlying neural Bayes estimators and discuss their implementation. In Section 3, we conduct extensive simulation studies that clearly demonstrate the utility of neural Bayes estimators. In Section 4, we use the spatial conditional extremes model to analyze sea-surface temperature in the Red Sea. In Section 5, we conclude and outline avenues for future research. A supplement is also available that contains more details and figures.

2 Methodology

In Section 2.1, we introduce neural Bayes estimators. In Section 2.2, we discuss the use of permutation-invariant neural networks for neural Bayes estimation from replicated data. In Section 2.3, we describe the general workflow for implementing neural Bayes estimators, and discuss some important practical considerations.

2.1 Neural Bayes Estimators

A parametric statistical model is a set of probability distributions on a sample space $\mathcal{S}$, where the probability distributions are parameterized via some p-dimensional parameter vector $\mathit{\theta }$ on a parameter space Θ (McCullagh 2002). Suppose that we have data, which we denote as Z, from one such distribution. Then, the goal of parameter point estimation is to come up with an estimate of the unknown $\mathit{\theta }$ from Z using an estimator, $$\widehat{\mathit{\theta }}:\mathcal{S}\to \Theta ,$$ which is a mapping from the sample space to the parameter space.

Estimators can be constructed within a decision-theoretic framework. Assume that the sample space is $\mathcal{S}={\mathbb{R}}^n$, and consider a nonnegative loss function, $L(\mathit{\theta },\widehat{\mathit{\theta }}(\mathit{Z}))$, which assesses an estimator $\widehat{\mathit{\theta }}(·)$ for a given $\mathit{\theta }$ and dataset $\mathit{Z}\sim f(\mathit{z}|\mathit{\theta })$, where $f(\mathit{z}|\mathit{\theta })$ is the probability density function of the data conditional on $\mathit{\theta }$. An estimator’s risk function is its loss averaged over all possible data realizations, (1) $$R(\mathit{\theta },\widehat{\mathit{\theta }}(·))\equiv {\int }_{\mathcal{S}}L(\mathit{\theta },\widehat{\mathit{\theta }}(\mathit{z}))f(\mathit{z}|\mathit{\theta })\mathrm{d}\mathit{z}.$$(1)

So-called Bayes estimators minimize a (weighted) average of (1) known as the Bayes risk, (2) $$r_{\Omega }(\widehat{\mathit{\theta }}(·))\equiv {\int }_{\Theta }R(\mathit{\theta },\widehat{\mathit{\theta }}(·))\mathrm{d}\Omega (\mathit{\theta }),$$(2) where $\Omega (·)$ is a prior measure for $\mathit{\theta }$. Note that in (1) and (2) we deviate slightly from classical notation and notate the estimator as $\widehat{\mathit{\theta }}(·)$, rather than just $\widehat{\mathit{\theta }}$, to stress that it is constructed as a function of the data.

Bayes estimators are theoretically attractive: for example, unique Bayes estimators are admissible and, under suitable regularity conditions and the squared-error loss, are consistent and asymptotically efficient (Lehmann and Casella 1998, chap. 5, thm. 2.4; chap. 6, thm. 8.3). Further, for a large class of prior distributions, every set of conditions that imply consistency of the maximum likelihood estimator also imply consistency of Bayes estimators (Strasser 1981). Unfortunately, however, Bayes estimators are typically unavailable in closed form for the complex models often encountered in practice. A way forward is to assume a flexible parametric model for $\widehat{\mathit{\theta }}(·)$, and to optimize the parameters within that model in order to approximate the Bayes estimator. Neural networks are ideal candidates, since they are universal function approximators (e.g., Hornik, Stinchcombe, and White 1989; Zhou 2018), and because they are fast to evaluate, usually involving only simple matrix-vector operations.

Let $\widehat{\mathit{\theta }}(\mathit{Z};\mathit{\gamma })$ denote a neural point estimator, that is, a neural network that returns a point estimate from data Z, where $\mathit{\gamma }$ contains the neural-network parameters (i.e., the so-called “weights” and “biases”). Bayes estimators may be approximated with $\widehat{\mathit{\theta }}(·;{\mathit{\gamma }}^{*})$ by solving the optimization problem, (3) $${\mathit{\gamma }}^{*}\equiv \underset{\mathit{\gamma }}{\text{arg} \min }{r}_{\Omega }(\widehat{\mathit{\theta }}(·;\mathit{\gamma })).$$(3)

Typically, $r_{\Omega }(·)$ defined in (2) cannot be directly evaluated, but it can be approximated using Monte Carlo methods. Specifically, given a set ϑ of K parameter vectors sampled from the prior $\Omega (·)$ and, for each $\mathit{\theta }\in \vartheta$, J samples from $f(\mathit{z}|\mathit{\theta })$ collected in the set ${\mathcal{Z}}_{\mathit{\theta }}$, (4) $$r_{\Omega }(\widehat{\mathit{\theta }}(·;\mathit{\gamma }))\approx \frac{1}{K}\sum _{\mathit{\theta }\in \vartheta }\frac{1}{J}\sum _{\mathit{Z}\in {\mathcal{Z}}_{\mathit{\theta }}}L(\mathit{\theta },\widehat{\mathit{\theta }}(\mathit{Z};\mathit{\gamma })).$$(4)

Note that (4) does not involve evaluation, or knowledge, of the likelihood function.

The surrogate objective function (4) can be straightforwardly minimized with respect to $\mathit{\gamma }$ using back-propagation and stochastic gradient descent with deep-learning software packages such as Flux (Innes 2018). For sufficiently expressive architectures, the point estimator targets a Bayes estimator with respect to $L(·,·)$ and $\Omega (·)$. We therefore call the fitted neural point estimator a neural Bayes estimator. Note that neural Bayes estimators, like all other neural networks, are function approximators that, at least in theory, can approximate continuous functions arbitrarily well. However, the discrepancy between $\widehat{\mathit{\theta }}(·;{\mathit{\gamma }}^{*})$ and the true Bayes estimator will depend on a number of factors, such as the neural-network architecture and the specific optimization procedure used to minimize (4). We provide more discussion on these practical considerations in Section 2.3.

Like Bayes estimators, neural Bayes estimators target a specific point summary of the posterior distribution. For instance, the 0–1, absolute-error, and squared-error loss functions lead to neural Bayes estimators that approximate the posterior mode, median, and mean, respectively. Further, posterior quantiles, which can be used to construct credible intervals, may be estimated by using the quantile loss function (e.g., Cressie 2023, eq. 7(7) $$\begin{array}{c}R(\mathit{\theta },\widehat{\mathit{\theta }}(·;\mathit{\gamma }))\equiv \sum _{m\in \mathcal{M}}\text{Pr}(M=m)\\ \hspace{1em}\times \left({\int }_{{\mathcal{S}}^m}L(\mathit{\theta },\widehat{\mathit{\theta }}({\mathit{z}}^{(m)};\mathit{\gamma }))f({\mathit{z}}^{(m)}|\mathit{\theta })\mathrm{d}{\mathit{z}}^{(m)}\right).\end{array}$$(7) ). This important link between neural point estimators and Bayes estimators is often overlooked in the literature.

2.2 Neural Bayes Estimators for Replicated Data

While neural Bayes estimators have been employed, sometimes inadvertently, in various applications, much less attention has been given to neural Bayes estimators for replicated data. In what follows, we denote a collection of m independent and identically distributed replicates (iid) ${\mathit{Z}}_1,\dots ,{\mathit{Z}}_m$ as ${\mathit{Z}}^{(m)}\equiv ({\mathit{Z}}_1^{\prime },\dots ,{\mathit{Z}}_m^{\prime })\prime$. In Section 2.2.1 we show how the so-called DeepSets framework is ideally placed to construct Bayes estimators for replicated data; in Section 2.2.2 we discuss approaches one can adopt to design neural point estimators that are approximately Bayes for any number of replicates; and in Section 2.2.3 we illustrate the neural Bayes estimator for replicated data on a relatively simple problem where the Bayes estimator is known in closed form.

2.2.1 The DeepSets Framework

Parameter estimation from replicated data is commonly required in statistical applications and, as discussed in Section 1, naïve approaches to constructing neural Bayes estimators for replicated data can lead to bias, an explosion in the number of neural-network parameters, and an inability to generalize to different sample sizes. We overcome these issues by leveraging an important property of Bayes estimators for replicated data, namely that they are invariant to permutations of the replicates. In Section S1 of the supplementary material, we give a formal proof of this property.

The permutation-invariance property of Bayes estimators for replicated data has largely been ignored to date in the context of neural point estimation. We propose enforcing our neural Bayes estimators to be permutation invariant by couching them in the DeepSets framework (Zaheer et al. 2017), a computationally convenient special case of Janossy pooling (Murphy et al. 2019). In this framework, the Bayes estimator is represented by (5) $$\begin{array}{c}\widehat{\mathit{\theta }}({\mathit{Z}}^{(m)};\mathit{\gamma })=\mathit{\varphi }(\mathit{T}({\mathit{Z}}^{(m)};{\mathit{\gamma }}_{\psi });{\mathit{\gamma }}_{\varphi }),\\ \mathit{T}({\mathit{Z}}^{(m)};{\mathit{\gamma }}_{\psi })=\mathit{a}(\{\mathit{\psi }({\mathit{Z}}_i;{\mathit{\gamma }}_{\psi }):i=1,\dots ,m\}),\end{array}$$(5) where $\mathit{\psi }:{\mathbb{R}}^n\to {\mathbb{R}}^q$ and $\mathit{\varphi }:{\mathbb{R}}^q\to {\mathbb{R}}^p$ are (deep) neural networks parametrized by ${\mathit{\gamma }}_{\psi }$ and ${\mathit{\gamma }}_{\varphi }$, respectively, $\mathit{\gamma }\equiv ({\mathit{\gamma }}_{\varphi }^{\prime },{\mathit{\gamma }}_{\psi }^{\prime })\prime$, and $\mathit{a}:{({\mathbb{R}}^q)}^m\to {\mathbb{R}}^q$ is a permutation-invariant (set) function. Figure 1 shows a schematic of this representation. Each element of $\mathit{a}(·)$, $a_j(·),j=1,\dots ,q$, can be chosen to be a simple aggregation function, such as elementwise addition, average, or maximum, but may also be parameterized (Soelch et al. 2019); in this article, we use the elementwise average. Beyond its attractive parsimonious form, (5) has appealing theoretical properties. For example, Han et al. (2022) show that functions of the form (5) can approximate any continuously differentiable permutation-invariant real-valued function for sufficiently large q, the dimension of $\mathit{T}(·)$; see also Wagstaff et al. (2019, 2022) for relevant discussion. The permutation-invariance property of Bayes estimators (Section S1) coupled with the representational capacity of (5) make DeepSets a principled and theoretically motivated framework for constructing neural Bayes estimators for replicated data.

Fig. 1 Schematic of the DeepSets representation. Each independent replicate ${\mathit{Z}}_i$ is transformed independently using the function $\mathit{\psi }(·)$. The set of transformed inputs are then aggregated elementwise using a permutation-invariant function, $\mathit{a}(·)$, yielding the summary statistic T. Finally, the summary statistic is mapped to parameter estimates estimates $\widehat{\mathit{\theta }}$ by the function $\mathit{\varphi }(·)$. Many classical estimators take this form; in this work, we use neural networks to model $\mathit{\psi }(·)$ and $\mathit{\varphi }(·)$.

Furthermore, the representation (5) is similar in form to many classical estimators when viewed as a nonlinear mapping $\mathit{\varphi }(·)$ of summary statistics $\mathit{T}(·)$. For example, best (i.e., minimum variance) unbiased estimators for exponential family models are of the form (5) where $\mathit{T}(·)$ is sufficient for $\mathit{\theta }$ (Casella and Berger 2001, chap. 7). This connection provides an additional rationale for adopting the structure in (5), and provides interpretability: the functions $\mathit{\psi }(·)$ and $\mathit{a}(·)$ together extract summary statistics from the data, while $\mathit{\varphi }(·)$ maps these learned summary statistics to parameter estimates.

Some summary statistics are available in closed form, simple to compute, and highly informative (e.g., sample quantiles). Denote these statistics as $\mathit{S}(·)$. One may choose to explicitly incorporate these statistics by making $\mathit{\varphi }(·)$ in (5) a function of both $\mathit{T}(·)$ (learned) and $\mathit{S}(·)$ (user-defined). The estimator remains permutation invariant provided that $\mathit{S}(·)$ is permutation invariant. Since $\mathit{T}(·)$ can theoretically approximate well any summary statistic as a continuous function of the data, the choice to include $\mathit{S}(·)$ in (5) is mainly a practical one that could be useful in certain applications.

2.2.2 Variable Sample Size

Estimators of the form (5) can be applied to datasets of arbitrary size. However, the Bayes estimator for replicated data is a function of the number of replicates (see, e.g., (8) in the illustrative example of Section 2.2.3). Therefore, a neural Bayes estimator of the form (5) trained on datasets with m replicates will generally not be Bayes for datasets containing $\tilde{m}\ne m$ replicates.

There are at least two approaches that could be adopted if one wishes to use a neural Bayes estimator of the form (5) with datasets containing an arbitrary number of replicates, m. First, one could train l neural Bayes estimators for different sample sizes, or groups thereof (e.g., a small-sample estimator and a large-sample estimator). Specifically, for sample-size change-points $m_1<m_2<\cdots <m_{l-1}$, one could construct a piecewise neural Bayes estimator, (6) $$\widehat{\mathit{\theta }}({\mathit{Z}}^{(m)};{\mathit{\gamma }}^{*})=\{\begin{array}{cc}\widehat{\mathit{\theta }}({\mathit{Z}}^{(m)};{\mathit{\gamma }}_{{\tilde{m}}_1}^{*})\hfill & m\le m_1,\hfill \\ \widehat{\mathit{\theta }}({\mathit{Z}}^{(m)};{\mathit{\gamma }}_{{\tilde{m}}_2}^{*})\hfill & m_1<m\le m_2,\hfill \\ \hspace{1em}⋮\hfill & \hfill \\ \widehat{\mathit{\theta }}({\mathit{Z}}^{(m)};{\mathit{\gamma }}_{{\tilde{m}}_l}^{*})\hfill & m>m_{l-1},\hfill \end{array}$$(6) where, here, ${\mathit{\gamma }}^{*}\equiv ({\mathit{\gamma }}_{{\tilde{m}}_1}^{*},\dots ,{\mathit{\gamma }}_{{\tilde{m}}_l}^{*})$, and where ${\mathit{\gamma }}_{\tilde{m}}^{*}$ are the neural-network parameters optimized for sample size $\tilde{m}$; $\widehat{\mathit{\theta }}(·;{\mathit{\gamma }}_{\tilde{m}}^{*})$ is then assumed to be approximately Bayes over the range of sample sizes for which it is then applied in (6). Typically, ${\tilde{m}}_1\le m_1,m_1<{\tilde{m}}_2\le m_2$, and so on. We find that this approach works well in practice, and that it is less computationally burdensome than it first appears when used in conjunction with pre-training (see Section 2.3). Note that the relative influence of the prior distribution diminishes as the sample size increases, so that only a single “large sample” estimator is needed. Alternatively, one could treat the sample size as a random variable, M, with support over a set of positive integers, $\mathcal{M}$, in which case (1) becomes (7) $$\begin{array}{c}R(\mathit{\theta },\widehat{\mathit{\theta }}(·;\mathit{\gamma }))\equiv \sum _{m\in \mathcal{M}}\text{Pr}(M=m)\\ \hspace{1em}\times \left({\int }_{{\mathcal{S}}^m}L(\mathit{\theta },\widehat{\mathit{\theta }}({\mathit{z}}^{(m)};\mathit{\gamma }))f({\mathit{z}}^{(m)}|\mathit{\theta })\mathrm{d}{\mathit{z}}^{(m)}\right).\end{array}$$(7)

This does not materially alter the workflow, except that one must sample the number of replicates before simulating data for (4). These two approaches can also be combined, so that each sub-estimator in (6) is trained using (7). We illustrate the importance in accounting for the dependence of the Bayes estimator on m in Section S2 of the supplementary material.

2.2.3 Illustrative Example

We now present a relatively simple example that demonstrates that the DeepSets representation can approximate well Bayes estimators for replicated data. Consider the problem of estimating θ from data $Z_1,\dots ,Z_m$ that are iid according to a uniform distribution on $[0,\theta ]$. Suppose that we choose a conjugate Pareto $(\alpha ,\beta )$ prior for θ with shape α = 4 and scale β = 1; that is, $\text{Pr}(\theta \le x)=1-{(x/\beta )}^{-\alpha },x\ge \beta$. Under the absolute-error loss, the unique Bayes estimator is the posterior median which, for this model, is available in closed form: for a Pareto $(\alpha ,\beta )$ prior distribution, the Bayes estimator given m independent replicates ${\mathit{Z}}^{(m)}\equiv (Z_1,\dots ,Z_m)\prime$ is (8) $${\widehat{\theta }}_{\text{Bayes}}({\mathit{Z}}^{(m)})=2^{\frac{1}{\alpha +m}}\max (Z_1,\dots ,Z_m,\beta ),$$(8) which is clearly invariant to permutations of the observations, and can be expressed as a scalar rendition of (5), where $\psi (z)\equiv z,a(\{·\})\equiv \max (\{·\})$, and $\varphi (t)\equiv 2^{1/(\alpha +m)}\max (t,\beta )$. Note that when $a(\{·\})$ is instead chosen to be the set average, alternative nonlinear functions $\psi (z)$ and $\varphi (t)$ exist such that the estimator (5) approximates (8) arbitrarily well by approximating the maximum as a log-sum-exp.

For this model, the “one-at-a-time” estimator, which applies the single-replicate Bayes estimator to each replicate independently and averages the resulting estimates, is given by (9) $${\widehat{\theta }}_0({\mathit{Z}}^{(m)})=\frac{1}{m}\sum _{i=1}^m{2}^{\frac{1}{\alpha +1}}\max (Z_i,\beta ).$$(9)

Although (9) is permutation invariant, it is clearly very different to (8). In Section S3 of the Supplementary Material we show that the one-at-a-time estimator for this model is not consistent, while the Bayes estimator is consistent. We compare our proposed estimator to the one-at-a-time estimator in this example and in the remainder of the article, to stress the importance of properly accounting for replication when constructing neural Bayes estimators.

We now proceed with the design of the neural Bayes estimator as though we had no knowledge of its closed-form expression. We construct our neural Bayes estimator to have the form (5), where we use three-layer neural networks to model $\mathit{\varphi }(·;\mathit{\gamma })$ and $\mathit{\psi }(·;\mathit{\gamma })$, and we let J = 1 and $K={10}^6$ in (4) when training the network. For simplicity, we consider a single value for m and fix m = 10. We compare the distributions of the estimators at $\theta =4/3$ and m = 10 by simulating 30,000 datasets of size m = 10 from a uniform distribution on $[0,4/3]$ and then applying the estimators to each of these simulated datasets. Figure 2 shows the kernel-smoothed distribution of the true Bayes estimator (8), the neural Bayes estimator, the one-at-a-time estimator (9), and the maximum likelihood estimator which, for this model, is given by $\max (Z_1,\dots ,Z_m)$. The distribution of our neural Bayes estimator is clearly very similar to that of the true Bayes estimator, while the one-at-a-time estimator is clearly biased. The kernel-smoothed distribution of the maximum likelihood estimator emphasizes that (neural) Bayes estimators are influenced by the prior distribution; in this case, the prior distribution leads to a Bayes estimator with a bimodal distribution due to a point mass in the lower endpoint of the support. This example also shows clearly that neural point estimators trained using (4) target the Bayes estimator, and not necessarily the maximum likelihood estimator; this notion is sometimes misconstrued in the literature.

Fig. 2 Kernel density approximations to the distribution of the Bayes estimator (green line), our neural Bayes estimator (red line), the maximum likelihood estimator (purple line), and the one-at-a-time estimator (orange line), for θ from $\text{Unif}(0,\theta )$ data, where $\theta =4/3$ (gray dashed line) and where the sample size is m = 10.

2.3 Implementation

Neural Bayes estimators are conceptually simple and can be used in a wide range of problems where other approaches are computationally infeasible. They also have marked practical appeal, as the general workflow for their construction is only loosely connected to the model being considered. The workflow for constructing a neural Bayes estimator is as follows:

1. Define the prior, $\Omega (·)$.

2. Choose a loss function, $L(·,·)$, typically the absolute-error or squared-error loss.

3. Design a suitable neural-network architecture for the neural point estimator $\widehat{\mathit{\theta }}(·;\mathit{\gamma })$.

4. Sample parameters from $\Omega (·)$ to form training/validation/test parameter sets.

5. Given the above parameter sets, simulate data from the model, to form training/validation/test datasets.

6. Train the neural network (i.e., estimate $\mathit{\gamma }$) by minimizing the loss function averaged over the training sets. That is, perform the optimization task (3) where the Bayes risk is approximated using (4) with the training sets. During training, monitor performance and convergence using the validation sets.

7. Assess the fitted neural Bayes estimator, $\widehat{\mathit{\theta }}(·;{\mathit{\gamma }}^{*})$, using the test set.

We elaborate on the steps of this workflow below. A crucial factor in the viability of any statistical method is the availability of software; to facilitate the adoption of neural Bayes estimators by statisticians, we accompany the article with user-friendly open-source software in the Julia (Bezanson et al. 2017) and R (R Core Team 2023) programming languages (available at https://github.com/msainsburydale/NeuralEstimators) that can be used in a wide range of settings.

2.3.1 Defining the Prior

Prior disributions are typically determined by the applied problem being tackled. However, the choice of prior has practical implications on the training phase of the neural network. For example, an informative prior with compact and narrow support reduces the volume of the parameter space that must be sampled from when evaluating (4). In this case, a good approximation of the Bayes estimator can typically be obtained with smaller values of K and J in (4) than those required under a diffuse prior. This consideration is particularly important when the number of parameters, p, is large, since the volume of the parameter space increases exponentially with p. On the other hand, if the neural Bayes estimator needs to be re-used for several applications, it might be preferable to employ prior distributions that are reasonably uninformative. Lenzi et al. (2023) suggest using likelihood-based estimates to elicit an informative prior from the data: however, this requires likelihood estimation to be feasible in the first place, which is often not the case in applications for which neural Bayes estimators are attractive.

2.3.2 Designing the Neural-Network Architecture

The main consideration when designing the neural-network architecture is the structure of the data. For example, if the data are gridded, a convolutional neural network (CNN) may be used, while a dense neural network (DNN; also known as a multi-layer perceptron, MLP) is more appropriate for unstructured data (for further examples, see Goodfellow, Bengio, and Courville 2016). When estimating parameters from replicated data using the representation (5), the architecture of $\mathit{\psi }(·)$ is dictated by the structure of the data (since it acts on the data directly), while $\mathit{\varphi }(·)$ is typically a DNN (since $\mathit{T}(·)$ is a fixed-dimensional vector). Note that the training cost of the neural Bayes estimator is only amortized if the structure of new datasets conforms with the chosen architecture (otherwise the neural Bayes estimator will need to be re-trained with a new architecture). Once the general neural-network class is identified, one must specify architectural hyperparameters, such as the number of layers and number of neurons in each layer. We found that our results were not overly-sensitive to the specific choice of hyperparameters. In practice, the network must be sufficiently large so that the universal approximation theorem applies, but not so large as to make training prohibitively expensive. We give details on the architectures used in our experiments in Section 3.

2.3.3 Simulating Parameters/Data and Training the Network

In standard applications of neural networks, the amount of training data is fixed. One of the biggest risks one faces when fitting a large neural network is that the amount of training data is too “small”, which can result in the neural network overfitting the training set (Goodfellow, Bengio, and Courville 2016, Section 5.2). Such a neural network is said to have a high generalization error. However, when constructing a neural Bayes estimator, we are able to simulate as much training data as needed. That is, we are able to set K and J in (4) as large as needed to avoid overfitting. The amount of training data needed would depend on the model, the number of parameters, the neural-network architecture, and the optimization algorithm that is used. Providing general guidance is therefore difficult; however, our experience is that K needs to be at least 10⁴–10⁵ for the neural network to have low generalization error. Note that J can be kept small (on the order of 10°–10¹) since data are simulated for every sampled parameter vector. One also has the option to simulate training data “on-the-fly”, in the sense that new training data are simulated continuously during training (Chan et al. 2018). This approach facilitates the use of large networks with a high representational capacity, since then the data used in stochastic-gradient-descent updates are always different from those in previous updates (see Section S4 of the supplementary material for more details and for an illustration). On-the-fly simulation also allows the data to be simulated “just-in-time,” in the sense that the data can be simulated from a small batch of parameters, used to train the estimator, and then removed from memory; this can reduce memory requirements when a large amount of training data are required to avoid overfitting. Chan et al. (2018) also continuously refresh the parameter vectors in the training set, which has similar benefits. Keeping these parameters fixed, however, allows computationally expensive terms, such as Cholesky factors, to be reused throughout training, which can substantially reduce the training time with some models (Gerber and Nychka 2021).

The parameters of a neural network trained for one task can be used as initial values for the parameters of another neural network intended for a slightly different task. This is known as pre-training (Goodfellow, Bengio, and Courville 2016, chap. 8). Pre-training is ideal when developing a piecewise neural Bayes estimator (6), whereby one may randomly initialize and train $\widehat{\mathit{\theta }}(·;{\mathit{\gamma }}_{{\tilde{m}}_1})$, use the optimized parameters ${\mathit{\gamma }}_{{\tilde{m}}_1}^{*}$ as initial values when training $\widehat{\mathit{\theta }}(·;{\mathit{\gamma }}_{{\tilde{m}}_2})$, where ${\tilde{m}}_2>{\tilde{m}}_1$, and so on. Since each estimator need only account for a larger sample size than that used by its predecessor, scant computational resources are needed to train subsequent estimators. This approach can be useful even when only a single large-sample estimator is needed: doing most of the learning with small, computationally cheap sample sizes can substantially reduce training time when compared to training a single large-sample estimator from scratch. We demonstrate the benefits of this strategy through an illustration in Section S5 of the supplementary material.

2.3.4 Assessing the Estimator

Once a neural Bayes estimator is trained, its performance needs to be assessed on unseen test data. Simulation-based (empirical) methods are ideal, since simulation is already required for constructing the estimator. Neural Bayes estimators are very fast to evaluate, and can therefore be applied to thousands of simulated datasets at almost no computational cost. One can therefore quickly and accurately assess the estimator with respect to any property of its sampling distribution (e.g., bias, variance, etc.).

3 Simulation Studies

We now conduct several simulation studies that clearly demonstrate the utility of neural Bayes estimators in increasingly complex settings. In Section 3.1, we outline the general setting. In Section 3.2 we estimate the parameters of a Gaussian process model with three unknown parameters. Since the likelihood function is available for this model, here we compare the efficiency of our neural Bayes estimator to that of a gold-standard likelihood-based estimator. In Section 3.3 we consider a spatial extremes setting and estimate the two parameters of Schlather’s max-stable model (Schlather 2002); the likelihood function is computationally intractable for this model, and we observe substantial improvements over the classical composite-likelihood approach. In Section 3.4 we consider the highly-parameterized conditional extremes model (Wadsworth and Tawn 2022). We provide simulation details and density functions for each model in Section S6 of the supplementary material. For ease of notation, we omit dependence on the neural-network parameters $\mathit{\gamma }$ in future references to the neural Bayes estimators, so that $\widehat{\mathit{\theta }}(·;\mathit{\gamma })$ is written simply as $\widehat{\mathit{\theta }}(·)$.

3.1 General Setting

Across the simulation studies we assume, for ease of exposition, that our processes are spatial. Our spatial domain of interest, $\mathcal{D}$, is $[0,16]\times [0,16]$, and we simulate data on a regular grid with unit-square cells, yielding ${16}^2=256$ observations per spatial field. CNNs are a natural choice for regularly spaced gridded data, and we therefore use a CNN architecture, summarized in Table S3 of the Supplementary Material. To implement the neural point estimators, we use the accompanying package NeuralEstimators that is written in Julia and which leverages the package Flux (Innes 2018). We conduct our studies using a workstation with an AMD EPYC 7402 3.00GHz CPU with 52 cores and 128 GB of CPU RAM, and an NVIDIA Quadro RTX 6000 GPU with 24 GB of GPU RAM. All results presented in the remainder of this article can be generated using reproducible code at https://github.com/msainsburydale/NeuralBayesEstimators.

We consider two neural point estimators which are both based on the architecture given in Table S3 of the supplementary material. The first estimator, ${\widehat{\mathit{\theta }}}_0(·)$, is the permutation-invariant one-at-a-time neural point estimator considered by Gerber and Nychka (2021). The second estimator, ${\widehat{\mathit{\theta }}}_{\text{DS}}(·)$, is the piecewise neural point estimator (6), where each sub-estimator employs the DeepSets representation (5) with $\mathit{\psi }(·)$ and $\mathit{\varphi }(·)$ constructed using the first four and last two rows of Table S3, respectively. Five sub-estimators are used, with training sample sizes ${\tilde{m}}_1=1,{\tilde{m}}_2=10,{\tilde{m}}_3=35,{\tilde{m}}_4=75$, and ${\tilde{m}}_5=150$, and with sample-size changepoints $m_1=1,m_2=20,m_3=50,$ and $m_4=100$.

We assume that the parameters are a priori independent and uniformly distributed on an interval that is parameter dependent. We train the neural point estimators under the absolute-error loss. We set K in (4) to 10,000 and 2000 for the training and validation parameter sets, respectively, and we keep these sets fixed during training. We construct the training and validation datasets by simulating J = 10 sets of m model realizations for each parameter vector in the training and validation parameter sets. During training, we fix the validation data, but simulate the training data on-the-fly; hence, in this article, we define an epoch as a pass through the training sets when doing stochastic gradient descent, after which the training data (i.e., the J = 10 datasets at each of the $10,000$ parameter samples) are refreshed. We cease training when the risk evaluated using the validation set has not decreased in five consecutive epochs.

We compare the neural Bayes estimators to the maximum a posteriori (MAP) estimator, (10) $${\widehat{\mathit{\theta }}}_{\text{MAP}}({\mathit{Z}}^{(m)})\equiv \underset{\mathit{\theta }\in \Theta }{\text{arg} \max } \sum _{i=1}^m\mathcal{l}(\mathit{\theta };{\mathit{Z}}_i)+\hspace{0.17em}\log \hspace{0.17em}p(\mathit{\theta }),$$(10) where $\mathcal{l}(\mathit{\theta };·)$ is the log-likelihood function for a single replicate and $p(\mathit{\theta })$ is the prior density. We solve (10) using the true parameters as initial values for the Nelder–Mead algorithm. We advantage the competitor MAP estimator, which minimizes (2) under the 0–1 loss, by assessing all estimators with respect to the 0–1 loss (we assign zero loss if an estimate is within 10% of the true value), with K = 500 test parameter vectors.

3.2 Gaussian Process Model

Spatial statistics is concerned with modeling data that are collected across space; reference texts include Cressie (1993), Banerjee, Carlin, and Gelfand (2004), and Diggle and Ribeiro (2007). Here, we consider a classical spatial model, the linear Gaussian-Gaussian model, (11) $$Z_{ij}=Y_i({\mathit{s}}_j)+{ϵ}_{ij},\hspace{1em}i=1,\dots ,m,j=1,\dots ,n,$$(11) where ${\mathit{Z}}_i\equiv (Z_{i1},\dots ,Z_{in})\prime$ are data observed at locations $\{{\mathit{s}}_1,\dots ,{\mathit{s}}_n\}$ on a spatial domain $\mathcal{D}$, $\{Y_i(·)\}$ are iid spatially-correlated mean-zero Gaussian processes, and ${ϵ}_{ij}\sim \mathrm{N}(0,{\sigma }_{ϵ}^2)$ is Gaussian white noise. The covariance function, $C(\mathit{s},\mathit{u})\equiv \text{cov}\mathrm{}(Y_i(\mathit{s}), Y_i(\mathit{u}))$, for $\mathit{s},\mathit{u}\in \mathcal{D}$ and $i=1,\dots ,m$, is the primary mechanism for capturing spatial dependence. Note that, since $\{Y_i(·)\}$ are iid, we have that $\text{cov}\mathrm{}(Y_i(\mathit{s}), Y_{i\prime }(\mathit{u}))=0$ for all $i\ne i\prime$ and $\mathit{s},\mathit{u}\in \mathcal{D}$. Here, we use the popular isotropic Matérn covariance function, (12) $$C(\mathit{s},\mathit{u})={\sigma }^2\frac{2^{1-\nu }}{\Gamma (\nu )}{\left(\frac{⏧\mathit{s}-\mathit{u}⏧}{\rho }\right)}^{\nu }{K}_{\nu }\mathrm{}\left(\frac{⏧\mathit{s}-\mathit{u}⏧}{\rho }\right),$$(12) where ${\sigma }^2$ is the marginal variance, $\Gamma (·)$ the gamma function, $K_{\nu }(·)$ is the Bessel function of the second kind of order ν, and $\rho >0$ and $\nu >0$ are the range and smoothness parameters, respectively. We follow Gerber and Nychka (2021) and fix ${\sigma }^2=1$. This leaves three unknown parameters that need to be estimated: $\mathit{\theta }\equiv ({\sigma }_{ϵ},\rho ,\nu )\prime$.

Fixing ν simplifies the estimation task, since the remaining parameters are well-identified from a single replicate; we consider this model, which was also considered by Gerber and Nychka (2021), in Section S7 of the supplementary material. Here, we consider the case where all three parameters are unknown. Estimation in this case is more challenging since ρ and ν are only weakly identifiable from a single replicate (Zhang 2004), but inference on ν is important since it controls the differentiability of the process (Stein 1999, chap. 2).

We use the priors ${\sigma }_{ϵ}\sim \text{Unif}(0.1,1)$, $\rho \sim \text{Unif}(2,10)$, and $\nu \sim \text{Unif}(0.5,3)$. The total training time for ${\widehat{\mathit{\theta }}}_0(·)$ is 15 min. The training time for ${\widehat{\mathit{\theta }}}_{\text{DS}}(·)$, despite it consisting of several sub-estimators and being trained with large sample sizes, is only slightly more, at 33 min, due to the pre-training strategy discussed in Section 2.3. This modest increase in training time is a small cost when accounting for the improvement in statistical efficiency, as shown in Figure 3. Clearly, ${\widehat{\mathit{\theta }}}_{\text{DS}}(·)$ substantially improves over ${\widehat{\mathit{\theta }}}_0(·)$, and performs similarly well to the MAP estimator. The neural point estimators only take 0.1 and 2 sec to yield estimates from all the test data when m = 1 and m = 150, respectively, while the MAP estimator takes between 600 and 800 sec.

Fig. 3 The risk with respect to the 0–1 loss against the number of replicates, m, evaluated using the parameter vectors in the test parameter set, for the estimators considered in Section 3.2. The estimators ${\widehat{\mathit{\theta }}}_0(·)$, ${\widehat{\mathit{\theta }}}_{\text{DS}}(·)$, and the MAP estimator are shown in green, orange, and purple, respectively.

Next, we compare the empirical joint distributions of the estimators for a single parameter vector. Figure 4 shows the true parameters (red cross) and corresponding estimates obtained by applying each estimator to 100 sets of m = 150 independent model realizations. The MAP estimator and ${\widehat{\mathit{\theta }}}_{\text{DS}}(·)$ are both approximately unbiased, and both are able to capture the correlations expected when using the parameterization (12); both are clearly much better than ${\widehat{\mathit{\theta }}}_0(·)$. Empirical joint distributions for additional parameter vectors are shown in Figure S11 of the supplementary material, and lead to similar conclusions to those drawn here. Overall, these results show that ${\widehat{\mathit{\theta }}}_{\text{DS}}(·)$ is a substantial improvement over the prior art, ${\widehat{\mathit{\theta }}}_0(·)$, and that ${\widehat{\mathit{\theta }}}_{\text{DS}}(·)$ is competitive with the likelihood-based estimator for this model.

Fig. 4 The empirical joint distribution of the estimators considered in Section 3.2 for a single parameter vector. The true parameters are shown in red, while estimates from ${\widehat{\mathit{\theta }}}_0(·)$, ${\widehat{\mathit{\theta }}}_{\text{DS}}(·)$, and the MAP estimator are shown in green, orange, and purple, respectively. Each estimate was obtained from a simulated dataset of size $m=150$.

3.3 Schlather’s Max-Stable Model

We now consider models used for spatial extremes, useful reviews for which are given by Davison and Huser (2015), Davison, Huser, and Thibaud (2019), and Huser and Wadsworth (2022). Max-stable processes are the cornerstone of spatial extreme-value analysis, being the only possible nondegenerate limits of properly renormalized pointwise block maxima of iid random fields. However, their practical use has been severely hampered due to the computational bottleneck in evaluating their likelihood function. They are thus natural models to consider in our experiments, and we here consider a fairly simple max-stable model with only two parameters. We consider Schlather’s max-stable model (Schlather 2002), given by (13) $$Z_{ij}=\underset{k\in \mathbb{N}}{⏧}{\zeta }_{ik}\max \{0,Y_{ik}({\mathit{s}}_j)\},\hspace{1em}i=1,\dots ,m,j=1,\dots ,n,$$(13) where $⏧$ denotes the maximum over the indexed terms, ${\mathit{Z}}_i\equiv (Z_{i1},\dots ,Z_{in})\prime$ are observed at locations $\{{\mathit{s}}_1,\dots ,{\mathit{s}}_n\}\subset \mathcal{D}$, $\{{\zeta }_{ik}:k\in \mathbb{N}\}$ for $i=1,\dots ,m$ are iid Poisson point processes on $(0,\infty )$ with intensity measure $\mathrm{d}\Lambda (\zeta )={\zeta }^{-2}\mathrm{d}\zeta$, and $\{Y_{ik}(·):i=1,\dots ,m,k\in \mathbb{N}\}$ are iid mean-zero Gaussian processes scaled so that $\mathbb{E}[\max \{0,Y_{ik}(·)\}]=1$. Here, we model each $Y_{ik}(·)$ using the Matérn covariance function (12), with ${\sigma }^2=1$. Hence, $\mathit{\theta }\equiv (\rho ,\nu )\prime$.

We use the same uniform priors for ρ and ν as in Section 3.2. Realizations from the present model, here expressed on unit Fréchet margins, tend to have highly varying magnitudes, and we reduce this variability by log-transforming our data to the unit Gumbel scale. The total training time for ${\widehat{\mathit{\theta }}}_0(·)$ and ${\widehat{\mathit{\theta }}}_{\text{DS}}(·)$ is 14 and 66 min, respectively.

As in Section 3.2, we assess the neural point estimators by comparing them to a likelihood-based estimator. For Schlather’s model (and other max-stable models in general), the full likelihood function is computationally intractable, since it involves a summation over the set of all possible partitions of the spatial locations (see, e.g., Padoan 2010; Huser et al. 2019, and the references therein). A popular substitute is the pairwise likelihood (PL) function, a composite likelihood formed by considering only pairs of observations; specifically, the pairwise log-likelihood function for the ith replicate is (14) $${\mathcal{l}}_{\text{PL}}(\mathit{\theta };{\mathit{z}}_i)\equiv \sum _{j=1}^{n-1}\sum _{j\prime =j+1}^n\hspace{0.17em}\log \hspace{0.17em}f(z_{ij},z_{ij\prime }|\mathit{\theta }),$$(14) where $f(·,·|\mathit{\theta })$ denotes the bivariate probability density function for pairs in ${\mathit{z}}_i$. Hence, in this subsection, we compare the neural Bayes estimators to the pairwise MAP (PMAP) estimator, that is, (10) with the full log-likelihood function $\mathcal{l}(\mathit{\theta };·)$ replaced by ${\mathcal{l}}_{\text{PL}}(\mathit{\theta };·)$. Often, both computational and statistical efficiency can be drastically improved by using only a subset of pairs that are within a fixed cutoff distance, d (see, e.g., Bevilacqua et al. 2012; Sang and Genton 2012). A line-search for d (see Figure S12 of the supplementary material for details) shows that, here, d = 3 units (used hereafter) provides good results.

The left and center panels of Figure 5 show the estimators’ risk against the number of independent replicates. For small samples, both neural point estimators improve over the PMAP estimator. For moderate-to-large samples, ${\widehat{\mathit{\theta }}}_0(·)$ hits a performance plateau, while ${\widehat{\mathit{\theta }}}_{\text{DS}}(·)$ continues to substantially outperform the PMAP estimator. The run time for the neural point estimators to estimate all test parameters scales linearly between 0.1 and 1.5 sec for m = 1 and m = 150, respectively, while the PMAP estimator takes between 750 and 1900 sec. The right panel of Figure 5 shows the empirical joint distribution of the estimators for a single parameter vector, where each estimate was obtained from m = 150 replicates. Again, ${\widehat{\mathit{\theta }}}_0(·)$ is strongly biased, while the PMAP estimator is unbiased but is less efficient than ${\widehat{\mathit{\theta }}}_{\text{DS}}(·)$. Empirical joint distributions for additional parameter vectors are shown in Figure S13 of the supplementary material, and lead to similar conclusions to those drawn here. Overall, the proposed estimator, ${\widehat{\mathit{\theta }}}_{\text{DS}}(·)$, is statistically and computationally superior to the likelihood-based technique for Schlather’s max-stable model.

Fig. 5 Diagnostic plots for the simulation study of Section 3.3. (Left and center) The risk with respect to the 0–1 loss against the number of replicates, m, evaluated using the parameter vectors in the test parameter set. (Right) True parameters (red) and corresponding estimates, each of which was obtained using a simulated dataset of size m = 150. In all panels, ${\widehat{\mathit{\theta }}}_0(·)$, ${\widehat{\mathit{\theta }}}_{\text{DS}}(·)$, and the pairwise MAP estimator are shown in green, orange, and purple, respectively.

3.4 Spatial Conditional Extremes Model

While max-stable processes are asymptotically justified for modeling spatial extremes defined as block maxima, they have strong limitations in practice (Huser and Wadsworth 2022). Beyond the intractability of their likelihood function in high dimensions, max-stable models have an overly rigid dependence structure, and in particular cannot capture a property known as asymptotic tail independence. Therefore, different model constructions, justified by alternative asymptotic conditions, have been recently proposed to circumvent the restrictions imposed by max-stability. In particular, the spatial conditional extremes model, first introduced by Wadsworth and Tawn (2022) as a spatial extension of the multivariate Heffernan and Tawn (2004) model, and subsequently studied and used in applications by, for example, Richards, Tawn, and Brown (2022) and Simpson, Opitz, and Wadsworth (2023), is especially appealing. This model has a flexible dependence structure capturing both asymptotic tail dependence and independence (Wadsworth and Tawn 2022) and leads to likelihood-based inference amenable to higher dimensions. Nevertheless, parameter estimation remains challenging because this model is complex and is typically highly parameterized. Here, we consider a version of the spatial conditional extremes model that involves eight dependence parameters in total.

Our formulation of the model is similar to that originally proposed by Wadsworth and Tawn (2022). Specifically, we model the process, expressed on Laplace margins, conditional on it exceeding a threshold, u, at a conditioning site, ${\mathit{s}}_0\in \mathcal{D}$, as (15) $$\begin{array}{c}{Z}_{ij}|Z_{i0}>u\stackrel{\mathrm{d}}{=}a({\mathit{h}}_j,Z_{i0})+b({\mathit{h}}_j,Z_{i0})Y_i({\mathit{s}}_j),\\ i=1,\dots ,m,j=1,\dots ,n,\end{array}$$(15) where “$\stackrel{\mathrm{d}}{=}$” denotes equality in distribution, $Z_{i0}$ is the datum at ${\mathit{s}}_0,{\mathit{h}}_j\equiv {\mathit{s}}_j-{\mathit{s}}_0$, and $Z_{i0}-u|Z_{i0}>u$ is a unit exponential random variable that is independent of the residual process, $Y_i(·)$, which we describe below. We model $a(·,·)$ and $b(·,·)$ using parametric forms proposed by Wadsworth and Tawn (2022), namely $$\begin{array}{cccc}a(\mathit{h},z)\hfill & =z\hspace{0.17em}\exp \hspace{0.17em}\{-{(⏧\mathit{h}⏧/\lambda )}^{\kappa }\},\hfill & \hfill & \hspace{1em}\lambda >0,\kappa >0,\hfill \\ b(\mathit{h},z)\hfill & =1+a{(\mathit{h},z)}^{\beta },\hfill & \hfill & \hspace{1em}\beta >0,\hfill \end{array}$$ where $\mathit{h}\equiv \mathit{s}-{\mathit{s}}_0$ for $\mathit{s}\in \mathcal{D}$ and $z\in \mathbb{R}$. We construct the residual process, $Y_i(·)$, by first defining ${\tilde{Y}}_i^{(0)}(·)\equiv {\tilde{Y}}_i(·)-{\tilde{Y}}_i({\mathit{s}}_0)$, with ${\tilde{Y}}_i(·)$a mean-zero Gaussian process with Matérn covariance function and unit marginal variance, and then marginally transforming it to the scale of a delta-Laplace (generalized Gaussian) distribution (Subbotin 1923) which, for $y\in \mathbb{R}$ and parameters $\mu \in \mathbb{R},\tau >0,\delta >0$, has density function $$f_S(y|\mu ,\tau ,\delta )=\frac{\delta }{2\tau \Gamma (1/\delta )}\exp \hspace{0.17em}\left(-{\left|\frac{y-\mu }{\tau }\right|}^{\delta }\right).$$

We model δ as decaying from 2 to 1 as the distance to ${\mathit{s}}_0$ increases; specifically, (16) $$\delta (\mathit{h})=1+\hspace{0.17em}\exp \hspace{0.17em}\left\{-{(\Vert \mathit{h}\Vert /{\delta }_1)}^2\right\},\hspace{1em}{\delta }_1>0.$$(16)

We defer to Wadsworth and Tawn (2022) for model justification and interpretation of model parameters. The threshold u is a modeling decision made to ensure that the data are sufficiently extreme and that, in a real-data setting, we have sufficiently many fields available for estimation: here, we set u to the 0.975 quantile of the unit-Laplace distribution. For simplicity, we consider ${\mathit{s}}_0$ fixed and in the center of $\mathcal{D}$.

We again use uniform priors, with $\kappa \sim \text{Unif}(1,2)$, $\lambda \sim \text{Unif}(2,5)$, $\beta \sim \text{Unif}(0.05,1)$, $\mu \sim \text{Unif}(-0.5,0.5)$, $\tau \sim \text{Unif}(0.3,0.9)$, ${\delta }_1\sim \text{Unif}(1.3,3)$, and with the same priors for ρ and ν as used in the preceding sections. We use the cube-root function as a variance-stabilizing transformation. The training time for ${\widehat{\mathit{\theta }}}_0(·)$ and ${\widehat{\mathit{\theta }}}_{\text{DS}}(·)$ is 22 and 43 min, respectively.

Figure 6 shows the risk against the sample size. Figure 7 shows the empirical joint distribution of the estimators for a single parameter vector (panels in the lower triangle) and simulations from the corresponding model (panels in the upper triangle). The estimator ${\widehat{\mathit{\theta }}}_{\text{DS}}(·)$ is approximately unbiased for all parameters, and captures the expected negative correlation between ρ and ν. Overall, the estimator ${\widehat{\mathit{\theta }}}_{\text{DS}}(·)$ is clearly appropriate for this highly parameterized model (unlike ${\widehat{\mathit{\theta }}}_0(·)$), which is an important result for this framework.

Fig. 6 The risk with respect to the 0–1 loss against the number of replicates, m, evaluated using the parameter vectors in the test parameter set, for the estimators considered in Section 3.4. The estimators ${\widehat{\mathit{\theta }}}_0(·)$ and ${\widehat{\mathit{\theta }}}_{\text{DS}}(·)$ are shown in green and orange, respectively.

Fig. 7 (Lower triangle) The empirical joint distribution of the estimators considered in Section 3.4 for a single parameter vector. The true parameters are shown in red, while estimates from ${\widehat{\mathit{\theta }}}_0(·)$ and ${\widehat{\mathit{\theta }}}_{\text{DS}}(·)$ are shown in green and orange, respectively. Each estimate was obtained from a simulated dataset of size m = 150. (Upper triangle) Simulations from the model.

4 Application to Red Sea Surface Temperature

We now apply our methodology to the analysis of sea-surface temperature data in the Red Sea, which have also been analyzed by Hazra and Huser (2021), Simpson and Wadsworth (2021), and Simpson, Opitz, and Wadsworth (2023), among others, and have been the subject of a competition in the prediction of extreme events (Huser 2021). The dataset we analyze comprises daily observations from the years 1985 to 2015, for 16,703 regularly spaced locations across the Red Sea; see Donlon et al. (2012) for further details. Following Simpson, Opitz, and Wadsworth (2023), we focus on a southern portion of the Red Sea, consider only the summer months to approximately eliminate the effects of seasonality, and retain only every third longitude and latitude value; this yields a dataset with 678 unique spatial locations that are regularly spaced but contained within an irregularly-shaped spatial domain, $\mathcal{D}$. To account for $\mathcal{D}$ lying away from the equator, we follow Simpson, Opitz, and Wadsworth (2023) in scaling the longitude and latitude so that each unit distance in the spatial domain corresponds to approximately 100 km.

To model these data, we use the spatial conditional extremes model described in Section 3.4. We transform our data to the Laplace scale and set the threshold, u in (15), to the 95th percentile of the transformed data. This yields 141 spatial fields for which the transformed datum at the conditioning site, ${\mathit{s}}_0$, here chosen to lie in the center of $\mathcal{D}$, is greater than u; that is, we estimate parameters based on m = 141 replicates of the spatial process. A randomly-selected sample of these extreme fields are shown in the upper panels of Figure S14 of the supplementary material.

The irregular shape of $\mathcal{D}$ means that CNNs are not directly applicable. We use the standard technique of padding empty regions with zeros, so that each field is a 29 × 37 rectangular array consisting of a data region and a padded region, as shown in Figure S15 of the supplementary material. We use the same prior distributions as in Section 3.4, but with those associated with range parameters scaled appropriately. Our architecture is the same as that given in Table S3, but with an additional convolutional layer that transforms the input array to dimension 16 × 16. We validate our neural Bayes estimator using the approach taken in Section 3 (figures omitted for brevity).

Once training is complete, we can compute parameter estimates for the observed data. Table 1 gives estimates and 95% confidence intervals for the estimates. These confidence intervals are obtained using the nonparametric bootstrap procedure described in Section S8 of the Supplementary Material, which accounts for temporal dependence between the spatial fields. Estimation from a single set of 141 fields takes only 0.008 sec, meaning that bootstrap confidence intervals can be obtained very quickly. This is clearly an advantage of using neural point estimators, as opposed to likelihood-based techniques for which uncertainty assessment in complex models is usually a computational burden.

Table 1 Parameter estimates and 95% bootstrap confidence intervals (provided via the 2.5 and 97.5 percentiles of the bootstrap distribution) for the Red Sea dataset of Section 4.

	κ	λ	β	ρ	ν	μ	τ	δ1
Estimate	1.00	2.74	0.24	1.34	0.78	0.09	0.56	1.22
2.5%	0.88	2.13	0.14	0.95	0.70	0.04	0.49	0.91
97.5%	1.19	3.79	0.40	1.85	0.87	0.14	0.65	1.63

Next, following Simpson, Opitz, and Wadsworth (2023), we separate $\mathcal{D}$ into 17 nonoverlapping regions, which are shown in the left panel of Figure 8. Given that the transformed datum at ${\mathit{s}}_0$ exceeds u, we estimate (and quantify the uncertainty in) the proportion of locations in each region that also exceed u, using both model-based and empirical methods. These are shown in the right panel of Figure 8, and indicate overall agreement between the model and the observed data, with some minor lack of fit at very short distances. The empirical estimates are not monotonically decreasing as a function of distance; this could be due to complex spatial dynamics or non-stationarity in the data, or it could simply be an artifact of sampling variability. In either case, the fit is reasonable, which suggests that our neural Bayes estimator has provided reasonable parameter estimates for this dataset.

Fig. 8 (Left) The spatial domain of interest for the Red Sea study of Section 4, with the conditioning site, ${\mathit{s}}_0$, shown in red, and the remaining locations separated into 17 regions; the region labels begin at 1 in the center of the domain and increase with distance from the center. (Right) The estimated proportion of locations for which the process exceeds u given that it is exceeded at ${\mathit{s}}_0$ (points) and corresponding 95% confidence intervals (vertical segments) using model-simulated datasets (blue) and bootstrap samples of the observed dataset (red).

5 Conclusion

Neural Bayes estimators are a class of likelihood-free estimators that approximate Bayes estimators using neural networks. Their connection to classical estimation theory is often under-appreciated, and this article serves to increase the awareness and adoption of this powerful estimation tool by statisticians. This article also proposes a principled way to construct neural Bayes estimators for replicated data via the DeepSets architecture, which has the same structure as that of well-known conventional estimators, such as best unbiased estimators with exponential family models. Using these estimators that are able to automatically learn suitable summary statistics from the data, we jointly estimate the range and smoothness parameters in a Gaussian process model and in Schlather’s max-stable model, and estimate parameters in the highly parameterized spatial conditional extremes model. These estimators are implemented with just a few lines of code, thanks to the package NeuralEstimators, which is released with this article and is available at https://github.com/msainsburydale/NeuralEstimators.

As with all estimation methods, neural Bayes estimators come with advantages and disadvantages that will either make them highly applicable, or impractical, depending on the application. At the time of writing we see five main drawbacks of neural Bayes estimators. First, it is unclear how one could verify that a trained neural Bayes estimator is “close” to the true Bayes estimator, and more theory is needed to provide guarantees and guidelines on how to choose the number of training samples in practice. This limitation can be somewhat mitigated by running empirical checks (as we do in our examples); these checks are not computationally costly given the amortized nature of the neural Bayes estimator. Second, the training cost is amortized only if the estimator is used repeatedly for the same estimation problem; in applications where estimation needs to be done just once, and the likelihood function is available and computationally tractable, classical likelihood-based methods supported by theoretical guarantees and extensive software availability might be more attractive. Third, one needs to be able to simulate relatively easily from the data generating process; although this sampling phase is parallelizable, this step could make the use of neural Bayes estimators with some models impractical. Fourth, the fast construction of neural Bayes estimators requires graphics processing units (GPUs); however, GPUs suitable for deep learning are now commonplace in high-end workstations, and this hardware is only needed for training, and not for evaluating the estimator once trained. Fifth, the volume of the parameter space increases exponentially fast with the number of parameters; we therefore expect training the neural Bayes estimator to be more difficult with highly-parameterized models, especially those with non-orthogonal parameterizations. On the other hand, neural point estimators have substantial advantages that will, in our opinion, make them the preferred and the de facto option in several applications in the near future. First, they are particularly useful when the likelihood function is unavailable or computationally intractable. Second, due to their amortized nature, they are ideal for settings in which the same statistical model must be fit repeatedly (e.g., when solving inverse problems with remote sensing data; see Cressie 2018), or in scenarios where accurate bootstrap-based uncertainty estimates are needed.

There are many avenues for future research. In this work, we have illustrated neural Bayes estimation with spatial data observed over a regular grid. Parameter estimation from irregular spatial data is an important problem, and one that we consider briefly with DNNs in Section S9 of the supplementary material. An attractive way forward in this regard is the use of graph neural networks (GNNs; Wu et al. 2021), which generalize the convolution operation to irregular data; their use in the context of parameter estimation is the subject of ongoing work. Ways to incorporate covariate information also need to be explored. A possible criticism of neural Bayes estimators is that they are not robust to model misspecification since neural networks are, generally, poor extrapolators. However, we did not find this to be the case in our work. In Section S10 of the supplementary material, we provide some empirical evidence that neural Bayes estimators can be used on data that are very different to those used during training, but further research on this topic is needed. Neural Bayes estimators of the form (5) are ideally placed for online learning problems; investigating their potential in this context is the subject of future work. Finally, while this article focuses on neural point estimation, there is a growing literature on neural approaches that approximate the full posterior distribution (Radev et al. 2022; Pacchiardi and Dutta 2022): this is also a promising avenue for future research.

Supplementary Materials

In Section S1, we give a proof for the permutation-invariance property of Bayes estimators for replicated data. In Section S2, we illustrate how neural Bayes estimators depend on the sample size. In Section S3, we show that the Bayes estimator is consistent in the example of Section 2.2.3, but that the one-at-a-time estimator is not. In Sections S4 and S5, we conduct experiments to illustrate the benefits of “on-the-fly” simulation and pre-training, respectively. In Section S6, we describe how we simulate from the models considered in the main text. In Section S7, we describe our analysis for a Gaussian process model with known smoothness parameter. In Section S8, we detail the bootstrap techniques used in Section 4. In Section S9, we apply our methodology to the analysis of sea-surface temperature data in the Red Sea, where the data are measured irregularly in space. In Section S10, we provide some empirical evidence that our neural Bayes estimators can be used on data that are different to those used during training. Finally, in Section S11, we provide additional exploratory and diagnostic plots for the experiments given in the main text.

Acknowledgments

The authors would like to thank Yi Cao for technical support. Thanks also go to Emma Simpson, Jennifer Wadsworth, and Thomas Opitz for providing access to the preprocessed Red Sea dataset. The authors are also grateful to Noel Cressie and Jonathan Rougier for providing comments on the manuscript. We are also grateful to the reviewers and the editors for their helpful comments and suggestions that improved the quality of the manuscript.

Disclosure Statement

The authors report that there are no competing interests to declare.

Additional information

Funding

Matthew Sainsbury-Dale’s and Andrew Zammit-Mangion’s research was supported by an Australian Research Council (ARC) Discovery Early Career Research Award, DE180100203. Matthew Sainsbury-Dale’s research was also supported by an Australian Government Research Training Program Scholarship and a 2022 Statistical Society of Australia (SSA) top-up scholarship. Raphaël Huser was partially supported by the King Abdullah University of Science and Technology (KAUST) Office of Sponsored Research (OSR) under Award No. OSR-CRG2020-4394. The authors would like to acknowledge Johan Barthelemy, NVIDIA, SMART Infrastructure Facility of the University of Wollongong, as well as a 2021 Major Equipment Grant from the University of Wollongong’s University Research Committee, that together have provided access to extensive GPU computing resources.