Search

Modeling spatial processes that exhibit both smooth and rough features poses a significant challenge. This is especially true in fields where complex physical variables are observed across spatial domains. Traditional spatial techniques, such as Gaussian processes (GPs), are ill-suited to capture sharp transitions and discontinuities in spatial fields. In this paper, we propose a new approach incorporating non-Gaussian processes (NGPs) into a hybrid model which identifies both smooth and rough components. Specifically, we model the rough process using scaled mixtures of Gaussian distributions in a Bayesian hierarchical model (BHM). Our motivation comes from the Community Earth System Model Large Ensemble (CESM-LE), where we seek to emulate climate sensitivity fields that exhibit complex spatial patterns, including abrupt transitions at ocean-land boundaries. We demonstrate that traditional GP models fail to capture such abrupt changes and that our proposed hybrid model, implemented through a full Gibbs sampler. This significantly improves model interpretability and accurate recovery of process parameters. Through a multi-factor simulation study, we evaluate the performance of several scaled mixtures designed to model the rough process. The results highlight the advantages of using these heavier tailed priors as a replacement to the Bayesian fused LASSO. One prior in particular, the normal Jeffrey's prior stands above the rest. We apply our model to the CESM-LE dataset, demonstrating its ability to better represent the mean function and its uncertainty in climate sensitivity fields. This work combines the strengths of GPs for smooth processes with the flexibility of NGPs for abrupt changes. We provide a computationally efficient Gibbs sampler and include additional strategies for accelerating Monte Carlo Markov Chain (MCMC) sampling., Comment: 10 pages, 8 figures, 1 table

Bias correction is a common pre-processing step applied to climate model data before it is used for further analysis. This article introduces an efficient adaptation of a well-established bias-correction method - quantile mapping - for global horizontal irradiance (GHI) that ensures corrected data is physically plausible through incorporating measurements of clearsky GHI. The proposed quantile mapping method is fit on reanalysis data to first bias correct for regional climate models (RCMs) and is tested on RCMs forced by general circulation models (GCMs) to understand existing biases directly from GCMs. Additionally, we adapt a functional analysis of variance methodology that analyzes sources of remaining biases after implementing the proposed quantile mapping method and considered biases by climate region. This analysis is applied to four sets of climate model output from NA-CORDEX and compared against data from the National Solar Radiation Database produced by the National Renewable Energy Lab., Comment: 28 pages, 15 figures

In geostatistics, traditional spatial models often rely on the Gaussian Process (GP) to fit stationary covariances to data. It is well known that this approach becomes computationally infeasible when dealing with large data volumes, necessitating the use of approximate methods. A powerful class of methods approximate the GP as a sum of basis functions with random coefficients. Although this technique offers computational efficiency, it does not inherently guarantee a stationary covariance. To mitigate this issue, the basis functions can be "normalized" to maintain a constant marginal variance, avoiding unwanted artifacts and edge effects. This allows for the fitting of nearly stationary models to large, potentially non-stationary datasets, providing a rigorous base to extend to more complex problems. Unfortunately, the process of normalizing these basis functions is computationally demanding. To address this, we introduce two fast and accurate algorithms to the normalization step, allowing for efficient prediction on fine grids. The practical value of these algorithms is showcased in the context of a spatial analysis on a large dataset, where significant computational speedups are achieved. While implementation and testing are done specifically within the LatticeKrig framework, these algorithms can be adapted to other basis function methods operating on regular grids., Comment: Version 3

This study introduces a stochastic model of COVID-19 transmission tailored to the Colorado School of Mines campus and evaluates surveillance testing strategies within a university context. Enhancing the conventional SEIR framework with stochastic transitions, our model accounts for the unique characteristics of disease spread in a residential college, including specific states for testing, quarantine, and isolation. Employing an approximate Bayesian computation (ABC) method for parameter estimation, we navigate the complexities inherent in stochastic models, enabling an accurate fit of the model with the campus case data. We then studied our model under the estimated parameters to evaluate the efficacy of different testing policies that could be implemented on a university campus. This framework not only advances understanding of COVID-19 dynamics on the Mines campus but serves as a blueprint for comparable settings, providing insights for informed strategies against infectious diseases., Comment: Under review

