In didiergirard/CGEMEV: Conditional Gibbs Energy Mean - Empirical Variance
knitr::opts_chunk$set(echo = TRUE)
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
- a simple function for simulating realizations of a stationary isotropic Gaussian process when the correlation belongs to the common Matérn family with smoothness index nu > 0
- tools for estimating, when nu is known and the level sigmaNOISE of the possible ``measurement errors'' is either known or unknown, the correlation range, also called ''decorrelation length'' (or its inverse that will be denoted theta) from one realization on a (possibly incomplete) lattice.
The variance of the field is simply estimated by the empirical variance that will be denoted bEV.
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
- Simulating one realization
- Plotting (and saving) several realizations
- Setting the uncomplete lattice
- Plotting data
- Computing (and plotting) the estimating function at log-equispaced ranges
- Estimating theta and the micro-ergodic parameter
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.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 <- 108
# 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.
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 iMac (Mid-2017) 4.2 GHz Core i7 (I7-7700K) :
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 with a iMac (Mid-2017) 4.2 GHz I7-7700K (using one core) :
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=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(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 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")
didiergirard/CGEMEV documentation built on Aug. 14, 2019, 10:08 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = TRUE)
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
- a simple function for simulating realizations of a stationary isotropic Gaussian process when the correlation belongs to the common Matérn family with smoothness index nu > 0
- tools for estimating, when nu is known and the level sigmaNOISE of the possible ``measurement errors'' is either known or unknown, the correlation range, also called ''decorrelation length'' (or its inverse that will be denoted theta) from one realization on a (possibly incomplete) lattice. The variance of the field is simply estimated by the empirical variance that will be denoted bEV.
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
- Simulating one realization
- Plotting (and saving) several realizations
- Setting the uncomplete lattice
- Plotting data
- Computing (and plotting) the estimating function at log-equispaced ranges
- Estimating theta and the micro-ergodic parameter
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.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 <- 108 # 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.
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 iMac (Mid-2017) 4.2 GHz Core i7 (I7-7700K) :
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 with a iMac (Mid-2017) 4.2 GHz I7-7700K (using one core) :
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=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(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 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")
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.