Re: Using the STk package for Kriging

[Julien Bect]

Le 14/11/2016 à 16:54, Francesco Potortì a écrit :

Francesco Potortì:

I am experimenting with Kriging interpolation and I am studying the STK
package.

I am surprised at the fact that the provided covariance functions are
Metérn, Gauss, and empirical.  The commonly used exponential and
spherical models are missing, like other common models.

Julien Bect:

I have added the exponential and spherical covariance functions to STK 

It will be available in STK 2.4.0, which I hope should be released 
soon... (within a month or two).

In the meantime, you can try the beta version that is available from 
SourceForge: <http://people.ku.edu/~gbohling/cpe940>.

Thanks a lot, my colleagues are working on it now!

Since I am at it, here are some comments on the STK package.

I have had some difficulty in using STK because the concepts and
terminology are different form what I have read in the literature.  The
main example is the introductory paper on Kriging by Geoff Bohling that
you find at <http://people.ku.edu/~gbohling/cpe940>.  But I have found
the same concepts and terminology used in other papers and in the
"Kriging" and "variogram" Wikipedia pages.

Hi Francesco,

Thanks for sharing these difficulties.

The reason is that "kriging" (a.k.a. Gaussian process (GP) modelling) is used in several domains : geostastistics, DACE, machine learning, Bayesian (nonparametric) statistics, numerical analysis (RBF), etc.  Therefore, there are sometimes several names given to to same thing, or (even worse from beginners) several equivalent mathematical descriptions for the same "object".

Possible solutions, as I see it :

 * A long term solution would be a manual with introductions to STK for people coming from various backgrounds (but this is not likely to happen soon).

 * A more realistic solution would be to add a "geostatistics" subdirectory in the example directory, where examples are given with the vocabulary and tools of geostatistics.

 * An even more flexible solution would be to add a page on the Octave wiki explaining STK to people coming from the field of geostatistics (which seems to be your case).  Since the wiki is open, perhaps can you help to make this happen ?

Specifically, your code does not mention the variogram and does not directly use parameters such as nugget, sill and range.

theoretical variogram : in STK we prefer the concept of covariance functions (actually, generalized covariance functions, i.e., conditionnally positive definite kernels).

empirical variogram : it is true that there is no function to plot an empirical variogram.  We don't provide variogram-based methods for parameter estimation, but restricted maximum likelihood (see Michael Stein's book, for instance, for a discussion of both types of approaches).  Also, empirical variograms are not commonly used outside geostatistics (e.g., in DACE or machine learning) since they are mostly useful for 2D or 3D isotropic data.  But none of this is a real excuse ;-)  and I agree that having some tools in STK related to empirical variograms would be a nice addition for some users.

nugget : this term is very specific to geostatistics -> the equivalent concept in STK is "noise variance"  (Note : nugget and noise variance are equivalent when there are no replicated measurements.)

sill : this is even more specific to geostatistics -> what you call the sill is equal in STK to the sum of the noise variance and the GP variance.

range : ah ah, this one is used !   See for instance "help stk_sphcov" ;-)

Since I am a beginner at this, I would appreciate if you could please help me check
that I have correctly understood how to use the STK code.

For example, when using the new spherical variogram model, these are the
parameters that I should use:

sill is the long-range variance
nugget is the variance at null distance
range is where the variance becomes practically equal to the sill

model.lognoisevariance = log(sill);
model.name = 'stk_spherical_iso';
model.param = [log(sill-nugget); -log(range)];

Is that correct?

No.  Based on what I said above, it should be:

model.lognoisevariance = log(nugget);

I think that the learning curve for the STK package would be much
smoother if the concepts of variogram, sill, nugget and range were
explicitely mentioned and an example were given on how to convert them
for use in the various models in STK.

I agree, this would make it easier for people with a geostatistical background to get started with STK.  See above: you can definitely contribute on this point.

By the way, if you have some interesting geostatistical datatsets that are (or can be) published under an appropriate licence (say, CC0) I am interested as well.

Please keep us informed of your progress and difficulties with STK.

@++
Julien