STACI: Spatio-Temporal Aleatoric Conformal Inference
Authors: Brandon R. Feng
David Keetae Park
Xihaier Luo
Arantxa Urdangarin
Shinjae Yoo
Brian Reich
Venue: NeurIPS 2025 poster ()
Keywords: Gaussian Process
Spatio-Temporal
Conformal Inference
Deep Learning
Primary area: probabilistic_methods

Abstract: Fitting Gaussian Processes (GPs) provides interpretable aleatoric uncertainty quantification for estimation of spatio-temporal fields. Spatio-temporal deep learning models, while scalable, typically assume a simplistic independent covariance matrix for the response, failing to capture the underlying correlation structure. However, spatio-temporal GPs suffer from issues of scalability and various forms of approximation bias resulting from restrictive assumptions of the covariance kernel function. We propose STACI, a novel framework consisting of a variational Bayesian neural network approximation of non-stationary spatio-temporal GP along with a novel spatio-temporal conformal inference algorithm. STACI is highly scalable, taking advantage of GPU training capabilities for neural network models, and provides statistically valid prediction intervals for uncertainty quantification. STACI outperforms competing GPs and deep methods in accurately approximating spatio-temporal processes and we show it easily scales to datasets with millions of observations.

STACI: Temporal Aleatoric Conformal Inference
Brandon R. Feng
David Keetae Park Brookhaven National Laboratory
North Carolina State University
Xihaier Luo
Arantxa Urdangarin
Brookhaven National Laboratory
University of the Basque Country
Shinjae Yoo
Brian J. Reich"
Brookhaven National Laboratory
North Carolina State University
Abstract
Fitting Gaussian Processes (GPs) provides interpretable alcatoric uncertainty quan- tification for estimation of spatio-temporal fields. Spatio-temporal deep learning models, while scalable, typically assume simplistic independent covariance matrix for the response, failing to capture the underlying correlation structure. However, spatio-temporal GPs suffer from issues of scalability and various forms of approximation bias resulting from restrictive assumptions of the covariance kernel function. We propose STACI, a novel framework consisting of variational Bayesian neural network approximation of non-stationary spatio-temporal GP along with novel ppatio-temporal conformal inference algorithm. STACI is highly scalable, taking advantage GPU training capabilities for neural network models, and provides statistically valid prediction intervals for uncertainty quantification. STACI outperforms competing GPs and deep methods in accurately approximating spatio-temporal processes and we show it easily scales datasets with millions of observations.
1 Introduction
Accurate estimation of spatio-temporal (ST) fields is complex task that has relevance in wide variety domains I 21 These domains range from, but are not limited to medical 3maging[3 4) remote sensing [S) 6) climate modeling [7 and video quality [9]. Analysis of ST data is difficult as there are often obstructions in the field causing missing data, the surface is rarely smooth across both space and time, and the volume of data requires efficient estimators to be used. These issues have resulted the need of flexible models providing both accurate prediction of surfaces and uncertainty quantification (UQ) providing interpretability of results.
GPs provide both a flexible predictive surface and precise UQ. Unfortunately, it is not scalable, with the computational cost of likelihood evaluation being cubic in number of locations. Sparse GPs [13 [4] trained from subset of inducing points, spectral GPs [15] that project the GP into a low-rank spectral domain, and nearest neighbor methods [16: 17 that assume a local dependent structure are three popular classes of approximate GPs that reduce this cost. In much of the approximate STGP literature. the covariance kernel assumed to have some known stationary and separable form to maintain
Gaussian Process (GP) regression heavily used in spatio-temporal statisties [10
12
"Corresponding author, bjreich @nesu.edu
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

However, this assumption is often simplistic and fails to hold in complex fields, such as environmental data [2] 22 23]. Non-stationary ST kernels have been derived, but they can be difficult to compute and can still be mis-specified, requiring more general
computational tractability 19; 20]
approaches to modeling non-stationarity [4
1 25
26 27
Deep learning methods have risen in popularity, providing scalable and flexible estimation of ST surfaces [28] 29 30). Implicit neural representations (INRs) are class of neural network architectures specifically developed to accurately model the non-linear surface of coordinate-based fields. They typically use periodic activation functions in multilayer perceptron (MLP) architecture map spatial (IR2) or ST coordinates (IR3) to the corresponding complex signal domain [31 32 33 34 [3]]
However while they are unparalleled in estimation flexibility, they are deterministic functions that do not inherently provide UQ and can overfit to the observed data.
Deep GP models attempt bridge the gap between neural network flexibility and the interpretability of regular GP models with easily computable measures of variance and subsequent prediction intervals [35] (3). The covariance kernel is modeled as the oulput of another GP, and this can be done to an arbitrary depth to approximate any non-stationary kernel. However, similar to regular GPs, there is trade-off of scalability as even a few layers can create computational bottlenecks. Recently, 37 developed highly scalable variant of ST deep GP, projecting separate spatial and temporal GP layers into the spectral domain before concatenating into one unified prediction output. They approximate the spectral process using probabilistic Bayesian neural network trained using backpropagation, making it highly scalable. However, as this a deep GP, the covariance parameters lose interpretability as they are used in multiple GP layers. This approximation is most similar to our work that we will now introduce.
The drawback to these approximations is the potential loss of nominal coverage of prediction intervals as the variance estimate is no longer exact. Conformal inference [38 39 40 is a model-free approach
to calculating valid prediction intervals that achieve frequentist coverage probability by assuming exchangeability of the data. Recent advancements in spatial conformal prediction E1 [2] motivate our extension to the ST setting, guaranteeing validity of our method prediction intervals.
We propose the STACI algorithm': two-stage approach consisting of combining Bayesian neural network approximation of a spectral non-stationary STGP with a ST conformal inference algorithm providing valid prediction intervals for accurate UQ. STACI provides:
1. Scalable, flexible interpolation for datasets with millions of ST locations.
2. Ability to both directly compute the underlying correlation structure and interpret estimated covariance parameters, allowing for prior choices to be user-friendly, with the ability to informative priors if desired.
3. Adaptable prediction interval lengths, clearly showing areas of high and low uncertainty while maintaining desired coverage properties and reducing the impact of both spectral and neural network approximations on UQ accuracy.
2 Preliminaries
Problem Statement. Assume the ST process Y(s, 1), indexed by location s E [0, 1]2 and time > can be decomposed as
Y(s, t) (s, 1) + 7 (s. f) c(s
for mean function μ4 mean-zero GP and noise c(s. t) iid NN 72). A common model [e.g.. [0
for the covariance of Z is the stationary Matérn function
2221
21-
0²
I(v) (V2vd) K, (V2vd).
(2)
for d² = - (t /r) and Bessel function K2. The covariance is defined by the variance 2², spatial range Ps. temporal range Pt and smoothness parameter v. For simplicity. assume the observed data has already been centered and scaled with the mean and variance estimated separately.
'Code and data: https//gitnub.comm/b5552444444

---

Challenges. Applying this GP framework I 2 directly to large, complex real-world ST datasets faces several significant hurdles:
1. Computational Scalability: Exact GP inference requires operations on the N N co- variance matrix (where N is the number of observations), entailing ((N9) computational complexity and (NN2) memory requirements. This quickly becomes prohibitive for datasets with thousands, let alone millions, of points [12]
2. Non-stationarity: The stationarity assumption, implying a covariance structure dependent only on ST lag (h., he), is often unrealistic. Environmental heterogeneities, varying physical dynamics, or boundary effects can cause the dependence structure (e.g., correlation range, variance) to change across space and time
3. UQ: While various approximation techniques are employed to address scalability (C1) and non-stationarity (C2), these approximations often break the theoretical guarantees of the exact GP model. Consequently, the resulting predictive uncertainties (e.g., variances, prediction intervals) may lack statistical validity, potentially failing to achieve the desired nominal coverage probability [10 38).
Spectral Methods for CI. To overcome the computational bottleneck (CI), prominent class of methods leverages the spectral representation of stationary GPs [10 12. Bochner's theorem guarantees that any stationary covariance function C(h,, ht) = Cov (s, t), Z(s + h,, + he)] is the Fourier transform of spectral measure. For the Matérn and many other processes, this measure has a density o² (w), where (w) is the normalized spectral density (f f(w)da 1) and o² C(O, 0) is the process variance. The covariance is recovered from the spectral density via the Wiener-Khinchin theorem:
C(h,, (1)) (1) = o² { cos(w h. (whht)f())) dw,
R3
where w = (w,, wi) represents the ST frequencies. The spectral density f(w) corresponding to the Matérn covariance (2 known to have the form of multivariate Student's t-distribution [TO]. This spectral view motivates computationally efficient approximations based on sampling frequencies w from (w), as detailed next. This spectral view motivates computationally efficient approximations using J random basis functions [45] such as Bayesian Random Fourier Features (BRFF) [46]. The
GP. Z. is approximated as:
(s, = Σ [oosww?ys wr, + sin(w",ss wrnt)byl
In the fully Bayesian BRFF approach, both frequencies and amplitudes are treated as random variables. Priors are chosen such that the resulting process Z, approximates Z: frequencies rid
rid sampled (wa.f ;) ~ f(w) (the Matérn spectral density), and amplitudes are aj, by ~ AN(0, o²/
Miller and Reich [46] demonstrated the computational efficiency and predictive performance of this using a fully Bayesian approach in spatial setting. However, their approach is unable to scale to large numbers of spatial, or ST observations, requiring a more computationally efficient method of utilizing the BRFF approximation.
Dimension Expansion for C2. To relax the often unrealistic stationarity assumption, dimension expansion methods propose that a non-stationary process ((s, 1) can be viewed as stationary in an augmented space [s, L(s, ()] 47 44). Here, L(s, () [L (s, ().... L, 4)] is mapping p-dimensional latent space, representing unobserved factors (e.g., local environmental conditions, dynamic regimes) that modulate the covariance structure. Given L. stationarity is recovered, and
a stationary kernel like Matérn (2) can be used with modified distance incorporating the latent variables:
= --1s-s'11 s' 1/1///22
1)2 + ((li.t))Lsssttee L(s', t()
(5)
The latent dimension p and the processes L(s, t) themselves are typically unknown and must be inferred from the data, adding complexity the modeling task.
Conformal Inference for C3. Addressing the challenge of obtaining statistically valid uncertainty estimates, especially when using approximations for scalability or non-stationarity, can be achieved
3

---

using conformal inference [38 48.. Consider observations Wi...., W.. where each W, (Y, X,) consists of a response Y; E IR and covariates X, E IRd. The core assumption of conformal inference that the sequence W....., Wa, Wn+ (including a new, unseen test point) is exchangeable. meaning their joint distribution is invariant under permutation. Conformal inference operates using non-conformity measure A, where dę A(W, W-i) quantifies how dissimilar W; is from the set W- (Wi...., Wn+1) (W.). The specific function A is user-defined; common choices involve residuals from a fitted model (e.g.. A(W, W-)) I (X,() where ğ fitted on W-1). To construct a prediction interval for Ya41 given X,+1 and training data (W.).. conformal considers hypothetical values y for Ya+1. For each y. computes the con-conformity score (1) = AW,, W-.). The plausibility of is measured by its conformal p-value:
n++1
1
p(u) = n+1 166,
(6)
where 6; for S n are typically computed using leave-one-out retraining or, more efficiently, using residuals held-out calibration set (split conformal inference). The 100(1 a3% prediction interval for Yn+1 comprises all values y deemed sufficiently plausible:
n R+1 {y E R p(y) > a).
Under only the exchangeability assumption, this interval guarantees P(Y E r n+1 2 38].
While spatial variants of conformal inference have been developed 42 it has not yet been well explored in the ST setting. The challenges remaining in CI, C2 and C3 motivate our STACI
methodology outlined in Figure !
3
Methodology
Approximate STGP architecture
Predict Posterior mean surface
Valid Uncertainty Quamification
Coendrates
Spece
INR
Approcizaaee CP
M
Estiman Protzrite with the distribution
of Appraacianee GPss
Contomaal Infersce
Figure I:
STACI algorithm pipeline. Data from coordinates across space and time are first fed into our approximate spatio-temporal GP architecture and trained using SVGD. A conformal inference step is then fit to provide customized prediction intervals for valid uncertainty quantification.
3.1 Bayesian Neural Network Architecture
The greatest drawbacks to the BRFF approximation in (4) are its reliance on assuming a stationary GP kernel, as well as limited scalability stemming from using full MCMC samplers for the frequencies and covariance parameters.
To resolve the first issue, we model Z in ①I using the non-stationary dimension-expansiii covariance function in (5). Given latent field L., the spectral approximation is
((s t) Σ cos [W','' + Weyt + wT, L(s, t)] a, + sin [wi s we,t tt L
L(s, t)] by,
(8)
with frequencies Wij (wT Wi.ji wf,.) To approximate Matérn covariance, model the fre- quencies with a multivariate t distribution with v degrees of freedom, location zero, and diagonal
iid scale matrix with diagonal elements (ps: Pss Pt, Pt. Pt,J), i.e., w ~ MVT,(0, D). The am- (ss
plitudes are modeled as d,, by äid No,mal a2 /J). For computational simplicity and identifiability concerns, we set Pr, the following theorem justifies the spectral approximation to the covariance function. The proof included in Appendix A.I.
Pt.J pr- Extending (46) to the nonstationary spatiotemporal case,

---

Theorem The prior mean of the spatiotemporal covariance function of the discrete process in s equals the Matérn correlation with distance defined as in 5 for all and the point-wise prior variance decreases at rate J.
Latent field L is modeled by an INR with three dimensional input, (s, (). and p-dimensional output, t). This neural network model for the latent space provides both flexibility and computational scalability for estimating the nonstationary covariance function.
(s, ... (s, (s.,).
We note that the form of ( (8) closely resembles a single hidden layer MLP where the first layer is passed through cosine and sine transforms with the bias terms removed and the hidden layer dimension J. Turning this into a Bayesian MLP with priors set to the aforementioned distributions related to (8) this is exactly the form of the BRFF. We then utilize a skip connection, appending the original location set and times to L(s, 1) resulting in the final architecture to estimate ((s, t). We place relatively uninformative priors the covariance parameters and utilize the Stein variational gradient descent (SVGD) algorithm developed by 49 to efficiently train our model and obtain variational posterior distributions.
3.1.1 Posterior Estimation
We use the SVGD variational inference (VI) framework to approximate the posterior distributions of neural network model. For some parameter set 9 = (01,..., 0p) R and data Y, VI approximates a target posterior, p(0 Y) = p(0)Y) using a simpler distribution q(0), found a predefined family Q {9)) by minimizing the KL divergence
q" (0) arg min ( (gllp)) (()) E, (eglog(((
Here, /)o|Y) is the un-normalized posterior distribution. VI turns posterior estimation into an optimization problem. allowing for greater scalability than sampling- baaseingg--sedddddd methods of posterior estimation [50; 52]. However, there can be a tradc-off in approximation accuracy. For example, pre-specifying the variational family can result in under-estimation of the posterior variance [22
SVGD is a VI algorithm that does not require specifying a variational posterior priori, making more generalized approach. For parameter set 0 = { .... 0,), SVGD initializes set of M independent particles (copies), {0} 0 O,,..., Os, that will be trained to approximate the posterior distribution through minimizing the KL divergence between these copies and the target posterior. Given data Y, priors (01),,.., π (0p)). joint prior (0)) = I π (0,), joint likelihood L(Y)0) and kernel function <(...) .), for each 0, (0). the update at iteration with step-size Eς is
oCt 0 + ccd(0))
(10)
with smooth optimal perturbation
M
(00) = M Σ [n(0) o + log L(Y)0))) + Vo 0)) ()))
(11)
As .M + x0, the distribution of {0} approaches the variational posterior (0) In (1I). the first term draws particles towards high probability areas of posterior p p(0) Y) while the second term is repulsive force that discourages particles from grouping in local modes of the posterior.
3.1.2 Neural Network Uncertainty Quantification
Using the SVGD training method gives M trained neural networks (.). Thus for observation i, the MAP estimate of Y(s,, 4) is
1 M
(s, [ = (s, (.) + c 14)
Σ [x(x,)
(12)
M j=1
with marginal variance
(IY = =V E[Z(s,,,,,
t.) E[Z(s;,
(13)
5

---

We use the MAP estimate from the trained neural networks, f², to estimate the nugget effect. The nugget variance is the random/mmeasurement eπror term- aleatoric uncertainty. Given desired type-1 error a, letting the mean estimate (2) be Y, and variance (3.122 be à? we are able to construct 100(1 a)% coverage level credible interval (CI):
C = (L; U:) = Y; + (1-0)//20
(14)
where (1-) is the value of the standard normal cumulative distribution function (CDF) the (1- (1-a)/2 percentile. However, to this point, we have made multiple approximations. First, is the approximation of the full STGP process with our trained Bayesian neural network. Second, is the approximation of the aleatoric nugget variance, 72. Thus, our estimated CI cannot be considered statistically valid and we could have over/under estimation of the true variance. The second conformal inference part of STACI serves to correct this issue.
3.2
Spatio-temporal Conformal Prediction
We extend the local spatial conformal prediction algorithm of [41 to a ST setting to provide valid prediction intervals for our approximate GP model. To simplify notation, for observation denote the response as Y; = Y (s;, 44), the coordinates as X; = (s, (4) and the pair as W, - (Y;, X:). The Bayesian computations in the previous section provide fitted values Y, and working standard errors oN As the data is non-stationary over space and time, is unreasonable to assume exchangeability between all n observations approximately exchangeable neighbors across space and time. The neighbors for prediction site x,+1 selected as the K closest training observations based squared distance
in Section 2 Thus, the first step is identify a local set of K
|||s; - S 1//p_ + [|t, - fn+ 1/
(15)
where the estimates of spatial and temporal range, Ps and Pr. are given by the neural network in the first step of STACI. We denote the indices of K neighbors and the prediction location + 1 as Nx+1 C (L, n + 1}.
we provisionally set Yn+1 = yn++ and use discrepancy measure 6, = =1y; Y5l/āj, the plausibility
is
plyn+1)
1
Σ JEN
160 2 5m+1
(16)
To identify the prediction interval, we search over the range of Y, among the K neighbors and take the interval as the set of Un+1 with plausibility at least a. Combining S3.T and 3.2 the full
STACI algorithm provides a scalable non-stationary STGP approximation with valid 100 a)%
prediction intervals during ST interpolation.
4
Application
4.1
Setup
4.1.1
Data Description
We evaluate STACI's performance on two ST datasets: one synthetic and one real. The first synthetic dataset is simulated mean sea surface height (MSS) data of the Arctic sea [37] based on historical satellite data from 3 different tracks. The data spans 10 days from March 1st to March 10th 2020 and has 1,158,505 total datapoints. The second real dataset is Aerosol Optical Depth (AOD) data captured using the Moderate Resolution Imaging Spectroradiometer (MODIS) on NASA's Terra satellite [53]. The data is spread on 1400 x 720 grid spanning the Earth's 's surface and we use daily data from March 2025, equating to 3,189,641 total observations.
The MSS dataset is randomly split into a 80% train, 10% validation, 10% test split. As most observations are seen on each day, this setting tests ST interpolation with small artifacts such as cloud cover creating patches of missing data. For the AOD dataset, we sample 10% of observations randomly per day to comprise the training set. The validation set all observations over the first days while the test set is all observations on the 20th day. This mimics task of full reconstruction

---

over field when there are limited sensors providing the ground truth data. We note that even the full AOD field is sparse, adding further difficulty to the reconstruction task. For both datasets, the spatial coordinates are scaled to [0, 1] x [0, 1] grid and the response is normalized to the training set mean and standard deviation.
4.1.2 Model Settings
We test STACI with 3 different ofate-of-thhear INR architectures to estimate the latent space: Residual Multilayer Perceptron (ResMLP) [34]. Fourier Feature Network with Positional encoding (FFNP) and FFN with Gaussian encoding (FFNG) [3). Each INR backbone has 5 layers with a layer width of 1,024. We set J 5,000 for the final hidden layer width of STACI, representing the number of random fourier features. We then use M = 10 network copies to train using SVGD for the initial Bayesian UQ. Finally, we determine D through cross validation on a random subset of 100 training locations for the number of neighbors to use for conformal inference.
As we are an approximate GP, we compare our algorithm to state-of-the-art scalable GP methods that provide UQ. The first method is sparse variational GP (SVGP) from [4], representing a case of using an approximate stationary GP. The second method doubly stochastic deep GP from [35], representing an alternate method of estimating non -on-stationary GP. The third method is GPSat from [27]. modeling the full stattitttonary field as a mixture of local, stationary GPs. The final method is deep random features (DRF) from [7], representing similar algorithm of turning the spectral representation of GP into Bayesian deep neural network. SVGP and DeepGP are trained through variational inference, while DRF's is trained through Bayesian optimization. DeepGP is set with 4 total layers with layer width 9 and trained with 10 models. GPSat is initialized with 1,225 expert locations across the spatio-temporal domain. DRF is set with 5 hidden layers of width 1,024 and bottleneck layers of width 128 and trained with 10 models, All models are trained on NVIDIA A- GPUs for 15 epochs (optimization iterations for DRF) with batch size 1,024. Conformal prediction interval calculation is parallelized over NVIDIA A-100 GPUs.
4.1.3 Performance Metrics
We measure the performance of the algorithms in terms of both estimation and UQ quality on the test set of both datasets. For estimation quality, we use root mean square error (RMSE), negative Gaussian log likelihood (NLL). For UQ quality, we use the continuous ranked probability score (CRPS) metric and provide coverage of prediction intervals, interval score and interval width based on = 0.05. Finally, we track time per training epoch to compare computational efficiency. The time for DRF is the time for one optimization iteration. Note that GPSat time is the total time to fit across all expert locations. We provide both Bayesian and conformal UQ performance for STACI. The conformal time represents time needed for the entire conformal step. For RMSE, NLL, CRPS, interval score and interval width, lower value is better. For coverage, the value closest 0.95 is deemed the best as estimators providing both over-coverage and under-coverage are considered to be inefficient.
4.2
Main Results
Table results for the MSS dataset. We see that STACI with the FFNP latent model has the lowest RMSE and NLL, indicating accurate estimation of the surface. We also see this the lowest interval score and interval width of the non-conformal methods while maintaining desired coverage. The conformal interval score is at at worst half that of its competitors, while also being under half as wide as most others. This indicates the individualized interval widths provide much more efficient UQ while achieving the desired coverage level. Total fit time for conformal is also manageable, taking under 8 minutes for the FFN latent models and around 10 minutes for the ResMLP model. We see that SVGP performs much worse than its counterparts in this setting, showcasing the need for inclusion of non-stationarity. GPSat provides comparable estimation metrics to STACI, however the total runtime is over 3 times longer.
Table 2 shows results for AOD dataset. STACI with the FENP latent model again has the lowest test set RMSE and interval UQ metrics. The conformal addition lowered the interval score by over half, showing the efficiency of this UQ method despite the difficult training and data setting. GPSat had better NLL and near-identical CRPS, albeit again with higher runtime. We note that Deep RF had
7

---

Table 1: Comparisons the MSS dataset. Best performances highlighted in bold and
underlining for top and second-best methods, respectively.
Model
RMSE
NLLL
CRPS
Coverage
Interval Score
Interval Width
Time (8)
STACI-ResMI. (Bayes)
0.422
-0.363
0.221
0.948
2.509
1.683
132
STACI-ResMLP (Conf.)
NA
NA
NA
0.958
0.514
0.500
652
STACI-FFNP (Bayes)
0.161
-1.331
0.086
0.951
0.923
0.663
116
STACI-FFNP (Conf.)
NA
NA
NA
0.958
0,514
0.500
446
STACI-FFNG (Bayes)
0.247
-0.907
0.129
0.948
1.458
0.992
105
STACI-FFNG (Conf.)
NA
NA
NA
0.958
0.514
0.500
431
Deep RF
0.203
-1.112
0.106
0.970
1.172
0.941
83
Deep
0.277
-0.782
0.145
0.957
1.633
1.171
GPSat
0.200
-1.188
0.102
0.974
1.140
0.939
7020
SVGP
0.447
-0.304
0.236
0.951
2.571
1.844
difficulties with this high data sparsity setting and required fixing the nugget variance parameter to achieve convergence.
Table 2: Comparisons on the AOD dataset. Best performances are highlighted in bold and underlining for top and second-best methods, respectively.
Model
RMSE
NLL
CRPS
Coverage
Interval Score
Interval Width
Time (s)
STACI-ResMLP (Bayes)
0.720
0.361
0.400
0.916
5.506
2.656
(Conf.)
NA
NA
NA
0.948
1.850
1.555
598
STACI-FFNP (Bayes)
0.560
0.042
0.295
0.944
3.944
2.148
STACI-FFNP (Conf.)
NA
NA
NA
0.948
1.850
1.555
410
STACI-FFNG (Bayes)
0.675
0.275
0.372
0.929
5.098
2.542
STACI-FFNG (Conf.)
NA
NA
NA
0.948
1.850
1.555
391
Deep RF
0.733
0.175
0.452
0.954
5.433
2.264
103
Deep GP
0.633
0.191
0.356
0.941
4.781
2,788
70
GPSat
0.653
-0.064
0.287
0.958
3,476
2.272
2745
SVGP
0.690
0.276
0.395
0.956
5.150
3,237
20
The difference in prediction quality and UQ between the models can be visualized in Figure 2 The first of Figure represents h the ground truth AOD values and the predicted means from the STACI Bayesian and conformal variants as well as the competitor models. The color scheme for this row goes from blue (lower values of opacity) red (higher values of opacity). The second row represents the widths of the prediction intervals from each of the methods assuming a desired 95% coverage. The darker purple represents narrower widths that moves lighter into yellow which represents the wider. Comparing to the ground truth, we see STACI is able to capture more of the high pollution area than the competitors mean surface prediction. Meanwhile, for the interval width representing UQ. we see that the conformal variant of STACI is the only one that has highly varied interval widths, while the UQ from most other methods is fairly stagnant in the middle purple and yellow. The UQ benefit of GPSat is shown here where the local GPs provide personalized UQ similar to the conformal variant of STACI. The high variance areas from conformal and GPSat also correspond to areas with deep red high spikes in pollution which is commonly smoothed over in these classes of models.

---

Ground Truth STACI-B
STACI-C
DRF
DeepGP
GPSat
SVGP
Figure 2: Predicted AOD surface values. Top row: predicted surface. Red indicates higher AOD values.; Bottom row: interval widths for Bayesian and conformal (STACI-C) uncertainty on AOD data. Darker shades denote narrower intervals.
4.3
Ablation Studies
Here, we perform two ablation studies for the AOD dataset to quantify model robustness. The first ablation study is impact of latent model dimension size estimation error. The second ablation is impact of sampling percentage in training set construction on estimation error.
Ablation Study: Latent dimension
Figure 3a shows how the RMSE of each latent model varies we change the latent dimension size from [32. 64, 128, 256]. We see that the two FFN models are fairly stable at around 0.58 for positional encoding and around 0.67 for Gaussian encoding. ResMLP error sees noticeable decrease increasing from latent dimension 128 to 256. This indicates the FFN latent model results are stable to choice of latent dimension and the reported results are optimal.
Ablation Study: Sampling Percentage
Figure F shows how the RMSE of each approximate GP model varies as we change the training set sampling percentage at each time-step from [5%, 10%, 25%]. For STACI, we use the FFNP latent model. We see generally a decrease in estimation error with greater sampling percentage. Interestingly, DRF error rises slightly at the highest sampling percentage, but remains below the 5% case. We speculate this is due fixing the noise parameter during Bayesian optimization. STACI still achieves the lowest error across settings.
RMSE vs Latent Dimension
RMSE vs Training Sample Percentage
0.85
1.0
Model
0.80
FFNP
STACI
FFNG
0.9
DRF
0.75
ResMLP
DeepGP
0.8
RMSS
0.70
RASE
0.65
0.7
0.60
0.6
0.55
0.5
0.50
150
Latent Dimension
Sample Percentage
(a) RMSE latent dimension size for FENP, FFNG, (b) RMSE vs. train set sampling percentage for STACI.
ResMLP.
DRF. DeepGP, and SVGP.
Figure 3: Ablation Studies for STACI

---

Conclusion
In this work, introduce the STACI algorithm, an approximate STGP combining high scalability and prediction accuracy with valid UQ. Compared to other deep and neural network based models, STACI is able to provide interpretable covariance parameters from simpler stationary kernels using the latent dimension expansion, along with the ability to recover the data ST correlation structure. limitation of the approach is the computational expense both speed and space of the SVGD training algorithm with larger numbers of models. With smaller model number, we are not able to fully explore the posterior space, and may converge to suboptimal estimates. This impacts the strong assumption of exchangeability of selected neighbors is impacted by our STACI's ability accurately estimate the relevant covariance parameters. Additionally, the conformal algorithm as constructed is fairly simplistic, assigning equal importance each selected neighbor. This can be optimized further with weighting proportional to ST distance. Nevertheless, STACI provides scalable aleatoric UQ and can be used in variety of ST interpolation tasks.
One possible alternative direction is to use data-driven basis functions, such as the leading terms in Karhunen-Loc Loève expansion, as opposed to trigonometric functions in the output layer. We also believe that there are multiple possible extensions to non-Gaussian data. For example, adding a link function, slightly changing the likelihood, should extend to distributions such as Bernoulli, or Poisson model spatio-temporal correlated binary and count data. We also believe an extension to model spatial point processes is a feasible adaptation of the model. Following the work of [54], the Poisson process intensity function can be modeled as the square utilization of both random Fourier features and our introduced approximation for efficient likelihood computation. These extensions would allow STACI to reach a wider audience and cover additional types of data.
GP. This would allow continued
Acknowledgments
This work was supported by the National Institutes of Health (R01ES031651-01), the National Science Foundation (DMS2152887), and the U.S. Department of Energy (DOE), Office of Sci- ence (SC), Advanced Scientific Computing Research program under award DE-SC-0012704 and KJ0401010//C147 and used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility using NERSC award NERSC DDR-ERCAP0030592.
10

---

References
[1] Noel Cressie and Christopher K Wikle. Statistics for spatio-temporal data. John Wiley & Sons,
2011.
[2] Senzhang Wang, Jiannong Cao, and S Yu Philip. Deep learning for spatio-temppoal data mining:
A survey. IEEE transactions knowledge and data engineering. 34(8):3681-3700, 2020.
[3] Kota Katanoda, Yasumasa Matsuda, and Morihiro Sugishita. A spatio-temporal regression model for the analysis of functional mri data. Neurolmage, 17(3):415-14223 2002.
[4] Nikola K Kasabov. Neucube: A spiking neural network architecture for mapping, learning and understanding of spatio-temporal brain data. Neural networks, 52:62-76, 2014.
[5] Qiang Zhang, Qiangqiang Yuan, Chao Zeng. Xinghua Li, and Yancong Wei. Missing data recon- struction in remote sensing image with a unified neural network. IEEE Transactions on Geoscience and Remote Sensing, 56(8):4274 4288, 2018.
deep convolutional
[6] Matthias Katzfuss and Noel Cressie. Spatio-temporal smoothing and em estimation for massive remote-sensing data sets. Journal of Time Series Analysis, 32(4):430-446, 2011.
[7] Markus Reichstein, Gustau Camps- Valls, Bjom Stevens, Martin Jung, Joachim Denzler, Nuno Carvalhais, and F Prabhat. Deep learning and process understanding for data-driven earth system science. Nature, 566(7743):195-244, 2019.
[8] Ashesh Chattopadhyay, Pedram Hassanzadeh, and Saba Pasha. Predicting clustered weather patterns: A test case for applications of convolutional neural networks to spatio-temporal climate data. Scientific reports, 10(1):1317, 2020.
[9] Xuelong Li, Qun Guo, and Xiaoqiang Lu. Spatiotemporal statistics for video quality assessment. IEEE Transactions on Image Processing, 25(7):33293334)) 2016
[10] Michael L Stein. Interpolation of spatial data: Some theory for Kriging. Springer Science & Business Media, 1999.
[11] Sudipto Banerjee, Alan E. Gelfand, Andrew O Finley, and Huiyan Sang. Gaussian predictive process models for large spatial data sets, Journal of the Royal Statistical Society Series B: Statistical Methodology, 70(4):825-84) 2008.
[12] Christopher KI Williams and Carl Edward Rasmussen. Gaussian processes for machine learning, volume 2. MIT press Cambridge, MA, 2006.
[13] Michalis Titsias. Variational learning of inducing variables in sparse gaussian processes. In Artificial intelligence and statistics, pages 567-574. PMLR, 2009.
[14] James Hensman, Nicolo Fusi, and Neil D Lawrence. Gaussian processes for big data, arXiv preprint arXiv: 1309.683,, 2013.
[15] Montserrat Fuentes and Brian Reich. Spectral domain. Handbook of Spatial Statistics, pages 57-77, 2010.
[16] Joseph Guinness. Gaussian process learning via fisher scoring of vecchia's approximation. Statistics and Computing, 31(3):25, 2021.
[17] Abhirup Datta, Sudipto Banerjee, Andrew O Finley, Nicholas AS Hamm, and Martijn Schaap. Nonseparable dynamic nearest neighbor gaussian process models for large spatio-temporal data with application to particulate matter analysis. The annals of applied statistics, 10(3):1286. 2016.
[18] Oliver Hamelijnck, William Wilkinson, Niki Loppi, Amo Solin, and Theodoros Damoulas. Spatio-temporal variational Gaussian processes. Advances in Neural Information Processing Systems, 34:23621-23633. 2021.

---

[19] Simo Sarkka, Arn Solin, and Jouni Hartikainen. Spatiotemporal learning via infinite- dimensional Bayesian filtering and smoothing: A look Gaussian process regression through Kalman filtering. IEFF Signal Processing Magazine, 30(4):51-61, 2013.
[20] Junpeng Zhang. Yue Ju, Biqiang Mu, Renxin Zhong, and Tianshi Chen. An efficient implemen- tation for spatial-temporal Gaussian process regression and its applications. Automatica, 147: 110679, 2023.
[21] Sahil Garg, Amarject Singh, and Fabio Ramos. Learning non-stationary space-time models for environmental monitoring. In Proceedings of the AAAl Conference on Artificial Intelligence. volume 26, pages 288-294, 2012.
[22] Emilio Porcu, Reinhard Furrer, and Douglas Nychka. 30 years of space-time covariance functions. Wiley Intendisciplinarr Reviews: Computational Statistics, 13(2)=e151,, 2021.
[23] Lyndsay Shand and Bo Li. Modeling nonstationarity in space and time. Biometrics, 73(3): 759-768, 2017.
[24] Christopher K Wikle. A kernel-based spectral model for non-gaussian spatio-temporal processes. Statistical Modelling. 2(4):299-314, 2002.
[25] Sami Remes, Markus Heinonen, and Samuel Kaski. Non-stationary spectral kernels, Advances in neural information processing systems, 30, 2017.
[26] Wanfang Chen, Marc G Genton, and Ying Sun. Space-time covariance structures and models. Annual Review of Statistics and Its Application, 1(1):19 2021.
[27] William Gregory, Ronald MacFachern, So Takao, Isobel R Lawrence, Carmen Nab, Mare Peter Deisenroth, and Michel Tsamados. Scalable interpolation of satellite altimetry data with probabilistic machine learning. Nature Communications, 15(1):7453, 2024.
[28] Pratik Nag, Ying Sun, and Brian Reich. Spatio-temporal deepkriging for interpolation and probabilistic forecasting. Spatial Statistics, 57:100773, 2023.
[29] Bing Yu, Haoteng Yin, and Zhanxing Zhu. Spatio-temporal graph convolutional networks: A deep learning framework for traffic forecasting, arXiv preprint 1709, 04875, 2017.
[30] Xihaier Luo, Wei Xu, Yihui Ren, Shinjae Yoo, and Balu Nadiga. Continuous field reconstruction from sparse observations with implicit neural networks. arXiv preprint arXiv:2401.11611 2024.
[31] Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein. Im- plicit neural representations with periodic activation functions. Advances in Neural Information Processing Systems, 33:7462-7473, 2020.
[32] Matthew Tancik, Pratul Srinivasan, Ben Mildenhall, Sara Middvvich-Keil, Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan Barron, and Ren Ng. Fourier features let networks learn high frequency functions in low dimensional domains. Advances in Neural Information Processing Systems, 33:7537-7547, 2020.
[33] Yiheng Xie, Towaki Takikawa, Shunsuke Saito, Or Litany, Shiqin Yan, Numair Khan, Federico Tombari, James Tompkin, Vincent Sitzmann, and Srinath Sridhar. Neural fields in visual computing and beyond. In Computer Graphics Forum, volume 41, pages 641-676. Wiley Online Library, 2022.
[34] Langwen Huang and Torsten Hoefler. Compressing multidimensional weather and climate data into neural networks. arXiv preprint 1rXiv:2210.122 /2538, 2022.
[35] Hugh Salimbeni and Marc Deisenroth. Doubly stochastic variational inference for deep Gaussian processes. Advances in Neural Information Processing Systems, 30, 2017.
[36] Andreas Damianou and Neil Lawrence. Deep gaussian processes. Artificial intelligence and statistics, pages 207-215. PMLR, 2013.
12

---

[37] Weibin Chen, Azhir Mahmood, Michel Tsamados, and So Takao. Deep random features for scalable interpolation of spatiotemporal data. In The Thirteenth International Conference on Learning Representations, 2024.
[38] Vladimir Vovk, Alexander Gammerman, and Glenn Shafer. Algorithmic learning in random world, volume 29. Springer, 2005.
[39] Victor Chernozhukov, Kaspar Wüthrich, and Zhu Yinchu. Exact and robust conformal inference methods for predictive machine learning with dependent data. In Conference On learning theory, pages 732-749. PMLR, 2018.
[40] Jing Lei, Max G'Sell, Alessandro Rinaldo, Ryan J Tibshirani, and Larry Wasserman. Distribution-Fiee predictive inference for regression. Journal of the American Statistical Associ- ation, 113(523):104 1094-1111, 2018.
[41] Huiying Mao, Ryan Martin, and Brian J Reich. Valid model-free spatial prediction. Journal of the American Statistical Association, 119(546):904 914, 2024.
[42] Hanyang Jiang and Yao Xie. Spatial conformal inference through localized quantile regression. arXiv preprint arXiv:2412.01098, 2024.
[43] Christopher K Wikle, Andrew Zammit-Mangion, and Noel Cressie. Spatio-temporal statistics with R. Chapman and Hall/CRC, 2019.
[44] Christopher J Paciorek and Mark J Schervish. Spatial modelling using a new class of nonsta- tionary covariance functions. Environmetrics, 17(5):483-500 506, 2006.
[45] Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. Advances in Neural Information Processing Systems, 20, 2007.
[46] Matthew J Miller and Brian J Reich. Bayesian spatial modeling using random Fourier frequen- cics. Spatial Statistics, 48:1005,8, 2022.
[47] Luke Bornn, Gavin Shaddick, and James V Zidek. Modeling nonstationary processes through dimension expansion. Journal of the American Statistical Association, 107(497):211-28,, 2012.
[48] Glenn Shafer and Vladimir Vovk. A tutorial on conformal prediction. Journal of Machine Learning Research, 9(3), 2008
[49] Qiang Liu and Dilin Wang. Stein variational gradient descent: general purpose Bayesian inference algorithm. Advances in Neural Information Processing Systems, 29, 2016.
[50] Alex Graves. Practical variational inference for neural networks Advances in Neural Informa- tion Processing Systems, 24, 2011.
[51] Matthew D Hoffman, David M Blei, Chong Wang, and John Paisley. Stochastic variational inference. The Journal of Machine Learning Research, 1(1):1303-1347, 2013.
[52] David M Blei, Alp Kucukelbir, and Jon D McAuliffe. Variational inference: A review for statisticians. Journal of the American Statistical Association, 1 12(558) 2017.
[53] Alexei Lyapustin, Yujie Wang, Sergey Korkin, and Dong Huang. MODIS collection 6 MAIAC algorithm. Atmospheric Measurement Techniques, 11100):744-5565, 2018
[54] Chris Lloyd, Tom Gunter, Michael Osborne, and Stephen Roberts. Variational inference for gaussian process modulated poisson processes. In International Conference Machine Learning, pages 1814-1822. PMLR, 2015.
13

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: We show through two different real spatio-temporal datasets and settings that STACI provides lower prediction and more accurate and interpretable uncertainty quantification. We include relevant theoretical outcomes in the paper and proofs are included in the Appendix.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation as long as is clear that these goals attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
We mention potential limitations of this work in the conclusion of the paper.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and prool's
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
14

---

Answer:
[Yes]
Justification: Assumptions and proofs of theorems are included the Appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: Full data description with how the training details apply to real scenarios are provided in the text. Train/test details along with model hy hyperparameters are included in both the main text and the Appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reprodueibility may be tricky in some cases, in which authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
15

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification: Code is included in the supplement of the material with relevant scripts to recreate results.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://niis.cc
aublic/guudees/Codeettbiiss ott.eeeeuuuiiiiiiooooooiiooo for more details.
While encourage the release of code and data, understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
/ wwtps/cccppppgguuuaaaeoeaec
wtt.eeeeeeooSSeeesSSSeoeee.. wsw.sssoiicgyyy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Data deseription with train/test split details are given in the text and further details are given in the Appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented the core the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: Coverage properties of the estimator are included in the results section of the paper.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
16

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
OK to report 1-sigma error bars, but one should state The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: GPUs used are mentioned in the main text and further details are given in the Appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't 't make it into the paper).
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics httppp////errrp...cc/paa/////h/p
Answer: [Yes]
Justification:
Paper conforms to the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: Method is created for the scalable pnatio-temporal and uncertainty quan- tification research community general. can provide positives such as higher climate prediction accuracy, better MRI reconstruction and provide the ability to quantify uncertainty noisy spatio-temporal fields. The method is not directly tied to any negative applications far as the authors aware.
17

---

Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: No such risks for misuse.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: Comparison works are properly cited in the main text and further details are provided the appendix.
Guidelines:
The answer NA means that the paper does not existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
18

---

For scraped data from a particular source (e.g., website). the copyright and terms of service of that souree should be provided.
assets are released, the license, copyright information, and terms use in the datasets.paperssiii paperswi thcode. www.coo//aaaaa datasets package should be provided. For popular datasets,
has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification:
Details of dataset, code and model are included the Appendix.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moodl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
Paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals equivalent for research with human subjects
Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Question:
Answer: [NA]
Justification: Paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
19

---

Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if' applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: LLM does not impact core methodology, scientific rigorousness, or originality of the research.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our LLM policy (https://neurd) hate:p///earrrgscccooofeeeeeeeee / for what should or should not be described.
20


Reviews:
[{"review_id": "AQ3ktXmCTk", "reviewer": "Reviewer_8XoG", "summary": "The paper consider the problem of inferring GP models using a variational Bayesian neural network approximation. A conformal inference algorithm is proposed and compared with three baselines including SVGP, deep RF, and deep GP.", "questions": "- How does the algorithm compare with baselines mentioned above, both theoretically and numerically?\n- Can any guarantees be provided for the approximation methods used in this algorithm? \n- Can any theoretical analysis on the improvement be provided?", "limitations": "yes", "rating": 3, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 1, "originality": 1, "strengths_and_weaknesses": "+ The idea of using a variational Bayesian neural network to approximate non-stationary spatio-temporal GPs makes sense. \n+ The proposed solution is clearly described.\n- The evaluation is quite limited and missed many important baselines, including both GP and more advanced Gaussian Cox process models, such as variational Bayesian approximation, mean-field approximation, path integral approximations of GP, and change of kernel in a reproducing kernel Hilbert space. Without a complete comparison, the benefit of the proposed method is hardly convincing.\n- The pros and cons of the proposed algorithm compared to these baselines isn't clear.\n- Can any guarantees be provided for the approximation methods used in this algorithm? Can any theoretical analysis on the improvement be provided?", "was_revised": true, "final_justification": "The paper considers GP and not more complex problems like non-Guassian or Gaussian Cox processes. Based on the responses provided by the authors, I'm raising the review score.", "rebuttal": "We would like to thank the reviewer 8XoG for taking the time to read through our paper and providing much constructive feedback. We appreciate the reviewer's acknowledgment of the clarity of our work and the logic of the method. Rebuttals to specific comments are provided.\n\n>The evaluation is quite limited and missed many important baselines, including both\nGP and more advanced Gaussian Cox process models, such as variational Bayesian\napproximation, mean-field approximation, path integral approximations of GP, and\nchange of kernel in a reproducing kernel Hilbert space. Without a complete comparison, the benefit of the proposed method is hardly convincing.\n\nWhile we may have been unclear, we believe we have compared our method to many of these approaches. It is computationally infeasible to fit a GP to the full dataset, and all of the methods we compare to are approximations to a GP. We covered the variational Bayesian baseline by including SVGP. This KL-divergence based method was also shown to have better performance than the mean-field approximation in the original SVGP paper, so we did not include this in our comparisons. We also believe the Gaussian Cox process model, while advanced, is not in the scope of this work as we are not measuring count data across space and time. However, a potential area of future work is extending the method to non-Gaussian data, like the Gaussian Cox process model does. Replacing the Gaussian likelihood with other distributions such as Bernoulli or Poisson to handle geo-referenced binary or count data should be fairly straightforward. It should also be possible to extend the method to more challenging cases such as inhomogeneous Poisson processes for spatial point pattern data. For example, the model of Lloyd et al (2015) takes the Poisson process intensity function as the square of a GP. Thus the likelihood with random Fourier features approximating the GP can be computed analytically in O(n) flops, permitting the same computational approach taken here to fit the non-stationary Fourier model to spatial point pattern data.\n\n>The pros and cons of the proposed algorithm compared to these baselines isn’t clear.\n\nWe explain the drawbacks of current baselines in paragraphs 2-4 of the introduction and apologize if this wasn't clear. Our method provides greater interpretability of covariance parameters than other deep GP methods and provides greater flexibility than the traditional GP. The conformal addition of our method is able to provide guaranteed coverage properties for prediction intervals, assuming the exchangeability condition holds. The drawback to the STACI method is the computational cost required for both SVGD training and conformal inference, however we do not believe it is prohibitive for most use cases.\n\n> Can any guarantees be provided for the approximation methods used in this algorithm? Can any theoretical analysis on the improvement be provided?\n\nWe agree that theoretical guarantees need to be made for the approximation and provided a theorem on the spectral covariance approximation in supplement section A. We show that the spectral approximation is centered on the Matérn covariance function with accuracy that increases with $J$, the number of basis functions, or the width of the final hidden layer in our method.\n\n**Reference**\nLloyd, Chris, et al. \"Variational inference for Gaussian process modulated Poisson processes.\" International Conference on Machine Learning. PMLR, 2015."}, {"review_id": "39KrFy5uKT", "reviewer": "Reviewer_WrPZ", "summary": "This paper proposes STACI, a novel two-stage framework that combines a Bayesian neural approximation of non-stationary spatio-temporal GP with a conformal inference algorithm for uncertainty quantification. The authors address three key challenges in spatio-temporal modeling: (1) scalability, (2) non-stationarity, and (3) valid uncertainty quantification. The model leverages spectral approximations and a dimension expansion technique to approximate non-stationary covariance functions. STACI also incorporates conformal prediction to correct potential under/over-confidence introduced by approximations.\n\nEmpirical evaluations on both synthetic and real-world datasets demonstrate that STACI achieves superior prediction accuracy and tighter, valid prediction intervals compared to deep GPs, sparse GPs, and Deep RF methods.", "questions": "1.\tWhile the proposed spectral approximation is theoretically justified, it would be helpful to understand the error bounds associated with the overall approximation. Specifically, how do factors such as the number of random features or SVGD particles affect the fidelity of the STGP approximation and the resulting prediction intervals?\n\n2.\tHave you tried to compare with the mixture model of local GPs in GPSat library?\n\n3.\tWhy not using Continuous Ranked Probability Score(CRPS) to evaluate the quality of predicted uncertainties?", "limitations": "yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strength:\n\n1.\tThe authors propose a novel, scalable method that provides valid uncertainty quantification for non-stationary spatio-temporal processes—an important and under-explored area of research.\n\n2.\tThe paper is well-written, and the methodology and preliminaries are clearly presented, making complex concepts easy to follow.\n\n3.\tThe authors provide thorough experimental details, including ablations, training procedures, and model architecture. The visualization and results presented are clear and informative.\n\nWeakness:\n\n1.\tAlthough the authors compared the proposed model with several GP-based models in the experiments, they did not include a comparison with DNN-based models.\n\n2.\tThere are some typos and glitches here and there. E.g., Eq 9. Duplicate “at” in line 269.", "was_revised": true, "final_justification": "I appreciate the authors' thoughtful rebuttal and their efforts to incorporate feedback. The inclusion of CRPS results and the discussion on the spectral approximation and SVGD behavior address several of my concerns. However, two key points remain partially or fully unaddressed: (1) the lack of empirical comparison to standard DNN-based models, despite the method’s deep learning components, and (2) the omission of a known strong baseline (GPSat), which the authors acknowledged but have not yet evaluated. For these reasons, I have chosen to retain my original score.", "rebuttal": "We'd like to thank the reviewer WrPZ for providing valuable comments for our work and directions for further improvement. We appreciate highlighting the novelty of our approach and clarity of the paper format. For the typos, we appreciate the reviewer for pointing these out and will fix these for the camera-ready version. We provide some rebuttals for individual comments here.\n\n> Although the authors compared the proposed model with several GP-based models in the experiments, they did not include a comparison with DNN-based models.\n\n While we did not include DNN-based models, the backbone models for our latent space are popular Implicit Neural Representation DNN-based models currently in use for coordinate-based data. DNN-based models typically do not provide inherent uncertainty quantification, which was a major focus of our work. \n\n> While the proposed spectral approximation is theoretically justified, it would be helpful to understand the error bounds associated with the overall approximation. Specifically, how do factors such as the number of random features or SVGD particles affect the fidelity of the STGP approximation and the resulting prediction intervals?\n\nWe agree that theoretical justification is important for practical use of this method and show that the spectral approximation is centered on the Matérn covariance function. The approximation accuracy increases with $J$, the number of basis functions, or the width of the final hidden layer in STACI. In Appendix B.4, we also show that increasing the number of SVGD particles does not have much of an effect on RMSE, but slighty decreases the NLL, indicating slightly better variance estimation at the cost of higher computational time.\n\n> Have you tried to compare with the mixture model of local GPs in GPSat library?\n\nWhile we did not compare with the mixture model from GPSat in our original study, we thank the reviewer for bringing it to our attention. From the Deep RF paper, we see it performs very well in this area, placing near the top in each metric across multiple datasets and showing more flexibility than traditional GPs. Unfortunately we could not implement it in time thus far, however we will make an effort to have results for the camera-ready version.\n\n> Why not using Continuous Ranked Probability Score (CRPS) to evaluate the quality of predicted uncertainties?\n\nGreat point to use CRPS as a UQ quality metric! We include the CRPS results here for both datasets with the non-conformal variant of STACI. We see that STACI with FFNP backbone is still the highest performer across these use cases.\n\nFor the MSS Data:\n\n| Model          | Mean CRPS |\n| -------------- | --------- |\n| STACI (FFNP)   | 0.0978    |\n| STACI (FFNG)   | 0.1327    |\n| STACI (ResMLP) | 0.2442    |\n| DRF            | 0.1062    |\n| DeepGP         | 0.1449    |\n| SVGP           | 0.2361    |\n\nAOD Data:\n\n| Model          | Mean CRPS |\n| -------------- | --------- |\n| STACI (FFNP)   | 0.2904    |\n| STACI (FFNG)   | 0.3685    |\n| STACI (ResMLP) | 0.3768    |\n| DRF            | 0.4522    |\n| DeepGP         | 0.3557    |\n| SVGP           | 0.3950    |"}, {"review_id": "oZsLC0I5GR", "reviewer": "Reviewer_w1e9", "summary": "The problem with Gaussian Process approximations is that they only work well in low data cases. In this work the authors merge a Bayesian neural network approximation with the newly emerged conformal inference algorithm. This allows them to both produce a good approximation of the spatio-temporal data whilst still being able to accurately predict uncertainty quantification. The authors compare their work (in terms of RMSE and negative Gaussian log likelihood) with three common approaches from the literature on two datasets.", "questions": "- Please provide a clear description of figure 1.\n\n- Please ensure all concepts discussed in the paper are clearly defined before used.\n\n- Figure 2 is quite hard to interpret. Please provide more guidance.\n\n- Further consideration for the conclusions would help. Especially what to do next.", "limitations": "yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "The paper is technically sound and covers the problem well. However, the clarity of the presentation could be improved to make the paper more approachable for those not closely working in this area. For example, there is a figure which provides an overview of the whole process, but this is not discussed or explained in the text; the concept of nugget variance is discussed but not defined.", "was_revised": true, "final_justification": "Having read the reviews of others and the responses to those reviews I see that my original review was on the generous side. However, the comments to my questions lead me to believe that my score is now fitting for the paper. I do feel that the comments by WrPZ about missing baselines is a good one and strongly encourage the authors to achieve this for the final paper - this shouldn't just be a best effort.", "rebuttal": "We appreciate reviewer w1e9 for taking the time to understand and provide feedback for our method. We appreciate the praise for soundness of our work and apologize for some areas that lack clarity. We provide some rebuttals for individual comments here.\n\n> Please provide a clear description of figure 1.\n\n  We apologize for lack of clarity and will add the following. Figure 1: STACI algorithm pipeline. Data from coordinates across space and time are first fed into our approximate spatio-temporal GP architecture and trained using SVGD. This allows for estimation of the latent space and prediction of the mean surface through the variational posterior. A conformal inference step is then fit to provide customized prediction intervals for valid uncertainty quantification.\n\n> Please ensure all concepts discussed in the paper are clearly defined before used.\n\nWe apologize for not defining the nugget variance beforehand. The nugget variance is the portion of the variance that is fully random and can't be explained by our spatio-temporal model. This can be thought of as the random/measurement error term, or in other words the amount of aleatoric uncertainty. We will add this into the paper for the camera-ready version.\n\n> Figure 2 is quite hard to interpret. Please provide more guidance.\n\nThank you for pointing this out!  We will add a more informative description akin to:\n\nThe first row of Figure 2 represents the ground truth AOD values and the predicted means from the STACI Bayesian and conformal variants as well as the competitor models. The color scheme for this row goes from blue (lower values of opacity) to red (higher values of opacity). The second row represents the widths of the prediction intervals from each of the methods assuming a desired 95\\% coverage. The darker purple represents narrower widths that moves lighter into yellow which represents the wider. We can see that the conformal variant of STACI is the only one that has highly varied interval widths, while the UQ from other methods is fairly stagnant in the middle of purple and yellow. The high variance areas from conformal also correspond to areas with deep red high spikes in pollution which is commonly smoothed over in these classes of models. \n\n> Further consideration for the conclusions would help. Especially what to do next.\n\nWe thank the reviewer for highlighting this and further steps are important to extend this work. We believe that there are multiple possible extensions to non-Gaussian data. For example, adding a link function, or slightly changing the likelihood, should extend to distributions such as Bernoulli, or Poisson to model spatio-temporal correlated binary and count data. We also believe an extension to model spatial point processes is a feasible adaptation of the model. Following the work of Lloyd et al. (2015), the Poisson intensity can be modeled as the square of a GP. This would allow us to continue utilizing random Fourier features for efficient likelihood computation. These extensions would allow STACI to reach a wider audience and cover additional types of data.\n\n\n**Reference**\nLloyd, Chris, et al. \"Variational inference for Gaussian process modulated Poisson processes.\" International Conference on Machine Learning. PMLR, 2015."}, {"review_id": "UojuOKwKh2", "reviewer": "Reviewer_AwKK", "summary": "This paper proposes a novel methodology for estimating non-stationary random fields using feature construction, BNNs, and conformal predictions. The methods is well validated on two good examples and properly compared to SOTA methods in the field with promising results.", "questions": "- How would one perform model selection in this framework; as opposed to marginal likelihood maximization in GPs? Could this be elaborated on?\n\n- The proposed method is, in some parts, complex. Which parts of the model could be simplified if the inference task at hand is perhaps simpler (for example not none-stationary)?\n\n- What is the relation between (4) and Karhunen-Loeve expansions? If there is one, could this link be expanded on to help readers more familiar with these tools understand the proposed method?", "limitations": "Visual explanations of the use of conformal predictions could be helpful for readers less familiar with this class of methods.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 3, "originality": 4, "strengths_and_weaknesses": "Strengths:\n- The paper is well written, clear, and properly exposited. Shortcomings of other methods, and the proposed method, are well explained.\n- The explanation on conformal predictions is succinct and clear.\n- The numerics are well done and compare to many relevant methods.\n- Overall the methods is quite original and tackles a challenging problem in spatio-temporal inference.\n\n\nWeaknesses:\n- In the numerics, the metrics for assessing performance are poorly explained, i.e.: Coverage, Interval, Score, Interval Width.\n Could these be elaborated on and explained how they are computed as they would help the reader.\n-SVGD might suffer from mode-collapse in high dimensions, would this have a negative effect on the proposed methodology at test time on harder problems?", "was_revised": true, "final_justification": "Thank you to the authors for their time and effort. I believe this is a strong, well executed paper that solve an important problem in an effective and principled manner. I suggest acceptance.", "rebuttal": "We thank reviewer AwKK for taking the time to read through our work and provide valuable comments. We appreciate the praise for the clarity of our paper and originality of the methodology. We provide some rebuttals for individual comments here.\n\n> In the numerics, the metrics for assessing performance are poorly explained, i.e.: Coverage, Interval, Score, Interval Width. Could these be elaborated on and explained how they are computed as they would help the reader.\n\nWe agree that descriptions of these metrics are helpful for readers and elaborated on these metrics further in the Supplementary Material Section B.3.\n\n> SVGD might suffer from mode-collapse in high dimensions, would this have a negative effect on the proposed methodology at test time on harder problems?\n\n   You raise a great point on the dangers of mode-collapse for difficult high-dimensional problems. While we have not found this to be a problem for our tested datasets, we acknowledge this could happen in other use cases and would affect the UQ from the Bayesian variant of STACI. However, we believe the conformal addition would still provide valid intervals as the exchangeability calculations are based on the mean estimates of the covariance parameters and not their predictive variances.\n\n> How would one perform model selection in this framework; as opposed to marginal likelihood maximization in GPs? Could this be elaborated on?\n\n  Model selection would be done through cross validation in our methodology for the model fitting and conformal hyperparameters. Priors for the covariance parameters used in the kernel function can also be adjusted for a more direct connection to marginal likelihood maximization in GPs.\n\n> The proposed method is, in some parts, complex. Which parts of the model could be simplified if the inference task at hand is perhaps simpler (for example not none-stationary)?\n\nWe agree that this can be a complex framework. In a simpler stationary case, the need for the INR backbone is eliminated and this becomes a single hidden layer Bayesian neural network for the first half of our framework. Removal of the temporal component further reduces the number of priors set and model parameters needed to be learned.\n\n> What is the relation between (4) and Karhunen-Loeve expansions? If there is one, could this link be expanded on to help readers more familiar with these tools understand the proposed method?\n\nThe KL transformation is an alternative to the Fourier transformation used in the spectral approximation. Both represent a stochastic process as a linear combination of basis functions. The KL transformation is a generalization of the principle components analysis (PCA) where the basis functions are orthonormal and chosen to maximize variance explained. A Fourier transformation restricts the basis functions to be trigonometric functions. Like PCA, application of the KL decomposition typically requires replication of the process to estimate the basis functions, whereas the Fourier approximation can be applied to a single realization of the space-time process by assuming a stationary covariance (in latent space, in our case)."}]

Author rebuttal:

Meta-review: Three out of four reviewers gave very positive evaluations of this submission. Reviewer WrPZ raised a valid concern regarding an important missing baseline, which should be added in the final version. Based on the above, I recommend acceptance of this work.
Decision: Accept (poster)