Automatic Structure Fitting

D. Renard, N. Desassis

21 August 2017

Abstract

RGeostats is a library of geostatistical procedures available in the R environment. It is based on the geostatistical library Geoslib written in C language. The specificity of RGeostats is that it can deal with infor- mation defined in a space of any dimension, and can process any number of variables simultaneously. This paper is meant to demonstrate the Automatic Structure Fitting procedure. It covers all the possibilities of this procedure, starting from case of a simply digitized model, to the more complex fit of a variogram map. In particular, it will cover the case of complex anisotropy, the 2-D and 3-D case.

Fitting a digitized variogram

This is the most simple illustration of the quality of the Automatic Model Fitting procedure.

Isotropic structure

The principle is to define an isotropic model composed of:

For the demonstration sake, this model is stored as an R object (called Exdemo_autofit1.model) which belongs to the RGeostats package that can be loaded by typing the command:

data(Exdemo_autofit1.model)

This model is sampled by calculating an experimental variogram with 100 lags of 0.1 each. Obviously the experimental variogram exactly represents the model (with no statistical fluctuation):

vario <- model.sample(Exdemo_autofit1.model,lag=0.1,nlag=100)

The experimental variogram is displayed in the following figure.

plot(vario,title="Isotropic Variogram to be Fitted")

We now wish to fit a model to this experimental variogram using the default parameters, i.e. a mixture of basic structures such as Nugget Effect, Gaussian, Cubic, Exponential and Spherical components:

model <- model.auto(vario,maxiter=100,title="Isotropic Model")

