Estimation

Multidimensional kriging estimation with noisy data observations can be performed using “krig”

    [d_est,d_var,lambda_sk,K_dd,k_du,inhood]=krig(pos_known,val_known,pos_est,V,options);
    
    
    ndata : number of data observations
    ndims : dimensions of data location (>=1)
    nest  : number of data locations to be estimated
    
    pos_known [ndata,ndims] : Locations of data observations
    val_known [ndata,1 or 2]  : col1 : Data value as measured at 'pos_known'
                                col2 : Data uncertainty as measured at
    'pos_known' (optional)
    pos_est   [1 ,ndims] : Location of data to be estimated
    V : Variogram model, e.g. '1 Sph(100)'
    options : kriging options.
    
  

Please note data uncertainty is listed as the second column of val_known. Data uncertainty is defined as the variance of Gaussian noise distribution associated to each data measurement. Specified in this manner, the noise on data observations is considered uncorrelated. If only one column is specified, data observations are treated as noise free.

Kriging Options

options, the last argument, is a Matlab structure that controls most aspect of the kriging, such as choosing the neighborhood, and kriging method. The names for most of these options are same as the names used to control GSTAT, see GSTAT from Matlab

Table 2.1. kriging options

option.[?]GSTAT equivdescriptionlink
sk_mean,meansk_meanSpecify the global mean, as used by simple kriging 
trendtrendKrig only the trend. 
maxmaxThe maximum number of data to use in the search neighborhood 
polytrenddThe polynomial order of the trend for each dimension 
xvalidxvalidIf specified as xvalid=1 cross validation on known data locations is performed 
isorange isorange=1 assumes no rotation in anisotropy. If chosen, each entry in the range section of a variogram corresponds to the range for a certain dimension. For example 'Sph(1,10,1,100)' corresponds to a range of 1 in the first dimension, 10 in the second dimension, 1 in the third dimension and 100 in the fourth dimension. If GSTAT is used for kriging in 1, 2 and 3 dimensions, the range selections are translates properly to GSTAT format using the isorange function. 


kriging methods

the Kriging methods: simple Kriging, ordinary Kriging and Kriging with a trend are all available using krig. The difference between these three methods is the way the trend is modelled. Simple Kriging assumed a constant and known mean. Ordinary kriging makes use of an varying unknown mean (that is estimated within the neighborhood). Kriging with a trend model the trend as a smoothly varying polynomial trend. Thus ordinary kriging is but a form of kriging with a trend,; a 0th order polynomial trend.

Simple kriging

By default ordinary Kriging in an exhaustive neighborhood (i.e. all data measurements are used all the time) is performed. This is identical to simple Kriging when options.mean=mean(pos_known(:,1))

Ordinary kriging

If the 'mean' is given as options.mean, Simple Kriging is performed.

Kriging with a trend / Universal kriging

For each dimension (direction) the order of the polynomial fit of the trend can be specified by options.polyfit

if the value of options.polyfit is an integer scalar, then the same polynomial fit is used for all dimensions. Thus the default linear trend in all directions is identical to specifying options.polyfit=1.

if the value of options.polyfit is an array of integers, each integer in the array must specify the order of the trend for each dimension. The length of the array must be equal to the size of the dimension of the observation. Thus, for 2D data observations one can specify a 2nd order polynomial trend in the first direction and a 0th order fit in the second direction as options.polyfit=[2 0].

kriging neighborhood

By default all data observations are always used. The kriging neighborhood denote the data that are actually used by the kriging systems.

A maximum number of data used by the kriging system is selected as options.max=10. The data locations closes to the point being estimated is retained in the data neighborhood.

See also “nhood” that controls the data neighborhood behavior.

Kriging examples

Based on the data observations below, a number of example are shown illustrating the use of kriging in mGstat

1D Kriging, no data uncertainty

In the following example 3 data measurements has been made : x(1)=0; x(5)=3; x(2)=2. Using a Spherical semivariogram model with a range of '0.2' and a sill of '1', the mean and variance of the distribution of the local probability density function at x(2) is found:

