Fitting a Powered Exponential Autocorrelation: Alternatives to Maximum Likelihood via Conjugate Gradient Linear Solvers

Initializing live version
Download to Desktop

Requires a Wolfram Notebook System

Interact on desktop, mobile and cloud with the free Wolfram Player or other Wolfram Language products.

The setting here is the same as described in [1], except that the autocorrelation function, which was the Matérn function, is now replaced by the powered exponential autocorrelation. For the definition of these correlations, see the documentation of the function VariogramModel.


Recall that the power can be any value in the half-open interval RowBox[{RowBox[{"(", RowBox[{"0", ",", "2"}]}], "]"}], but we restrict ourselves to as a rational fraction . Here RowBox[{"θ", ">", "0"}] is still denoted as the inverse-range parameter since SuperscriptBox["θ", RowBox[{"-", "1"}]] is the "range" (or "scale") in the autocorrelation .

As in [1], assume that one realization of the considered process is observed over the interval RowBox[{RowBox[{"(", RowBox[{"0", ",", "1"}]}], "]"}] at in the regular grid (thus is the size of each dataset). As for the Matérn (1) case, such a dataset is also simulated via fast Fourier transforms (FFTs).

The exponent is assumed known and the same three methods as in the Matérn (1) case are implemented here for estimating θ and the variance SuperscriptBox["τ", "2"], namely: Cressie's weighted least-squares method [2], the "hybrid" Zhang–Zimmerman method [3] and the recently-proposed CGEM-EV method [4], here in the case with no measurement error. For each of the three methods, the estimate of the so-called "microergodic coefficient," which is now (in the case , this coincides with the well-known expression of for the Matérn (1/2) case since the general expression for the Matérn (ν) case is ), is naturally defined as the product of the squared estimate of by the power of the estimate of θ. Recent mathematical studies [5, 6] have established that can indeed be easily estimated from one dense observation of over a finite interval (a so-called "infill" asymptotic framework). However, the separate estimation of θ or τ is more difficult.


Contributed by: Didier A. Girard (August 2022)
(CNRS-LJK and Université Grenoble Alpes)
Open content licensed under CC BY-NC-SA



See further details and references in [1].

The computation of , where , still uses a conjugate-gradient (CG) solver preconditioned by the classical factored sparse approximate inverse (FSAI) method, each product by being obtained via FFTs from the standard embedding of in a circulant matrix.


[1] D. A. Girard, "Estimating a Centered Matérn (1) Process: Three Alternatives to Maximum Likelihood via Conjugate Gradient Linear Solvers" from the Wolfram Demonstrations Project—A Wolfram Web Resource.

[2] N. Cressie, "Fitting Variogram Models by Weighted Least Squares," Journal of the International Association for Mathematical Geology, 17(5), 1985 pp. 563–586. doi:10.1007/BF01032109.

[3] H. Zhang and D. L. Zimmerman, "Hybrid Estimation of Semivariogram Parameters," Mathematical Geology, 39(2), 2007 pp. 247–260. doi:10.1007/s11004-006-9070-8.

[4] D. A. Girard, "Asymptotic Near-Efficiency of the 'Gibbs-Energy and Empirical-Variance' Estimating Functions for Fitting Matérn Models to a Dense (Noisy) Series."

[5] E. Anderes, "On the Consistent Separation of Scale and Variance for Gaussian Random Fields," The Annals of Statistics, 38(2), 2010 pp. 870–893. doi:10.1214/09-AOS725.

[6] W.-Li. Loh, "Estimating the Smoothness of a Gaussian Random Field from Irregularly Spaced Data via Higher-Order Quadratic Variations," The Annals of Statistics, 43(6), 2015 pp. 2766–2794. doi:10.1214/15-AOS1365.

Feedback (field required)
Email (field required) Name
Occupation Organization
Note: Your message & contact information may be shared with the author of any specific Demonstration for which you give feedback.