Cokriging (Goovaerts 1998), a multivariate extension of univariate kriging, is widely used for analysis and prediction of multivariate spatial fields in multivariate spatial statistics. Prediction (with uncertainty quantification) at new unobserved sites using the information gained from the observed locations is one of the common objectives of this strategy. For Gaussian random fields with known covariance structure, cokriging is the best linear unbiased predictor. However, in reality, this requirement of Gaussianity and properly stated covariance is rarely met. Because of this, modeling heavy-tailed or skewed data with complicated covariance structure requires a more adaptable prediction methodology.
Various approaches have been proposed to model non-Gaussian spatial data such as scale mixing Gaussian random fields (Fonseca and Steel 2011), multiple indicator kriging (Journal and Alabert 1989), skew-Gaussian processes (Zhang and El-Shaarawi 2010), copula-based multiple indicator kriging (Agarwal, Sun, and Wang 2021) and Bayesian nonparametrics (Reich and Fuentes 2015). To address the nonstationary behavior of the random fields over a large region, nonstationary Matérn covariance models (Cressie and Huang 1999; Fuentes 2001; Nychka, Wikle, and Royle 2002; Paciorek and Schervish 2003) have been introduced. Trans-Gaussian random fields CONTACT Pratik Nag pnag@uow.edu.au School of Mathematics and Applied Statistics, University of Wollongong, Wollongong NSW 2500, Australia. Supplementary materials for this article are available online. Please go to www.tandfonline.com/r/TECH. (Cressie 1993; De Oliveira, Fokianos, and Kedem 2002) find nonlinear transformations which enables application of Gaussian processes on the transformed data. However, marginal transformation to normality may not confirm joint normality at multiple locations and may change important properties of the variable (Changyong et al. 2014). One common drawback of many of the approaches is that there is no straight forward implementation of these methods in the bivariate setting and also they may not be optimal for more general spatial datasets. Another issue is that most of these models are based on Gaussian processes which uses cholesky decomposition of the N × N covariance matrix which requires O(N 3 ) time and O(N 2 ) memory complexity where N is the number of spatial locations.
Recently, Deep neural network (DNN) based algorithms have proved to be the most powerful prediction methodologies in the field of computer vision and natural language processing (LeCun, Bengio, and Hinton 2015). DNNs can handle more complex functions and in theory they can be applied to approximate any function which is appealing for modeling large-scale nonstationary and non-Gaussian spatial processes. DNNs are also computationally efficient and thus can be applied to large datasets. The computation time for training can be massively reduced by using GPUs (Najafabadi et al. 2015). Due to this broad applicability of the neural networks, several researchers are trying to incorporate DNNs in the spatial problems (Wikle and Zammit-Mangion 2023). Cracknell and Reading (2014) included spatial coordinates as features for DNNs. Wang, Guan, and Reich (2019) proposed a nearest neighbor neural network approach for Geostatistical modeling. Zammit-Mangion et al. (2021) used neural networks to estimate the warping function which transforms the spatial domain to fit stationary and isotropic covariance structure. Convolutional neural networks (CNNs) (Krizhevsky, Sutskever, and Hinton 2012) can capture the spatial dependence successfully, but require a huge amount of data on regular grid for model training. However, in environmental statistics scenarios one of the main objective is to give spatial interpolation at unobserved locations for irregular grid datasets which is infeasible for CNN modeling. Most of these DNN based methods are developed for univariate data and only concentrate on point predictions and ignore the prediction interval estimations. Recently, there has been few works which proposes density function estimation of the predictive distribution using neural networks, for example Li, Bondell, and Reich (2019) proposed a discretized density function approach and predict the predictive distribution probabilities by training a neural network classifier. Neal (2012) and Posch, Steinbrener, and Pilz (2019) applied Bayesian inference methodologies to neural networks to predict uncertainties via posterior distributions. But none of these methods are directly applicable to spatial applications.
To address these drawbacks Chen et al. ( 2022) introduced a spatially dependent deep neural network structure called DeepKriging, for univariate spatial prediction. They used basis functions to capture spatial dependence and use them as features to fit the model. They also provided an approach for uncertainty quantification by an histogram based approach. However, the histograms require thresholds which are mostly deterministic and choice of the thresholds effects the results drastically.
The objective of this study is to establish a nonparametric statistical framework for conducting bivariate spatial prediction and estimating prediction uncertainties. Our proposed approach, a bivariate extension of DeepKriging (BDK), aims to address these limitations. We introduce a spatially dependent neural network by incorporating an embedding layer of spatial coordinates using basis functions. Additionally, we propose a bootstrap and ensemble neural network-based prediction interval for prediction uncertainty. BDK is nonparametric, making it more versatile and applicable to non-Gaussian, nonstationary, and even categorical prediction problems. Through simulations, we demonstrate that this approach yields comparable results to cokriging when the underlying process is Gaussian, and it surpasses traditional statistical methodologies in non-Gaussian and nonstationary scenarios. Moreover, our proposed approach significantly reduces computation time compared to traditional methods.
While architecturally BDK is a straightforward extension of univariate DeepKriging, integrating this methodology with traditional statistical methods remains a challenging task. In this article, we illustrate that under certain conditions, the linear model of coregionalization (LMC) represents a special case of BDK. Additionally, we introduce a novel methodology for computing prediction intervals using bootstrap and ensemble neural networks. Additionally, the execution of BDK is computationally twenty times faster than traditional statistical approaches, effectively overcoming the bottleneck associated with applying this method to large datasets.
The remainder of our article is structured as follows. Section 2 outlines the proposed methodology. Section 3 compares the performance of the proposed method with traditional approaches. Lastly, Section 4 applies the BDK method to Saudi Arabian wind data. More detailed information on the exploratory data analysis of the wind dataset and further results on simulation studies can be found in the supplementary materials.
be a bivariate spatial process, and {Z(s 1 ), Z(s 2 ), . . . , Z(s N )} be the realization of the process at N spatial locations, where
A classical spatial model assumes
where (s) = { 1 (s), 2 (s)} T , known as the nugget effect, is a bivariate white-noise process with variances σ 2 1 (s) and σ 2 2 (s), respectively.
Given observations Z(s), one of the main objectives of spatial prediction is to find the best predictor of the true process Y(s 0 ) at the unobserved location s 0 . The optimal predictor with parameter set θ can be defined by minimizing the expected value of the loss function, that is,
where
Ŷopt (s 0 |Z N ) is the optimal predictor given Z N = {Z(s 1 ), Z(s 2 ), . . . , Z(s N )}. The function R(•) is the risk function, which is the approximation of the loss L.
Bivariate kriging, also known as Cokriging, is the best linear unbiased predictor (Stein and Corsten 1991;Goovaerts 1998) for bivariate spatial prediction. The spatial process Y(s) is typically modeled as
where γ (s) = {γ 1 (s), γ 2 (s)} is a spatially-dependent zeromean bivariate random process with cross-covariance function
The mean structure μ(s) can be modeled as μ(s) = X(s) T β, where
1 vector of coefficients corresponding to the covariates. Let Z vec = {Z 1 (s 1 ), . . . , Z 1 (s N ), Z 2 (s 1 ), . . . , Z 2 (s N )} T be a 2N × 1 vector, C is the crosscovariance matrix of order 2N × 2N, and
, the kriging prediction at an unobserved location s 0 is
where β = VX T C -1 Z vec , the generalized least square estimator of β, X 0 = X(s 0 ) is the 2 × (p 1 + p 2 ) matrix of covariates at location s 0 , C 0 is the 2N × 1 cross-covariance matrix between the response at location s 0 and the observed locations. The choice of cross-covariance functions in cokriging prediction plays a vital role in capturing the spatial dependence among locations and also among variables for multivariate data.
Among existing models, the multivariate Matérn crosscovariance function (Apanasovich, Genton, and Sun 2012) is one of the most popular classes of Matérn covariance function. The parsimonious Matérn covariance function can be defined as
where d(s i , s j ) = s is j is the Euclidean distance between the locations, M is the Matérn correlation function (Cressie and Huang 1999), ν is the smoothness, and α is the range parameter. Constraints on these parameters have been derived in Apanasovich, Genton, and Sun (2012).
Another example of a cross-covariance function is the linear model of coregionalization (LMC); (Genton and Kleiber 2015). In this approach, a multivariate random field is represented as a linear combination of r independent random fields with potentially different spatial correlation functions. For r = 2, we can write the process γ (s) in terms of the latent processes as γ (s) = AU(s), where A is the 2 × 2 matrix of weights and U(s) = {U 1 (s), U 2 (s)} T , cov(U 1 , U 2 ) = 0. Following this architecture, the LMC covariance function of a bivariate random field can be written as
i,j=1 is a 2 × r full-rank matrix. Although cokriging is the best linear unbiased predictor, it suffers from several limitations. First, maximum likelihood estimation (MLE) for Gaussian processes requires Cholesky factorization of the covariance matrix, which is computationally expensive, making it infeasible to implement traditional cokriging for large datasets. Moreover, many real-life datasets are non-Gaussian, and cokriging is not an optimal predictor in those scenarios.
DeepKriging, proposed by Chen et al. (2022), addresses these issues for univariate processes by rephrasing the problem as regression and introducing a spatially dependent neural network structure for spatial prediction. They have used the radial basis functions to embed the spatial locations into a vector of weights and passed them to the DNNs as inputs allowing them to capture the spatial dependence of the process. Similar to Chen et al.
(2022), we propose to use spatial basis functions for the bivariate processes to construct the embedding layer instead of passing the coordinates of the observations directly to the NN. This idea is motivated by the multivariate Karhunen-Loéve expansion (Daw et al. 2022) of a multivariate spatial random field, that is, a bivariate spatial process γ (s) can be represented as
where w b,u 's are independent random variables and φ b,u (s)'s are the pairwise orthonormal basis functions corresponding to variable u, u = 1, 2. Therefore, a bivariate process can be approximated by
where the weights w b,u can be estimated by minimizing the risk function R. Hence, the spatial problem can be viewed as multioutput linear regression (Borchani et al. 2015) with these basis functions φ b,u as covariates.
Among many existing basis functions, such as spline basis functions (Wahba 1990), wavelet basis functions (Vidakovic 2009), and radial basis functions (Hastie, Tibshirani, and Friedman 2001), we have chosen the multi-resolution compactly supported Wendland radial basis function (Nychka et al. 2015) defined via w(d) = (1-d) 6 3 (35d 2 + 18d + 3)1{0 ≤ d ≤ 1}. The spatial basis functions are then defined as φ i (s) = w( su i /θ ) with θ as the bandwidth parameter and anchor points (spatial locations) {u 1 , u 2 , . . ., u K }. We take {u 1 , u 2 , . . ., u K } to be a square grid of locations covering the spatial domain and θ to be 2.5 times the distance between adjacent anchor points. By selecting various K values, one can produce distinct sets of basis functions. This will enable us to record both the long-and short-range spatial dependencies. We must combine these basis functions into a single vector in order to send them to the neural network. See Chen et al. (2022) for detailed discussions on the basis function embedding. Therefore, these co-ordinates embedding can now be used as an input to the DNNs along with covariates X vec (s) = {X 1 (s) T , X 2 (s) T } T which will better capture the spatial features of the data. Note that under our assumption of basis functions, φ .,1 = φ .,2 for location s, hence, we can reduce our basis function set to φ(s
We have used a multi-output deep neural network structure to build the bivariate DeepKriging framework. We define X φ (s) as the vector of inputs containing the embedded vector of basis functions φ(s) and the covariates X vec (s). Then X φ (s) = (φ(s) T , X vec (s) T ) T . Taking X φ (s) as the inputs to a neural network with L layers and M l nodes in layer l = 1, . . ., L DNN can be specified as,
...
where Figure 1 gives a graphical illustration of the spatially dependent neural network structure.
The Bivariate DeepKriging (BDK) is a multivariate extension of the univariate DeepKriging. Here we propose to approximate the optimal predictor Ŷopt (s 0 ) with the help of deep neural network structure, that is; Ŷopt (s 0 ) ≈ f NN (X φ (s 0 )) We define the optimal neural network based predictor as
where
. and
Here σ 2 u is unknown and can be estimated through the sample variance of the uth variable. The final predictor of the neural network would be the bivariate prediction vector f opt NN (X φ (s 0 )) where Theorem. For a co-located bivariate spatial process, assuming that the latent variables are constructed with the same sets of basis functions, the Linear Model of Co-regionalization (LMC) represents a special case of BDK.
Proof. For the bivariate case, the construction of LMC (Genton and Kleiber 2015) provides insights into such a link. If Y(s) is a zero-mean process then it can simply be represented as the underlying bivariate spatial process γ (s) which again can be written as
where A = a 11 a 12 a 21 a 22 and U(s
Now, using the Karhunen-Loève Theorem (Adler 2010), we can approximate them as
We then write
. Now, as both the random variables U 1 (s) and U 2 (s) are colocated, we can choose the same set of basis functions for them, resulting in
which is a multioutput regression with φ(s) as covariates and γ (s) as the response, which is same as ( 4) with L = 1, b 1 = 0, and a linear activation function.
The utilization of neural networks for uncertainty quantification has garnered significant attention in recent years. Various prominent techniques have emerged (Khosravi et al. 2011; Nourani, Paknezhad, and Tanaka 2021) to enhance the estimation of prediction intervals. For instance, Cannon and Whitfield (2002) and Jeong and Kim (2005) employed ensemble DNNs for forming prediction intervals, while Srivastav, Sudheer, and Chaubey (2007) and Boucher, Perreault, and Anctil (2009) introduced bootstrap-based uncertainty quantification for DNNs. Additionally, Kasiviswanathan and Sudheer (2016) proposed a Bayesian framework for uncertainty quantification with DNNs.
The univariate DeepKriging model computed prediction intervals using a histogram-based approximation to the predictive distribution. Nonetheless, this method exhibits several limitations. First, the accuracy of the prediction intervals heavily relies on the choice of bins for the histogram. The selection of cut points in the histogram, which are randomly chosen from a uniform distribution, may impact the approach's performance. Moreover, the computational time increases when training a large number of ensemble models, and the ensemble approach suggested by this model might fail to capture the variation in the data. In this section, we propose an uncertainty quantification method to address these limitations.
For the bivariate spatial prediction problem, using ( 1) and ( 4), the prediction at an unobserved location s 0 can be expressed as Ẑ(s 0 ) = f opt NN (X φ (s 0 )) + (s 0 ). By employing ensembles, we can generate B replications of Ẑ(s 0 ) at s 0 . Consequently, the prediction can be articulated as
Employing the multidimensional Central Limit Theorem (Borovkov and Rogozin 1965) yields Ẑ(s 0 ) B ∼ N 2 (Z(s 0 ), 1 B Z(s 0 ) ). The prediction interval is thus given as
where, σ 2 (Z u (s 0 )) is the variance term associated with Z u (s 0 ), for u = 1, 2, which can be decomposed as σ 2 (Z u (s 0 )) = var(Y u (s 0 )) + var( u (s 0 )), assuming independence between Y u (s 0 ) and u (s 0 ). From the outcomes of B ensembles, we can express var(Y u (s 0 )) as
The noise variance var( u (s 0 )) can then be represented as
An empirical approximation of the above is
Consequently, based on the above formulation, the prediction interval at the α level of significance for variable u at location s 0 can be given as,
where z (1-α/2),df represents the 1α/2 quantile of the Normal distribution.
In connection with the preceding explanation, we now provide a comprehensive definition of the algorithm for calculating the prediction interval. We partition our training dataset D into two subsets, D 1 and D 2 . From D 1 , we extract a sample D 11 . Initially, we train a neural network model with L layers using the dataset D 11 , and then proceed to draw B bootstrap samples from D 1 to train B DNN models. For these B bootstrap fits, we keep the weights fixed for the first L 0 layers, allowing only the new models to be trained based on the last L -L 0 layers. This approach is adopted to control the degrees of freedom of the t-distribution in (9). By exclusively re-training the last L -L 0 layers, we ensure that the total number of trainable parameters p in each bootstrap computation remains smaller than N. Following this procedure, we compute Ẑu (s 0 ) and Var(Y u (s 0 )) based on the training data D 1 . Since we know E u (s 0 ) = 0, we have
u (s 0 ) as defined in ( 8) is infeasible as we do not have Z u (s 0 ). Hence, we estimate r 2 u (s 0 ) through the second part of the training dataset D 2 . To do this we use the nearest neighbor approach. We obtain a set D 20 = {s k : s k ∈ D 2 } with cardinality G of locations closest to s 0 and take average of the corresponding r 2 u (s k )'s, that is,
such that s g 's are the nearest G locations to s 0 . Algorithm 1 provides a step-by-step explanation of our prediction interval construction method as discussed in previous paragraph.
Intervals Algorithm Split D into D 1 and D 2 equally. Further split D 1 into D 11 and D 12 . Train a deep neural network (DNN) of L layers on D 11 . Take B random samples {D 1 1 , D 2 1 , . . ., D B 1 } from D 1 . for i ← 1 to B do Fix the weights of the first L 0 layers of the DNN and train the last L -L 0 layers on D i 1 . Train the DNN on D i 1 and store the result. Denote it as f opt NN (X φ (s)) i . end for for location s k in D 2 do Calculate Ẑu (s k ) B (6) and Var(Y u (s k )) (7). Calculate r 2 u (s k ) (8). end for For test location s 0 , obtain the set D 20 = {s k : s k ∈ D 2 } of the nearest G locations from s 0 . Calculate Ẑu (s 0 ) B (6), var(Y u (s k )) (7), and r2 u (s 0 ) (10). Calculate the prediction interval as defined in (9).
Where u stands for the uth variable, u = 1, 2.
To better understand the computational benefits of DNN over traditional kriging we can look at the time complexity of the DNN with kriging. DeepKriging involves matrix multiplication in several layers to give us the resulting output. Whereas, kriging involves a 2N × 2N matrix inversion. Note that, time complexity of multiplication of one m × n and n × p matrix is O(mnp). So a single layer l of the neural network with minibatch (a minibatch (Hinton, Srivastava, and Swersky 2012) is a randomly selected portion of the data used for neural network training.) size of b, M l-1 input nodes and M l output nodes will have the time complexity of O(bM l-1 M l ). Faster training is made possible by defining a minibatch. Hence, a neural network with L layers will have time complexity O( L l=1 bM l-1 M l ). On the other hand, time complexity of Kriging is O(8N 3 ). Hence, for large N deep kriging with adequate number of layers and nodes is more computationally efficient than traditional kriging.
We have conducted several experiments through Gaussian, non-Gaussian and nonstationary simulations to give a comprehen-sive comparison of our proposed method with the parsimonious bivariate Matérn (CMK) as well as the linear model of coregionalization (CLK). In order to compare our results we have computed the Root Mean Square Prediction Error for each variable as the validation metric. For N locations this can be given as
Our model performed at par with CMK and CLK in the Gaussian scenario. Details of the study has been provided in the supplimentary materials. In this section we mainly concentrate on non-Gaussian and nonstationary scenarios where BDK clearly outperforms the traditional statistical techniques.
Here we consider two simulation settings. The first one is bivariate non-Gaussian with covariates. We simulated the model from the following framework
where μ(s) is the mean function and γ (s) follows a non-Gaussian distribution. We define μ(s) = {μ 1 (s), μ 2 (s)} as a complex function of covariates given as:
where the covariates x i (s), i = 1, . . . , 5 were generated from multivariate Gaussian distributions with Matérn covariance with varying parameters. γ (s) was simulated from a multivariate Gaussian distribution with parsimonious Matérn covariance (2) with parameters σ 1 = 0.7, σ 2 = 0.8, ν 1 = 0.3, ν 2 = 0.6, α 2 1 = 0.05, α 2 2 = 0.1, R = 0.5. We generated 100 replicates of the random field at 1200 irregular locations on [0, 1] × [0, 1] by transforming the Gaussian fields from the previous simulation using the Tukey-g and h transformation (Xu and Genton 2017) defined as
with g = 0.8, h = 0.5 for variable 1 and g = -0.8, h = 0.5 for variable 2.
Figure 2 illustrates the boxplot of MSPE for the two variables across 100 replicates. The superiority of BDK over other conventional approaches is evident from the boxplots. It's worth noting that due to parameter configurations and data sample paths, optimization in CMK and CLK may present challenges. In contrast, BDK is devoid of such limitations and can be readily implemented on any dataset.
The second simulation explores nonstationary datasets. We simulate nonstationary random fields using the Wendland radial basis functions as discussed in Section 2.2. We apply a nonlinear transformation to the basis functions to obtain the nonstationary fields, defined as
where parameters {a, b, c, d, e} ∼ Uniform(-2.5, 2.5).
Figure 3 displays the RMSPE over 100 replicates. It is evident that BDK significantly outperforms the other comparative methods.
Several studies (Ho et al. 2001;Zhao et al. 2008) have proposed various approaches for quantifying the quality of the prediction interval. We have chosen here two specific measures for comparison of the prediction bands, namely prediction interval coverage probability (PICP) and mean prediction interval width (MPIW). In mathematical notations
where (L u (s), U u (s)) are the lower and upper prediction bounds for variable u at location s and n is the number test samples.
In this section we evaluate our proposed approach of the prediction interval computation with the parametric prediction intervals obtained from cokriging models. We have computed the prediction interval validation metrics (11) for the comparative models using the same set of simulations as defined in Section 3.1. In all of the following accessments we have computed the 95% prediction bound for predictions in the test set.
Table 1 shows the results for all the simulation scenarios. It can be seen that for the Gaussian simulations where CMK true is optimal the BDK model gives comparable performance. For non-Gaussian and nonstationary scenarios clearly the BDK outperforms the other methods by attaining the 95% prediction with smaller width. Note that, some of the simulation scenarios While Section 3.1 demonstrates how the nonlinearity in BDK aids in capturing the nonstationarity in the data, it's equally important to showcase the improvement in prediction intervals offered by the proposed BDK model. To investigate this, we generate nonstationary data with a nonstationary mean function. The simulation details are provided in supplementary materials.
Figure 4 displays the prediction intervals for BDK and the cokriging models for the nonstationary simulation. The prediction intervals for BDK adapt spatially. Specifically, BDK offers wider prediction bounds for regions with higher variance, whereas it provides narrower prediction bounds for regions with lower variability.
Based on the same simulation setting we investigate the computational time of BDK compared to CMK with different number of observations N. We have fully optimized the model and have measured the computation time for the whole process in each case. We have also compared CMK with both of our proposed methods for point prediction as well as interval prediction.
Figure 5 shows the comparison results. As it can be seen BDK always outperforms CMK in computational time hence, it is much more scalable when the sample size increases. For example, when N = 6400, which is the largest sample size we have considered, it takes more than 23 hr (83,484 sec) to train the cokriging model, which makes it computational infeasible for larger N. However, for the same data size, BDK with prediction interval only costs 7.85 min (471 sec) without GPU acceleration and 3.58 min (215 sec) with a Tesla P100 GPU. We used the GPU from Ibex a shared computing platform at KAUST. For a more powerful GPU, the computational cost will be further reduced.
In recent years, Saudi Arabia is seeking to reduce its reliance on fossil fuels for energy demand by investing in renewable energy sources. To achieve this goal, they have introduced several milestones for their upcoming smart cities. For King Abdullah City for Atomic and Renewable Energy (KA-CARE, 2012) they are planning for 54 GW of renewable energy portfolio by 2032, of which 9 GW will come from wind power. The Saudi Vision 2030 (2016) aims at achieving 9.5 GW of renewable energy by 2023. NEOM, an acronym for "New Future" and "New Enterprise Operating Model, " is the upcoming mega-city project for Saudi Arabia which is planned to use only renewable power sources (wind and solar). For instance, wind turbines are widely used throughout the world to produce wind energy. However, installing such turbines necessitates precise measurements of the local wind direction and speed. To achieve this goal, an accurate wind speed interpolation model is crucial.
The majority of previous publications modeled wind data using geostatistical methods. Kriging was employed by Wang et al. (2020) to model wind fields. Large wind datasets were modeled for interpolation by Salvaña et al. (2021) using the ExaGeoStat software, which uses bivariate parsimonious Matérn covariance structure. However, these techniques always presumptively assume that the data are stationary and Gaussian. In this article we propose a novel methodology to model the bivariate wind speed which can then be used for spatial prediction and get high resolution maps of wind fields in any particular region. For our application we have also provided a high resolution wind interpolation for the NEOM region.
We have fitted BDK on the wind data discussed in supplementary materials. We have used three layers of basis functions with K = 100, 361, 1369, respectively. The neural network consists of L = 7 layers with M 1 = M 2 = M 3 = M 4 = 100 for first 4 layers, M 5 = 50, M 6 = 50 and M 7 = 2 nodes in the final layer. By using a sample from a uniform distribution, the model weights were initialized, and a learning rate of 0.01 was selected. The dataset consists of 506,771 locations. However, due to computational issues, the cokriging models could only handle a small subsample of the data. For this reason, we have randomly chosen 147,000 locations for training of the cokriging model. On the other hand, the BDK can easily handle the full data. Hence, we have fitted two different sets of models with deep neural networks, first with 147,000 (BDK 147000 ) and second with 450,000 (BDK 450000 ) locations and we randomly chosen 1000 locations from the rest for testing. To make the cokriging model more competitive and to reduce the computational burden, we have splitted the data into 100 sub-regions and fitted the model on each sub-region independently.
Table 2 shows that BDK outperforms CMK in MSPE.
Next we look at the 95% prediction intervals provided by BDK and CMK. From Table 3 we can see that even though the MPIW for CMK is smaller than BDK, it fails drastically containing the true values within its bounds.
We have also captured the total computation time using BDK and CMK. For CMK total computation time was 2.18 days (188,352 sec) where as for BDK it took 16.81 min (1009 sec) for point prediction and 55.61 min (3337 sec) for interval prediction.
We have given a high resolution interpolation (from 5km × 5km to 1km × 1km) for the region near NEOM, an upcoming smart city in Saudi Arabia. This downscaling can improve our understanding of the wind patterns. Sanchez Gomez and Lundquist (2020) explored how the change in wind direction affect wind turbine performance.
Hence, understanding the exact wind patterns in a high resolution field will help in wind energy setup in the area. Figure 6 gives a visual representation of the spatial interpolation carried out using BDK.
In this work, we have extended the univariate DeepKriging to bivariate scenario for spatial prediction without any parametric assumption of the underlying distribution of the data. We have also introduced a nonparametric approach for prediction interval computation for this deep learning based spatial modeling architecture. Our model is generally compatible with nonstationarity, nonlinear relationships, and non-Gaussian data.
In our implementation we have chosen radial basis functions for spatial embedding. However for long range spatial dependence and circular covariance structures choice of other basis functions such as smoothing spline basis functions (Wahba 1990), wavelet basis functions (Vidakovic 2009) can be explored.
The loss function used to BDK is a weighted MSE where weights are proportional to the variance of the variables. More complex weights can be incorporated to capture complex terrain behaviors. In our study we have not considered topographical elevation. This can be in incorporated by using relevant covariates or by creating an adaptive loss function.
An interesting avenue of future work for prediction intervals is to extend the conformal prediction methods of Mao, Martin, and Reich (2024) to the deepKriging case. This method is model free and thus would avoid the assumption that the prediction intervals are a symmetric t-distribution.