In [6]:
if graphics:
ax = mydb.plot(nameColor="Pb")
ax.decoration(title="Data Set")
import numpy as np
import pandas as pd
import sys
import os
import matplotlib.pyplot as plt
import gstlearn as gl
import gstlearn.plot as gp
import gstlearn.document as gdoc
Global variables
verbose = True
graphics = True
gl.OptCst.define(gl.ECst.NTCOL,6)
gdoc.setNoScroll()
The data are stored in a CSV format in the file called Pollution.dat
filepath = gdoc.loadData("Pollution", "Pollution.dat")
mydb = gl.Db.createFromCSV(filepath,gl.CSVformat())
mydb.setLocators(["X","Y"],gl.ELoc.X)
mydb.setLocator("Zn",gl.ELoc.Z)
if verbose:
dbfmt = gl.DbStringFormat()
dbfmt.setFlags(flag_resume=True, flag_vars=True, flag_extend=True)
mydb.display(dbfmt)
Data Base Characteristics ========================= Data Base Summary ----------------- File is organized as a set of isolated points Space dimension = 2 Number of Columns = 5 Total number of samples = 102 Data Base Extension ------------------- Coor #1 - Min = 109.850 - Max = 143.010 - Ext = 33.16 Coor #2 - Min = 483.660 - Max = 513.040 - Ext = 29.38 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = X - Locator = x1 Column = 2 - Name = Y - Locator = x2 Column = 3 - Name = Zn - Locator = z1 Column = 4 - Name = Pb - Locator = p1
Accessing to the variable names
print("List of all variable names =",mydb.getAllNames())
List of all variable names = ('rank', 'X', 'Y', 'Zn', 'Pb')
Extracting the vector containing the Zn variable in order to perform a selection
tabZn = mydb.getColumn('Zn') # equivalent to mydb["Zn"] or mydb[3]
selZn = tabZn < 20
mydb.addSelection(selZn,'sel')
mydb.setLocator('Pb',gl.ELoc.Z)
dbfmt = gl.DbStringFormat.createFromFlags(flag_stats=True)
if verbose:
mydb.display(dbfmt)
Data Base Characteristics ========================= Data Base Summary ----------------- File is organized as a set of isolated points Space dimension = 2 Number of Columns = 6 Total number of samples = 102 Number of active samples = 99 Data Base Statistics -------------------- 1 - Name rank - Locator NA Nb of data = 102 Nb of active values = 99 Minimum value = 1.000 Maximum value = 102.000 Mean value = 51.808 Standard Deviation = 29.114 Variance = 847.650 2 - Name X - Locator x1 Nb of data = 102 Nb of active values = 99 Minimum value = 109.850 Maximum value = 143.010 Mean value = 119.924 Standard Deviation = 6.549 Variance = 42.890 3 - Name Y - Locator x2 Nb of data = 102 Nb of active values = 99 Minimum value = 483.660 Maximum value = 513.040 Mean value = 498.622 Standard Deviation = 8.000 Variance = 64.001 4 - Name Zn - Locator NA Nb of data = 102 Nb of active values = 99 Minimum value = 1.090 Maximum value = 12.100 Mean value = 2.881 Standard Deviation = 1.667 Variance = 2.779 5 - Name Pb - Locator z1 Nb of data = 102 Nb of active values = 99 Minimum value = 3.000 Maximum value = 12.700 Mean value = 5.658 Standard Deviation = 1.697 Variance = 2.881 6 - Name sel - Locator sel Nb of data = 102 Nb of active values = 99 Minimum value = 1.000 Maximum value = 1.000 Mean value = 1.000 Standard Deviation = 0.000 Variance = 0.000 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = X - Locator = x1 Column = 2 - Name = Y - Locator = x2 Column = 3 - Name = Zn - Locator = NA Column = 4 - Name = Pb - Locator = z1 Column = 5 - Name = sel - Locator = sel
Display my Data (with samples represented by color and size)
if graphics:
ax = mydb.plot(nameColor="Pb")
ax.decoration(title="Data Set")
We first define the geometry of the variogram calculations
myVarioParamOmni = gl.VarioParam()
mydir = gl.DirParam.create(npas=10, dpas=1.)
myVarioParamOmni.addDir(mydir)
We use the variogram definition in order to calculate the variogram cloud.
dbcloud = gl.db_variogram_cloud(mydb, myVarioParamOmni)
We recall that the Variogram cloud is calculated by filling an underlying grid where each cell is painted according to the number of pairs at the given distance and given variability. Representing the variogram cloud
if graphics:
ax = dbcloud.plot("Cloud*")
ax.decoration(title="Variogram Cloud")
Calculating the experimental omni-directional variogram
myVarioOmni = gl.Vario(myVarioParamOmni)
err = myVarioOmni.compute(mydb, gl.ECalcVario.VARIOGRAM)
if verbose:
myVarioOmni.display()
Variogram characteristics
=========================
Number of variable(s) = 1
Number of direction(s) = 1
Space dimension = 2
Variance-Covariance Matrix 2.881
Direction #1
------------
Number of lags = 10
Direction coefficients = 1.000 0.000
Direction angles (degrees) = 0.000 0.000
Tolerance on direction = 90.000 (degrees)
Calculation lag = 1.000
Tolerance on distance = 50.000 (Percent of the lag value)
For variable 1
Rank Npairs Distance Value
0 3.000 0.389 0.462
1 123.000 1.081 1.495
2 183.000 2.038 1.620
3 205.000 3.006 2.526
4 231.000 4.013 2.240
5 229.000 5.036 2.524
6 198.000 5.962 2.396
7 187.000 7.000 2.708
8 204.000 7.996 2.772
9 184.000 8.990 2.868
The variogram is represented graphically for a quick check
if graphics:
ax = myVarioOmni.plot()
ax.decoration(title="Omni-directional Variogram for Pb")
Calculate a variogram in several directions
myvarioParam = gl.VarioParam()
mydirs = gl.DirParam.createMultiple(4, 10, 1.)
myvarioParam.addMultiDirs(mydirs)
myvario = gl.Vario(myvarioParam)
myvario.compute(mydb,gl.ECalcVario.VARIOGRAM)
if verbose:
myvario.display()
Variogram characteristics
=========================
Number of variable(s) = 1
Number of direction(s) = 4
Space dimension = 2
Variance-Covariance Matrix 2.881
Direction #1
------------
Number of lags = 10
Direction coefficients = 1.000 0.000
Direction angles (degrees) = 0.000 0.000
Tolerance on direction = 22.500 (degrees)
Calculation lag = 1.000
Tolerance on distance = 50.000 (Percent of the lag value)
For variable 1
Rank Npairs Distance Value
0 1.000 0.410 0.180
1 29.000 1.094 1.634
2 47.000 2.079 1.415
3 53.000 3.003 2.824
4 63.000 3.999 2.348
5 66.000 5.035 2.319
6 60.000 5.978 3.115
7 52.000 7.045 2.746
8 52.000 8.020 3.927
9 37.000 8.980 2.554
Direction #2
------------
Number of lags = 10
Direction coefficients = 0.707 0.707
Direction angles (degrees) = 45.000 0.000
Tolerance on direction = 22.500 (degrees)
Calculation lag = 1.000
Tolerance on distance = 50.000 (Percent of the lag value)
For variable 1
Rank Npairs Distance Value
0 1.000 0.344 0.080
1 31.000 1.051 1.113
2 50.000 1.960 1.890
3 62.000 2.999 2.443
4 58.000 4.014 2.701
5 51.000 5.016 2.702
6 36.000 5.999 1.833
7 37.000 7.015 2.130
8 50.000 7.997 2.060
9 53.000 8.995 2.381
Direction #3
------------
Number of lags = 10
Direction coefficients = 0.000 1.000
Direction angles (degrees) = 90.000 0.000
Tolerance on direction = 22.500 (degrees)
Calculation lag = 1.000
Tolerance on distance = 50.000 (Percent of the lag value)
For variable 1
Rank Npairs Distance Value
1 32.000 1.149 1.631
2 39.000 2.080 1.670
3 39.000 2.979 2.511
4 48.000 4.012 2.120
5 51.000 5.029 3.055
6 47.000 5.939 2.856
7 49.000 6.965 2.386
8 42.000 7.952 2.708
9 41.000 9.018 2.320
Direction #4
------------
Number of lags = 10
Direction coefficients = -0.707 0.707
Direction angles (degrees) = 135.000 0.000
Tolerance on direction = 22.500 (degrees)
Calculation lag = 1.000
Tolerance on distance = 50.000 (Percent of the lag value)
For variable 1
Rank Npairs Distance Value
0 1.000 0.411 1.125
1 31.000 1.028 1.606
2 47.000 2.044 1.496
3 51.000 3.040 2.330
4 62.000 4.028 1.791
5 61.000 5.058 2.155
6 55.000 5.939 1.587
7 49.000 6.975 3.425
8 60.000 8.004 2.408
9 53.000 8.972 3.996
if graphics:
ax = myvario.plot(idir=-1)
ax.decoration(title="Multi-Directional Variogram of Pb")
Calculating the Variogram Map
myvmap = gl.db_vmap_compute(mydb,gl.ECalcVario.VARIOGRAM,[20,20])
if verbose:
myvmap.display()
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 5 Total number of samples = 1681 Grid characteristics: --------------------- Origin : -33.160 -29.380 Mesh : 1.658 1.469 Number : 41 41 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = x1 - Locator = x1 Column = 2 - Name = x2 - Locator = x2 Column = 3 - Name = VMAP.Pb.Var - Locator = z1 Column = 4 - Name = VMAP.Pb.Nb - Locator = NA
if graphics:
ax = myvmap.plot("*Var")
ax.decoration(title="Variogram Map")
Fitting a Model. We call the Automatic Fitting procedure providing the list of covariance functions to be tested.
mymodel = gl.Model.createFromDb(mydb)
err = mymodel.fit(myvario,[gl.ECov.EXPONENTIAL,gl.ECov.SPHERICAL])
Visualizing the resulting model, overlaid on the experimental variogram
if graphics:
ax = gp.varmod(myvario,mymodel)
ax.decoration(title="Model for Pb")
We can impose some constraints on the parameters during the fit. For instance here, we impose an equality constraint on the range (range = 1).
myModelConstrained = gl.Model.createFromDb(mydb)
constr = gl.Constraints()
paramid = gl.CovParamId(0,0,gl.EConsElem.RANGE,0,0)
constr.addItem(gl.ConsItem(paramid,gl.EConsType.EQUAL,1.))
err = myModelConstrained.fit(myVarioOmni,[gl.ECov.SPHERICAL],constr)
if (err > 0): print("Error while fitting the model")
myModelConstrained
Model characteristics ===================== Space dimension = 2 Number of variable(s) = 1 Number of basic structure(s) = 1 Number of drift function(s) = 0 Number of drift equation(s) = 0 Covariance Part --------------- Spherical - Sill = 2.101 - Range = 1.000 Total Sill = 2.101
We can impose inequality constraints by using EConsType.LOWER or EConsType.UPPER.
mymodel.setDriftIRF()
if verbose:
mymodel.display()
Model characteristics
=====================
Space dimension = 2
Number of variable(s) = 1
Number of basic structure(s) = 2
Number of drift function(s) = 1
Number of drift equation(s) = 1
Covariance Part
---------------
Exponential
- Sill = 1.035
- Ranges = 1.786 0.366
- Theo. Ranges = 0.596 0.122
- Angles = 45.023 0.000
- Rotation Matrix
[, 0] [, 1]
[ 0,] 0.707 -0.707
[ 1,] 0.707 0.707
Spherical
- Sill = 1.621
- Ranges = 7.051 5.132
- Angles = 136.897 0.000
- Rotation Matrix
[, 0] [, 1]
[ 0,] -0.730 -0.683
[ 1,] 0.683 -0.730
Total Sill = 2.656
Drift Part
----------
Universality_Condition
We initiate a Neigborhood (Moving with a small number of samples for Demonstration)
myneigh = gl.NeighMoving.create(flag_xvalid=False,nmaxi=6,radius=10)
if verbose:
myneigh.display()
Moving Neighborhood =================== Minimum number of samples = 1 Maximum number of samples = 6 Maximum horizontal distance = 10
We must first create a Grid which covers the area of interest
mygrid = gl.DbGrid.createCoveringDb(mydb,[],[0.5,0.5],[],[2,2])
if verbose:
mygrid.display()
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 2 Total number of samples = 5100 Grid characteristics: --------------------- Origin : 107.850 481.660 Mesh : 0.500 0.500 Number : 75 68 Variables --------- Column = 0 - Name = x1 - Locator = x1 Column = 1 - Name = x2 - Locator = x2
We can now test the neighborhood characteristics for each node of the previously defined grid.
err = gl.test_neigh(mydb,mygrid,mymodel,myneigh)
if (err > 0): print("Error while running test_neigh")
if verbose:
mygrid.display()
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 7 Total number of samples = 5100 Grid characteristics: --------------------- Origin : 107.850 481.660 Mesh : 0.500 0.500 Number : 75 68 Variables --------- Column = 0 - Name = x1 - Locator = x1 Column = 1 - Name = x2 - Locator = x2 Column = 2 - Name = Neigh.Pb.Number - Locator = NA Column = 3 - Name = Neigh.Pb.MaxDist - Locator = NA Column = 4 - Name = Neigh.Pb.MinDist - Locator = NA Column = 5 - Name = Neigh.Pb.NbNESect - Locator = NA Column = 6 - Name = Neigh.Pb.NbCESect - Locator = z1
We can visualize some of the newly created variables, such as:
if graphics:
ax = mygrid.plot("Neigh*Number")
ax.decoration(title="Number of Samples per Neighborhood")
if graphics:
ax=mygrid.plot("Neigh*MaxDist")
ax.decoration(title="Maximum Distance per Neighborhood")
We can now process the cross-validation step
err = gl.xvalid(mydb,mymodel,myneigh)
if (err > 0): print("Error while running xvalid")
if verbose:
mydb.display()
Data Base Characteristics ========================= Data Base Summary ----------------- File is organized as a set of isolated points Space dimension = 2 Number of Columns = 8 Total number of samples = 102 Number of active samples = 99 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = X - Locator = x1 Column = 2 - Name = Y - Locator = x2 Column = 3 - Name = Zn - Locator = NA Column = 4 - Name = Pb - Locator = NA Column = 5 - Name = sel - Locator = sel Column = 6 - Name = Xvalid.Pb.esterr - Locator = z1 Column = 7 - Name = Xvalid.Pb.stderr - Locator = NA
if graphics:
gp.histogram(mydb,name="Xvalid.Pb.stderr",bins=50,range=(-4,3))
We now perform the Estimation by Ordinary Kriging. The Neighborhood is changed into a Unique Neighborhood.
mydb.setLocator("Pb",gl.ELoc.Z)
myneigh = gl.NeighUnique.create()
err = gl.kriging(mydb,mygrid,mymodel,myneigh)
if (err > 0): print("Error while running kriging")
if verbose:
mygrid.display(dbfmt)
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 9 Total number of samples = 5100 Grid characteristics: --------------------- Origin : 107.850 481.660 Mesh : 0.500 0.500 Number : 75 68 Data Base Statistics -------------------- 1 - Name x1 - Locator x1 Nb of data = 5100 Nb of active values = 5100 Minimum value = 107.850 Maximum value = 144.850 Mean value = 126.350 Standard Deviation = 10.824 Variance = 117.167 2 - Name x2 - Locator x2 Nb of data = 5100 Nb of active values = 5100 Minimum value = 481.660 Maximum value = 515.160 Mean value = 498.410 Standard Deviation = 9.814 Variance = 96.313 3 - Name Neigh.Pb.Number - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 1.000 Maximum value = 6.000 Mean value = 5.381 Standard Deviation = 1.538 Variance = 2.366 4 - Name Neigh.Pb.MaxDist - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 0.935 Maximum value = 9.999 Mean value = 5.598 Standard Deviation = 2.497 Variance = 6.233 5 - Name Neigh.Pb.MinDist - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 0.014 Maximum value = 9.978 Mean value = 3.484 Standard Deviation = 2.514 Variance = 6.320 6 - Name Neigh.Pb.NbNESect - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 1.000 Maximum value = 1.000 Mean value = 1.000 Standard Deviation = 0.000 Variance = 0.000 7 - Name Neigh.Pb.NbCESect - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 0.000 Maximum value = 0.000 Mean value = 0.000 Standard Deviation = 0.000 Variance = 0.000 8 - Name Kriging.Pb.estim - Locator z1 Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.412 Maximum value = 11.456 Mean value = 6.123 Standard Deviation = 0.648 Variance = 0.419 9 - Name Kriging.Pb.stdev - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 0.248 Maximum value = 1.659 Mean value = 1.541 Standard Deviation = 0.159 Variance = 0.025 Variables --------- Column = 0 - Name = x1 - Locator = x1 Column = 1 - Name = x2 - Locator = x2 Column = 2 - Name = Neigh.Pb.Number - Locator = NA Column = 3 - Name = Neigh.Pb.MaxDist - Locator = NA Column = 4 - Name = Neigh.Pb.MinDist - Locator = NA Column = 5 - Name = Neigh.Pb.NbNESect - Locator = NA Column = 6 - Name = Neigh.Pb.NbCESect - Locator = NA Column = 7 - Name = Kriging.Pb.estim - Locator = z1 Column = 8 - Name = Kriging.Pb.stdev - Locator = NA
Visualizing the results
if graphics:
ax = mygrid.plot("Kriging.Pb.estim")
ax = mydb.plot("Pb")
ax.decoration(title="Estimate of Pb")
if graphics:
ax = mygrid.plot("Kriging.Pb.stdev")
ax = mydb.plot("Pb")
ax.decoration(title="St. Deviation of Pb")
We must first transform the Data into Gaussian
myanamPb = gl.AnamHermite(30)
myanamPb.fitFromLocator(mydb)
if verbose:
myanamPb.display()
Hermitian Anamorphosis
----------------------
Minimum absolute value for Y = -2.7
Maximum absolute value for Y = 2.6
Minimum absolute value for Z = 3.0029
Maximum absolute value for Z = 12.9777
Minimum practical value for Y = -2.7
Maximum practical value for Y = 2.6
Minimum practical value for Z = 3.0029
Maximum practical value for Z = 12.9777
Mean = 5.65758
Variance = 2.86296
Number of Hermite polynomials = 30
Normalized coefficients for Hermite polynomials (punctual variable)
[, 0] [, 1] [, 2] [, 3] [, 4] [, 5] [, 6]
[ 0,] 5.658 -1.625 0.440 -0.069 -0.017 0.082 -0.061
[ 7,] 0.001 0.036 -0.044 0.004 0.047 -0.030 -0.029
[ 14,] 0.037 0.007 -0.031 0.010 0.018 -0.019 -0.003
[ 21,] 0.019 -0.010 -0.014 0.019 0.006 -0.023 0.004
[ 28,] 0.022 -0.013
We can produce the Gaussian Anamorphosis graphically within its definition domain.
if graphics:
gp.anam(myanamPb)
The next step consists in translating the target variable ('Pb') into its Gaussian transform
mydb.setLocator("Pb",gl.ELoc.Z)
err = myanamPb.rawToGaussianByLocator(mydb)
if verbose:
mydb.display()
Data Base Characteristics ========================= Data Base Summary ----------------- File is organized as a set of isolated points Space dimension = 2 Number of Columns = 9 Total number of samples = 102 Number of active samples = 99 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = X - Locator = x1 Column = 2 - Name = Y - Locator = x2 Column = 3 - Name = Zn - Locator = NA Column = 4 - Name = Pb - Locator = NA Column = 5 - Name = sel - Locator = sel Column = 6 - Name = Xvalid.Pb.esterr - Locator = NA Column = 7 - Name = Xvalid.Pb.stderr - Locator = NA Column = 8 - Name = Y.Pb - Locator = z1
We quickly calculate experimental (omni-directional) variograms of the gaussian variable (locater Z) using the already defined directions
myvarioParam = gl.VarioParam()
mydir = gl.DirParam(10,1.)
myvarioParam.addDir(mydir)
myVario = gl.Vario(myvarioParam)
err = myvario.compute(mydb,gl.ECalcVario.VARIOGRAM)
We fit the model by automatic fit (with the constraints that the total sill be equal to 1).
mymodelG = gl.Model.createFromDb(mydb)
err = mymodelG.fit(myvario,[gl.ECov.EXPONENTIAL])
if graphics:
ax = gp.varmod(myvario,mymodelG)
ax.decoration(title="Model for Gaussian Pb")
We perform a set of 10 conditional simulations using the Turning Bands Method.
err = gl.simtub(mydb,mygrid,mymodel,myneigh,10)
if verbose:
mygrid.display()
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 19 Total number of samples = 5100 Grid characteristics: --------------------- Origin : 107.850 481.660 Mesh : 0.500 0.500 Number : 75 68 Variables --------- Column = 0 - Name = x1 - Locator = x1 Column = 1 - Name = x2 - Locator = x2 Column = 2 - Name = Neigh.Pb.Number - Locator = NA Column = 3 - Name = Neigh.Pb.MaxDist - Locator = NA Column = 4 - Name = Neigh.Pb.MinDist - Locator = NA Column = 5 - Name = Neigh.Pb.NbNESect - Locator = NA Column = 6 - Name = Neigh.Pb.NbCESect - Locator = NA Column = 7 - Name = Kriging.Pb.estim - Locator = NA Column = 8 - Name = Kriging.Pb.stdev - Locator = NA Column = 9 - Name = Simu.Y.Pb.1 - Locator = z1 Column = 10 - Name = Simu.Y.Pb.2 - Locator = z2 Column = 11 - Name = Simu.Y.Pb.3 - Locator = z3 Column = 12 - Name = Simu.Y.Pb.4 - Locator = z4 Column = 13 - Name = Simu.Y.Pb.5 - Locator = z5 Column = 14 - Name = Simu.Y.Pb.6 - Locator = z6 Column = 15 - Name = Simu.Y.Pb.7 - Locator = z7 Column = 16 - Name = Simu.Y.Pb.8 - Locator = z8 Column = 17 - Name = Simu.Y.Pb.9 - Locator = z9 Column = 18 - Name = Simu.Y.Pb.10 - Locator = z10
Some statistics on the Conditional simulations in Gaussian scale
mygrid.deleteColumn("Stats.*")
mygrid.statisticsBySample(["Simu.Y.*"],[gl.EStatOption.MINI,
gl.EStatOption.MAXI,
gl.EStatOption.MEAN,
gl.EStatOption.STDV],
True)
if verbose:
dbfmt = gl.DbStringFormat()
dbfmt.setFlags(flag_resume=True, flag_stats=True)
mygrid.display(dbfmt)
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 23 Total number of samples = 5100 Grid characteristics: --------------------- Origin : 107.850 481.660 Mesh : 0.500 0.500 Number : 75 68 Data Base Statistics -------------------- 1 - Name x1 - Locator x1 Nb of data = 5100 Nb of active values = 5100 Minimum value = 107.850 Maximum value = 144.850 Mean value = 126.350 Standard Deviation = 10.824 Variance = 117.167 2 - Name x2 - Locator x2 Nb of data = 5100 Nb of active values = 5100 Minimum value = 481.660 Maximum value = 515.160 Mean value = 498.410 Standard Deviation = 9.814 Variance = 96.313 3 - Name Neigh.Pb.Number - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 1.000 Maximum value = 6.000 Mean value = 5.381 Standard Deviation = 1.538 Variance = 2.366 4 - Name Neigh.Pb.MaxDist - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 0.935 Maximum value = 9.999 Mean value = 5.598 Standard Deviation = 2.497 Variance = 6.233 5 - Name Neigh.Pb.MinDist - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 0.014 Maximum value = 9.978 Mean value = 3.484 Standard Deviation = 2.514 Variance = 6.320 6 - Name Neigh.Pb.NbNESect - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 1.000 Maximum value = 1.000 Mean value = 1.000 Standard Deviation = 0.000 Variance = 0.000 7 - Name Neigh.Pb.NbCESect - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 0.000 Maximum value = 0.000 Mean value = 0.000 Standard Deviation = 0.000 Variance = 0.000 8 - Name Kriging.Pb.estim - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.412 Maximum value = 11.456 Mean value = 6.123 Standard Deviation = 0.648 Variance = 0.419 9 - Name Kriging.Pb.stdev - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 0.248 Maximum value = 1.659 Mean value = 1.541 Standard Deviation = 0.159 Variance = 0.025 10 - Name Simu.Y.Pb.1 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -5.862 Maximum value = 7.142 Mean value = 0.484 Standard Deviation = 1.573 Variance = 2.474 11 - Name Simu.Y.Pb.2 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -4.793 Maximum value = 5.789 Mean value = 0.142 Standard Deviation = 1.568 Variance = 2.459 12 - Name Simu.Y.Pb.3 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -4.882 Maximum value = 6.844 Mean value = 0.277 Standard Deviation = 1.603 Variance = 2.568 13 - Name Simu.Y.Pb.4 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -5.957 Maximum value = 5.443 Mean value = 0.117 Standard Deviation = 1.614 Variance = 2.604 14 - Name Simu.Y.Pb.5 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -5.515 Maximum value = 4.505 Mean value = 0.097 Standard Deviation = 1.432 Variance = 2.052 15 - Name Simu.Y.Pb.6 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -6.198 Maximum value = 5.125 Mean value = 0.214 Standard Deviation = 1.530 Variance = 2.340 16 - Name Simu.Y.Pb.7 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -5.859 Maximum value = 5.812 Mean value = 0.079 Standard Deviation = 1.560 Variance = 2.434 17 - Name Simu.Y.Pb.8 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -4.978 Maximum value = 5.598 Mean value = 0.271 Standard Deviation = 1.589 Variance = 2.526 18 - Name Simu.Y.Pb.9 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -5.340 Maximum value = 6.919 Mean value = 0.683 Standard Deviation = 1.603 Variance = 2.569 19 - Name Simu.Y.Pb.10 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -5.433 Maximum value = 5.800 Mean value = 0.157 Standard Deviation = 1.617 Variance = 2.614 20 - Name Stats.MINI - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -6.198 Maximum value = 1.423 Mean value = -2.101 Standard Deviation = 0.964 Variance = 0.929 21 - Name Stats.MAXI - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -0.843 Maximum value = 7.142 Mean value = 2.624 Standard Deviation = 1.021 Variance = 1.043 22 - Name Stats.MEAN - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -1.983 Maximum value = 2.417 Mean value = 0.252 Standard Deviation = 0.594 Variance = 0.353 23 - Name Stats.STDV - Locator z1 Nb of data = 5100 Nb of active values = 5100 Minimum value = 0.220 Maximum value = 2.788 Mean value = 1.419 Standard Deviation = 0.362 Variance = 0.131 Variables --------- Column = 0 - Name = x1 - Locator = x1 Column = 1 - Name = x2 - Locator = x2 Column = 2 - Name = Neigh.Pb.Number - Locator = NA Column = 3 - Name = Neigh.Pb.MaxDist - Locator = NA Column = 4 - Name = Neigh.Pb.MinDist - Locator = NA Column = 5 - Name = Neigh.Pb.NbNESect - Locator = NA Column = 6 - Name = Neigh.Pb.NbCESect - Locator = NA Column = 7 - Name = Kriging.Pb.estim - Locator = NA Column = 8 - Name = Kriging.Pb.stdev - Locator = NA Column = 9 - Name = Simu.Y.Pb.1 - Locator = NA Column = 10 - Name = Simu.Y.Pb.2 - Locator = NA Column = 11 - Name = Simu.Y.Pb.3 - Locator = NA Column = 12 - Name = Simu.Y.Pb.4 - Locator = NA Column = 13 - Name = Simu.Y.Pb.5 - Locator = NA Column = 14 - Name = Simu.Y.Pb.6 - Locator = NA Column = 15 - Name = Simu.Y.Pb.7 - Locator = NA Column = 16 - Name = Simu.Y.Pb.8 - Locator = NA Column = 17 - Name = Simu.Y.Pb.9 - Locator = NA Column = 18 - Name = Simu.Y.Pb.10 - Locator = NA Column = 19 - Name = Stats.MINI - Locator = NA Column = 20 - Name = Stats.MAXI - Locator = NA Column = 21 - Name = Stats.MEAN - Locator = NA Column = 22 - Name = Stats.STDV - Locator = z1
We visualize a conditional simulation in Gaussian scale
if graphics:
ax = mygrid.plot("Simu.Y.Pb.1")
ax = mydb.plot("Pb")
ax.decoration(title="One Simulation of Pb in Gaussian Scale")
We turn the Gaussian conditional simulations into Raw scale (using the Anamorphosis back transform) and get rid of the Gaussian conditional simulations.
myanamPb.gaussianToRaw(mygrid,"Simu.Y.*")
mygrid.deleteColumn("Simu.Y.*")
mygrid.deleteColumn("Stats.*")
if verbose:
mygrid.display()
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 19 Total number of samples = 5100 Grid characteristics: --------------------- Origin : 107.850 481.660 Mesh : 0.500 0.500 Number : 75 68 Variables --------- Column = 0 - Name = x1 - Locator = x1 Column = 1 - Name = x2 - Locator = x2 Column = 2 - Name = Neigh.Pb.Number - Locator = NA Column = 3 - Name = Neigh.Pb.MaxDist - Locator = NA Column = 4 - Name = Neigh.Pb.MinDist - Locator = NA Column = 5 - Name = Neigh.Pb.NbNESect - Locator = NA Column = 6 - Name = Neigh.Pb.NbCESect - Locator = NA Column = 7 - Name = Kriging.Pb.estim - Locator = NA Column = 8 - Name = Kriging.Pb.stdev - Locator = NA Column = 9 - Name = Z.Simu.Y.Pb.1 - Locator = z1 Column = 10 - Name = Z.Simu.Y.Pb.2 - Locator = z2 Column = 11 - Name = Z.Simu.Y.Pb.3 - Locator = z3 Column = 12 - Name = Z.Simu.Y.Pb.4 - Locator = z4 Column = 13 - Name = Z.Simu.Y.Pb.5 - Locator = z5 Column = 14 - Name = Z.Simu.Y.Pb.6 - Locator = z6 Column = 15 - Name = Z.Simu.Y.Pb.7 - Locator = z7 Column = 16 - Name = Z.Simu.Y.Pb.8 - Locator = z8 Column = 17 - Name = Z.Simu.Y.Pb.9 - Locator = z9 Column = 18 - Name = Z.Simu.Y.Pb.10 - Locator = z10
We calculate some statistics on the Conditional Simulations in Raw scale.
mygrid.deleteColumn("Stats.*")
mygrid.statisticsBySample(["Z.Simu.*"],[gl.EStatOption.MINI,
gl.EStatOption.MAXI,
gl.EStatOption.MEAN,
gl.EStatOption.STDV],
True)
if verbose:
dbfmt = gl.DbStringFormat()
dbfmt.setFlags(flag_resume=True, flag_stats=True)
mygrid.display(dbfmt)
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 23 Total number of samples = 5100 Grid characteristics: --------------------- Origin : 107.850 481.660 Mesh : 0.500 0.500 Number : 75 68 Data Base Statistics -------------------- 1 - Name x1 - Locator x1 Nb of data = 5100 Nb of active values = 5100 Minimum value = 107.850 Maximum value = 144.850 Mean value = 126.350 Standard Deviation = 10.824 Variance = 117.167 2 - Name x2 - Locator x2 Nb of data = 5100 Nb of active values = 5100 Minimum value = 481.660 Maximum value = 515.160 Mean value = 498.410 Standard Deviation = 9.814 Variance = 96.313 3 - Name Neigh.Pb.Number - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 1.000 Maximum value = 6.000 Mean value = 5.381 Standard Deviation = 1.538 Variance = 2.366 4 - Name Neigh.Pb.MaxDist - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 0.935 Maximum value = 9.999 Mean value = 5.598 Standard Deviation = 2.497 Variance = 6.233 5 - Name Neigh.Pb.MinDist - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 0.014 Maximum value = 9.978 Mean value = 3.484 Standard Deviation = 2.514 Variance = 6.320 6 - Name Neigh.Pb.NbNESect - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 1.000 Maximum value = 1.000 Mean value = 1.000 Standard Deviation = 0.000 Variance = 0.000 7 - Name Neigh.Pb.NbCESect - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 0.000 Maximum value = 0.000 Mean value = 0.000 Standard Deviation = 0.000 Variance = 0.000 8 - Name Kriging.Pb.estim - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.412 Maximum value = 11.456 Mean value = 6.123 Standard Deviation = 0.648 Variance = 0.419 9 - Name Kriging.Pb.stdev - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 0.248 Maximum value = 1.659 Mean value = 1.541 Standard Deviation = 0.159 Variance = 0.025 10 - Name Z.Simu.Y.Pb.1 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.871 Standard Deviation = 2.959 Variance = 8.757 11 - Name Z.Simu.Y.Pb.2 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.270 Standard Deviation = 2.785 Variance = 7.759 12 - Name Z.Simu.Y.Pb.3 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.532 Standard Deviation = 2.904 Variance = 8.434 13 - Name Z.Simu.Y.Pb.4 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.264 Standard Deviation = 2.835 Variance = 8.036 14 - Name Z.Simu.Y.Pb.5 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.125 Standard Deviation = 2.548 Variance = 6.492 15 - Name Z.Simu.Y.Pb.6 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.397 Standard Deviation = 2.731 Variance = 7.456 16 - Name Z.Simu.Y.Pb.7 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.174 Standard Deviation = 2.709 Variance = 7.340 17 - Name Z.Simu.Y.Pb.8 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.504 Standard Deviation = 2.871 Variance = 8.245 18 - Name Z.Simu.Y.Pb.9 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 7.261 Standard Deviation = 3.079 Variance = 9.481 19 - Name Z.Simu.Y.Pb.10 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.340 Standard Deviation = 2.852 Variance = 8.133 20 - Name Stats.MINI - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 8.083 Mean value = 3.412 Standard Deviation = 0.570 Variance = 0.324 21 - Name Stats.MAXI - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 4.282 Maximum value = 12.978 Mean value = 11.377 Standard Deviation = 2.109 Variance = 4.446 22 - Name Stats.MEAN - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.438 Maximum value = 11.597 Mean value = 6.474 Standard Deviation = 1.063 Variance = 1.131 23 - Name Stats.STDV - Locator z1 Nb of data = 5100 Nb of active values = 5100 Minimum value = 0.349 Maximum value = 4.357 Mean value = 2.539 Standard Deviation = 0.738 Variance = 0.544 Variables --------- Column = 0 - Name = x1 - Locator = x1 Column = 1 - Name = x2 - Locator = x2 Column = 2 - Name = Neigh.Pb.Number - Locator = NA Column = 3 - Name = Neigh.Pb.MaxDist - Locator = NA Column = 4 - Name = Neigh.Pb.MinDist - Locator = NA Column = 5 - Name = Neigh.Pb.NbNESect - Locator = NA Column = 6 - Name = Neigh.Pb.NbCESect - Locator = NA Column = 7 - Name = Kriging.Pb.estim - Locator = NA Column = 8 - Name = Kriging.Pb.stdev - Locator = NA Column = 9 - Name = Z.Simu.Y.Pb.1 - Locator = NA Column = 10 - Name = Z.Simu.Y.Pb.2 - Locator = NA Column = 11 - Name = Z.Simu.Y.Pb.3 - Locator = NA Column = 12 - Name = Z.Simu.Y.Pb.4 - Locator = NA Column = 13 - Name = Z.Simu.Y.Pb.5 - Locator = NA Column = 14 - Name = Z.Simu.Y.Pb.6 - Locator = NA Column = 15 - Name = Z.Simu.Y.Pb.7 - Locator = NA Column = 16 - Name = Z.Simu.Y.Pb.8 - Locator = NA Column = 17 - Name = Z.Simu.Y.Pb.9 - Locator = NA Column = 18 - Name = Z.Simu.Y.Pb.10 - Locator = NA Column = 19 - Name = Stats.MINI - Locator = NA Column = 20 - Name = Stats.MAXI - Locator = NA Column = 21 - Name = Stats.MEAN - Locator = NA Column = 22 - Name = Stats.STDV - Locator = z1
We visualize a Conditional Simulation in Raw Scale
if graphics:
ax = mygrid.plot("Z.Simu.Y.Pb.1")
ax = mydb.plot("Pb")
ax.decoration(title="One simulation of Pb in Raw Scale")
Let us now average the conditional simulations in order to have a comparison with the estimation by kriging.
mygrid.deleteColumn("Stats.*")
mygrid.statisticsBySample(["Z.Simu.*"],[gl.EStatOption.MEAN],True)
if verbose:
dbfmt = gl.DbStringFormat()
dbfmt.setFlags(flag_resume=True, flag_stats=True)
mygrid.display(dbfmt)
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 20 Total number of samples = 5100 Grid characteristics: --------------------- Origin : 107.850 481.660 Mesh : 0.500 0.500 Number : 75 68 Data Base Statistics -------------------- 1 - Name x1 - Locator x1 Nb of data = 5100 Nb of active values = 5100 Minimum value = 107.850 Maximum value = 144.850 Mean value = 126.350 Standard Deviation = 10.824 Variance = 117.167 2 - Name x2 - Locator x2 Nb of data = 5100 Nb of active values = 5100 Minimum value = 481.660 Maximum value = 515.160 Mean value = 498.410 Standard Deviation = 9.814 Variance = 96.313 3 - Name Neigh.Pb.Number - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 1.000 Maximum value = 6.000 Mean value = 5.381 Standard Deviation = 1.538 Variance = 2.366 4 - Name Neigh.Pb.MaxDist - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 0.935 Maximum value = 9.999 Mean value = 5.598 Standard Deviation = 2.497 Variance = 6.233 5 - Name Neigh.Pb.MinDist - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 0.014 Maximum value = 9.978 Mean value = 3.484 Standard Deviation = 2.514 Variance = 6.320 6 - Name Neigh.Pb.NbNESect - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 1.000 Maximum value = 1.000 Mean value = 1.000 Standard Deviation = 0.000 Variance = 0.000 7 - Name Neigh.Pb.NbCESect - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 0.000 Maximum value = 0.000 Mean value = 0.000 Standard Deviation = 0.000 Variance = 0.000 8 - Name Kriging.Pb.estim - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.412 Maximum value = 11.456 Mean value = 6.123 Standard Deviation = 0.648 Variance = 0.419 9 - Name Kriging.Pb.stdev - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 0.248 Maximum value = 1.659 Mean value = 1.541 Standard Deviation = 0.159 Variance = 0.025 10 - Name Z.Simu.Y.Pb.1 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.871 Standard Deviation = 2.959 Variance = 8.757 11 - Name Z.Simu.Y.Pb.2 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.270 Standard Deviation = 2.785 Variance = 7.759 12 - Name Z.Simu.Y.Pb.3 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.532 Standard Deviation = 2.904 Variance = 8.434 13 - Name Z.Simu.Y.Pb.4 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.264 Standard Deviation = 2.835 Variance = 8.036 14 - Name Z.Simu.Y.Pb.5 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.125 Standard Deviation = 2.548 Variance = 6.492 15 - Name Z.Simu.Y.Pb.6 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.397 Standard Deviation = 2.731 Variance = 7.456 16 - Name Z.Simu.Y.Pb.7 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.174 Standard Deviation = 2.709 Variance = 7.340 17 - Name Z.Simu.Y.Pb.8 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.504 Standard Deviation = 2.871 Variance = 8.245 18 - Name Z.Simu.Y.Pb.9 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 7.261 Standard Deviation = 3.079 Variance = 9.481 19 - Name Z.Simu.Y.Pb.10 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.340 Standard Deviation = 2.852 Variance = 8.133 20 - Name Stats.MEAN - Locator z1 Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.438 Maximum value = 11.597 Mean value = 6.474 Standard Deviation = 1.063 Variance = 1.131 Variables --------- Column = 0 - Name = x1 - Locator = x1 Column = 1 - Name = x2 - Locator = x2 Column = 2 - Name = Neigh.Pb.Number - Locator = NA Column = 3 - Name = Neigh.Pb.MaxDist - Locator = NA Column = 4 - Name = Neigh.Pb.MinDist - Locator = NA Column = 5 - Name = Neigh.Pb.NbNESect - Locator = NA Column = 6 - Name = Neigh.Pb.NbCESect - Locator = NA Column = 7 - Name = Kriging.Pb.estim - Locator = NA Column = 8 - Name = Kriging.Pb.stdev - Locator = NA Column = 9 - Name = Z.Simu.Y.Pb.1 - Locator = NA Column = 10 - Name = Z.Simu.Y.Pb.2 - Locator = NA Column = 11 - Name = Z.Simu.Y.Pb.3 - Locator = NA Column = 12 - Name = Z.Simu.Y.Pb.4 - Locator = NA Column = 13 - Name = Z.Simu.Y.Pb.5 - Locator = NA Column = 14 - Name = Z.Simu.Y.Pb.6 - Locator = NA Column = 15 - Name = Z.Simu.Y.Pb.7 - Locator = NA Column = 16 - Name = Z.Simu.Y.Pb.8 - Locator = NA Column = 17 - Name = Z.Simu.Y.Pb.9 - Locator = NA Column = 18 - Name = Z.Simu.Y.Pb.10 - Locator = NA Column = 19 - Name = Stats.MEAN - Locator = z1
Displaying the average of the Conditional Simulations
if graphics:
ax = mygrid.plot("Stats*MEAN")
ax = mydb.plot("Pb")
ax.decoration(title="Mean of Pb simulations")
The Gaussian transform of the Pb variable has already been calculated. It suffices to perform the Gaussian transform of the Zn variable
mydb.setLocator("Zn",gl.ELoc.Z)
myanamZn = gl.AnamHermite(30)
myanamZn.fitFromLocator(mydb)
if verbose:
myanamZn.display()
Hermitian Anamorphosis
----------------------
Minimum absolute value for Y = -2.5
Maximum absolute value for Y = 2.6
Minimum absolute value for Z = 1.1469
Maximum absolute value for Z = 12.1276
Minimum practical value for Y = -2.5
Maximum practical value for Y = 2.6
Minimum practical value for Z = 1.1469
Maximum practical value for Z = 12.1276
Mean = 2.88061
Variance = 2.76263
Number of Hermite polynomials = 30
Normalized coefficients for Hermite polynomials (punctual variable)
[, 0] [, 1] [, 2] [, 3] [, 4] [, 5] [, 6]
[ 0,] 2.881 -1.277 0.877 -0.447 -0.095 0.294 -0.121
[ 7,] -0.087 0.134 -0.029 -0.087 0.069 0.034 -0.065
[ 14,] 0.005 0.044 -0.026 -0.020 0.034 0.001 -0.033
[ 21,] 0.010 0.027 -0.016 -0.019 0.016 0.012 -0.014
[ 28,] -0.005 0.011
if graphics:
gp.anam(myanamZn)
We convert the raw data into its Gaussian equivalent
mydb.setLocator("Zn",gl.ELoc.Z)
err = myanamZn.rawToGaussianByLocator(mydb)
if verbose:
mydb.display()
Data Base Characteristics ========================= Data Base Summary ----------------- File is organized as a set of isolated points Space dimension = 2 Number of Columns = 10 Total number of samples = 102 Number of active samples = 99 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = X - Locator = x1 Column = 2 - Name = Y - Locator = x2 Column = 3 - Name = Zn - Locator = NA Column = 4 - Name = Pb - Locator = NA Column = 5 - Name = sel - Locator = sel Column = 6 - Name = Xvalid.Pb.esterr - Locator = NA Column = 7 - Name = Xvalid.Pb.stderr - Locator = NA Column = 8 - Name = Y.Pb - Locator = NA Column = 9 - Name = Y.Zn - Locator = z1
We now perform the multivariate variogram calculation
mydb.setLocators(["Y.Pb","Y.Zn"],gl.ELoc.Z)
myvario = gl.Vario(myvarioParam)
err = myvario.compute(mydb,gl.ECalcVario.VARIOGRAM)
mymodelM = gl.Model.createFromDb(mydb)
err = mymodelM.fit(myvario,[gl.ECov.EXPONENTIAL])
if graphics:
axs = gp.varmod(myvario,mymodelM)
gp.decoration(axs,title="Multivariate Model")
gp.geometry(axs,dims=[5,5])
We perform 10 bivariate conditional simulations (deleting the previous monovariate simulation outcomes first for better legibility)
mygrid.deleteColumn("Z.Simu*")
err = gl.simtub(mydb,mygrid,mymodelM,myneigh,nbsimu=10)
if verbose:
mygrid.display()
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 30 Total number of samples = 5100 Grid characteristics: --------------------- Origin : 107.850 481.660 Mesh : 0.500 0.500 Number : 75 68 Variables --------- Column = 0 - Name = x1 - Locator = x1 Column = 1 - Name = x2 - Locator = x2 Column = 2 - Name = Neigh.Pb.Number - Locator = NA Column = 3 - Name = Neigh.Pb.MaxDist - Locator = NA Column = 4 - Name = Neigh.Pb.MinDist - Locator = NA Column = 5 - Name = Neigh.Pb.NbNESect - Locator = NA Column = 6 - Name = Neigh.Pb.NbCESect - Locator = NA Column = 7 - Name = Kriging.Pb.estim - Locator = NA Column = 8 - Name = Kriging.Pb.stdev - Locator = NA Column = 9 - Name = Stats.MEAN - Locator = NA Column = 10 - Name = Simu.Y.Pb.1 - Locator = z1 Column = 11 - Name = Simu.Y.Pb.2 - Locator = z2 Column = 12 - Name = Simu.Y.Pb.3 - Locator = z3 Column = 13 - Name = Simu.Y.Pb.4 - Locator = z4 Column = 14 - Name = Simu.Y.Pb.5 - Locator = z5 Column = 15 - Name = Simu.Y.Pb.6 - Locator = z6 Column = 16 - Name = Simu.Y.Pb.7 - Locator = z7 Column = 17 - Name = Simu.Y.Pb.8 - Locator = z8 Column = 18 - Name = Simu.Y.Pb.9 - Locator = z9 Column = 19 - Name = Simu.Y.Pb.10 - Locator = z10 Column = 20 - Name = Simu.Y.Zn.1 - Locator = z11 Column = 21 - Name = Simu.Y.Zn.2 - Locator = z12 Column = 22 - Name = Simu.Y.Zn.3 - Locator = z13 Column = 23 - Name = Simu.Y.Zn.4 - Locator = z14 Column = 24 - Name = Simu.Y.Zn.5 - Locator = z15 Column = 25 - Name = Simu.Y.Zn.6 - Locator = z16 Column = 26 - Name = Simu.Y.Zn.7 - Locator = z17 Column = 27 - Name = Simu.Y.Zn.8 - Locator = z18 Column = 28 - Name = Simu.Y.Zn.9 - Locator = z19 Column = 29 - Name = Simu.Y.Zn.10 - Locator = z20
We back-transform each set of simulation outcomes using its own Gaussian Anamorphosis function. Finally we delete the Gaussian variables and ask for the statistics on the simulated variables in the Raw Scale.
err = myanamZn.gaussianToRaw(mygrid,"Simu.Y.Zn*")
err = myanamPb.gaussianToRaw(mygrid,"Simu.Y.Pb*")
mygrid.deleteColumn("Simu.Y*")
mygrid.deleteColumn("Stats.*")
mygrid.statisticsBySample(["Z.Simu.*"],[gl.EStatOption.MINI,
gl.EStatOption.MAXI,
gl.EStatOption.MEAN,
gl.EStatOption.STDV],
True)
if verbose:
dbfmt = gl.DbStringFormat()
dbfmt.setFlags(flag_resume=True, flag_stats=True)
mygrid.display(dbfmt)
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 33 Total number of samples = 5100 Grid characteristics: --------------------- Origin : 107.850 481.660 Mesh : 0.500 0.500 Number : 75 68 Data Base Statistics -------------------- 1 - Name x1 - Locator x1 Nb of data = 5100 Nb of active values = 5100 Minimum value = 107.850 Maximum value = 144.850 Mean value = 126.350 Standard Deviation = 10.824 Variance = 117.167 2 - Name x2 - Locator x2 Nb of data = 5100 Nb of active values = 5100 Minimum value = 481.660 Maximum value = 515.160 Mean value = 498.410 Standard Deviation = 9.814 Variance = 96.313 3 - Name Neigh.Pb.Number - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 1.000 Maximum value = 6.000 Mean value = 5.381 Standard Deviation = 1.538 Variance = 2.366 4 - Name Neigh.Pb.MaxDist - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 0.935 Maximum value = 9.999 Mean value = 5.598 Standard Deviation = 2.497 Variance = 6.233 5 - Name Neigh.Pb.MinDist - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 0.014 Maximum value = 9.978 Mean value = 3.484 Standard Deviation = 2.514 Variance = 6.320 6 - Name Neigh.Pb.NbNESect - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 1.000 Maximum value = 1.000 Mean value = 1.000 Standard Deviation = 0.000 Variance = 0.000 7 - Name Neigh.Pb.NbCESect - Locator NA Nb of data = 5100 Nb of active values = 4596 Minimum value = 0.000 Maximum value = 0.000 Mean value = 0.000 Standard Deviation = 0.000 Variance = 0.000 8 - Name Kriging.Pb.estim - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.412 Maximum value = 11.456 Mean value = 6.123 Standard Deviation = 0.648 Variance = 0.419 9 - Name Kriging.Pb.stdev - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 0.248 Maximum value = 1.659 Mean value = 1.541 Standard Deviation = 0.159 Variance = 0.025 10 - Name Z.Simu.Y.Zn.1 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 12.128 Mean value = 2.805 Standard Deviation = 1.580 Variance = 2.498 11 - Name Z.Simu.Y.Zn.2 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 12.128 Mean value = 2.733 Standard Deviation = 1.408 Variance = 1.983 12 - Name Z.Simu.Y.Zn.3 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 12.128 Mean value = 2.781 Standard Deviation = 1.566 Variance = 2.453 13 - Name Z.Simu.Y.Zn.4 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 12.128 Mean value = 2.907 Standard Deviation = 1.630 Variance = 2.658 14 - Name Z.Simu.Y.Zn.5 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 12.128 Mean value = 2.900 Standard Deviation = 1.551 Variance = 2.405 15 - Name Z.Simu.Y.Zn.6 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 12.128 Mean value = 2.919 Standard Deviation = 1.629 Variance = 2.654 16 - Name Z.Simu.Y.Zn.7 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 12.128 Mean value = 2.840 Standard Deviation = 1.500 Variance = 2.251 17 - Name Z.Simu.Y.Zn.8 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 12.128 Mean value = 2.851 Standard Deviation = 1.607 Variance = 2.581 18 - Name Z.Simu.Y.Zn.9 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 12.128 Mean value = 2.760 Standard Deviation = 1.395 Variance = 1.947 19 - Name Z.Simu.Y.Zn.10 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 12.128 Mean value = 2.934 Standard Deviation = 1.701 Variance = 2.895 20 - Name Z.Simu.Y.Pb.1 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.745 Standard Deviation = 1.770 Variance = 3.133 21 - Name Z.Simu.Y.Pb.2 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.728 Standard Deviation = 1.769 Variance = 3.130 22 - Name Z.Simu.Y.Pb.3 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.822 Standard Deviation = 1.811 Variance = 3.281 23 - Name Z.Simu.Y.Pb.4 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.634 Standard Deviation = 1.715 Variance = 2.940 24 - Name Z.Simu.Y.Pb.5 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.639 Standard Deviation = 1.642 Variance = 2.695 25 - Name Z.Simu.Y.Pb.6 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.742 Standard Deviation = 1.714 Variance = 2.939 26 - Name Z.Simu.Y.Pb.7 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.742 Standard Deviation = 1.853 Variance = 3.435 27 - Name Z.Simu.Y.Pb.8 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.699 Standard Deviation = 1.751 Variance = 3.065 28 - Name Z.Simu.Y.Pb.9 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.726 Standard Deviation = 1.784 Variance = 3.184 29 - Name Z.Simu.Y.Pb.10 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.784 Standard Deviation = 1.837 Variance = 3.374 30 - Name Stats.MINI - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 4.252 Mean value = 1.839 Standard Deviation = 0.287 Variance = 0.082 31 - Name Stats.MAXI - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.300 Maximum value = 12.978 Mean value = 9.189 Standard Deviation = 2.051 Variance = 4.206 32 - Name Stats.MEAN - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 2.386 Maximum value = 9.627 Mean value = 4.285 Standard Deviation = 0.491 Variance = 0.241 33 - Name Stats.STDV - Locator z1 Nb of data = 5100 Nb of active values = 5100 Minimum value = 0.498 Maximum value = 4.139 Mean value = 2.093 Standard Deviation = 0.485 Variance = 0.235 Variables --------- Column = 0 - Name = x1 - Locator = x1 Column = 1 - Name = x2 - Locator = x2 Column = 2 - Name = Neigh.Pb.Number - Locator = NA Column = 3 - Name = Neigh.Pb.MaxDist - Locator = NA Column = 4 - Name = Neigh.Pb.MinDist - Locator = NA Column = 5 - Name = Neigh.Pb.NbNESect - Locator = NA Column = 6 - Name = Neigh.Pb.NbCESect - Locator = NA Column = 7 - Name = Kriging.Pb.estim - Locator = NA Column = 8 - Name = Kriging.Pb.stdev - Locator = NA Column = 9 - Name = Z.Simu.Y.Zn.1 - Locator = NA Column = 10 - Name = Z.Simu.Y.Zn.2 - Locator = NA Column = 11 - Name = Z.Simu.Y.Zn.3 - Locator = NA Column = 12 - Name = Z.Simu.Y.Zn.4 - Locator = NA Column = 13 - Name = Z.Simu.Y.Zn.5 - Locator = NA Column = 14 - Name = Z.Simu.Y.Zn.6 - Locator = NA Column = 15 - Name = Z.Simu.Y.Zn.7 - Locator = NA Column = 16 - Name = Z.Simu.Y.Zn.8 - Locator = NA Column = 17 - Name = Z.Simu.Y.Zn.9 - Locator = NA Column = 18 - Name = Z.Simu.Y.Zn.10 - Locator = NA Column = 19 - Name = Z.Simu.Y.Pb.1 - Locator = NA Column = 20 - Name = Z.Simu.Y.Pb.2 - Locator = NA Column = 21 - Name = Z.Simu.Y.Pb.3 - Locator = NA Column = 22 - Name = Z.Simu.Y.Pb.4 - Locator = NA Column = 23 - Name = Z.Simu.Y.Pb.5 - Locator = NA Column = 24 - Name = Z.Simu.Y.Pb.6 - Locator = NA Column = 25 - Name = Z.Simu.Y.Pb.7 - Locator = NA Column = 26 - Name = Z.Simu.Y.Pb.8 - Locator = NA Column = 27 - Name = Z.Simu.Y.Pb.9 - Locator = NA Column = 28 - Name = Z.Simu.Y.Pb.10 - Locator = NA Column = 29 - Name = Stats.MINI - Locator = NA Column = 30 - Name = Stats.MAXI - Locator = NA Column = 31 - Name = Stats.MEAN - Locator = NA Column = 32 - Name = Stats.STDV - Locator = z1
We compare the initial variable 'Pb' with a set of disjoint intervals. The 'Pb' values varying from 3 to 12.7, we consider three classes:
We first build the indicators for each class
limits = gl.Limits([np.nan, 4., 6., np.nan])
if verbose:
limits.display()
Bound( 1 ) : ] -Inf ; 4 [ Bound( 2 ) : [ 4 ; 6 [ Bound( 3 ) : [ 6 ; +Inf [
We apply the set of limits previously defined in order to transform the input variable into Indicators of the different classes.
err = limits.toIndicator(mydb,"Pb")
if verbose:
mydb.display()
Data Base Characteristics ========================= Data Base Summary ----------------- File is organized as a set of isolated points Space dimension = 2 Number of Columns = 13 Total number of samples = 102 Number of active samples = 99 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = X - Locator = x1 Column = 2 - Name = Y - Locator = x2 Column = 3 - Name = Zn - Locator = NA Column = 4 - Name = Pb - Locator = NA Column = 5 - Name = sel - Locator = sel Column = 6 - Name = Xvalid.Pb.esterr - Locator = NA Column = 7 - Name = Xvalid.Pb.stderr - Locator = NA Column = 8 - Name = Y.Pb - Locator = NA Column = 9 - Name = Y.Zn - Locator = NA Column = 10 - Name = Indicator.Pb.Class.1 - Locator = z1 Column = 11 - Name = Indicator.Pb.Class.2 - Locator = z2 Column = 12 - Name = Indicator.Pb.Class.3 - Locator = z3
We calculate the variogram of the Indicators for future use
myvarioindParam = gl.VarioParam()
myvarioindParam.addDir(mydir)
myvarioInd = gl.Vario(myvarioindParam)
err = myvarioInd.compute(mydb,gl.ECalcVario.VARIOGRAM)
if verbose:
myvarioInd.display()
Variogram characteristics
=========================
Number of variable(s) = 3
Number of direction(s) = 1
Space dimension = 2
Variance-Covariance Matrix
[, 0] [, 1] [, 2]
[ 0,] 0.107 -0.062 -0.044
[ 1,] -0.062 0.250 -0.187
[ 2,] -0.044 -0.187 0.231
Direction #1
------------
Number of lags = 10
Direction coefficients = 1.000 0.000
Direction angles (degrees) = 0.000 0.000
Tolerance on direction = 90.000 (degrees)
Calculation lag = 1.000
Tolerance on distance = 50.000 (Percent of the lag value)
For variable 1
Rank Npairs Distance Value
0 3.000 0.389 0.000
1 123.000 1.081 0.081
2 183.000 2.038 0.126
3 205.000 3.006 0.156
4 231.000 4.013 0.132
5 229.000 5.036 0.159
6 198.000 5.962 0.152
7 187.000 7.000 0.107
8 204.000 7.996 0.096
9 184.000 8.990 0.068
For variables 2 and 1
Rank Npairs Distance Value
0 3.000 0.389 0.000
1 123.000 1.081 -0.065
2 183.000 2.038 -0.077
3 205.000 3.006 -0.085
4 231.000 4.013 -0.093
5 229.000 5.036 -0.085
6 198.000 5.962 -0.061
7 187.000 7.000 -0.045
8 204.000 7.996 -0.042
9 184.000 8.990 -0.038
For variable 2
Rank Npairs Distance Value
0 3.000 0.389 0.167
1 123.000 1.081 0.199
2 183.000 2.038 0.221
3 205.000 3.006 0.251
4 231.000 4.013 0.292
5 229.000 5.036 0.258
6 198.000 5.962 0.237
7 187.000 7.000 0.254
8 204.000 7.996 0.228
9 184.000 8.990 0.234
For variables 3 and 1
Rank Npairs Distance Value
0 3.000 0.389 0.000
1 123.000 1.081 -0.016
2 183.000 2.038 -0.049
3 205.000 3.006 -0.071
4 231.000 4.013 -0.039
5 229.000 5.036 -0.074
6 198.000 5.962 -0.091
7 187.000 7.000 -0.061
8 204.000 7.996 -0.054
9 184.000 8.990 -0.030
For variables 3 and 2
Rank Npairs Distance Value
0 3.000 0.389 -0.167
1 123.000 1.081 -0.134
2 183.000 2.038 -0.145
3 205.000 3.006 -0.166
4 231.000 4.013 -0.199
5 229.000 5.036 -0.172
6 198.000 5.962 -0.177
7 187.000 7.000 -0.209
8 204.000 7.996 -0.186
9 184.000 8.990 -0.196
For variable 3
Rank Npairs Distance Value
0 3.000 0.389 0.167
1 123.000 1.081 0.150
2 183.000 2.038 0.194
3 205.000 3.006 0.237
4 231.000 4.013 0.238
5 229.000 5.036 0.247
6 198.000 5.962 0.268
7 187.000 7.000 0.270
8 204.000 7.996 0.240
9 184.000 8.990 0.226
ax = gp.varmod(myvarioInd)
gp.geometry(ax,dims=[4,4])
Then we build a categorical variable which gives the index of the class to which each sample belongs
err = limits.toCategory(mydb,"Pb")
if verbose:
dbfmt = gl.DbStringFormat()
dbfmt.setFlags(flag_stats=True)
dbfmt.setNames(["Category*"])
dbfmt.setMode(2)
mydb.display(dbfmt)
Data Base Characteristics ========================= Data Base Summary ----------------- File is organized as a set of isolated points Space dimension = 2 Number of Columns = 14 Total number of samples = 102 Number of active samples = 99 Data Base Statistics -------------------- 14 - Name Category.Pb - Locator z1 Nb of data = 102 Nb of active values = 99 Class 1 = 12 ( 12.121%) Class 2 = 51 ( 51.515%) Class 3 = 36 ( 36.364%) Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = X - Locator = x1 Column = 2 - Name = Y - Locator = x2 Column = 3 - Name = Zn - Locator = NA Column = 4 - Name = Pb - Locator = NA Column = 5 - Name = sel - Locator = sel Column = 6 - Name = Xvalid.Pb.esterr - Locator = NA Column = 7 - Name = Xvalid.Pb.stderr - Locator = NA Column = 8 - Name = Y.Pb - Locator = NA Column = 9 - Name = Y.Zn - Locator = NA Column = 10 - Name = Indicator.Pb.Class.1 - Locator = NA Column = 11 - Name = Indicator.Pb.Class.2 - Locator = NA Column = 12 - Name = Indicator.Pb.Class.3 - Locator = NA Column = 13 - Name = Category.Pb - Locator = z1