Many industrial and engineering processes monitored as times series have smooth trends that indicate normal behavior and occasionally anomalous patterns that can indicate a problem. This kind of behavior can be modeled by a smooth trend, such as a spline or Gaussian process, and a disruption based on a sparser representation. Our approach is to expand the process signal into two sets of basis functions: one set uses L2 penalties on the coefficients, and the other set uses L1 penalties to control sparsity. From a frequentist perspective, this results in a hybrid smoother that combines cubic smoothing splines and the LASSO. As a Bayesian hierarchical model (BHM), this is equivalent to priors giving a Gaussian process and a Laplace distribution for anomaly coefficients. For the hybrid smoother, we propose two new ways of determining the penalty parameters that use effective degrees of freedom and contrast this with the BHM that uses loosely informative inverse gamma priors. Several reformulations are used to make sampling the BHM posterior more efficient, including some novel features in orthogonalizing and regularizing the model basis functions. This methodology is motivated by a substantive application, offline monitoring of a water treatment process for municipal water filtration. We also test the robustness of these methods with a Monte Carlo study designed to inspect a range trended time series under an array of conditions and compare this new approach to multiple existing modern methods. Both the hybrid smoother and the full BHM give comparable results with small false positive and false negative rates. Besides being successful in the water treatment application, this work can be easily extended to other Gaussian process models and other features that represent process disruptions in offline data., Comment: 32 pages, 17 figures, 2 table

The heavy-tailed behavior of the generalized extreme-value distribution makes it a popular choice for modeling extreme events such as floods, droughts, heatwaves, wildfires, etc. However, estimating the distribution's parameters using conventional maximum likelihood methods can be computationally intensive, even for moderate-sized datasets. To overcome this limitation, we propose a computationally efficient, likelihood-free estimation method utilizing a neural network. Through an extensive simulation study, we demonstrate that the proposed neural network-based method provides Generalized Extreme Value (GEV) distribution parameter estimates with comparable accuracy to the conventional maximum likelihood method but with a significant computational speedup. To account for estimation uncertainty, we utilize parametric bootstrapping, which is inherent in the trained network. Finally, we apply this method to 1000-year annual maximum temperature data from the Community Climate System Model version 3 (CCSM3) across North America for three atmospheric concentrations: 289 ppm $\mathrm{CO}_2$ (pre-industrial), 700 ppm $\mathrm{CO}_2$ (future conditions), and 1400 ppm $\mathrm{CO}_2$, and compare the results with those obtained using the maximum likelihood approach., Comment: 19 pages, 6 figures

Initial steps in statistical downscaling involve being able to compare observed data from regional climate models (RCMs). This prediction requires (1) regridding RCM output from their native grids and at differing spatial resolutions to a common grid in order to be comparable to observed data and (2) bias correcting RCM data, via quantile mapping, for example, for future modeling and analysis. The uncertainty associated with (1) is not always considered for downstream operations in (2). This work examines this uncertainty, which is not often made available to the user of a regridded data product. This analysis is applied to RCM solar radiation data from the NA-CORDEX data archive and observed data from the National Solar Radiation Database housed at the National Renewable Energy Lab. A case study of the mentioned methods over California is presented., Comment: 16 pages, 5 figures, submitted to: Advances in Statistical Climatology, Meteorology and Oceanography

Computing an ensemble of random fields using conditional simulation is an ideal method for retrieving accurate estimates of a field conditioned on available data and for quantifying the uncertainty of these realizations. Methods for generating random realizations, however, are computationally demanding, especially when the estimates are conditioned on numerous observed data and for large domains. In this article, a \textit{new}, \textit{approximate} conditional simulation approach is applied that builds on \textit{circulant embedding} (CE), a fast method for simulating stationary Gaussian processes. The standard CE is restricted to simulating stationary Gaussian processes (possibly anisotropic) on regularly spaced grids. In this work we explore two possible algorithms, namely local Kriging and nearest neighbor Kriging, that extend CE for irregularly spaced data points. We establish the accuracy of these methods to be suitable for practical inference and the speedup in computation allows for generating conditional fields close to an interactive time frame. The methods are motivated by the U.S. Geological Survey's software \textit{ShakeMap}, which provides near real-time maps of shaking intensity after the occurrence of a significant earthquake. An example for the 2019 event in Ridgecrest, California is used to illustrate our method.

Gaussian processes (GPs) are a popular model for spatially referenced data and allow descriptive statements, predictions at new locations, and simulation of new fields. Often a few parameters are sufficient to parameterize the covariance function, and maximum likelihood (ML) methods can be used to estimate these parameters from data. ML methods, however, are computationally demanding. For example, in the case of local likelihood estimation, even fitting covariance models on modest size windows can overwhelm typical computational resources for data analysis. This limitation motivates the idea of using neural network (NN) methods to approximate ML estimates. We train NNs to take moderate size spatial fields or variograms as input and return the range and noise-to-signal covariance parameters. Once trained, the NNs provide estimates with a similar accuracy compared to ML estimation and at a speedup by a factor of 100 or more. Although we focus on a specific covariance estimation problem motivated by a climate science application, this work can be easily extended to other, more complex, spatial problems and provides a proof-of-concept for this use of machine learning in computational statistics., Comment: 12 pages, 5 figures

Surface registration, the task of aligning several multidimensional point sets, is a necessary task in many scientific fields. In this work, a novel statistical approach is developed to solve the problem of nonrigid registration. While the application of an affine transformation results in rigid registration, using a general nonlinear function to achieve nonrigid registration is necessary when the point sets require deformations that change over space. The use of a local likelihood-based approach using windowed Gaussian processes provides a flexible way to accurately estimate the nonrigid deformation. This strategy also makes registration of massive data sets feasible by splitting the data into many subsets. The estimation results yield spatially-varying local rigid registration parameters. Gaussian process surface models are then fit to the parameter fields, allowing prediction of the transformation parameters at unestimated locations, specifically at observation locations in the unregistered data set. Applying these transformations results in a global, nonrigid registration. A penalty on the transformation parameters is included in the likelihood objective function. Combined with smoothing of the local estimates from the surface models, the nonrigid registration model can prevent the problem of overfitting. The efficacy of the nonrigid registration method is tested in two simulation studies, varying the number of windows and number of points, as well as the type of deformation. The nonrigid method is applied to a pair of massive remote sensing elevation data sets exhibiting complex geological terrain, with improved accuracy and uncertainty quantification in a cross validation study versus two rigid registration methods.

Modeling data with non-stationary covariance structure is important to represent heterogeneity in geophysical and other environmental spatial processes. In this work, we investigate a multistage approach to modeling non-stationary covariances that is efficient for large data sets. First, we use likelihood estimation in local, moving windows to infer spatially varying covariance parameters. These surfaces of covariance parameters can then be encoded into a global covariance model specifying the second-order structure for the complete spatial domain. The resulting global model allows for efficient simulation and prediction. We investigate the non-stationary spatial autoregressive (SAR) model related to Gaussian Markov random field (GMRF) methods, which is amenable to plug in local estimates and practical for large data sets. In addition we use a simulation study to establish the accuracy of local Mat\'ern parameter estimation as a reliable technique when replicate fields are available and small local windows are exploited to reduce computation. This multistage modeling approach is implemented on a non-stationary climate model output data set with the goal of emulating the variation in the numerical model ensemble using a Gaussian process.

In many practical applications, spatial data are often collected at areal levels (i.e., block data) and the inferences and predictions about the variable at points or blocks different from those at which it has been observed typically depend on integrals of the underlying continuous spatial process. In this paper we describe a method based on Fourier transform by which multiple integrals of covariance functions over irregular data regions may be numerically approximated with the same level of accuracy to traditional methods, but at a greatly reduced computational expense., Comment: 14 pages, 7 figures, 7 tables. Submitted to Stat (Wiley Online Library) in December 2019

Gaussian process (GP) models are widely used to analyze spatially referenced data and to predict values at locations without observations. In contrast to many algorithmic procedures, GP models are based on a statistical framework, which enables uncertainty quantification of the model structure and predictions. Both the evaluation of the likelihood and the prediction involve solving linear systems. Hence, the computational costs are large and limit the amount of data that can be handled. While there are many approximation strategies that lower the computational cost of GP models, they often provide only sub-optimal support for the parallel computing capabilities of current (high-performance) computing environments. We aim at bridging this gap with a parameter estimation and prediction method that is designed to be parallelizable. More precisely, we divide the spatial domain into overlapping subsets and use cross-validation (CV) to estimate the covariance parameters in parallel. We present simulation studies, which assess the accuracy of the parameter estimates and predictions. Moreover, we show that our implementation has good weak and strong parallel scaling properties. For illustration, we fit an exponential covariance model to a scientifically relevant canopy height dataset with 5 million observations. Using 512 processor cores in parallel brings the evaluation time of one covariance parameter configuration to less than 1.5 minutes. The parallel CV method can be easily extended to include approximate likelihood methods, multivariate and spatio-temporal data, as well as non-stationary covariance models., Comment: 11 pages, 6 Figs. Under review

