Conditional expectation
gstlearn Team
February 2023
rm(list=ls())
library(gstlearn)
library(ggplot2)
library(ggpubr)
library(ggnewscale)## Data points
fileNF = loadData("Scotland", "Scotland_Temperatures.NF")
dat = Db_createFromNF(fileNF)
## Target grid
fileNF = loadData("Scotland", "Scotland_Elevations.NF")
grid = DbGrid_createFromNF(fileNF)\frametitle{Histogram of the raw variable (temperature)}
p = ggplot()
p = p + plot.hist(dat,"January*")
p = p + plot.decoration(title="Temperatures")
ggPrint(p)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
\frametitle{Gaussian scores}
anam = AnamHermite(30)
err = anam$fitFromLocator(dat)
err = anam$rawToGaussian(dat, "January_temp")
anam##
## Hermitian Anamorphosis
## ----------------------
## Minimum absolute value for Y = -2.8
## Maximum absolute value for Y = 2.7
## Minimum absolute value for Z = 0.62599
## Maximum absolute value for Z = 5.24756
## Minimum practical value for Y = -2.8
## Maximum practical value for Y = 2.7
## Minimum practical value for Z = 0.62599
## Maximum practical value for Z = 5.24756
## Mean = 2.81457
## Variance = 1.01677
## Number of Hermite polynomials = 30
## Normalized coefficients for Hermite polynomials (punctual variable)
## [, 0] [, 1] [, 2] [, 3] [, 4] [, 5] [, 6]
## [ 0,] 2.815 -1.003 0.010 0.067 0.005 0.030 -0.007
## [ 7,] -0.035 0.009 0.027 -0.011 -0.019 0.014 0.013
## [ 14,] -0.017 -0.008 0.019 0.004 -0.020 -0.001 0.020
## [ 21,] -0.002 -0.018 0.004 0.016 -0.005 -0.014 0.006
## [ 28,] 0.011 -0.005\frametitle{Gaussian scores}
p = ggplot()
p = p + plot.XY(dat["Y.January_temp"], dat["January_temp"], flagLine=FALSE, flagPoint=TRUE)
p = p + plot.decoration(xlab="Gaussian", ylab="Raw")
ggPrint(p)
\frametitle{Gaussian scores}
p = ggplot()
p = p + plot.hist(dat,"Y.January*")
p = p + plot.decoration(title="Temperatures (Gaussian scale)")
ggPrint(p)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
\frametitle{Kriging of the Gaussian scores}
We calculate the experimental directional variogram of the gaussian scores and fit the Model (with the constraints that sill should be 1)
varioparam = VarioParam_createMultiple(ndir=2, npas=40, dpas=10)
vario_gauss2dir = Vario_create(varioparam)
err = vario_gauss2dir$compute(dat)
fitmodgauss = Model()
err = fitmodgauss$fit(vario_gauss2dir,
types=ECov_fromKeys(c("NUGGET", "SPHERICAL","CUBIC")),
constraints = Constraints(1))ggplot() + plot.varmod(vario_gauss2dir, fitmodgauss)
\frametitle{Kriging of the Gaussian scores}
neighU = NeighUnique_create()
err = kriging(dat, grid, fitmodgauss, neighU)p = ggDefaultGeographic()
p = p + plot(grid,"*estim")
p = p + plot(dat)
p = p + plot.decoration(title="Kriging of Gaussian scores")
ggPrint(p)
\frametitle{Kriging of the Gaussian scores}\framesubtitle{Standard Deviation}
p = ggDefaultGeographic()
p = p + plot(grid,"*stdev")
p = p + plot(dat, flagCst=TRUE)
p = p + plot.decoration(title="St. Dev. of Gaussian scores")
ggPrint(p)
\frametitle{Conditional expectation}
Use the Turning Bands method with 1000 simulations
selectivity = Selectivity_createByKeys(c("Z"), flag_est=TRUE, flag_std=TRUE)
err = ConditionalExpectation(grid, anam, selectivity, "K*.estim", "K*.stdev",
nbsimu=100,
namconv=NamingConvention("CE",FALSE,TRUE,FALSE))\frametitle{Conditional expectation}
p = ggDefaultGeographic()
p = p + plot(grid,"CE*estim")
p = p + plot(dat)
p = p + plot.decoration(title = "Conditional Expectation")
ggPrint(p)