The previous command automatically fits the optimal model (which will be stored in the object called model, and plots the model together with the experimental variogram, as shown in the next figure. As the fit is perfect, we cannot distinguish the two curves.

When comparing the characteristics of the fitted model to the ones of the digitized one, we can check that the procedure has performed perfectly.

Anisotropic structure

The same work-flow is performed on an model in the 2-D space with nested structures and complex anisotropies:

data(Exdemo_autofit2.model)
vario <- model.sample(Exdemo_autofit2.model,lag=0.1,nlag=100,
      dirvect=c(0,45,90,135))
plot(vario,title="Anisotropic Variogram to be Fitted")

The model is sampled in the 4 main directions (0, 45, 90 and 135 degrees) and represented in the next figure.

model <- model.auto(vario,maxiter=500,
                    title="Anisotropic Model")

print(model)
## 
## 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
## ---------------
## - Cubic
##   Sill        =      1.621
##   Anisotropy  =      1.000     0.391
##   Ranges      =     13.400     5.235
##   Rotation    =
##                 [,1]      [,2]
##       [1,]     0.712    -0.702
##       [2,]     0.702     0.712
##   Angle       =    44.615 degrees
## - Exponential
##   Sill        =      1.990
##   Anisotropy  =      1.000     0.546
##   Ranges      =      5.786     3.161
##   Theo. Ranges=      1.931     1.055
##   Rotation    =
##                 [,1]      [,2]
##       [1,]     0.774    -0.633
##       [2,]     0.633     0.774
##   Angle       =    39.245 degrees
## Total Sill    =      3.611
## 
## Drift Part
## ----------
## Universality Condition

Note that the characteristics of the model are found by the Automatic Model Fitting even if the four directions along which the experimental variograms have been calculated do not match main orientation of the anisotropy ellipsoid.

Fitting a real experimental variogram

In this second part, we perform the Automatic Structure Fitting from real experimental variograms.

The 2-D case

We define a 2-D model first composed of:

We calculate one simulation outcome (using the turning band method) of the model on a fine 2-D square grid with 100 nodes of 0.1 along each axis. The experimental variograms calculated from this grid in the four main directions is displayed in the next figure:

data(Exdemo_autofit3.vario)
plot(Exdemo_autofit3.vario,
    title="Variogram in 2-D space to be Fitted")

This variogram (computed in 4 directions) is fitted using the Automatic Model Fitting procedure. The resulting model is displayed together with the experimental variograms in the next figure:

model <- model.auto(Exdemo_autofit3.vario,maxiter=200,
        title="Model in 2-D space")

model
## 
## Model characteristics
## =====================
## Space dimension              = 2
## Number of variable(s)        = 1
## Number of basic structure(s) = 3
## Number of drift function(s)  = 1
## Number of drift equation(s)  = 1
## 
## Covariance Part
## ---------------
## - Nugget Effect
##   Sill        =      1.186
## - Gaussian
##   Sill        =      0.979
##   Anisotropy  =      1.000     0.882
##   Ranges      =      3.580     3.157
##   Theo. Ranges=      2.068     1.824
##   Rotation    =
##                 [,1]      [,2]
##       [1,]     0.382    -0.924
##       [2,]     0.924     0.382
##   Angle       =    67.559 degrees
## - Spherical
##   Sill        =      5.511
##   Anisotropy  =      0.048     1.000
##   Ranges      =     35.310   742.877
##   Rotation    =
##                 [,1]      [,2]
##       [1,]    -0.612    -0.791
##       [2,]     0.791    -0.612
##   Angle       =   127.754 degrees
## Total Sill    =      7.676
## 
## Drift Part
## ----------
## Universality Condition

The fit is graphically acceptable, even if its characteristics are quite different from the ones of the initial model:

The 3-D case

We perform the same type of operation in the 3-D space. The model is composed of:

data(Exdemo_autofit4.vario)
plot(Exdemo_autofit4.vario,
     title="Variogram in 3-D space to be Fitted")

After 1000 iterations, the model perfectly reproduces the experimental var- iogram as shown in the following figure:

model <- model.auto(Exdemo_autofit4.vario,maxiter=100,
        title="Model in 3-D space")

The characteristics of the fitted model are:

Fitting in the multivariate case

This example deals with a 2-D case with a large number of variables (18). The experimental variogram already has been calculated and stored as an R object that can be downloaded by simply typing:

data(Exdemo_autofit5.vario)

The interested user can visualize the experimental variogram. However, due to the large number of variables, there would be 18*19/2=171 simple and cross- variograms, each one represented on a very small portion of the figure. This is the reason why, by default, only the simple and cross-variograms of the first five variables are visualized.

plot(Exdemo_autofit5.vario)

The automatic model fitting procedure is performed using the defaulted five basic structures (Nugget Effect, Gaussian, Cubic, Exponential and Spherical components) and having the Goulard option switched on.

model <- model.auto(Exdemo_autofit5.vario,maxiter=100,
                    draw=FALSE)
## 
## Convergence not reached after 100 iterations (4 parameters)

As there is only one calculation direction, we are not looking for any anisotropy: each model is characterized by its range (except the Nugget Effect) and the ma- trix of sills (of dimension 18*18). Using the Goulard option reduces the number of parameters to be fitted within the Foxleg algorithm down to 4 parameters (the ranges of each structure). All five basic structures are kept in the final model which is represented in the next figure:

plot(Exdemo_autofit5.vario,varcols=seq(1,5))
plot(model,Exdemo_autofit5.vario,varcols=seq(1,5),add=T)

Breaking the ambiguity

This example has been requested by a reviewer of the paper presenting the method. We use a regular 2-D square grid with 100 nodes and a mesh of 1m: the grid covers 100m by 100m. We simulate the outcome of a random function which follows a structure composed of three nested isotropic spherical components:

The variogram is calculated along the four main directions (two axes and two bisectors) of the grid which gives the following graphic representation. Note that, although the model is isotropic, the different experimental variograms shows some statistical fluctuations that could be interpreted as anisotropy.

data(Exdemo_autofit6.vario)
plot(Exdemo_autofit6.vario,
    title="Directional variograms on the grid")

The first attempt consists in running the automatic model fitting procedure using the defaulted five basic structures (Nugget Effect, Gaussian, Cubic, Exponential and Spherical components), leading to the following result composed of anisotropic exponential and spherical components. The fit is reasonably correct although the convergence has not been reached after 1000 iterations.

model <- model.auto(Exdemo_autofit6.vario,maxiter=100,
        title="Directional Model on the grid")
## 
## Convergence not reached after 100 iterations (12 parameters)

The next trial is performed by defining the set of possible basic structures as a combinaison of three spherical components, as follows:

As expected, the resulting model is composed of three nested anisotropic spherical structures:

model <- model.auto(Exdemo_autofit6.vario,maxiter=100,
      struct=melem.name(c(3,3,3)),title="Fitting with 3 sphericals")

Fitting a Variogram Map

The principle is to create a square grid with 200 nodes and a mesh of 1m:

grid <- db.create(nx=c(200,200))

The model is stored as an R object and can be loaded using the statement:

data(Exdemo_autofit7.model)

A simulation outcome is calculated over the grid using the Turning Bands simulation technique (with 1000 bands):

grid  <- simtub(,grid,Exdemo_autofit7.model)

Plotting the simulated outcome produces the following figure:

plot(grid,title="Simulation Outcome")

The next step consists in calculating the variogram map on the previous outcome: the extension of the variogram map is fixed to 61 by 61 nodes (i.e. the radius is equal to 30 nodes).

vmap  <- vmap.grid(grid,nx=30,ny=30,flag.fft=TRUE)

which leads to the following variogram map representation:

plot(vmap,title="Experimental Vmap")

The Automatic Variogram Map Fitting procedure is performed by defining the set of possible basic structures to a combinaison of a nugget effect and two spherical components:

model <- vmap.auto(vmap,maxiter=100,
                   struct=melem.name(c(1,3,3)))

Although the convergence has not been reached after 1000 iterations, the graphic result is quite convincing, where the initial experimental variogram map is displayed together with the theoretical variogram map obtained with the fitted model: