Generalized Covariance Function Research Articles

Linking physics and spatial statistics: A new family of Boltzmann-Gibbs random fields

We investigate a connection between spatial statistics and statistical physics to obtain new covariance functions with direct physical interpretation for spatial random fields. These covariance functions are based on the exponential Boltzmann-Gibbs representation and use an energy functional to represent interactions between the values of the random field at different points in space. This formulation results in closed-form generalized covariance functions, which display infinite variance in Euclidean spaces of dimension larger than one. We propose regularization schemes in real and reciprocal (spectral) space that lead to well-behaved covariance structures. The real-space regularization parameter allows a continuous interpolation between the Boltzmann-Gibbs covariance and the exponential covariance. We also propose discretized approximations on regular grids, and we show that they represent reparametrized versions of the well-known Besag and Leroux lattice models. We then discuss parameter estimation and spatial prediction for the regularized Boltzmann-Gibbs covariance model in two dimensions. We recommend using the pairwise difference likelihood that combines satisfactory estimation performance and good scalability with many observation points. The predictive performance of the regularized covariance function is assessed by means of cross-validation statistics. Irregularly-spaced samples from the Walker Lake dataset are used, and spatial prediction is conducted by means of ordinary kriging. The regularized Boltzmann-Gibbs covariance yields improved predictive performance compared to the exponential covariance model.

A Bayesian nonparametric approach to option pricing

Accurately modeling the implied volatility surface is of great importance to option pricing, trading and hedging. In this paper, we investigate the use of a Bayesian nonparametric approach to fit and forecast the implied volatility surface with observed market data. More specifically, we explore Gaussian Processes with different kernel functions characterizing general covariance functions. We also obtain posterior distributions of the implied volatility and build confidence intervals for the predictions to assess potential model uncertainty. We apply our approach to market data on the S&P 500 index option market in 2018, analyzing 322,983 options. Our results suggest that the Bayesian approach is a powerful alternative to existing parametric pricing models

Easy-But-Effective Domain Sub-Similarity Learning for Transfer Regression

Transfer covariance function, which can model domain similarity and adaptively control the knowledge transfer across domains, is widely used in transfer learning. In this paper, we concentrate on Gaussian process ( <i>GP</i> ) models using a transfer covariance function for regression problems in a black-box learning scenario. Precisely, we investigate a family of rather general transfer covariance functions, <inline-formula><tex-math notation="LaTeX">${T}_{*}$</tex-math></inline-formula> , that can model the heterogeneous sub-similarities of domains through multiple kernel learning. A necessary and sufficient condition to obtain valid <i>GP</i> s using <inline-formula><tex-math notation="LaTeX">${T}_{*}$</tex-math></inline-formula> ( <inline-formula><tex-math notation="LaTeX">$GP_{T_{*}}$</tex-math></inline-formula> ) for any data is given. This condition becomes specially handy for practical applications as (i) it enables semantic interpretations of the sub-similarities and (ii) it can readily be used for model learning. In particular, we propose a computationally inexpensive model learning rule that can explicitly capture different sub-similarities of domains. We propose two instantiations of <inline-formula><tex-math notation="LaTeX">$GP_{T_{*}}$</tex-math></inline-formula> , one with a set of predefined constant base kernels and one with a set of learnable parametric base kernels. Extensive experiments on 36 synthetic transfer tasks and 12 real-world transfer tasks demonstrate the effectiveness of <inline-formula><tex-math notation="LaTeX">$GP_{T_{*}}$</tex-math></inline-formula> on the sub-similarity capture and the transfer performance.

Composite likelihood estimation for a Gaussian process under fixed domain asymptotics

We study the problem of estimating the covariance parameters of a one-dimensional Gaussian process with exponential covariance function under fixed-domain asymptotics. We show that the weighted pairwise maximum likelihood estimator of the microergodic parameter can be consistent or inconsistent. This depends on the range of admissible parameter values in the likelihood optimization. On the other hand, the weighted pairwise conditional maximum likelihood estimator is always consistent. Both estimators are also asymptotically Gaussian when they are consistent. Their asymptotic variances are larger or strictly larger than that of the maximum likelihood estimator. A simulation study is presented in order to compare the finite sample behavior of the pairwise likelihood estimators with their asymptotic distributions. For more general covariance functions, an additional inconsistency result is provided, for the weighted pairwise maximum likelihood estimator of a variance parameter.