One of the commonly used approaches to modeling extremes is the peaks-over-threshold (POT) method. The POT method models exceedances over a threshold that is sufficiently high or low so that the exceedance has approximately a generalized Pareto distribution (GPD). This method requires the selection of a threshold that might affect the estimates. Here we propose an alternative method, the Log-Histospline (LHSpline), to explore modeling the tail behavior and the remainder of the density in one step using the full range of the data. LHSpline applies a smoothing spline model to a finely binned histogram of the log transformed data to estimate its log density. By construction, a LHSpline estimation is constrained to have polynomial tail behavior, a feature commonly observed in daily rainfall observations. We illustrate the LHSpline method by analyzing precipitation data collected in Houston, Texas., Comment: 32 pages, 13 figures, 2 table

In spatial statistics, a common method for prediction over a Gaussian random field (GRF) is maximum likelihood estimation combined with kriging. For massive data sets, kriging is computationally intensive, both in terms of CPU time and memory, and so fixed rank kriging has been proposed as a solution. The method however still involves operations on large matrices, so we develop an alteration to this method by utilizing the approximations made in fixed rank kriging combined with restricted maximum likelihood estimation and sparse matrix methodology. Experiments show that our methodology can provide additional gains in computational efficiency over fixed-rank kriging without loss of accuracy in prediction. The methodology is applied to climate data archived by the United States National Climate Data Center, with very good results., Comment: Sankhya, Series A, accepted for publication

Geophysical and other natural processes often exhibit non-stationary covariances and this feature is important to take into account for statistical models that attempt to emulate the physical process. A convolution-based model is used to represent non-stationary Gaussian processes that allows for variation in the correlation range and vari- ance of the process across space. Application of this model has two steps: windowed estimates of the covariance function under the as- sumption of local stationary and encoding the local estimates into a single spatial process model that allows for efficient simulation. Specifically we give evidence to show that non-stationary covariance functions based on the Mat`ern family can be reproduced by the Lat- ticeKrig model, a flexible, multi-resolution representation of Gaussian processes. We propose to fit locally stationary models based on the Mat`ern covariance and then assemble these estimates into a single, global LatticeKrig model. One advantage of the LatticeKrig model is that it is efficient for simulating non-stationary fields even at 105 locations. This work is motivated by the interest in emulating spatial fields derived from numerical model simulations such as Earth system models. We successfully apply these ideas to emulate fields that de- scribe the uncertainty in the pattern scaling of mean summer (JJA) surface temperature from a series of climate model experiments. This example is significant because it emulates tens of thousands of loca- tions, typical in geophysical model fields, and leverages embarrassing parallel computation to speed up the local covariance fitting, Comment: 32 pages total, 10 figures

The Gaussian process is an indispensable tool for spatial data analysts. The onset of the "big data" era, however, has lead to the traditional Gaussian process being computationally infeasible for modern spatial data. As such, various alternatives to the full Gaussian process that are more amenable to handling big spatial data have been proposed. These modern methods often exploit low rank structures and/or multi-core and multi-threaded computing environments to facilitate computation. This study provides, first, an introductory overview of several methods for analyzing large spatial data. Second, this study describes the results of a predictive competition among the described methods as implemented by different groups with strong expertise in the methodology. Specifically, each research group was provided with two training datasets (one simulated and one observed) along with a set of prediction locations. Each group then wrote their own implementation of their method to produce predictions at the given location and each which was subsequently run on a common computing environment. The methods were then compared in terms of various predictive diagnostics. Supplementary materials regarding implementation details of the methods and code are available for this article online.

The fundamental equations that model turbulent flow do not provide much insight into the size and shape of observed turbulent structures. We investigate the efficient and accurate representation of structures in two-dimensional turbulence by applying statistical models directly to the simulated vorticity field. Rather than extract the coherent portion of the image from the background variation, as in the classical signal-plus-noise model, we present a model for individual vortices using the non-decimated discrete wavelet transform. A template image, supplied by the user, provides the features to be extracted from the vorticity field. By transforming the vortex template into the wavelet domain, specific characteristics present in the template, such as size and symmetry, are broken down into components associated with spatial frequencies. Multivariate multiple linear regression is used to fit the vortex template to the vorticity field in the wavelet domain. Since all levels of the template decomposition may be used to model each level in the field decomposition, the resulting model need not be identical to the template. Application to a vortex census algorithm that records quantities of interest (such as size, peak amplitude, circulation, etc.) as the vorticity field evolves is given. The multiresolution census algorithm extracts coherent structures of all shapes and sizes in simulated vorticity fields and is able to reproduce known physical scaling laws when processing a set of voriticity fields that evolve over time.