Conditional Expectation¶

In [1]:
import gstlearn as gl
import gstlearn.plot as gp
import gstlearn.document as gdoc
import matplotlib.pyplot as plt
import numpy as np
import os

gdoc.setNoScroll()
In [2]:
ndim = 2
gl.defineDefaultSpace(gl.ESpaceType.RN, ndim);
In [3]:
## Load observations
temp_nf = gdoc.loadData("Scotland", "Scotland_Temperatures.NF")
dat = gl.Db.createFromNF(temp_nf)

## Load grid
elev_nf = gdoc.loadData("Scotland", "Scotland_Elevations.NF")
grid = gl.DbGrid.createFromNF(elev_nf)
In [4]:
unique_neigh = gl.NeighUnique.create()

Histogram of the raw variable (Temperature)

In [5]:
ax = gp.histogram(dat, name="January*", bins=20)
ax.decoration(title="Temperatures")

Gaussian scores

In [6]:
anam = gl.AnamHermite(30)
err = anam.fitFromLocator(dat)
err = anam.rawToGaussian(dat, "January_temp")
In [7]:
anam.display()
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
 

Plot the Gaussian scores

In [8]:
ax = gp.sortedcurve(tabx=dat["Y.January_temp"], taby=dat["January_temp"])
ax.decoration(xlabel="Gaussian",ylabel="Raw")

Draw the histogram of the Gaussian transformed values

In [9]:
ax = gp.histogram(dat, name="Y.January*", bins=20)
ax.decoration(title="Temperatures (Gaussian scale)")

We calculate the experimental directional variogram of the gaussian scores and fit the Model (with the constraints that sill should be 1)

In [10]:
varioparam = gl.VarioParam.createMultiple(ndir=2, npas=40, dpas=10)
vario_gauss2dir = gl.Vario.create(varioparam)
err = vario_gauss2dir.compute(dat)

fitmodgauss = gl.Model()
err = fitmodgauss.fit(vario_gauss2dir, 
                types=[gl.ECov.NUGGET, gl.ECov.SPHERICAL, gl.ECov.CUBIC],
                constraints = gl.Constraints(1))
In [11]:
ax = gp.varmod(vario_gauss2dir, fitmodgauss)
In [12]:
neighU = gl.NeighUnique.create()

Kriging of Gaussian scores

In [13]:
err = gl.kriging(dat, grid, fitmodgauss, neighU)
In [14]:
gp.setDefaultGeographic(dims=[8,8])
ax = grid.plot("*estim")
ax = dat.plot()
ax.decoration(title="Kriging of Gaussian scores")
In [15]:
ax = grid.plot("*stdev")
ax = dat.plot(flagCst=True)
ax.decoration(title="St. Dev. of Gaussian scores")

Conditional Expectation¶

We use the Monte Carlo method with 1000 outcomes.

In [16]:
selectivity = gl.Selectivity.createByKeys(["Z"], [], flag_est=True, flag_std=True)
err = gl.ConditionalExpectation(grid, anam, selectivity, "K*.estim", "K*.stdev", nbsimu=100,
                                namconv=gl.NamingConvention("CE",False,True,False))

Display of the Conditional Expectation

In [17]:
ax = grid.plot("CE*estim")
ax = dat.plot()
ax.decoration(title = "Conditional Expectation")

Display of the Conditional Standard Deviation

In [18]:
ax = grid.plot("CE*stdev")
ax = dat.plot(flagCst=True)
ax.decoration(title="Conditional Standard Deviation")

Conditional Probability below 0

In [19]:
selectivity = gl.Selectivity.createByKeys(["PROP"], zcuts=[0],flag_est=True, flag_std=True)
err = gl.ConditionalExpectation(grid, anam, selectivity, "K*.estim", "K*.stdev",
                                namconv=gl.NamingConvention("CE",False,True,False))
In [20]:
ax = grid.plot("CE.Proba*estim")
ax = dat.plot()
ax.decoration(title = "Conditional Probability below 0")

Conditional Probability above 1

In [21]:
selectivity = gl.Selectivity.createByKeys(["T"], zcuts=[1],flag_est=True, flag_std=True)
err = gl.ConditionalExpectation(grid, anam, selectivity, "K*.estim", "K*.stdev",
                                namconv=gl.NamingConvention("CE",False,True,False))
In [22]:
ax = grid.plot("CE.T*estim-1")
ax = dat.plot()
ax.decoration(title = "Conditional Probability above 1")
In [23]:
ax = grid.plot("CE.T*stdev-1")
ax = dat.plot(flagCst=True)
ax.decoration(title = "Conditional probability (Standard Deviation)")