On the Choice of Momentum Control Variables and Covariance Modeling for Mesoscale Data Assimilation

Abstract For mesoscale variational data assimilation with high-resolution observations, there has been an issue concerning the choice of momentum control variables and related covariance modeling. This paper addresses the theoretical aspect of this issue. First, relationships between background error covariance functions for differently chosen momentum control variables are derived, and different choices of momentum control variables are proven to be theoretically equivalent in the sense that they lead to the same optimally analyzed incremental wind field in the limit of infinitely high spatial resolution provided their error covariance functions satisfy the derived relationships. It is then shown that when the velocity potential χ and streamfunction ψ are used as momentum control variables with their background error autocovariance functions modeled by single-Gaussian functions, the derived velocity autocovariance functions contain significant negative sidelobes. These negative sidelobes can represent background wind error structures associated with baroclinic waves on the synoptic scale but become unrepresentative on the mesoscale. To reduce or remove these negative sidelobes for mesoscale variational data assimilation, Gaussian functions are used with two types of modifications to model the velocity covariance functions in consistency with the assumed homogeneity and isotropy in variational data assimilation. In this case, the random (χ, ψ) background error fields have no classically valid homogeneous and isotropic covariance functions, but generalized (χ, ψ) covariance functions can be derived from the modified velocity covariance functions for choosing (χ, ψ) as momentum control variables. Mathematical properties of generalized covariance functions are explored with physical interpretations. Their important implications are discussed for mesoscale data assimilation.

Simulation of intrinsic random fields of order k with a continuous spectral algorithm

Intrinsic random fields of order k, defined as random fields whose high-order increments (generalized increments of order k) are second-order stationary, are used in spatial statistics to model regionalized variables exhibiting spatial trends, a feature that is common in earth and environmental sciences applications. A continuous spectral algorithm is proposed to simulate such random fields in a d-dimensional Euclidean space, with given generalized covariance structure and with Gaussian generalized increments of order k. The only condition needed to run the algorithm is to know the spectral measure associated with the generalized covariance function (case of a scalar random field) or with the matrix of generalized direct and cross-covariances (case of a vector random field). The algorithm is applied to synthetic examples to simulate intrinsic random fields with power generalized direct and cross-covariances, as well as an intrinsic random field with power and spline generalized direct covariances and Matern generalized cross-covariance.

Spatiotemporal interpolation and forecast of irradiance data using Kriging

Solar power variability is a concern to grid operators as unanticipated changes in PV plant power output can strain the electric grid. The main cause of solar variability is clouds passing over PV modules. However, geographic diversity across a region leads to a reduction in the cloud-induced variability. In this paper, spatiotemporal correlations of irradiance data are analyzed and spatial and spatiotemporal ordinary Kriging methods are applied to model irradiation at an arbitrary point based on the given time series of irradiation at some observed locations. The correlations among the irradiances at observed locations are modeled by general parametric covariance functions. Besides the isotropic covariance function (which is independent of direction), a new non-separable anisotropic parametric covariance function is proposed to model the transient clouds. Also, a new approach is proposed to estimate the spatial and temporal decorrelation distances analytically using the applied parametric covariance functions, which reduce the size of the computations without loss in accuracy (parameter shrinkage). The analysis has been performed and the Kriging method is first validated by using two spatially and temporally resolved artificial irradiance datasets generated from Large Eddy Simulation. Then, the spatiotemporal Kriging method is applied on real irradiance and output power data in California (Sacramento and San Diego areas) where the cloud motion had to be estimated during the process using cross-correlation method (CCM). Results confirm that the anisotropic model is most accurate with an average normalized root mean squared error (nRMSE) of 7.92% representing a 66% relative improvement over the persistence model.

Sampling image segmentations for uncertainty quantification.

