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

This paper introduces two innovative convex optimization formulations to simultaneously optimize the H2/Hinf observer gain and sensing precision, and guarantee a specified estimation error bound for nonlinear systems in LPV form. Applied to the design of an onboard celestial navigation system for cislunar operations, these formulations demonstrate the ability to maintain accurate spacecraft positioning with minimal measurements and theoretical performance guarantees by design.

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

Global and regional climate model projections are useful for gauging future patterns of climate variables, including solar radiation, but data from these models is often too coarse to assess local impacts. Within the context of solar radiation, the changing climate may have an effect on photovoltaic (PV) production, especially as the PV industry moves to extend plant lifetimes to 50 years. Predicting PV production while taking into account a changing climate requires data at a resolution that is useful for building PV plants. Although temporal and spatial downscaling of solar radiation data is widely studied, we present a novel method to downscale solar radiation data from daily averages to hourly profiles, while maintaining spatial correlation of parameters characterizing the diurnal profile of solar radiation. The method focuses on the use of a diurnal template which can be shifted and scaled according to the time or year and location and the use of thin plate splines for spatial downscaling. This analysis is applied to data from the National Solar Radiation Database housed at the National Renewable Energy Lab and a case study of the mentioned methods over several sub-regions of continental United States is presented., Comment: 35 pages, 14 figures

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

In this paper, we address the issue of quantifying maximum actuator degradation in linear time-invariant dynamical systems. We present a new unified framework for computing the state-feedback controller gain that meets a user-defined closed-loop performance criterion while also maximizing actuator degradation. This degradation is modeled as a first-order filter with additive noise. Our approach involves two novel convex optimization formulations that concurrently determine the controller gain, maximize actuator degradation, and maintain the desired closed-loop performance in both the $H_2$ and $H_{\infty}$ system norms. The results are limited to open-loop stable systems. We demonstrate the application of our results through the design of a full-state feedback controller for a model representing the longitudinal motion of the F-16 aircraft., Comment: 6 pages, 6 figures

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

Deriving the strength and direction of the three-dimensional (3D) magnetic field in the solar atmosphere is fundamental for understanding its dynamics. Volume information on the magnetic field mostly relies on coupling 3D reconstruction methods with photospheric and/or chromospheric surface vector magnetic fields. Infrared coronal polarimetry could provide additional information to better constrain magnetic field reconstructions. However, combining such data with reconstruction methods is challenging, e.g., because of the optical-thinness of the solar corona and the lack and limitations of stereoscopic polarimetry. To address these issues, we introduce the Data-Optimized Coronal Field Model (DOCFM) framework, a model-data fitting approach that combines a parametrized 3D generative model, e.g., a magnetic field extrapolation or a magnetohydrodynamic model, with forward modeling of coronal data. We test it with a parametrized flux rope insertion method and infrared coronal polarimetry where synthetic observations are created from a known "ground truth" physical state. We show that this framework allows us to accurately retrieve the ground truth 3D magnetic field of a set of force-free field solutions from the flux rope insertion method. In observational studies, the DOCFM will provide a means to force the solutions derived with different reconstruction methods to satisfy additional, common, coronal constraints. The DOCFM framework therefore opens new perspectives for the exploitation of coronal polarimetry in magnetic field reconstructions and for developing new techniques to more reliably infer the 3D magnetic fields that trigger solar flares and coronal mass ejections., Comment: 14 pages, 6 figures; Accepted for publication in ApJ

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 Coronal Multichannel Polarimeter (CoMP) routinely performs coronal polarimetric measurements using the Fe XIII 10747 $\AA$ and 10798 $\AA$ lines, which are sensitive to the coronal magnetic field. However, inverting such polarimetric measurements into magnetic field data is a difficult task because the corona is optically thin at these wavelengths and the observed signal is therefore the integrated emission of all the plasma along the line of sight. To overcome this difficulty, we take on a new approach that combines a parameterized 3D magnetic field model with forward modeling of the polarization signal. For that purpose, we develop a new, fast and efficient, optimization method for model-data fitting: the Radial-basis-functions Optimization Approximation Method (ROAM). Model-data fitting is achieved by optimizing a user-specified log-likelihood function that quantifies the differences between the observed polarization signal and its synthetic/predicted analogue. Speed and efficiency are obtained by combining sparse evaluation of the magnetic model with radial-basis-function (RBF) decomposition of the log-likelihood function. The RBF decomposition provides an analytical expression for the log-likelihood function that is used to inexpensively estimate the set of parameter values optimizing it. We test and validate ROAM on a synthetic test bed of a coronal magnetic flux rope and show that it performs well with a significantly sparse sample of the parameter space. We conclude that our optimization method is well-suited for fast and efficient model-data fitting and can be exploited for converting coronal polarimetric measurements, such as the ones provided by CoMP, into coronal magnetic field data., Comment: 23 pages, 12 figures, accepted in Frontiers in Astronomy and Space Sciences

Discussion of "A statistical analysis of multiple temperature proxies: Are reconstructions of surface temperatures over the last 1000 years reliable?" by B.B. McShane and A.J. Wyner [arXiv:1104.4002], Comment: Published in at http://dx.doi.org/10.1214/10-AOAS398K the Annals of Applied Statistics (http://www.imstat.org/aoas/) by the Institute of Mathematical Statistics (http://www.imstat.org)

Nonlinear regression is a useful statistical tool, relating observed data and a nonlinear function of unknown parameters. When the parameter-dependent nonlinear function is computationally intensive, a straightforward regression analysis by maximum likelihood is not feasible. The method presented in this paper proposes to construct a faster running surrogate for such a computationally intensive nonlinear function, and to use it in a related nonlinear statistical model that accounts for the uncertainty associated with this surrogate. A pivotal quantity in the Earth's climate system is the climate sensitivity: the change in global temperature due to doubling of atmospheric $\mathrm{CO}_2$ concentrations. This, along with other climate parameters, are estimated by applying the statistical method developed in this paper, where the computationally intensive nonlinear function is the MIT 2D climate model., Comment: Published in at http://dx.doi.org/10.1214/08-AOAS210 the Annals of Applied Statistics (http://www.imstat.org/aoas/) by the Institute of Mathematical Statistics (http://www.imstat.org)