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