Kriging is an optimal, linear estimation method that produces an UNBIASED, MINIMUM ESTIMATION ERROR VARIANCE estimate of a function of an unknown random field Z(u), The form of the estimator is,

> P:=Sum(lambda[i]*Z(u[i]),i=1...n);

The unbiasedness condition leads to an optimal set of weights that must meet the following restriction,

> Sum(lambda[i],i=1...n)=1;

The optimal set of weights is obtained by minimizing the estimation error variance. The minimization problem leads to a set of (n+1) linear equations in lambda that can be solved using standard linear algebra procedures.

The form of these equations depends on the kind of function being estimated. If the function is an integral of the unknown random field Z(u), as it would be in estimating the Mean Areal Precipitation as,

> MAP:=1/A*Int(Z(u),u=A..A);

then the set of equations is of the form,

> Sum(lambda[j]*cov(u[i]-u[j]),j=1..n)+mu=1/A*Int(cov(u[i]-u),u=A..A); Sum(lambda[i],i=1...n)=1;

If the function is the value of the random field Z(u) at a point , as it would be in interpolating Z(u), then the set of equations is of the form,

> Sum(lambda[j]*cov(u[i]-u[j]),j=1..n)+mu=cov(u[i]-u); Sum(lambda[i],i=1...n)=1;restart;

These sets of linear equations can be recast into matrix form such that,

where [A] is the matrix of coefficients, [l] is the vector of weights, which includes the Lagrange multiplier m, and [B] is the independent or constant vector. Note that this vector depends on the tipe of function being estimated.

Following is an example of the application of Kriging in the estimation of the two kinds of functions indicated above.

Let Z(u) represent the precipitation along a straight line, as a function of position u on the line.
a) The total length of the line, LL, is 3.
b) There are three raingauges located at coordinates u1 = 0, u2 = 1 and u3 = 3.
c) The total rainfall measurements at the three gauges are: Z(u1) = 3, Z(u2) = 8, and Z(u3) = 1.
d) The covariance function of Z(u) is:

> cov:=(h)->sigma^2*exp(-abs(h)/L);

In the covariance equation above let the standard deviation, s,be equal to 1. L is the characteristic length scale of the covariance function.

1) Use Kriging to estimate the value of the mean total precipitation over the line. This problem is equivalent to obtaining an estimate of the integral of an unknown function, Z(u), but of which a limited number of values are known (the observations). Thus, Kriging can be used as an alternative to other commonly used methods for numerical approximation of definite integrals.
 
 

> with(linalg);with(student);

> sigma:=1; L:=0.5;LL:=3;

Define the vector of observations Z(u), and the location, u, of the gauges. Observe that although only three gauges are defined, an additional element is included which will include the location of the unknown value and the corresponding estimate of that value.

> z:=matrix(4,1,[3,8,1,0]);u:=matrix(4,1,[0,1,3,0]);

Define the matrix of coefficients, [A], and obtain its inverse,

> A:=matrix(4,4,[[evalf(cov(u[1,1]-u[1,1])), evalf(cov(u[1,1]-u[2,1])), evalf(cov(u[1,1]-u[3,1])), 1], [evalf(cov(u[2,1]-u[1,1])), evalf(cov(u[2,1]-u[2,1])), evalf(cov(u[2,1]-u[3,1])), 1], [evalf(cov(u[3,1]-u[1,1])), evalf(cov(u[3,1]-u[2,1])), evalf(cov(u[3,1]-u[3,1])), 1], [1,1,1,0]]);

> AINV:=inverse(A);

Define the independent vector [B1] for the case of mean areal precipitation estimation.

> B1:=matrix(4,1,[evalf(int(cov(x-u[1,1]), x=0..LL)/LL), evalf(int(cov(x-u[2,1]), x=0..LL)/LL), evalf(int(cov(x-u[3,1]),x=0..LL)/LL), 1]);

Compute the set of optimal weights for mean areal precipitation estimation.

> lambda:=multiply(AINV,B1);

Compute the mean areal precipitation along the entire line,

> zave := multiply(transpose(lambda),z);

Compute the variance of the estimation error,

> variance := evalf(simpson(evalf(simpson(cov(x-y),y=0..LL,4)),x=0..LL,4))/LL^2 - 2/LL*add(lambda[k,1]*evalf(int(cov(x-u[k,1]),x=0..LL)),k=1..3) + sum(lambda[k,1]*sum(lambda[l,1]*cov(abs(u[k,1]-u[l,1])),l=1..3),k=1..3);

2) Use Kriging to interpolate the function Z(u) along the line from 0 to 3. Plot both the interpolated function Z(u) and the corresponding estimation error variance.

Define procedure to do optimal interpolation using Kriging.

> krig := proc() local x,u,B0,lambda,sigma,zo: fopen(kriging,WRITE,TEXT): for x from 0 by 0.1 while x <= 3 do u:=matrix(4,1,[0,1,3,x]);
B0:=matrix(4,1,[evalf(cov(u[4,1]-u[1,1])), evalf(cov(u[4,1]-u[2,1])), evalf(cov(u[4,1]-u[3,1])), 1]);
lambda:=multiply(AINV,B0): zo:=evalf(multiply(transpose(lambda),z));
sigma:=cov(0)-2*sum(lambda[k,1]*cov(abs(u[4,1]-u[k,1])), k=1..3)+sum(lambda[k,1]*sum(lambda[l,1]*cov(abs(u[k,1]-u[l,1])),l=1..3),k=1..3):
fprintf(kriging,`%g,%g,%g,\n`,x,zo[1,1],evalf(sigma)): od: fclose(kriging): end;

> krig();