# Ordinary kriging stationary case

I am trying to understand ordinary kriging.

Say I have 3 elevation measurements: Z1, Z2, and Z3 taken at X positions: X1, X2 and X3.

I am also assuming some semivariogram: g(h) and that the process is stationary.

From what I understand the kriging weights in this case should be ?:

$L = M^{-1}*Gs$

Where:

$L=\ \left( \begin{array}{ccc} l_1 \\ l_2 \\ l_3 \\ \mu\end{array} \right)$

$\ M = \left( \begin{array}{ccc} g(0) && g(X2-X1) && g(X3-X1) && 1 \\ g(X2-X1) && g(0) && g(X3-X2) && 1 \\ g(X3-X1) && g(X3-X2) && g(0) && 1 \\ 1 && 1 && 1 && 0 \end{array} \right)$

$Gs=\ \left( \begin{array}{ccc} g(|X1 - Xs|) \\ g(|X2 - Xs|) \\ g(|X3 - Xs|) \\ 1 \end{array} \right)$

And the elevation Zs at x position Xs should be estimated by ?:

$\ Z(Xs)=\sum_{i=1}^3 l_i*Zi$

Will the inverse of the matrix pose a problem? I found a Matlab script on the internet that was using "Moore-Penrose pseudoinverse" instead.

The reason I am asking this question is that I coded a Matlab script on this which drew random X and Z; and it seemed to work just fine most times, but every now and then it returned a "crazy" estimate for Z(Xs).

So I am just trying to understand if I have misunderstood the theory, if the inverse of the matrix sometimes become "ill posed" or if I have just tripped up in my Matlab coding.

Here is my Matlab script:

% kriging.m
n = 3;
X = 1000*rand(1,n);
Z = 1000 + 100*randn(1,n);
stem(X,Z,'b');
hold on;

a0 = 300; % range
c0 =   4; % scale

M = ones(n+1,n+1);
for r = 1 : n
for c = 1 : n
Xr = X(r);
Xc = X(c);
h = abs(Xr - Xc);
M(r,c) = c0*(1-exp( -(h/a0)^2 ));
end
end
M(n+1,n+1) = 0

MI = inv(M)

Xs = 500;
hs = abs(X - Xs);
Gs = c0*(1-exp( -(hs/a0).^2 ));
Gs = [Gs 1]'

L = MI*Gs;
L = L(1:n)

Zs = Z*L

stem(Xs,Zs,'r')


And here is an example of an crazy result:

X =

117.4177  296.6759  318.7783

Z =

1.0e+003 *

0.9679    1.0012    0.6971

M =

0    1.2010    1.4508    1.0000
1.2010         0    0.0217    1.0000
1.4508    0.0217         0    1.0000
1.0000    1.0000    1.0000         0

MI =

-0.8331    5.2214   -4.3883    1.0956
5.2214  -55.8175   50.5960   -6.3666
-4.3883   50.5960  -46.2077    6.2710
1.0956   -6.3666    6.2710   -1.4516

Gs =

3.2134
1.4732
1.2229
1.0000

L =

0.7441
-9.9423
10.1982

Zs =

-2.1255e+003


-
If you are using a valid semivariogram, you should never run into these difficulties. However, for certain data configurations and certain variograms you can get close to singular matrices. (I have never seen kriging software that uses a pseudoinverse--that would be appropriate only when data fall on a regular grid or when all data are used at once; directly solving the system $ML=G_s$ is more stable and far more efficient in most applications.) What semivariogram are you using that causes problems? BTW, did you really mean to put absolute values in the expression for $G_s$ but not in $M$? –  whuber Jun 1 '12 at 15:37
Hi whuber. Thanks for the help! I am using a Gaussian semivariogram. –  Andy Jun 1 '12 at 16:10
Gaussians can produce such strange results, Andy: because they and their derivatives are zero at the origin, they impose extreme smoothness on the kriged estimates. In the bottom picture, fitting a smooth spline (as a proxy for the kriging) clearly would extrapolate high values in the 150-250 range. Such extrapolation requires weights beyond the 0-1 range. The extremely low value you have plotted at 500 similarly looks consistent with the data and variogram. Why not experiment first with an exponential variogram (and add a tiny nugget to it)? –  whuber Jun 1 '12 at 17:36
Thanks. I understand now that the result value is not so crazy after all. Besides I want to do interpolation, not extrapolation so this case is not so important. –  Andy Jun 2 '12 at 10:07
The kriging equations are set up to guarantee that the matrix $M$ is always of full rank, Michael. It is an error to use a semivariogram $g$ that would produce a singular $M$. This issue is likely tangential to understanding OK and the Matlab script. –  whuber Jun 1 '12 at 16:08