Conditional expectation

rm(list=ls())
library(gstlearn)

## Data points
fileNF = loadData("Scotland", "Scotland_Temperatures.NF")
dat = Db_createFromNF(fileNF)

## Target grid
fileNF = loadData("Scotland", "Scotland_Elevations.NF")
grid = DbGrid_createFromNF(fileNF)

---

p = plot.init()
p = p + plot.hist(dat,"January*")
p = p + plot.decoration(title="Temperatures")
plot.end(p)

## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

---

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

---

p = plot.init()
p = p + plot.XY(dat["Y.January_temp"], dat["January_temp"], flagLine=FALSE, flagPoint=TRUE)
p = p + plot.decoration(xlab="Gaussian", ylab="Raw")
plot.end(p)

---

p = plot.init()
p = p + plot.hist(dat,"Y.January*")
p = p + plot.decoration(title="Temperatures (Gaussian scale)")
plot.end(p)

## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

---

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, nlag=40, dlag=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))

plot.init() + plot.varmod(vario_gauss2dir, fitmodgauss)

---

neighU = NeighUnique_create()

err = kriging(dat, grid, fitmodgauss, neighU)

p = plot.init()
p = p + plot(grid,"*estim")
p = p + plot(dat)
p = p + plot.decoration(title="Kriging of Gaussian scores")
plot.end(p)

---

p = plot.init()
p = p + plot(grid,"*stdev")
p = p + plot(dat, flagCst=TRUE)
p = p + plot.decoration(title="St. Dev. of Gaussian scores")
plot.end(p)

---

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))

---

p = plot.init()
p = p + plot(grid,"CE*estim")
p = p + plot(dat)
p = p + plot.decoration(title = "Conditional Expectation")
plot.end(p)

---

p = plot.init()
p = p + plot(grid, "CE*stdev")
p = p + plot(dat, flagCst=TRUE)
p = p + plot.decoration(title="Conditional Standard Deviation")
plot.end(p)

---

selectivity = Selectivity_createByKeys(c("PROP"), zcuts=c(0),
                                       flag_est=TRUE, flag_std=TRUE)
err = ConditionalExpectation(grid, anam, selectivity, 
                             "K*.estim", "K*.stdev",
                             namconv=NamingConvention("CE",FALSE,TRUE,FALSE))

p = plot.init()
p = p + plot(grid,"CE.Proba*estim")
p = p + plot(dat)
p = p + plot.decoration(title = "Conditional Probability below 0")
plot.end(p)

---

selectivity = Selectivity_createByKeys(c("T"), zcuts=c(1),
                                       flag_est=TRUE, flag_std=TRUE)
err = ConditionalExpectation(grid, anam, selectivity, "K*.estim", "K*.stdev",
                             namconv=NamingConvention("CE",FALSE,TRUE,FALSE))

---

p = plot.init()
p = p + plot(grid,"CE.T*estim-1")
p = p + plot(dat)
p = p + plot.decoration(title = "Conditional Probability above 1")
plot.end(p)

---

p = plot.init()
p = p + plot(grid, "CE.T*stdev-1")
p = p + plot(dat, flagCst=TRUE)
p = p + plot.decoration(title = "Conditional probability (Standard Deviation)")
plot.end(p)