{'fraction_non_alphanumeric': 0.032774285635859574, 'fraction_numerical': 0.011391397436249348, 'mean_word_length': 4.780104712041885, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Neural Networks for Parameter Estimation in Intractable Models

We propose to use deep learning to estimate parameters in statistical models when standard likelihood estimation methods are computationally infeasible. We show how to estimate parameters from max-stable processes, where inference is exceptionally challenging even with small datasets but simulation is straightforward. We use data from model simulations as input and train deep neural networks to learn statistical parameters. Our neural-network-based method provides a competitive alternative to current approaches, as demonstrated by considerable accuracy and computational time improvements. It serves as a proof of concept for deep learning in statistical parameter estimation and can be extended to other estimation problems.


Introduction
The large datasets that are increasingly available because of advances in data storage and sensor technology inevitably display complex dependencies, bringing new opportunities and new challenges to statistical modeling and prediction. Environmental processes give an important example for which datasets can have complicated interactions across multiple scales of spatial and temporal variability. Specifically, inference for parametric statistical models designed for non-Gaussian data that account for spatiotemporal dependencies are computationally challenging. For instance, consider models for multivariate extremes, where data are usually sampled at a large number of spatial locations but large computational complexities limit their usefulness. In order to circumvent this bottleneck, model approximations that avoid full likelihood computation have been developed. One such approximation is based on composite likelihood, usually constructed from pairs of observations (Tawn, 1988;Padoan et al., 2010;Ribatet et al., 2012). However, even though the maximum pairwise likelihood estimator is strongly consistent and asymptotically normal under mild regularity conditions (Xu and Reid, 2011), no obvious strategy exists for selecting the composite likelihood terms optimally, and the loss in efficiency makes inference for high-dimensional data still challenging (Varin et al., 2011;Huser et al., 2016).
As an alternative to composite likelihood inference in models for spatial extremes, Stephenson and Tawn (2005) suggested augmenting the componentwise block maxima data with their occurrence times and showed that partitioning the locations based on whether maxima coincided simplifies the likelihood function substantially. To avoid the bias that emerges from fixing these partitions, Huser et al. (2019) constructed a stochastic expectation-maximization algorithm for exact maximum likelihood estimation of max-stable processes by treating the partitions as latent variables. Similarly, in a Bayesian setting, Thibaud et al. (2016)

developed a Markov chain
Monte Carlo algorithm for max-stable processes conditioning on the partition values. A possible problem with these approaches is that they assume independence in time, and how much bias this temporal dependence might cause is unclear. Moreover, likelihood-based techniques are still 1 hindered by the computational bottleneck of evaluating multivariate max-stable cumulative distribution functions. A key aspect of models for dependence is that often simulation is fast and tractable, although likelihood computation is burdensome. Taking advantage of that, Erhardt and Smith (2012) employed a simulation-based approach using approximate Bayesian computation (ABC) for estimating parameters of max-stable models. In its most basic form, ABC first samples a set of parameters from the prior distribution, which are then used to simulate data under the assumed statistical model. For each new dataset, a procedure of evaluating discrepancies between the simulations and the observed data must be repeated. As a result, neither of the previously mentioned approaches are well suited for large datasets. Therefore, developing new inference methods for computationally challenging likelihoods remains an open challenge, as is the case of max-stable processes and many other applications.
In this work we propose an approach to estimate the parameters of intractable statistical models using deep neural networks (NNs). Deep learning algorithms utilize NNs as multiple interconnected layers to predict or classify complex maps between datasets, and in recent years they have shown remarkable success at predicting high-dimensional nonlinear processes. Literature on deep NN techniques applied directly in the context of parameter estimation is scarce. Within the ABC framework, Creel (2017) proposed to use informative statistics to train a deep NN upon which classical or Bayesian indirect inference may be based. Jiang et al. (2017) predicted summary statistics from simulated data, which were then used as an estimate of the posterior mean of an ABC procedure. Recently, Radev et al. (2020) proposed to perform Bayesian inference using informative summary statistics from the data as the input to invertible NNs. However, these papers provide no general guidance on how to construct suitable summary statistics for the networks in practice. Convolutional neural networks (CNNs) have been successfully applied to represent covariance structures of Gaussian processes by using sensor measurements (Liu et al., 2018), time-series data (Cremanns and Roos, 2017), and spatial data (Gerber and Nychka, 2020). Rudi et al. (2021) proposed to solve a deterministic inverse problem using dense and convolutional layers in NNs, where parameters from a nonlinear system of ordinary differential equations are estimated from simulated data. To our knowledge, the present work is the first study in which deep NN methodologies are utilized with the goal of performing inference based solely on simulated data and parameters.
Similar to ABC, we perform statistical inference based on simulations rather than direct evaluations of computationally expensive or intractable likelihood functions. A crucial difference from ABC is that instead of measuring similarity as a function of informative summary statistics of the simulated data with the actual observations, our approach learns to perform the mapping directly between data and parameters. Therefore, our method can be seen as a form of automatic ABC, where the data compression is done inside the NN without the difficulties of picking potentially highly sensitive metrics. Our new construction uses simulated data from the statistical model as the input to a deep NN. The output consists of the parameters used to simulate the input data. We leverage the feasibility to simulate arbitrarily large amounts of training data in order to ensure that the deep NNs approximate the actual parameters as well as possible. After the deep NN is trained, the observed data are used as the input in the testing data, returning the model parameters of interest as output. We show that the deep NN can estimate parameters of max-stable models for spatial extremes with higher accuracy than traditional (approximate) likelihood approaches can, with a considerable speed-up in computations. As the number of parameters increases, a substantial decrease in the computational efficiency is observed for ABC methods, and the same will likely occur with our proposed framework. This work holds promising results for parameter estimation and serves as a proof of concept on how machine learning can be used for statistical inference when the traditional techniques fail.
We demonstrate our framework on a problem of determining parameters from max-stable models. Spatial data in extremes can be viewed as a particular case of 2D images, making NNs with convolutional layers a promising candidate due to their success in learning features of complex images. In this work we initially carry out simulated examples for Brown-Resnick 3 (Kabluchko et al., 2009) and Schalther's (Schlather, 2002) models and compare the accuracy of our approach with the standard pairwise likelihood technique. We show that CNN architectures are fast to fit and can provide less biased and less variable predictions than the typical approaches provide. Then, using 10-day blocks of reanalysis temperature data in the midwestern United States, we illustrate the performance and benefits of the proposed framework in terms of accuracy and computational savings over the standard approach. Another advantage of our approach over pairwise likelihood estimation methods is that no detailed knowledge in the field (e.g., extremes) is necessary since the target statistical computation is done inside the deep NN. Therefore, our deep NN-based framework does not require any assumptions about the data generation process.
Instead, it is designed to have a flexible and computationally appealing form, with the network's weights learned by training, which allows accurate representation of a large class of functions.
Although our focus is on spatial extreme modeling, the new approach can be used for other inference problems as long as simulations from the probability model are possible.
The remainder of this paper is organized as follows. Section 2 provides an overview of the proposed approach and briefly introduces CNNs. In Section 3 we recall the definition of maxstable processes and inference based on pairwise likelihoods and conduct a simulation study to assess the performance of the estimators. In Section 4 we illustrate our method on a dataset of temperature minima, and we conclude with a discussion in Section 5.

Methodology
This section outlines our new framework for estimating parameters in statistical models by using deep NNs and gives a brief overview of the type of networks used in our examples.

Parameter estimation framework
In what follows, we will use the terminology parameters to refer to the parameters in statistical models of interest; the term NN weights denotes the parameters of NNs. As in classical statistical inference, our goal is to use observed data to infer the underlying characteristics of a stochastic process by estimating parameters of an assumed probability model. Suppose Y is a random variable with values in a measurable space and probability density function p with respect to a reference measure µ (often the Lebesgue measure). Consider that the probability density function p(·; θ) depends on a finite-dimensional parameter set θ ∈ Θ ⊂ R K , where Θ denotes all possible values of θ and K ∈ N. The statistical problem is to recover the unknown θ given y = (y 1 , . . . , y n ) , where y i are independent copies from p(·; θ). Evaluating p(·; θ) at observed data samples y gives a real-valued function p(y; θ), which is called the likelihood function. In most cases, this likelihood function is explicitly known and can be evaluated either analytically or numerically for pairs (y; θ). To estimate θ, one commonly uses the principle of maximum likelihood, which assumes that the most reasonable values are those for which the density of the observed sample is largest. The valueθ =θ(y) that maximizes the log-likelihood function is called the maximum likelihood estimate (MLE): θ = arg max θ∈Θ p(y; θ). (1) Often, however, no closed-form solution to the maximization problem is known or available, and the MLE can be found only via numerical optimization. In even more extreme cases, such as many non-Gaussian models for dependencies, writing down the full likelihood is impossible, and advanced numerical optimization strategies are impractical even at moderate dimensions.
We propose to tackle this problem using deep NNs, a type of machine learning that exploits a connected multilayer set of models to learn complex input-output maps. The idea is to use samples of parameters and data pairs to learn the function p(·; θ). In other words, our goal is to find a useful approximation to p(·; θ) that underlies the predictive relationship between θ and y. Our estimation technique is a deep NN taking y as the independent input data and θ as the dependent output data F : y → θ;θ = arg min d{y, F w (θ)}, where F w is a function that depends on a set of l 1 deep NN weights and l 2 NN biases, in other words, w = (w 1 , . . . w l 1 , b 1 , . . . , b l 2 ) . These weights are determined by training, whereas d in (2) represents a measure of distance between y and F w (θ). The function F w is analogous to the arg max operation, but it approximates the relation between y and θ by means of a deep NN instead of maximizing the likelihood function. In our approach we seek the function F for predicting θ given values of the input y. This problem requires a criterion for choosing F, namely, a loss function that penalizes errors in prediction. For regression problems, the mean squared error (MSE) loss is the most common and convenient choice (Friedman et al., 2001) and is given by An optimization algorithm uses the gradients, usually calculated via automatic differentiation (Paszke et al., 2017), of the loss function in (3) applied to the output of the deep NN and data.
The idea is that the deep NN finds the set of weights w such that the fit gets as close to the observed points as possible based on realizations of the output, where proximity is measured by the MSE(w). For a linear model, a simple closed-form solution to the minimization problem of (3) exists; otherwise, the solution often requires iterative methods. Several numerical optimizers built around batch gradient descent methods are implemented to perform this task (e.g., see Kingma and Ba (2014)).
The problem of optimizing F w involves generating training data for the deep NN. For this purpose we produce training samples of parameters and data pairs (θ train configurations covering the parameter domain of interest Θ. Different approaches can be used to find such configurations. One possible method is to obtain crude estimates of θ from the observed data y (e.g., using approximate maximum likelihood methods) and simulate (θ train in a large enough neighborhood of those. Then, after the parameter space Θ has been defined and samples of (θ train j ) J j=1 have been generated, these can be used to simulate corresponding data samples (y train j ) J j=1 . The number of training samples depends on the difficulty of the problem and directly affects the quality of the estimator (Juba and Le, 2019). However, we find that appropriate transformations of the input and output greatly reduce the number of samples needed for decent accuracy. Once the deep NN has been trained, it can be used to accomplish our goal of performing statistical inference, namely, returning the probability model's parameters based on new observed data.

Introduction to convolutional neural networks
Here we give details on the specific form assumed for F w in (2) responsible for learning the mapping between simulated data (input) and the statistical model parameters (output) of the deep NN. The choice of the deep neural network will likely depend on the type of input and output. In the following sections we analyze spatial 2D data; therefore, we focus on CNNs, a class of deep NNs most commonly applied to visual imagery that accounts for local dependencies in signals. We consider F w to be a composition of m nonlinear functions such that , where each f (·) is a function of at least one hidden layer. The general structure of a CNN alternates convolution layers followed by pooling layers, with the last layers being fully connected (dense). The convolution operation enables performing weighted averaging of inputs such that the network learns filters that are activated when specific spatial patterns in the input are detected. Formally, given input images y and a kernel k, the discrete convolution operator is given by where s 1 and s 2 are pixels of the image. The sums are finite in practice since the kernel K has compact support and depends on the number of pixels. Each kernel K is a matrix of trainable weights (also known as filters) applied to small regions of the input image, and the same kernel is translated across the dimensions of the image. Several independently trained 7 kernels can be applied when multiple filters are used in a convolutional layer. Depending on the values of the kernel weights, one can get different properties associated with the image, such as edges of objects. We refer to Pinaya et al. (2020) for more information on CNNs and to https://tensorflow.org/api_docs for instructions on the TensorFlow implementation of convolutional layers.

Proof of concept for max-stable processes
Max-stable distributions are commonly used for studying extreme events recorded in space and time. As a proof of concept, we exemplify the methodology from Section 2 by estimating parameters of max-stable processes. These processes are well known for having full likelihoods that are computationally intractable even in moderate dimensions. However, simulation is manageable.
With this example, we aim to demonstrate that our method can accurately recover the parameters of a model with intractable likelihood by learning the mapping from raw data. We use here the algorithm provided by Schlather (2002) for simulation of max-stable models implemented in the SpatialExtremes R-package.
In Section 3.1 we briefly outline max-stable distributions and processes and the most commonly used approach for estimating parameters in such models. Sections 3.2 and 3.3 outline the parameter estimation setup and the CNN used in the estimation procedure, respectively. In Sections 3.4 and 3.5 the performance of the CNN estimator is assessed and compared with a benchmark through a simulation study using two popular max-stable models.

Definition and the pairwise likelihood approach
Consider a sequence of independent and identically distributed (i.i.d.) stochastic processes . ., where s ∈ S ⊂ R d is a spatial site. A max-stable process Z is the limit process of normalized pointwise block-maximum process (with block size m), provided that sequences of 8 functions a m (s) > 0 and b m (s) exist: in other words, that pointwise maxima of independent copies of Z(s) remain in the same family up to a location and scale functions. Each marginal distribution of Z(s) is a generalized extreme value (GEV) distribution with location, scale, and shape parameters that potentially depend on the spatial locations (De Haan and Ferreira, 2007).
We now introduce a representation that will serve as the basis for constructing a parametric max-stable process in the following sections. Let {ξ i } i≥1 be the points of a Poisson process on (0, ∞) with intensity dΛ(ξ) = ξ −2 dξ, and let W 1 (s), W 2 (s), . . . be independent copies of a nonnegative stochastic process such that W (s) ≥ 0 with mean equal to one. The processes W i and the points of the Poisson process {ξ i } i≥1 are assumed to be independent. Then is a max-stable process with unit Fréchet margins: P{Z(s ≤ z)} = exp(−1/z), z > 0. Suitable choices for W (·) yield to different max-stable processes (Schlather, 2002). In the next sections we will consider two different choices: Brown-Resnick (Kabluchko et al., 2009) and Schalther's model (Schlather, 2002).
• The Brown-Resnick model arises when (5), where i (s) are independent copies of a standard stationary Gaussian process such that (0) = 0 almost surely and with semivariogram γ(h) = ( h /λ) ν , where h is the spatial separation, λ > 0 is the range, and ν ∈ (0, 2] is a smoothness parameter. Brown-Resnick processes are popular in practice because of their flexibility compared with other choices and the ability to generalize several stationary max-stable models, such as the geometric Gaussian model . For a precise justification of its practical performance, see Thibaud and Opitz (2015).
• The second model is the characterization proposed in Schlather (2002), with W i (s) = √ 2πmax{0, i (s)}, where 1 (s), 2 (s), . . . are independent copies of a stationary Gaussian process with unit variance and correlation function ρ that is a powered exponential: Schlather (2002) showed that one can obtain exact simulations of max-stable models on a finite sampling region if the support of the random Poisson process ξ (see (5)) is either included in a ball b(o; r) or is uniformly bounded by a constant C. For Brown-Resnick and Schlather's models, these conditions are not satisfied since Gaussian processes are not uniformly bounded. However, is sufficiently small, the simulation procedure is still accurate.
It follows from (5) that the joint cumulative distribution of Z(s) at a finite collection of sites where the exponent function V (z 1 , . . . , z D ) = E[max{W (s 1 )/z 1 , . . . , W (s D )/z D }], satisfies homogeneity and marginal constraints (De Haan et al., 1984). By differentiating the distribution (6) with respect to its variables z 1 , . . . , z D , one can deduce the corresponding density function, or the full likelihood as where P denotes the collection of all partitions π = {π 1 , . . . , π L } of {1, . . . , D} and V π l denotes the partial derivative of the function V with respect to the variables indexed by the set π l .
The sum in Equation (7) is taken over the set of all possible partitions, which equals the Bell number of order D. This leads to an explosion of terms even for a moderate D, and the full likelihood quickly becomes intractable. Castruccio et al. (2016) concluded that with modern computer technologies, full likelihood inference is still limited to D ≤ 12. In the simple bivariate case, however, the density is computationally tractable, which explains the common practice of adopting the pairwise likelihood in place of the full likelihood (see, e.g., Padoan et al. (2010); Davis et al. (2013); Shang et al. (2015)). The corresponding log weighted pairwise likelihood for model (6) is given by where φ ∈ Φ ⊂ R p is the vector of unknown parameters and α j 1 ,j 2 ≥ 0 denotes the likelihood weight assigned to the pair {j 1 , j 2 }. We suppress the dependence of V 1 , V 2 , and V 12 on φ for notation simplicity. The weights are often set to one only for pairs of locations up to a certain distance: h = s 1 − s 2 < δ, where δ > 0 is an appropriate cut-off distance chosen, for instance, via cross-validation (Ribatet, 2008). The block maximum recorded at the jth station is denoted by z j , where P = {(j 1 , j 2 ) : 1 ≤ j 1 < j 2 ≤ D}. For multiple independent replicates of the z j 's, the maximum pairwise likelihood estimatorφ in (8) is strongly consistent and asymptotically Gaussian. The asymptotic variance for such misspecified likelihood estimators is of the sandwich form (Padoan et al., 2010). However, calculating this sandwich variance is time-consuming, since it requires working with four points at a time. The variability ofφ may also be assessed by using block bootstrap techniques.

Proposed parameter estimate setup
For simplicity, we describe below our proposed estimation approach for simulated data on uniform margins, such that we estimate only the dependence parameters. We simulate I independent replicates (θ i , y i ) I i=1 (also known as testing data), where θ i = (λ test i , ν test

Pairwise-likelihood estimatorθ
The models are fitted to data by using the pairwise likelihood (see (8)) through the R-function fitmaxstab. This function performs optmization using another R-function, optim, and we set the method to L-BFGS-B. After a small cross-validation study, we find that including pairs that are more than three units apart worsens the estimation, besides increasing the fitting time. Therefore, we combine only likelihood contributions from pairs that are at most three units apart with equal weights α i 1 ,i 2 = 1, and the remaining weights are set to zero. We initialize the parameters by giving the optimizer multiple random starting pairs around the actual values. Then we run the full optimization from the five pairs with the highest pairwise likelihood among the random pairs as the starting point. We show the results using the pair with the pairwise likelihood maximizer as initial values.

CNN architecture
We use a sequential CNN taking simulated images as input and mapping it onto two scalar values. For each combination, which we call a scenario, we simulate 50 i.i.d. realizations of the respective Brown-Resnick process, which will be used to assess the uncertainty of the estimators. Examples of realizations covering the simulated parameters are illustrated in Figure 1 Figure 2 (a) and might be explained by the low efficiency of using only pairs of variables.
Moreover, we consider one replicate for each test data, and increasing this number could make the estimation more stable for the pairwise. However, investigating higher-order composite likelihood methods is outside this paper's scope, and we refer the reader to Genton et al. (2011) for details on this approach.  Figure 1 (b) shows examples of simulations from the Schlather's model. This model produces data that are more noisy and presents a more challenging scenario when compared with 15 the Brown-Resnick process presented in the preceding section (see Figure 1 (a)).

Results for Schlather's model
Each scatterplot in Figure 2 (b) displays estimated values of log(λ) versus log{ν/(2 − ν)} for 50 replicates and a different scenario where the estimation is performed by using either the CNN (green) or pairwise likelihood (red). Compared with the analyses for Brown-Resnick processes (see Figure 2 (a)), both methods produce estimates with less variability, and the bias from the pairwise is less prominent here. Exceptions are scenarios with small ranges for which the pairwise likelihood underestimates the range in a few datasets (see the first column of Figure 2 (b)). The The improved pattern from the CNN over the pairwise likelihood estimator is confirmed in Table 2, with RMSE and MAE values that are around 60% to 80% smaller. The only case where the pairwise likelihood beats the CNN is in terms of mean biases for the smoothness parameter, which is primarily seen for small ranges, namely, when the spatial dependence is relatively weak (see Figure 2 (b)). In such cases the CNN estimates of smoothness seem to have a distribution that is not strongly dependent on the actual smoothness. The computations were carried out on  Table 2: RMSE, MAE, and mean bias from the Brown-Resnick and Schlather's models using the CNN and the pairwise likelihood approaches. The two numbers in each column for the first three rows represent scores for estimating range and smoothness, respectively.

Application to temperature data
In this section we discuss an application of our method to a surface temperature reanalysis dataset, which we describe in Section 4.1. Section 4.2 outlines the estimation framework and highlights the additional steps needed to estimate parameters from real data compared with the simulation study in Section 3. In Section 4.3 we analyze the results.

NLDAS data and preprocessing
Reanalysis temperature data are extracted from the North American Land Data Assimilation System (NLDAS-2) (Mitchell et al., 2004;Xia et al., 2012b,a) and are freely available at https: //ldas.gsfc.nasa.gov/nldas/v2/models. NLDAS data is quality-controlled and available at a grid spacing of 12 km and the hourly temporal resolution. Hourly data are extracted between 1991 and 2019.
We select locations in the midwestern United States on a 107 × 55 grid-scale and perform spatial downsampling by selecting 25 equally spaced points in the latitude and longitude directions so that images are the same size as in the simulation study (see Section 3). To decrease the effect of seasonality, we restrict the data to 60 days over April-May, although there is still large variability within this period. We consider daily minima, and the models are fit to the negative of this quantity over six 10-day blocks, resulting in 174 datasets of size 25 × 25. Figure 3 shows images of minimum temperature data from the six 10-day data blocks over April-May (rows 1-6) at four equally spaced years between the 29 years available (rows 1-4). As expected, temperature minimas are lower at the beginning of April.
Although the spectral representation of max-stable processes in (5) relies on Fréchet margins, in practice the observations are rarely a unit Fréchet distribution. We follow the usual approach of fitting a GEV distribution before transforming the data to the unit Fréchet scale. We model marginal distributions with GEVs using componentwise minima over six 10-day blocks at each site separately. We estimate a single scale and shape parameter for the considered April-May period. In contrast, the location parameter is specific to each of the six 10-day blocks to better meet the identically distributed assumption. Quantile-quantile plots (not shown) suggest that a GEV reasonably approximates marginal distributions. Images of the estimated locations, scale, and shape are displayed in Figure 5. Location parameters for each of the 10-day blocks gradually increase as the weather gets warmer. Shape parameters are primarily negative, indicating a bounded upper tail. Negative shape parameter estimates are commonly found when fitting GEV distributions of temperature extremes; see, for example, Kharin and Zwiers (2000); Huang et al. (2016). Exceptions are close to the Great Lakes (see mid-upper right corner), where the positive shape parameter implies a heavy tail. We then use these estimated parameters from the GEV to transform the data to the unit Fréchet scale.

Implementation details of the estimation framework
In Section 3 we generated the training set for the CNN based on the data and parameters from the testing set, which were simulated and were therefore known. In this section the parameters of the testing set are unknown. Consequently, we first need to define a range of parameter values that reflect the assumed max-stable distribution of the unit Fréchet transformed data.
The diagram in Figure 4 illustrates the several steps needed to obtain parameter estimates of max-stable processes using the temperature data. Notice that in the simulated setup described in Section 3, only steps 4, 5, 6, and 7 were needed. Different approaches can be taken to choose the region to sample the parameters from. We present one possible method, where the permissible parameter region from fitting a pairwise likelihood to the data is expanded by simulating data in a broader grid with uniform probabilities.
The first task is to fit Brown-Resnick models to each of the 174 transformed data images separately by using the pairwise likelihood estimator (8), for which pairs of locations are at most three units apart. Figure 6 displays histograms of the pairwise likelihood (red) estimated range (left) and smoothness (right). Whereas the range parameter varies from approximately 12% to 50% of the total size of the domain [0, 20], the estimated smoothness is often relatively large. Next, we train the CNN with simulated training data of size 2,000 generated by sampling both parameters independently from uniform distributions. The lower and upper bounds of the uniform are obtained by perturbing the pairwise likelihood estimates based on their standard deviations: where a ∧ b denotes the maxmium between a and b,Θ PL = (Θ PL 1 , . . . ,Θ PL 174 ) ,Θ PL min = min(Θ PL 1 , . . . ,Θ PL 174 ) andΘ PL max = max(Θ PL 1 , . . . ,Θ PL 174 ) for Θ = λ, ν. Using the range and smoothness values generated from (9), we simulate Brown-Resnick processes with spatial dimension D = 25 2 . We then train a CNN model using the logarithm transformed data from these simulations as inputs and the variance stabilizing transformation on the outputs (see Section 3.3 for details). The CNN model configuration is summarized in Table 3. After the model is trained, we predict the output based on input images of logarithm transformed data. Histograms of estimated parameters at their original scale from the CNN are displayed in Figure 6 (green).
Compared with the pairwise likelihood estimates (red), the range parameter (left) estimated from the CNN is usually larger and more uniform across datasets. The smoothness (right) is closer to normally distributed with a higher probability of values close to one.

Results
We detail here the results from the CNN and pairwise likelihood approaches with the data and implementation described in the preceding sections. For clarity, the results are presented on the Fréchet transformed data scale.
The left panel of Figure 7 compares the empirical F-madogram estimates (black dots) of the extremal coefficients between pairs of sites, plotted against the spatial distance h = s i − s j , with respect to their pairwise likelihood (red curve) and CNN (green curve) model counterparts.
The empirical F-madogram is calculated from 174 datasets at the 25 2 locations, and estimates of the binned F-madogram for the two models are obtained by using 100 bins. The pairwise likelihood fit appears to underestimate extremal dependence at distances greater than five (red curve is generally above the back dots). In contrast, the CNN provides a more realistic estimation of spatial dependence, matching empirical and estimated extremal coefficients.
The uncertainty from both model fits can be assessed from the right panels of Figure  only the CNN can capture pairs with high extremal dependencies (e.g., θ(h) < 1.5), a feature that is more evident at larger spatial separation (the green area is below the red in the rightmost plot).
The quantile-quantile plots in Figure 8 compare the observed and predicted minima (left),

Discussion
In this work we have proposed a new approach to estimate parameters in statistical models based on deep neural networks. As a proof of concept, we have tested the methodology based on the highly challenging problem of inference for max-stable distributions and processes. We illustrated the benefit of using deep NN models over the classical pairwise likelihood approach through simulated and real-world data. Results from our simulation study for the popular Brown-Resnick and Schlather's models showed better performances (more negligible bias and variance) of our method compared with the current approach with speedups in computations.
Unlike pairwise likelihood, the proposed deep learning approach relies on simulations of the full model and does not involve approximations of spatial dependence.
We follow the paradigm that if it is possible to simulate from a model, one can perform inference. Our approach uses data from simulations as input to a deep NN and learns to predict statistical parameters. Similar to the well-known ABC, the idea is that the deep NN becomes a stand-in for compressing data and recognizing model parameters. However, a key difference from ABC is that our method automatically maps data and parameters inside the deep NN, eliminating the problem of choosing summary statistics and metrics to compare data. Another difference is that in ABC the entire estimation procedure must be run again from scratch for each individual dataset, whereas in the proposed framework the estimation is performed upfront in a single training phase, followed by a cheaper prediction phase able to handle multiple testing sets. Moreover, our new approach benefits from recent advances in algorithms for conventional deep NNs, which can learn features automatically from highly complex images while being fast and easy to apply. Problems with our approach will emerge when there are many parameters in the model, up to a point where generating training data for high-dimensional parameter spaces becomes practically intractable.
In this work we have focused on modeling multivariate extremes and followed a proof of concept approach. We can use the same idea presented here for inferring parameters in other statistical models. In particular, non-Gaussian models for dependencies tend to be hard to estimate by using classical approaches. Extensions for future work include nonhomogeneous Poisson processes (Moller and Waagepetersen, 2003), epidemiology models (Lawson, 2008), discrete stochastic population dynamics models (Wood, 2010), and random set models.
In summary, we show that deep NN models are easy to construct, are fast to fit, and show promising results for performing statistical inference, especially when the usual techniques fail.