pos_known=[1;5;10];
val_known=[0;3;2]; % adding some uncertainty
V='1 Sph(.2)';      % Select variogram model
pos_est=[2]';
[d_est,d_var]=krig(pos_known,val_known,pos_est,V)


d_est =

    1.6667


d_var =

    1.3333

from: mGstat/examples/mgstat_examples/krig_ex1.m

Kriging a series of points - noise free data

pos_known=[1;5;10];
val_known=[0;3;2]; % 
V='1 Sph(.2)';      % Select variogram model
pos_est=[0:.1:10]';
[d_est,d_var]=krig(pos_known,val_known,pos_est,V);
plot(pos_est,d_est,'k',pos_known,val_known,'ro')
print -dpng krigex2

Kriging a series of point - SK, OK, Ktrend

pos_known=[1;5;10]; %
val_known=[0;3;2];  % 
V='1 Sph(.2)';      % Select variogram model
pos_est=[0:.1:10]';
[d_est_ok,d_var_ok]=krig(pos_known,val_known,pos_est,V);
options.mean=2;
[d_est_sk,d_var_sk]=krig(pos_known,val_known,pos_est,V,options);
options.trend=1;
[d_est_kt,d_var_kt]=krig(pos_known,val_known,pos_est,V,options);
plot(pos_est,[d_est_sk,d_est_ok,d_est_kt],'-',pos_known,val_known,'ro')
legend('SK','OK','KT','Data')
print -dpng krigex3

Kriging a series of point - SK, OK, Ktrend

rand('seed',1)
ndata=10;
pos_known=rand(ndata,1)*10;
val_known=randn(ndata,1); % 
V='1 Sph(.2)';      % Select variogram model
pos_est=[0:.1:10]';
clear options;
[d_est_ok,d_var_ok]=krig(pos_known,val_known,pos_est,V,options);
options.mean=2;
[d_est_sk,d_var_sk]=krig(pos_known,val_known,pos_est,V,options);
options.trend=1;
[d_est_kt,d_var_kt]=krig(pos_known,val_known,pos_est,V,options);
plot(pos_est,[d_est_sk,d_est_ok,d_est_kt],'-',pos_known,val_known,'k*')
legend('SK','OK','KT','Data')
print -dpng krigex4

Kriging a series of point - SK, OK, Ktrend - Neighborhood

rand('seed',1)
ndata=10;
pos_known=rand(ndata,1)*10;
val_known=randn(ndata,1); % 
V='1 Sph(.2)';      % Select variogram model
pos_est=[0:.1:10]';
clear options;options.max=4;
[d_est_ok,d_var_ok]=krig(pos_known,val_known,pos_est,V,options);
options.mean=2;
[d_est_sk,d_var_sk]=krig(pos_known,val_known,pos_est,V,options);
options.trend=1;
[d_est_kt,d_var_kt]=krig(pos_known,val_known,pos_est,V,options);
plot(pos_est,[d_est_sk,d_est_ok,d_est_kt],'-',pos_known,val_known,'k*')
legend('SK','OK','KT','Data')
print -dpng krigex5

Kriging a series of point - SK, OK, Ktrend - Neighborhood, noisy data

rand('seed',1)
ndata=30;
pos_known=rand(ndata,1)*10;
val_known=randn(ndata,1); % 
val_var=zeros(ndata,1)+.1; % 
V='1 Sph(.1)';      % Select variogram model
pos_est=[-2:.1:12]';
clear options;options.max=4;
[d_est_ok,d_var_ok]=krig(pos_known,[val_known val_var],pos_est,V,options);
options.mean=mean(val_known);
[d_est_sk,d_var_sk]=krig(pos_known,[val_known val_var],pos_est,V,options);
options.trend=1;
[d_est_kt,d_var_kt]=krig(pos_known,[val_known val_var],pos_est,V,options);
plot(pos_est,[d_est_sk,d_est_ok,d_est_kt],'-',pos_known,val_known,'k*')
legend('SK','OK','KT','Data')
print -dpng krigex6