knitr::opts_chunk$set(echo = TRUE)
The CGEMEV
R package (dependent on the fields, spam, SpatialVx and RCPP packages) provides
These tools implement the estimation method proposed in Girard, D.A., 2016, Asymptotic near-efficiency of the Gibbs-energy and empirical-variance estimating functions for fitting Matérn models I: Densely sampled processes. Statistics and Probability Letters 110, 191-197, and in the preprint https://hal.archives-ouvertes.fr/hal-00515832 .
In this first version of CGEMEV
, the region of missing observations is defined as a union of disks.
Three main functions are used: gaussian.matern()
, simulate()
and fsai11Precond.GEevalOnThetaGrid()
, fsaiThreePrecond.fastRandzedCGEMEVbisectionLogScaleSearch()
and a fifth function grid.domain()
is required to precompute preconditioning sparse matrices.
These R-fonctions can be applied to a quite large grid even on a laptop (for example 512x512, provided the ''extension factor'' required for simulation, see below, is not too big). Indeed quite fast computation of the quadratic form which occurs in the estimating equation is possible by using a conjugate-gradient (CG) solver preconditioned by a classical factored sparse approximate inverse (FSAI) preconditioning, since the matrix-vector product, required in each CG iteration, can be obtained via FFT from the standard embedding of the correlation matrix in a circulant matrix.
In this first version of the CGEMEV
package, for simplicity, we restrict the spatial domain to be the unit square (0,1)X(0,1). For the example here, the simulations and the choice of observed sites are done using a 128x128 grid which partitions this domain. We consider the example nu =0.5
(that this the widely used exponential correlation).
Let us choose the correlation range such that its inverse, denoted $\theta$, satisfies $\sqrt{2 nu} \theta^{-1} =0.166667$.
The resulting correlation function can then be considered as one with an ‘’effective range’‘ equal to $0.166667$ (the formulae for Matern correlations often use this quantity denoted $\rho$, see https://en.wikipedia.org/wiki/Mat%C3%A9rn_covariance_function ).
library(CGEMEV) library(pdist) n1grid <- 128 # gaussian matern creation nu <-0.5 effectiveRange <- 0.166667 gm <- gaussian.matern(grid.size=n1grid,smoothness=nu, range=effectiveRange/sqrt(2*nu),factor=2)
NB: in the previous setting, factor=2
specifies the required extension factor (assumed, for simplicity, to be an integer) of the observation domain. Indeed for this example the choice factor=1
would entail (when calling simulate(gm)
) the message “FFT of covariance has negative values” which means that generating a realization via the classical embedding method (which doubles each length of the considered rectangular domain) would not work.
This is simply:
set.seed(321) # so that it is reproducible # simulate(gm)
We can plot (and save), the previous realization and, for example, 5 further realizations:
fullLattice.sixZs<- array(NA,c(n1grid*n1grid,6)) set.panel(2,3) plot(gm) fullLattice.sixZs[,1]<-gm$look[1:gm$n1,1:gm$n1] ut <- system.time( for (indexReplcitate in 2:6){ set.seed(320+indexReplcitate) simulate(gm) plot(gm) fullLattice.sixZs[,indexReplcitate]<-gm$look[1:gm$n1,1:gm$n1] })
The following timing (in seconds) is for a iMac (Mid-2017) 4.2 GHz Core i7 (I7-7700K) :
ut # for the simulation of 5 realizations :
Let us now define the regions (actually 5 disks) where the observations will be missing, and precompute the preconditioning matrix:
# md=missing.domains listOfHalfDiameters <- c(1,0.4,0.2,0.1,0.05)*sqrt(.1/ pi) listOfCenters <- t(matrix(c(c(0.5+sqrt(.1/ pi),0.5+sqrt(.1/ pi)),c(0.1+sqrt(.1/ pi), 0.05+sqrt(.1/ pi)),c(0.8,0.25),c(0.2,0.8),c(0.3,0.7)),2,length(listOfHalfDiameters))) ex1.md <- list( list(center=listOfCenters[1,],radius= listOfHalfDiameters[1]), list(center=listOfCenters[2,],radius= listOfHalfDiameters[2]), list(center=listOfCenters[3,],radius= listOfHalfDiameters[3]), list(center=listOfCenters[4,],radius= listOfHalfDiameters[4]), list(center=listOfCenters[5,],radius= listOfHalfDiameters[5]) ) # gd=grid.domain print(system.time(ex1WithN1eq128And5missindDisks.gd <- grid.domain(missing.domains=ex1.md,grid.size=n1grid, smoothness=nu)))
Size of the data set:
(nObs <- dim(ex1WithN1eq128And5missindDisks.gd$sparseG)[1])
Average width of the sparse preconditioner:
length(ex1WithN1eq128And5missindDisks.gd$sparseG) / dim(ex1WithN1eq128And5missindDisks.gd$sparseG)[1]
Each dataset is made up of the observations of one realization of the previous type at the following points:
xFull<- (1./n1grid)*matrix( c(rep(1: n1grid, n1grid),rep(1: n1grid,each= n1grid)), ncol=2) x <- xFull[!ex1WithN1eq128And5missindDisks.gd$missing.sites,] plot(x,asp=1, xlim=c(0,1), ylim=c(0,1), pch=".")
Choose a grid of candidates for the range parameter (more precisely, for the inverse-range parameter denoted theta
) at which the estimating function is computed:
(thetaTrue <- 1/gm$range)
candidateThetas1DGrid <- thetaTrue * 10**seq(-0.8,0.8,,15)
Consider the first one of the above realizations, and the naive variance estimator bEV:
# only observed outside the disks: #z <- gm$look[1:gm$n1,1:gm$n1][!ex1WithN1eq128And5missindDisks.gd$missing.sites] indexReplcitate <- 1 z <- fullLattice.sixZs[,indexReplcitate][!ex1WithN1eq128And5missindDisks.gd$missing.sites] # (bEV <- mean(z**2))
For this dataset z
, let us give the whole output of the function fsai11Precond.GEevalOnThetaGrid()
:
(out <- fsai11Precond.GEevalOnThetaGridNEW(z,candidateThetas1DGrid,nu=gm$smoothness, grid.domain=ex1WithN1eq128And5missindDisks.gd,tolPGC=1e-04) )
Let us repeat this computation for the five next data sets obtained from the above realizations, and plot the results:
set.panel(2,3) {plot(candidateThetas1DGrid, out$values, type="l", col=1,lty= "dashed",log="xy") abline(h= bEV, lty= "dotted" ) } ut <- system.time( for (indexReplcitate in 2:6){ z <- fullLattice.sixZs[,indexReplcitate][!ex1WithN1eq128And5missindDisks.gd$missing.sites] bEV <- mean(z**2) out <- fsai11Precond.GEevalOnThetaGridNEW(z, candidateThetas1DGrid, gm$smoothness, ex1WithN1eq128And5missindDisks.gd ,tolPGC=1e-04) # plot(candidateThetas1DGrid, out$values, type="l", col=1, lty= "dashed",log="xy") abline(h= bEV, lty= "dotted") })
Timing with a iMac (Mid-2017) 4.2 GHz I7-7700K (using one core) :
ut # for computing the estimating equation for 5 realizations :
For solving the CGEMEV estimating equation in theta, a specific bisection method has been elaborated, since it is useful to use "historical" results (more precisely the starting point of the PCG iterations at a given theta is chosen as the solution from the linear-solve associated with the previously considered theta); indeed a classic use of the uniroot
R-fonction would be much slower. Let us analyse N=100 datasets Yk, k=1,...,100 each one being Yk = Zk +epsilonk , where epsilonk is a field of white normal noise; they are simulated such that the true SNR (i.e. var(Z)/var(noise)) $= 33^2$. Only nu and the variance of the noise (or nugget-effect variance, $=1$) are assumed known.
sigma2Process<- 33^2 sigma2Noise <- 1. bTrue <- sigma2Process/sigma2Noise nbReplicates <- 100 bHatEV<-matrix(NA,nbReplicates) cHatCGEMEV<-matrix(NA, nbReplicates) thetaHatCGEMEV<-matrix(NA, nbReplicates) nCGiterationsMaxForY<-matrix(NA, nbReplicates) # ut <- system.time( for (indexReplcitate in 1: nbReplicates){ set.seed(520+indexReplcitate) simulate(gm) z <- gm$look[1:gm$n1,1:gm$n1][!ex1WithN1eq128And5missindDisks.gd$missing.sites] y <- sqrt(bTrue)* z + c(rnorm(nObs)) # bEV <- ( mean(y**2)-sigma2Noise ) / sigma2Noise if (bEV <0) { break} w <- c(rnorm(nObs)) out <- fsaiThreePrecond.fastRandzedCGEMEVbisectionLogScaleSearch( sigma2Noise,y, w, gm$smoothness, ex1WithN1eq128And5missindDisks.gd ,tolPGC=1e-04, 0.2,50, tolBis=1e-05) # thetaHatCGEMEV[indexReplcitate]<- (out$root) cHatCGEMEV[indexReplcitate]<- bEV*(out$root)**(2*gm$smoothness) nCGiterationsMaxForY[indexReplcitate]<- max(out$niterCGiterationsHistory) })
Timing for these 100 replicates, with a iMac (Mid-2017) 4.2 GHz I7-7700K (using one core) :
ut
Summary of the results:
print(c("log10 of true range =",log(sqrt(2*nu)*gm$range,10)))
summary(log(sqrt(2*nu)/thetaHatCGEMEV,10))
cTrue <- bTrue*(1/gm$range)**(2*gm$smoothness) summary(cHatCGEMEV/cTrue)
sd(cHatCGEMEV/cTrue)
This sd can be compared to the CR lower bound for the case known theta and no-noise
sqrt(2/length(z))
Let us plot an estimate of the density of the CGEMEV estimates of the effective range, and of the density of the relative errors in the CGEMEV estimates of the microergodic-parameter cTrue:
set.panel(2,2)
## plot window will lay out plots in a 2 by 2 matrix
den <- density(log(sqrt(2*nu)/thetaHatCGEMEV,10)) den$x <- 10**(den$x) plot(den, log="x",main="") title("CGEMEVestimates of the range (=0.166)") plot(density((cHatCGEMEV-cTrue)/cTrue),main="") title(" N errors (cHatCGEMEV-cTrue)/cTrue")
In a first stage, the variance of the noise (or nugget-effect), whose true value $=1$, is estimated by extrapolating at zero the classical semi-variogram (using maximum lag set to 3); this estimate is stored in sigma2NoiseVEZ
(this method was proposed in Matthias Katzfuss and
Noel Cressie. Bayesian hierarchical spatio‐temporal smoothing for very large datasets. Environmetrics. February 2012. pp 94-107). The second stage is a simple "plugin" of sigma2NoiseVEZ
in CGEM-EV in place of the true variance :
################################################################# ## definition of sig2noiseByVariogExtrapolAtZeroNew ############# ################################################################# library(SpatialVx) #help(structurogram.matrix) sig2noiseByVariogExtrapolAtZeroNew <- function(y,n1grid,n2grid, maxLag=5) {# missingPositions have zero value tmp <- array(y,c(n1grid,n2grid)) # outWithMissing<-structurogram.matrix(tmp,zero.out=TRUE,R=maxLag) variog<-outWithMissing$vgram lag<-outWithMissing$d lag2<-lag^2 quadratic.model <-lm(variog ~ lag + lag2) # sig2noiseVEZ <-quadratic.model$coefficients[1] # sig2noiseVEZ } ## end of sig2noiseByVariogExtrapolAtZero ################ ##########################################################
nbReplicates <-100 failureOfVEZ<-matrix(FALSE,nbReplicates)
sigma2Process<- 33^2 sigma2Noise <- 1. bTrue <- sigma2Process/sigma2Noise bHatEVsigma2NoiseVEZ<-matrix(NA,nbReplicates) cHatCGEMEVsigma2NoiseVEZ<-matrix(NA, nbReplicates) thetaHatCGEMEVsigma2NoiseVEZ<-matrix(NA, nbReplicates) nCGiterationsMaxForYsigma2NoiseVEZ<-matrix(NA, nbReplicates) sigma2NoiseVEZ<-matrix(NA, nbReplicates) # maxLag <-3 ut <- system.time( for (indexReplcitate in 1: nbReplicates){ set.seed(520+indexReplcitate) simulate(gm) z <- gm$look[1:gm$n1,1:gm$n1][!ex1WithN1eq128And5missindDisks.gd$missing.sites] y <- sqrt(bTrue)* z + c(rnorm(nObs)) # compute bHatEVsigma2NoiseVEZ Yrepl<-rep(NA,gm$n1*gm$n1) Yrepl[!ex1WithN1eq128And5missindDisks.gd$missing.sites] <- y minOfimageMinus1 <- min(Yrepl,na.rm=TRUE)-1 YreplOffset<-Yrepl-rep(minOfimageMinus1,gm$n1*gm$n1) YreplOffsetWithZeroInsteadNA <- YreplOffset YreplOffsetWithZeroInsteadNA[is.na(Yrepl)] <-0 estimation <- sig2noiseByVariogExtrapolAtZeroNew(YreplOffsetWithZeroInsteadNA,gm$n1,gm$n1,maxLag=maxLag) # lowerBound <- 10**(-12)*mean(y**2) sigma2NoiseVEZ[indexReplcitate] <- max(c(lowerBound,estimation),na.rm=TRUE) if (sigma2NoiseVEZ[indexReplcitate]==lowerBound) failureOfVEZ[indexReplcitate] <- TRUE # else {bEV <- ( mean(y**2)-sigma2NoiseVEZ[indexReplcitate] )/ sigma2NoiseVEZ[indexReplcitate] if (bEV <0) { break} w <- c(rnorm(nObs)) out <- fsaiThreePrecond.fastRandzedCGEMEVbisectionLogScaleSearch( sigma2NoiseVEZ[indexReplcitate] ,y, w, gm$smoothness, ex1WithN1eq128And5missindDisks.gd ,tolPGC=1e-04, 0.2,50, tolBis=1e-05) # # thetaHatCGEMEVsigma2NoiseVEZ[indexReplcitate] <- (out$root) cHatCGEMEVsigma2NoiseVEZ[indexReplcitate]<- ( mean(y**2)-sigma2NoiseVEZ[indexReplcitate] ) * (out$root)**(2*gm$smoothness) # nCGiterationsMaxForYsigma2NoiseVEZ[indexReplcitate]<- max(out$niterCGiterationsHistory) } })
Number of non-failure of VEZ
(numberOfNonFailureOfVEZ <- sum(failureOfVEZ==FALSE))
Timing for these replicates (including the extrapolations of each empirical varogriam), using a iMac (Mid-2017) 4.2 GHz I7-7700K (using one core)
ut
Summary of the results:
print(c("log10 of true range =",log(sqrt(2*nu)*gm$range,10)))
summary(log(sqrt(2*nu)/thetaHatCGEMEVsigma2NoiseVEZ,10))
cTrue <- bTrue*(1/gm$range)**(2*gm$smoothness) summary(cHatCGEMEVsigma2NoiseVEZ/cTrue)
sd(cHatCGEMEVsigma2NoiseVEZ/cTrue,na.rm=TRUE)
This sd can be compared to the CR lower bound for the case known theta and no-noise
sqrt(2/length(z))
Let us plot an estimate of the density of the CGEM-EV-VEZ estimates of the effective range, and of the density of the relative errors in the CGEM-EV-VEZ estimates of the microergodic-parameter cTrue:
set.panel(2,2)
den <- density(log(sqrt(2*nu)/thetaHatCGEMEVsigma2NoiseVEZ,10),na.rm=TRUE) den$x <- 10**(den$x) plot(den, log="x",main="") title("CGEMEVVEZestimates of the range (=0.166)") plot(density((cHatCGEMEVsigma2NoiseVEZ-cTrue)/cTrue,na.rm=TRUE),main="") title(" N errors (cHatCGEMEVVEZ-cTrue)/cTrue")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.