In this paper, we introduce a method to automatically produce plausible image segmentation samples from a single expert segmentation. A probability distribution of image segmentation boundaries is defined as a Gaussian process, which leads to segmentations which are spatially coherent and consistent with the presence of salient borders in the image. The proposed approach is computationally efficient, and generates visually plausible samples. The variability between the samples is mainly governed by a parameter which may be correlated with a simple Dice's coefficient, or easily set by the user from the definition of probable regions of interest. The method is extended to the case of several neighboring structures, but also to account for under or over segmentation, and the presence of excluded regions. We also detail a method to sample segmentations with more general non-stationary covariance functions which relies on supervoxels. Furthermore, we compare the generated segmentation samples with several manual clinical segmentations of a brain tumor. Finally, we show how this approach can have useful applications in the field of uncertainty quantification, and an illustration is provided in radiotherapy planning, where segmentation sampling is applied to both the clinical target volume and the organs at risk.

Simulation of Intrinsic Random Fields of Order $$k$$ k with Gaussian Generalized Increments by Gibbs Sampling

This work pertains to the simulation of an intrinsic random field of order $$k$$ with a given generalized covariance function and multivariate Gaussian generalized increments. An iterative algorithm based on the Gibbs sampler is proposed to simulate such a random field at a finite set of locations, resulting in a sequence of random vectors that converges in distribution to a random vector with the desired distribution. The algorithm is tested on synthetic case studies to experimentally assess its rate of convergence, showing that few iterations are sufficient for convergence to take place. The sequence of random vectors also proves to be strongly mixing, allowing the generation of as many independent realizations as desired with a single run of the algorithm. Another interesting property of this algorithm is its versatility, insofar as it can be adapted to construct realizations conditioned to pre-existing data and can be used for any number and configuration of the target locations and any generalized covariance model.

On a class of space–time intrinsic random functions

Power law generalized covariance functions provide a simple model for describing the local behavior of an isotropic random field. This work seeks to extend this class of covariance functions to spatial-temporal processes for which the degree of smoothness in space and in time may differ while maintaining other desirable properties for the covariance functions, including the availability of explicit convergent and asymptotic series expansions.

A Resampling-Based Stochastic Approximation Method for Analysis of Large Geostatistical Data

The Gaussian geostatistical model has been widely used in modeling of spatial data. However, it is challenging to computationally implement this method because it requires the inversion of a large covariance matrix, particularly when there is a large number of observations. This article proposes a resampling-based stochastic approximation method to address this challenge. At each iteration of the proposed method, a small subsample is drawn from the full dataset, and then the current estimate of the parameters is updated accordingly under the framework of stochastic approximation. Since the proposed method makes use of only a small proportion of the data at each iteration, it avoids inverting large covariance matrices and thus is scalable to large datasets. The proposed method also leads to a general parameter estimation approach, maximum mean log-likelihood estimation, which includes the popular maximum (log)-likelihood estimation (MLE) approach as a special case and is expected to play an important role in analyzing large datasets. Under mild conditions, it is shown that the estimator resulting from the proposed method converges in probability to a set of parameter values of equivalent Gaussian probability measures, and that the estimator is asymptotically normally distributed. To the best of the authors’ knowledge, the present study is the first one on asymptotic normality under infill asymptotics for general covariance functions. The proposed method is illustrated with large datasets, both simulated and real. Supplementary materials for this article are available online.

Efficient methods for large‐scale linear inversion using a geostatistical approach

In geophysical inverse problems, such as estimating the unknown parameter field from noisy observations of dependent quantities, e.g., hydraulic conductivity from head observations, stochastic Bayesian and geostatistical approaches are frequently used. To obtain best estimates and conditional realizations it is required to perform several matrix‐matrix computations involving the covariance matrix of the discretized field of the parameters. In realistic three‐dimensional fields that are finely discretized, these operations as performed in conventional algorithms become extremely expensive and even prohibitive in terms of memory and computational requirements. Using Hierarchical Matrices, we show how to reduce the complexity of forming approximate matrix‐vector products involving the Covariance matrices in log linear complexity for an arbitrary distribution of points and a wide variety of generalized covariance functions. The resulting system of equations is solved iteratively using a matrix‐free Krylov subspace approach. Furthermore, we show how to generate unconditional realizations using an approximation to the square root of the covariance matrix using Chebyshev matrix polynomials and use the above to generate conditional realizations. We demonstrate the efficiency of our method on a few standard test problems, such as interpolation from noisy observations and contaminant source identification.

