title: "Example5" output: html_document: keep_md: true


knitr::opts_chunk$set(echo = TRUE, warning=FALSE,message=FALSE,fig.path = "README_figs/README-")

Fitting a Matérn covariance to a (possibly incomplete) lattice observation

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.

Contents

Setting the probabilistic model

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 108x108 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.0833 (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 <- 108
# gaussian matern creation
nu <-0.5
effectiveRange <- 0.0833333
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.

Simulating one realization

This is simply:

set.seed(321)  # so that it is reproducible #
simulate(gm)

Plotting (and saving) several realizations

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 MacBookPro 13inch 2018 - Intel I5 2.3GHz - ram8GB :

ut   # for the simulation of 5 realizations :

Setting the uncomplete lattice

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(ex1WithN1eq108And5missindDisks.gd <- grid.domain(missing.domains=ex1.md,grid.size=n1grid,
            smoothness=nu)))

Size of the data set:

(nObs <- dim(ex1WithN1eq108And5missindDisks.gd$sparseG)[1])

Average width of the sparse preconditioner:

length(ex1WithN1eq108And5missindDisks.gd$sparseG) / dim(ex1WithN1eq108And5missindDisks.gd$sparseG)[1]

Plotting observation sites

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[!ex1WithN1eq108And5missindDisks.gd$missing.sites,]
plot(x,asp=1, xlim=c(0,1), ylim=c(0,1), pch=".")

No-noise case : Computing (and plotting) the estimating function at log-equispaced ranges

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][!ex1WithN1eq108And5missindDisks.gd$missing.sites]
indexReplcitate <- 1
z <- fullLattice.sixZs[,indexReplcitate][!ex1WithN1eq108And5missindDisks.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=ex1WithN1eq108And5missindDisks.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][!ex1WithN1eq108And5missindDisks.gd$missing.sites]
  bEV  <- mean(z**2)
  out <-     fsai11Precond.GEevalOnThetaGridNEW(z, candidateThetas1DGrid,
          gm$smoothness, ex1WithN1eq108And5missindDisks.gd ,tolPGC=1e-04)
#  
  plot(candidateThetas1DGrid, out$values, type="l",
                 col=1, lty= "dashed",log="xy")
    abline(h= bEV, lty= "dotted")
})

Timing for a MacBookPro 13inch 2018 - Intel I5 2.3GHz - ram8GB

ut   # for computing the estimating equation for 5 realizations :

Estimating theta and the micro-ergodic parameter from noisy observations with (unknown) SNR $= 33^2$, known sigmaNOISE

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=50 datasets Yk, k=1,...,50 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(720+indexReplcitate)
  simulate(gm)
  z <- gm$look[1:gm$n1,1:gm$n1][!ex1WithN1eq108And5missindDisks.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, ex1WithN1eq108And5missindDisks.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 50 replicates, using a MacBookPro 13inch 2018 - Intel I5 2.3GHz - ram8GB :

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.083)")
plot(density((cHatCGEMEV-cTrue)/cTrue),main="")
title(" N errors (cHatCGEMEV-cTrue)/cTrue")

Estimating theta and the micro-ergodic parameter from noisy observations with (unknown) SNR $= 33^2$, unknown sigmaNOISE: CGEM-EV-VEZ method

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). THe second stage is a simple "plugin" of this variance 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 ################
##########################################################
failureOfVEZ<-matrix(FALSE,nbReplicates)
sigma2Process<- 33^2
sigma2Noise <- 1.
bTrue <- sigma2Process/sigma2Noise

nbReplicates <-100
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(720+indexReplcitate)
  simulate(gm)
  z <- gm$look[1:gm$n1,1:gm$n1][!ex1WithN1eq108And5missindDisks.gd$missing.sites]
  y <-  sqrt(bTrue)* z +  c(rnorm(nObs))
# compute bHatEVsigma2NoiseVEZ
  Yrepl<-rep(NA,gm$n1*gm$n1)
      Yrepl[!ex1WithN1eq108And5missindDisks.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, ex1WithN1eq108And5missindDisks.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 MacBookPro 13inch 2018 - Intel I5 2.3GHz - ram8GB :

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.083)")
plot(density((cHatCGEMEVsigma2NoiseVEZ-cTrue)/cTrue,na.rm=TRUE),main="")
title(" N errors (cHatCGEMEVVEZ-cTrue)/cTrue")


didiergirard/CGEMEV documentation built on Aug. 14, 2019, 10:08 a.m.