title: "README" output: html_document: keep_md: 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)
## Loading required package: Rcpp
##
## Attaching package: 'CGEMEV'
## The following object is masked from 'package:graphics':
##
## grid
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)
## Loading required package: spam
## Loading required package: dotCall64
## Loading required package: grid
## Spam version 2.2-0 (2018-06-19) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
##
## Attaching package: 'spam'
## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve
## Loading required package: maps
## See www.image.ucar.edu/~nychka/Fields for
## a vignette and other supplements.
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 window will lay out plots in a 2 by 3 matrix
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 :
## user system elapsed
## 0.982 0.019 1.015
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)))
## user system elapsed
## 4.379 0.304 4.504
Size of the data set:
(nObs <- dim(ex1WithN1eq128And5missindDisks.gd$sparseG)[1])
## [1] 14431
Average width of the sparse preconditioner:
length(ex1WithN1eq128And5missindDisks.gd$sparseG) / dim(ex1WithN1eq128And5missindDisks.gd$sparseG)[1]
## [1] 72.06645
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)
## [1] 5.999988
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))
## [1] 1.535706
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)
)
## $values
## [,1]
## [1,] 6.1235805
## [2,] 4.7069899
## [3,] 3.6182450
## [4,] 2.7815080
## [5,] 2.1385032
## [6,] 1.6444457
## [7,] 1.2649374
## [8,] 0.9735590
## [9,] 0.7500493
## [10,] 0.5789036
## [11,] 0.4483112
## [12,] 0.3493698
## [13,] 0.2755327
## [14,] 0.2222467
## [15,] 0.1868012
##
## $niterForY
## [1] 16 14 14 13 13 14 11 9 7 6 6 8 11 17 20
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 window will lay out plots in a 2 by 3 matrix
{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 :
## user system elapsed
## 10.950 0.775 11.753
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
## user system elapsed
## 361.480 25.167 387.098
Summary of the results:
print(c("log10 of true range =",log(sqrt(2*nu)*gm$range,10)))
## [1] "log10 of true range =" "-0.778150381795548"
summary(log(sqrt(2*nu)/thetaHatCGEMEV,10))
## V1
## Min. :-1.0545
## 1st Qu.:-0.8525
## Median :-0.7886
## Mean :-0.7817
## 3rd Qu.:-0.7186
## Max. :-0.4458
cTrue <- bTrue*(1/gm$range)**(2*gm$smoothness)
summary(cHatCGEMEV/cTrue)
## V1
## Min. :0.9658
## 1st Qu.:0.9932
## Median :1.0007
## Mean :0.9996
## 3rd Qu.:1.0068
## Max. :1.0273
sd(cHatCGEMEV/cTrue)
## [1] 0.01195334
This sd can be compared to the CR lower bound for the case known theta and no-noise
sqrt(2/length(z))
## [1] 0.01177245
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
## 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 :
## Loading required package: spatstat
## Loading required package: spatstat.data
## Loading required package: nlme
## Loading required package: rpart
##
## spatstat 1.56-1 (nickname: 'Invisible Friend')
## For an introduction to spatstat, type 'beginner'
##
## Note: R version 3.5.1 (2018-07-02) is more than 9 months old; we strongly recommend upgrading to the latest version
## Loading required package: smoothie
## Loading required package: smatr
## Loading required package: turboEM
## Loading required package: doParallel
## Loading required package: foreach
## Loading required package: iterators
## Loading required package: parallel
## Loading required package: numDeriv
## Loading required package: quantreg
## Loading required package: SparseM
##
## Attaching package: 'SparseM'
## The following object is masked from 'package:spam':
##
## backsolve
## The following object is masked from 'package:base':
##
## backsolve
##
## Attaching package: 'turboEM'
## The following objects are masked from 'package:numDeriv':
##
## grad, hessian
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))
## [1] 78
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
## user system elapsed
## 295.924 20.114 316.445
Summary of the results:
print(c("log10 of true range =",log(sqrt(2*nu)*gm$range,10)))
## [1] "log10 of true range =" "-0.778150381795548"
summary(log(sqrt(2*nu)/thetaHatCGEMEVsigma2NoiseVEZ,10))
## V1
## Min. :-1.0652
## 1st Qu.:-0.8435
## Median :-0.7859
## Mean :-0.7769
## 3rd Qu.:-0.7150
## Max. :-0.4500
## NA's :22
cTrue <- bTrue*(1/gm$range)**(2*gm$smoothness)
summary(cHatCGEMEVsigma2NoiseVEZ/cTrue)
## V1
## Min. :0.9692
## 1st Qu.:0.9973
## Median :1.0124
## Mean :1.0169
## 3rd Qu.:1.0230
## Max. :1.1698
## NA's :22
sd(cHatCGEMEVsigma2NoiseVEZ/cTrue,na.rm=TRUE)
## [1] 0.03453066
This sd can be compared to the CR lower bound for the case known theta and no-noise
sqrt(2/length(z))
## [1] 0.01177245
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:
## plot window will lay out plots in a 2 by 2 matrix
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.