Application of FFT-based Algorithms for Large-Scale Universal Kriging Problems

Looking at kriging problems with huge numbers of estimation points and measurements, computational power and storage capacities often pose heavy limitations to the maximum manageable problem size. In the past, a list of FFT-based algorithms for matrix operations have been developed. They allow extremely fast convolution, superposition and inversion of covariance matrices under certain conditions. If adequately used in kriging problems, these algorithms lead to drastic speedup and reductions in storage requirements without changing the kriging estimator. However, they require second-order stationary covariance functions, estimation on regular grids, and the measurements must also form a regular grid. In this study, we show how to alleviate these rather heavy and many times unrealistic restrictions. Stationarity can be generalized to intrinsicity and beyond, if decomposing kriging problems into the sum of a stationary problem and a formally decoupled regression task. We use universal kriging, because it covers arbitrary forms of unknown drift and all cases of generalized covariance functions. Even more general, we use an extension to uncertain rather than unknown drift coefficients. The sampling locations may now be irregular, but must form a subset of the estimation grid. Finally, we present asymptotically exact but fast approximations to the estimation variance and point out application to conditional simulation, cokriging and sequential kriging. The drastic gain in computational and storage efficiency is demonstrated in test cases. Especially high-resolution and data-rich fields such as rainfall interpolation from radar measurements or seismic or other geophysical inversion can benefit from these improvements.

Generalized analysis of spatial variation in yield monitor data

Australian lupin and cotton yield monitor data were analysed using spatial models from the Matérn class of spatial covariance functions. Despite difficulties with the spatial disposition of the data, the analysis supports the statistical model in which the variation is a linear combination of white noise and the de Wijs process. The de Wijs process, also called the logarithmic covariance function, is a generalized covariance function that is conformally invariant and suggests that there is variation at all spatial scales. The present work also indicates that anisotropy and convolution are properties of yield monitor data and that it is hard to distinguish the two. The degree and causes of anisotropy require further investigation. Fitting this model is relatively easy for small, precision-agriculture datasets and open source software is available to this end. Comparing the de Wijs model with more general models in the Matérn class is computationally intensive for precision-agriculture datasets.

An Application of Bayesian Inverse Methods to Vertical Deconvolution of Hydraulic Conductivity in a Heterogeneous Aquifer at Oak Ridge National Laboratory

A Bayesian inverse method is applied to two electromagnetic flowmeter tests conducted in fractured weathered shale at Oak Ridge National Laboratory. Traditional deconvolution of flowmeter tests is also performed using a deterministic first-difference approach; furthermore, ordinary kriging was applied on the first-difference results to provide an additional method yielding the best estimate and confidence intervals. Depth-averaged bulk hydraulic conductivity information was available from previous testing. The three methods deconvolute the vertical profile of lateral hydraulic conductivity. A linear generalized covariance function combined with a zoning approach was used to describe structure. Nonnegativity was enforced by using a power transformation. Data screening prior to calculations was critical to obtaining reasonable results, and the quantified uncertainty estimates obtained by the inverse method led to the discovery of questionable data at the end of the process. The best estimates obtained using the inverse method and kriging compared favorably with first-difference confirmatory calculations, and all three methods were consistent with the geology at the site.

Local stationarity and simulation of self-affine intrinsic random functions

Gaussian intrinsic random functions with power law generalized covariance functions, which in one dimension essentially correspond to fractional and integrated fractional Brownian motions, form a class of self-affine models for random fields with a wide range of smoothness properties. These random fields are nonstationary, but appropriately filtered versions of them are stationary. This work proves that most such random functions are locally stationary in a certain well-defined sense. This result yields an efficient and exact method of simulating all fractional and integrated fractional Brownian motions.

