Non stationary covariances¶
This tutorial aims to show how to handle non stationary covariances for kriging, covariance evaluation, and SPDE approach.
For the covariance based approach (used in the gstlearn function kriging) the covariance formulas can be obtain in this monograph of M. Stein.
For the SPDE approach, a detailed material can be obtained in the Phd thesis of M. Pereira.
In gstlearn, there are two ways to specify the non stationary covariance :
- by providing maps of parameters through a Db which covers the domain
- by providing a function (AFunctional structure of gstlearn)
We will explain the two approaches.
We first create a reference model.
import gstlearn as gl
import numpy as np
model =gl.Model.createFromParam(gl.ECov.MATERN,range=10,param=1, sill = 2)
We extract the CovAniso structure
cova = model.getCova(0)
cova
Matern (Third Parameter = 1) - Sill = 2.000 - Range = 10.000 - Theo. Range = 2.887
We start by the first approach (non-stationarity defined through a Db).
We define a Db of parameters named grid1 for which we fill some columns of values (arbitrary in this tutorial).
grid1 = gl.DbGrid.create([10,10])
grid1["param1"] = 2 * grid1["x1"] +12
grid1["param2"] = 2 * grid1["x2"] +5
Then we can make the different parameters of the covariance non stationary by first attaching the Db and then providing the column name of the non stationary parameter. For instance:
cova.attachNoStatDb(grid1)
cova.makeAngleNoStatDb("param1")
cova.makeRangeNoStatDb("param2",1) #Range parameter in the second direction
cova
Matern (Third Parameter = 1) - Sill = 2.000 - Range = 10.000 - Theo. Range = 2.887 Non-Stationary Parameters ------------------------- 1 - Range : IDir=2 2 - Angle : IdAngle=1
Another way of doing this, is to specify the Db when you give the specification of the parameter. In that case, you don't need to explicitely attach a Db with attachDbNoStat.
model =gl.Model.createFromParam(gl.ECov.MATERN,range=10,param=1, sill = 2)
cova = model.getCova(0)
cova.makeAngleNoStatDb("param1",db = grid1)
cova.makeRangeNoStatDb("param2",1)
cova
Matern (Third Parameter = 1) - Sill = 2.000 - Range = 10.000 - Theo. Range = 2.887 Non-Stationary Parameters ------------------------- 1 - Range : IDir=2 2 - Angle : IdAngle=1
In that case, grid1 is recorded as the reference Db for the next parameters.
You can even have different Db for different parameters.
grid2 = gl.DbGrid.create([10,10])
grid2["paramnew"] = 3 * grid2["x2"] +24
cova.makeSillNoStatDb("paramnew",db = grid2)
cova
Matern (Third Parameter = 1) - Sill = 2.000 - Range = 10.000 - Theo. Range = 2.887 Non-Stationary Parameters ------------------------- 1 - Sill : Ivar=0 Jvar=0 2 - Range : IDir=2 3 - Angle : IdAngle=1
But the reference Db is the first one which has been attached (by attachNoStatDb or by the first parameter specification)
cova.makeRangeNoStatDb("paramnew")
No variable name corresponding to your criterion You have to specified a name of a column of the reference Db
Of course, you have to provide a Db to be allowed to give a specification:
model =gl.Model.createFromParam(gl.ECov.MATERN,range=10,param=1, sill = 2)
cova = model.getCova(0)
cova.makeAngleNoStatDb("param1")
cova
You have to define a Db (with attachNoStatDb or by specifying a Db here)
Matern (Third Parameter = 1) - Sill = 2.000 - Range = 10.000 - Theo. Range = 2.887
If you give several specifications for a given parameter, you have a warning message and the last specification is kept.
model =gl.Model.createFromParam(gl.ECov.MATERN,range=10,param=1, sill = 2)
cova = model.getCova(0)
cova.makeAngleNoStatDb("param1",db = grid1)
cova.makeAngleNoStatDb("param2")
Warning, this non stationarity was already specified. It has been replaced with the new specifications.
You can specify the non stationarity by using the scale instead of the practical range.
model =gl.Model.createFromParam(gl.ECov.MATERN,range=10,param=1, sill = 2)
cova = model.getCova(0)
cova.makeScaleNoStatDb("param1",db = grid1)
cova
Matern (Third Parameter = 1) - Sill = 2.000 - Range = 10.000 - Theo. Range = 2.887 Non-Stationary Parameters ------------------------- 1 - Scale : IDir=1
Same king of message if you specify the range after having specified the scale (for the same direction).
model =gl.Model.createFromParam(gl.ECov.MATERN,range=10,param=1, sill = 2)
cova = model.getCova(0)
cova.makeScaleNoStatDb("param1",db = grid1)
cova.makeRangeNoStatDb("param1",db = grid1)
cova
Warning, you gave a non-stationary specification for the range but it was already given for the scale. The new specification has replaced the previous one.
Matern (Third Parameter = 1) - Sill = 2.000 - Range = 10.000 - Theo. Range = 2.887 Non-Stationary Parameters ------------------------- 1 - Range : IDir=1
You can't work in scale for one dimension and in range for another one:
model =gl.Model.createFromParam(gl.ECov.MATERN,range=10,param=1, sill = 2)
cova = model.getCova(0)
cova.makeScaleNoStatDb("param1",0,db = grid1)
cova.makeRangeNoStatDb("param1",1,db = grid1)
cova
You try to specify non stationarities for range whereas you had already specified one for the scale in another dimension. It is invalid
Matern (Third Parameter = 1) - Sill = 2.000 - Range = 10.000 - Theo. Range = 2.887 Non-Stationary Parameters ------------------------- 1 - Scale : IDir=1
You can make back a parameter stationary
model =gl.Model.createFromParam(gl.ECov.MATERN,range=10,param=1, sill = 2)
cova = model.getCova(0)
cova.makeScaleNoStatDb("param1",0,db = grid1)
cova.makeScaleStationary(0)
cova
Matern (Third Parameter = 1) - Sill = 2.000 - Range = 10.000 - Theo. Range = 2.887
If you make a range stationary whereas the non-stationarity was previously defined for scale, the parameter is made stationary as well.
model =gl.Model.createFromParam(gl.ECov.MATERN,range=10,param=1, sill = 2)
cova = model.getCova(0)
cova.makeScaleNoStatDb("param1",0,db = grid1)
cova.makeRangeStationary(0)
cova
Matern (Third Parameter = 1) - Sill = 2.000 - Range = 10.000 - Theo. Range = 2.887
You can clear all the non-stationarity specifications.
model =gl.Model.createFromParam(gl.ECov.MATERN,range=10,param=1, sill = 2)
cova = model.getCova(0)
cova.makeScaleNoStatDb("param1",0,db = grid1)
cova.makeAngleNoStatDb("param1",0,db = grid1)
cova.makeStationary()
cova
Matern (Third Parameter = 1) - Sill = 2.000 - Range = 10.000 - Theo. Range = 2.887
The non-stationarity is transfered when you clone your CovAniso
model =gl.Model.createFromParam(gl.ECov.MATERN,range=10,param=1, sill = 2)
cova = model.getCova(0)
cova.makeScaleNoStatDb("param1",0,db = grid1)
cova.makeAngleNoStatDb("param1",0,db = grid1)
cova2 = cova.clone()
cova2
Matern (Third Parameter = 1) - Sill = 2.000 - Range = 10.000 - Theo. Range = 2.887 Non-Stationary Parameters ------------------------- 1 - Angle : IdAngle=1 2 - Scale : IDir=1
You can specify non-stationarities for the terms of the sill matrix in the multivariate setting:
model =gl.Model.createFromParam(gl.ECov.MATERN,range=10,param=1, sills = np.array([2.,1.,1.,3.]))
cova = model.getCova(0)
cova.makeSillNoStatDb("param1",0,1,db = grid1)
cova
Matern (Third Parameter = 1)
- Sill matrix:
[, 0] [, 1]
[ 0,] 2.000 1.000
[ 1,] 1.000 3.000
- Range = 10.000
- Theo. Range = 2.887
Non-Stationary Parameters
-------------------------
1 - Sill : Ivar=0 Jvar=1
Non stationarity by a function¶
To understand how to play with a function, let's have a look on a part of the AFunctional C++ class of the gstlearn library.
We can see that the class is built by using an integer (the space dimension), and that there is mainly to methods. The first one is pure virtual, that means that the users who want to create a class which inherits from AFunctional have to implement this method. The aim of this method is to give the value at a given position (pos) specified through a vector of dimension _ndim.
Here is an example of class in python to inherit from the AFunctional class.
class MyFunction(gl.AFunctional):
def __init__(self,ndim):
super(MyFunction,self).__init__(ndim)
pass
def getFunctionValue(self, coords):
return coords[0] + 1
You can then instantiate an object from the class MyFunction and use it to define non stationarity:
myfunc = MyFunction(2)
model =gl.Model.createFromParam(gl.ECov.MATERN,range=10,param=1, sills = np.array([2.,1.,1.,3.]))
cova = model.getCova(0)
cova.makeSillNoStatFunctional(myfunc,0,1)
cova
Matern (Third Parameter = 1)
- Sill matrix:
[, 0] [, 1]
[ 0,] 2.000 1.000
[ 1,] 1.000 3.000
- Range = 10.000
- Theo. Range = 2.887
Non-Stationary Parameters
-------------------------
1 - Sill : Ivar=0 Jvar=1
Functional
You can even mix different types of non-stationarities (with AFunctional and with NoStatDb):
model =gl.Model.createFromParam(gl.ECov.MATERN,range=10,param=1, sills = np.array([2.,1.,1.,3.]))
cova = model.getCova(0)
cova.makeSillNoStatDb("param1",0,0,grid1)
cova.makeSillNoStatFunctional(myfunc,0,1)
cova.makeSillNoStatDb("param1",1,1,grid1)
cova
Matern (Third Parameter = 1)
- Sill matrix:
[, 0] [, 1]
[ 0,] 2.000 1.000
[ 1,] 1.000 3.000
- Range = 10.000
- Theo. Range = 2.887
Non-Stationary Parameters
-------------------------
1 - Sill : Ivar=1 Jvar=1
2 - Sill : Ivar=0 Jvar=1
Functional
3 - Sill : Ivar=0 Jvar=0