The geostatistician takes another approach, called krigging. Instead of 
solving a global optimization problem the geostatistician
solves a series of small inversion problems.
They find all of the points around the point they
are trying to estimate.  They  extract the covariance
at the distance between known points $C_{ij}$ and
the covariance at the distance between known points and
estimation point $C_{i0}$.
They then solve for the weights   ($w_1 ... w_n$) by solving the set
of equations applied by
\beq 
\left[
\begin{array}{cccc} 
C_{11}  &...& C_{1n} & 1 \\
.  &...& . & . \\
.  &...& . & . \\
.  &...& . & . \\
C_{n1}  &...& C_{nn} &  1 \\
1 &...& 1 & 0 
\end{array}
\right]
\left[
\begin{array}{c}
w_1 \\
. \\
. \\
. \\
w_n \\
\mu 
\end{array}
\right]
=
\left[
\begin{array}{c}
C_{10} \\
. \\
. \\
. \\
C_{n0} \\
1 
\end{array} \label{eq:krig}
\right] .
\eeq
To guarantee that the matrix is invertible geostatisticians approximate
the covariance function
 through a linear combination of a limited set of 
functions.
Each function is described by:
\begin{description} 
\item [{\it range}], the distance at which the covariance function is
essentially gone to 0,
\item [{\it anisotropy}], the amount the covariance function is asymmetric, and
\item [{\it orientation}], the major orientation axis of the covariance function.
\end{description}
These parameters are used to describe either a spherical, exponential,
Gaussian, or power model
that are guaranteed to produce a positive definite covariance matrix.

The left panel of Figure~\ref{fig:geostat-mis} shows the geostatistical
approximation of the model covariance for our missing data problem.
This approximation does a good job characterizing the primary dip
of the covariance function, but
does not  accurately describe the {\it range} of the covariance function.
The right panel of Figure~\ref{fig:geostat-mis} shows the result of applying 
krigging
to the same missing data problem.  The result is fairly comparable to the
geophysics result. Instead of adding a second dip at every location we have
only forced the primary dip. 
Because this approximation adds insufficient {\it range} for the
covariance function our answer returns to the local average
between the two wells with the greates separation.
Overall the result is as unsatisfactory as the geophysics approach to the
problem.
\plot{geostat-mis}{width=6.0in,height=3.0in}
{Left panel shows the approximated covariance matrix, right
panel is the result of solving (\ref{eq:krig}) at every
unknown model point.}