Using AVHRR images for spatial prediction of clay content in the lower Namoi Valley of eastern Australia

Soil information is needed at the regional scale to enable planning of land utilization in accordance with its capacity. Because existing soil maps are inadequate in Australia to meet this demand, there is the need to develop models that could be used to improve soil maps at this scale for aggregation up to the national or continental scale. The most efficient and cheapest means of achieving this is by using remotely sensed data in multivariate spatial prediction models. This study therefore examines the soil spectral properties as depicted by the National Oceanic and Atmospheric Administration (NOAA) advanced very high resolution radiometer (AVHRR) data, with the main aim of developing prediction models for improved mapping of, for example, topsoil % clay in the Lower Namoi Valley of New South Wales (NSW). The paper compares several prediction models: multiple linear regression (MLR) using an external training set (MLR-ETS), interpolation by MLR — MLR-INT, kriging based on a generalized covariance function of order 1 (IRF-1), and a mixed model of MLR and ordinary kriging, termed as regression/kriging (RK). Comparison was based on an independent validation set ( N=40), using the root mean square error (RMSE) of prediction. The MLR-ETS performed very poorly (RMSE=35.0%), due to some degree of contrariety of the external training data with the data from the prediction area. The RK is superior to all the methods in predicting the topsoil % clay with RMSE of 10.2%. This performance by RK, compared with MLR-ITS (RMSE=13.3%) and IRF-1 (RMSE=12.6%), is quite remarkable at the regional scale of consideration. The correlation coefficient ( ρ) between the actual and predicted % clay confirms the order of prediction performance. Isarithmic maps of the topsoil % clay, as predicted by the best two methods, indicate that IRF-1 over-smoothed the predicted clay values as it removed both the low and high spikes in the spatial distribution of the % clay. However, the map of predicted topsoil clay by RK, which incorporated the AVHRR bands and indices in the prediction model reflects most of the local variability. Thus this study has demonstrated that a few soil-sampled sites, combined with the AVHRR data, could be adequate for regional soil inventory of good quality and known precision using RK. The basic tenet of this study can be extended to any situation where ancillary attributes have relatively high correlation with soil variables.

Approximation of harmonic covariance functions on the sphere by non-harmonic locally supported functions

Three methods to construct positive definite functions with compact support for the approximation of general geophysical harmonic covariance functions are presented. The theoretical background is given and simulations carried out, for three types of covariance functions associated with the determination of the anomalous gravity potential from gravity anomalies. The results are compared with those of the finite covariance function of Sanso and Schuh (1987).

Generalized covariance functions associated with the Laplace Equation and Their use in interpolation and inverse problems

This work examines which generalized covariance function when used in the stochastic approach produces the flattest possible estimate of an unknown function that is consistent with the data. Such an estimate is the plainest possible continuous function, thus in a sense eliminating details that are irrelevant or unsupported by data. The answer is found from the solution of the following variational problem: Determine the function that reproduces the data, has the smallest gradient (in the square norm sense), and has a gradient that vanishes at large distances from the observations. The generalized covariance functions are shown to be the Green's functions for the free‐space Laplace equation: the linear distance, in one dimension; the logarithmic distance in two dimensions; and the inverse distance in three dimensions. It is demonstrated that they are appropriate covariance functions for intrinsic random fields, a modification is proposed to facilitate numerical implementation, and a couple of examples are presented to illustrate the applicability of the methodology.

GCINFE: A computer program for inference of polynomial generalized covariance functions

GCINFE, an ANSI FORTRAN-77 coded computer program for statistical inference of polynomial generalized covariance parameters, is examined in this paper. The generalized covariance function describes a nonstationary random field within the framework of the theory of intrinsic random functions of order k. The program utilizes four different estimation methods: minimum norm, minimum variance unbiased quadratic, minimum variance unbiased quadratic iterative and maximum likelihood. An important feature of these methods is that they provide a measure of uncertainty of the estimated parameters as defined by the standard deviation of the estimation error. The sampling distribution statistics of the four different methods evaluated through the Monte Carlo simulation technique are also presented for an example of the generalized covariance model.