next up previous contents
Next: Dispersion Variance Up: Variances and Regularization Previous: Variances and Regularization   Contents

Estimation Error, Estimation Variance

Every estimation method involves an estimation error, arising from the simple fact that the quantity to be estimated generally differs from its estimator $/hat{z}$, thus implying an error of estimation $z - /hat{z}$.

When the mean grade $z_v({/boma x}_i)$ of a vertical bore-hole through the middle of a block $V$ is used to estimate the true mean grade $z_V({/boma x}_i)$ of the block, the error involved is $r({/boma x}_i) = z_v({/boma x}_i) - z_V({/boma x}_i)$. Just as the regionalized variable $z({/boma x})$ has been interpreted as a particular realization of the random function $Z({/boma x})$, so the error $r({/boma x}_i)$ also appears as a particular realization of the random variable $R({/boma x}_i)= Z_v({/boma x}_i) - Z_V({/boma x}_i)$ at the point ${/boma x}_i$.

Suppose now that the entire deposit is divided into blocks of equal size $V$, and that each block is intersected by a vertical bore-hole passing through its center. If the deposit is an homogeneous mineralization, i.e., if the random function $Z({/boma x})$ can be considered as stationary, then the error $R({/boma x})$ is also stationary and any, two errors $r({/boma x_i})$ and $r({/boma x_j})$ can be considered as two different realizations of the same stationary random function $R({/boma x}) = Z_v({/boma x})-Z_V({/boma x})$. If the histogram of experimental errors is available in a control zone, then it will be possible, because of the stationary assumption, to infer the complete distribution function of $R({/boma x})$. Even if such a histogram is not available, it will still be possible to calculate the stationary expectation $m_E = E/{R({/boma x})/}$ and variance $/sigma^2_E = Var/{R({/boma x})/}$ of the distribution function of the error.

The particular error $r({/boma x}_i) = z_v({/boma x}_i) - z_V({/boma x}_i)$ involved when estimating the block $V({/boma x}_i)$ remains unknown, but the mean and variance of the errors (or the complete distribution function if it is known) will provide a measure of the quality of the estimation.

In case of the assumption of second-order stationarity of the random function $R({/boma x})$ we have

(i)
the mathematical expectation,

/begin{displaymath}E/{R({/boma x})/} =m_E = constant, / / /forall {/boma x}/ / ;/end{displaymath}

(ii)
a variance, called the ``estimation variance'',

/begin{displaymath}Var/{R({/boma x})/} = E/{[R({/boma x})]^2/} - m^2_E = /sigma_E^2 =
constant, / / /forall {/boma x}/ / ./end{displaymath}

The expectation $m_E$ characterizes the mean error, while the variance $/sigma^2_E$ is a dispersion index of the error. Thus, a good estimation procedure must be such that it ensures
(i)
a mean error close to zero, this property of the estimator is known as unbiasedness;
(ii)
a dispersion of errors very concentrated around this mean zero value, this second criterion being expressed by a small value of the estimation variance $/sigma_E^2$.

At the estimation stage, the type of the distribution function of the errors is unknown in most cases. But, since the two most important characteristics of this function-its expectation and variance-can be calculated, we shall refer to a standard two-parameter ($m_E$ and $/sigma^2_E$) function which will provide an order of magnitude of the required confidence interval.

Among all the two-parameter distribution functions, the one most often used to characterize an error is the normal distribution. The main justification for using this distribution is that it is the one most often observed in practice, particularly in mining practice.

In mining applications, the error distribution functions are generally symmetric with a slightly more pronounced mode and larger tails than a normal distribution with the same expectation and variance (cf. Figure 5.1). Thus, in relation to the normal distribution, there are more small errors (in the region of $m_E = 0$) and more large errors (in the tails of the distribution). However, the classic confidence interval $[m_{E} /pm 2
/sigma_{E}]$ contains approximately 95% of the observed errors. The 95% normal confidence interval is then a good criterion for judging the quality of a geostatistical estimation, but it should always, as far as possible, be verified experimentally through test or control zones.

Figure 5.1: Distribution of Estimation Errors.
/begin{figure}/begin{center}
/mbox
{/beginpicture
/setcoordinatesystem units <1c...
... /
/endpicture}
/end{center}/index{Verteilung der Sch/uml atzfehler}/end{figure}

Until now we have denoted an estimated value by $/hat{z}(/boma{x}_{i})$. How can an estimator look like, which produces such estimated values of a particular realization. In general, it must be a function of $n$ random variables, the sample variables, namely

/begin{displaymath}/hat{Z}=f(Z(/boma{x}_{1}),/ldots,Z(/boma{x}_{n}))/ / / ./end{displaymath}

The function $f$ has to fulfill certain conditions, e.g.

(i)
unbiasedness:

/begin{displaymath}E(Z_{V}-/hat{Z}) = 0/ / / ;/end{displaymath}

(ii)
it should be reasonably simple, so as to permit the calculation of the estimation variance

/begin{displaymath}/sigma^2_E=E(Z_{V}-/hat{Z})^{2}=E(Z_{V}^{2})+E(/hat{Z}^{2})-2E(Z_{V}/hat{Z}) / / ./end{displaymath}

The simplest proposal for an estimator is a linear function:

/begin{displaymath}/hat{Z}=/sum_{i=1}^{n} /lambda_{i}Z(/boma{x}_{i})/ / / ./end{displaymath}

Consider first the second-order stationarity of $Z({/boma x})$ and the simple case of the estimation of the arithmetic mean $z_K$ of the $K$ unknown values $/{z({/boma x}_k), k=1,/ldots,K/}$:

/begin{displaymath}z_{K}=/frac{1}{K} /sum_{k=1}^{K}z(/boma{x}_{k})/ / ./end{displaymath}

The (simplest) linear estimator $/hat{z}_{K}$ is the arithmetic mean of the $n$ available data values $/{z({/boma x}_i), i=1,/ldots,n/}$:

/begin{displaymath}/hat{z}_{K}=/frac{1}{n} /sum_{i=1}^{n}z(/boma{x}_{i})/ / / ./end{displaymath}

These values $z(/boma{x}_{i})$ are interpreted as realizations of the random variables $Z(/boma{x}_{1}),/ldots,Z(/boma{x}_{n})$ which have the same distribution as $Z$. The unbiasedness condition is fulfilled because

/begin{displaymath}EZ_{K}=/frac{1}{K} /sum_k EZ(/boma{x}_{k})=m /; /; /mbox{and} /; /;
E /hat{Z}_{K}=/frac{1}{n} /sum_i EZ(/boma{x}_{i})=m / / ,/end{displaymath}

which implies

/begin{displaymath}E(Z_{K}-/hat{Z}_{K})=0/ / / ./end{displaymath}

The estimation variance is written as

/begin{displaymath}/sigma_E^2=E(Z_{K}-/hat{Z}_{K})^{2}=EZ_{K}^{2}+E
/hat{Z}_{K}^{2}-2EZ_{K}/hat{Z}_K/end{displaymath}

with

/begin{displaymath}/begin{array}{lcl}
EZ_{K}^{2} & = & /frac{1}{K^2} E/{/sum_k /...
...m_{k'} [C(/boma{x}_{k}-/boma{x}_{k'})+m^{2}]/ / / ,
/end{array}/end{displaymath}

where $C(/boma{x}_{k}-/boma{x}_{k'})$ only depends on the difference $/boma{x}_{k}-/boma{x}_{k'}$. Similar things hold for $/hat{Z}_K$,

/begin{displaymath}E/hat{Z}_{K}^{2}=/frac{1}{n^2} /sum_i /sum_j
[C(/boma{x}_{i}-/boma{x}_{j})+m^{2}]/ / / ,/end{displaymath}

and further,

/begin{displaymath}EZ_{K}/hat{Z}_{K}=/frac{1}{nK} /sum_k /sum_i
[C(/boma{x}_{k}-/boma{x}_{i})+m^{2}]/ / / ./end{displaymath}

$m^{2}$ is eliminated, leaving

/begin{displaymath}/sigma_E^2=E(Z_{k}-/hat{Z}_{k})^{2}/end{displaymath}


/begin{displaymath}=/frac{1}{K^2}/sum_k /sum_{k'} C(/boma{x}_{k}-/boma{x}_{k'}) ...
...2 /frac{1}{nK}
/sum_k /sum_i C(/boma{x}_k - /boma{x}_i)/ / / ./end{displaymath}

Let us denote the 3 averages by $/bar{C}((K),(K)),$ $/bar{C}((n),(n))$ respectively $/bar{C}((K),(n)).$ These are mean values of the covariance $C({/boma h})$ when one extremity of the vector ${/boma h}$ varies through one set and the other end varies through the other one. Therefore we obtain

/begin{displaymath}/sigma_E^2=/bar{C}((K),(K)) + /bar{C}((n),(n)) - 2/bar{C}((K),(n))/ /
/ ./end{displaymath}

Here we have estimated the average of $K$ points. The generalization to the continuous case is straightforward.

$K$ points $/rightarrow$ volume $V(/boma{x})$
$n$ points $/rightarrow$ volume $v(/boma{x}')$,
where $/boma{x}$ resp. $/boma{x}'$ denote the centers of the volumes. Then we get

/begin{displaymath}z_{K} /rightarrow z_{V}(/boma{x})=/frac{1}{V} /int_{V(/boma{x})}
z(/boma{y})d/boma{y}/ / / ./end{displaymath}


/begin{displaymath}/hat{z}_{K} /rightarrow /hat{z}_{v}(/boma{x}')=/frac{1}{v}
/int_{v(/boma{x}')} z(/boma{y})d/boma{y}/ / / ./end{displaymath}

$z_{V}(/boma{x})$ and $/hat{z}_{v}(/boma{x}')$ are again interpreted as realizations of corresponding random variables, and the estimation variance is defined by

/begin{displaymath}/sigma_E^2=E(Z_{V}(/boma{x})-/hat{Z}_{v}(/boma{x}'))^{2}/end{displaymath}


/begin{displaymath}=/frac{1}{V^2} /int_{V (/boma{x})} /int_{V( /boma{x})} C(/bom...
.../int_{v(/boma{x}')} C(/boma{y}-/boma{y'})d /boma{y} d /boma{y'}/end{displaymath}


/begin{displaymath}- /frac{2}{Vv} /int_{V(/boma{x})} /int_{v(/boma{x}')} C(/boma{y}-/boma{y}')
d/boma{y}'d/boma{y}/ / / ./end{displaymath}

Simple expressions for the mean values lead to the formula
/begin{displaymath}
/sigma_E^2=/bar{C}(V,V)+/bar{C}(v,v)-2 /bar{C}(V,v)/ / / ,
/end{displaymath} (5.1)

where $/bar{C}$ stands for averaged covariances between respective volumes.

If the covariance exists, $/gamma$ also exists, and because of

/begin{displaymath}C(/boma{h})=C(/boma{0}) - /gamma (/boma{h})/end{displaymath}

we have
/begin{displaymath}
/sigma_E^2=2 /bar{/gamma}(V,v) - /bar{/gamma}(V,V) - /bar{/gamma}(v,v)/
/ / .
/end{displaymath} (5.2)

$/bar{/gamma}(V,v)$ stands for the average value, when one end of $/boma{h}$ varies in $V$ and one end in $v$.

Remark 1: The estimation variance of $V$ by $v$ is sometimes referred to as the variance of extending the grade of $v$ to $V$ or simply the extension variance of $v$ to $V$ and is then denoted by $/sigma_E^2(v/V).$

Remark 2: The two formulae for $/sigma_E^2$ are completely general whatever the domains $V$ and $v$. In particular, the domains need not necessarily be connex (single continuous regions); the domain $V$ to be estimated may be, for example, two distinct blocks $V = V_1 + V_2$ and the information domain $v$ may consist of several samples $v = /sum_i S_i$, some of which may be situated within the domain $V$.

Remark 3: The variogram $2/gamma({/boma h})$ can itself be interpreted as the elementary estimation variance of a variable $Z({/boma x})$ by another variable $Z({/boma x}+{/boma h})$ at a distance $h$ from $Z({/boma x})$:

/begin{displaymath}/sigma_E^2 = E(Z(/boma{x}+/boma{h})-Z(/boma{x}))^{2} = 2 /gamma(/boma{h})/end{displaymath}

Remark 4: The quality of the estimation of $V$ by $v$ is a function of the following.

(i)
Because $/gamma(/vert /boma{h}/vert)$ in general increases with $/vert /boma{h}/vert$, $/bar{/gamma}(V,V)$ will also increase with the size of $V$, and therefore, $/sigma_E^2$ will decrease. This means that $V$ is easier to be estimated if it is larger. Moreover $/sigma_E^2$ depends on the form of $V$.
(ii)
The distances between $V$ and $v$ enter directly in the computation of $/sigma_E^2$.
(iii)
The form of $v$ is important: As larger $v$ as smaller $/sigma_E^2$. However, even we have the same size of $v$, its form also plays an important role. For example, the configuration

will better estimate $V$ than
units <1.0cm,1.0cm> x from -2.0 to 3.0, y from -1.0 to 3.0 0.0 -0.05 1.6 -0.05 / 0.1 1.65 -1.3 1.65 / 0.9 1.65 2.2 1.65 / corners at -0.8 -0.45 and 0.0 0.35 corners at 0.1 1.2 and 0.9 2.0 $V$ at -0.4 -0.1 $V$ at 0.5 1.6 360 degrees from 1.65 0.4 center at 1.6 0.4 360 degrees from 1.65 -0.55 center at 1.6 -0.5 $v_{1}$ [l] <0.15cm,0.0cm> at 1.65 0.4 $v_{2}$ [l] <0.15cm,0.0cm> at 1.65 -0.55 $d$ at 0.8 0.3 360 degrees from -1.3 1.65 center at -1.35 1.65 360 degrees from 2.3 1.65 center at 2.25 1.65 $v_{1}$ [t] <0.0cm,-0.15cm> at -1.35 1.65 $v_{2}$ [t] <0.0cm,-0.15cm> at 2.25 1.65 $d$ at 1.55 2.1 $d$ at -0.6 2.1

This intuitive notion of the importance of the sample configuration, which is formalized in geostatistics through the term $/bar{/gamma}(v, v)$, is ignored by the more usual estimation methods such as weighting by inverse distances or their squares.

(iv)
The quality of the estimation should depend on the structural characteristics (anisotropies, degree of regularity) of the regionalized phenomenon under study. The semi-variogram $/gamma$ may depend on the direction $a$ of the vector $h$, i.e., the semi-variogram will be a function of distance and direction, $/gamma(/vert{/boma h}/vert, {/boma /alpha})$. In horizontally stratified mineralizations, for example, the variation in grade values will be much more continuous in a horizontal direction than they will be vertically, and, by taking the anisotropic semi-variogram into account, more weight would be given to sample $v_1$ in estimating block $V$ in the same stratum than to sample $v_2$ in a different stratum.



units <.8cm,1.0cm> x from 0.0 to 16.8, y from 0.0 to 2.0 0.0 0.0 5.65 0.0 / 0.0 1.0 5.65 1.0 / 0.0 1.9 5.65 1.9 / 10.95 0.0 16.8 0.0 / 10.95 1.0 16.8 1.0 / 10.95 1.9 16.8 1.9 / corners at 0.25 1.15 and 1.15 1.75 corners at 11.2 1.15 and 12.1 1.75 $V$ at 0.7 1.45 $V$ at 11.65 1.45 360 degrees from 4.4 1.45 center at 4.35 1.45 360 degrees from 15.55 0.5 center at 15.5 0.5 $v$ [l] <0.15cm,0.0cm> at 4.4 1.45 $v$ [l] <0.15cm,0.0cm> at 15.55 0.5 is at 8.3 1.6 different at 8.3 1.0 of at 8.3 0.4

This fundamental geological notion of the spatial continuity of a mineralization as a function of direction is quite often overlooked or simply ignored in the application of many estimation methods such as that of weighting by inverse distances.

Remark 5: The above formulas also cover the particular case of estimating the mean value $Z_V$ of a block $V$ by a linear combination $/hat{Z}$ of $n$ available data values $Z({/boma x}_i)$ taken at the points $/{{/boma
x_i},i=1,/ldots,n/}$:

/begin{displaymath}Z_{V}=/frac{1}{V} /int_{V(/boma{x})} Z(/boma{y})dy, /; /;
/hat{Z}_{V}=/sum_{i=1}^{n} /lambda_{i}Z(/boma{x}_{i})/ / / ./end{displaymath}

The non-bias condition entails that $E/{/bar{Z}_V/} = m$, i.e.

/begin{displaymath}E(/hat{Z}_{V})= /sum_i /lambda_i EZ(/boma{x}_{i})= m /sum
/lambda_{i}=E(Z_{V})=m/ / / ,/end{displaymath}

i.e.

/begin{displaymath}/sum_i /lambda_{i}=1/ / / ./end{displaymath}

The estimation variance can then be calculated by (proof later)

/begin{displaymath}/sigma_E^2 = E/{Z_{V}-/hat{Z}_{V}/}^{2}= 2 /sum_i /lambda_{i}...
...j /lambda_{i} /lambda_j /gamma
(/boma{x}_{i}-/boma{x}_j)/ / / ,/end{displaymath}

where $/bar{/gamma}({/boma x_i}, V)$ denotes the mean value of $/gamma({/boma h})$ when one extremity of the distance vector $h$ is fixed at point $/boma{x}_i$ and the other extremity independently describes volume $V$.

This expresses the estimation variance $/sigma_E^2$ as a linear function of the $n$ weights $/lambda_i$. The estimation procedure known as kriging determines the optimal set of weights $/{/lambda_i,
i = 1, /ldots, n/}$, i.e., the weights $/lambda_i$ which minimize the variance $/sigma_E^2$ subject to the non-bias condition $/sum_i/lambda_i = 1$. Hence, kriging appears as the best linear unbiased estimator (in short a ``BLUE'' estimator).

Example 5.1: (Journel and Huijbregt, 1978[11]): Block kriging at Chuquicamata. Chuquicamata is an open-pit porphyry-copper deposit in Chile where short-term planning is based on data from blast-holes. The mean copper grade of each block $V(20 m /times 20 m /times 13 m)$ is estimated by kriging from the blast-hole data available from the nine neighboring blocks which have already been mined out, as shown on Figure 5.2, i.e., each block $V$ is assigned a weighted average of the blast-hole values, the weights being chosen to ensure the minimum estimation variance. The kriged estimator of block $V$ will be written $/hat{Z}_V$. Next, the block $V$ is put into production and blast-hole samples are taken from it, an average of six being taken from each block. Experimental observation has shown that the arithmetic mean of these six holes, $Z_V$, can be taken as the true grade of block $V$.

Figure 5.2: Short-term Kriging at Chuquicamata.
/begin{figure}/begin{center}
/mbox
{/beginpicture
/setcoordinatesystem units <1c...
...335 2.41668
/put {5} at 10.08335 1.76668
}
/endpicture}
/end{center}/end{figure}

A test zone of homogeneous mineralization which comprised 397 blocks gave the histogram of the observed errors $(Z_{V}-/hat{Z}_{V})$ (Figure 5.3).

Figure 5.3: Histogram of Kriging Errors.
/begin{figure}/begin{center}
/mbox
{/beginpicture
/setcoordinatesystem units <1c...
...0.75 5.35
1.5 3.3
2.65 0.8
3.9 0.1 /
/endpicture}
/end{center}/end{figure}

The arithmetic mean of the 397 errors was $+0.01/% Cu$ and, given that the mean copper grade over the test zone was $2/% Cu$, it is obvious that the non-bias condition has been met.

The experimental variance of the 397 errors was $0.273 (/% Cu)^2$ while the predicted theoretical estimation variance calculated was $0.261 (/% Cu)^2$. Therefore, the estimation variance was predicted with the very acceptable relative precision of 4%: (0.273-0.261)/0.273 = 0.04.

A normal distribution with zero mean and a variance equal to that of the experimentally observed errors ( $/sigma_E^2 = 0.273$) is also shown on Figure 5.3. It fits reasonably well the experimental distribution of errors with the following remarks.

(i)
The normal distribution underestimates the proportion of both small errors (between $/pm 0.5/% Cu$) and large errors (greater than $1.5/% Cu$).
(ii)
The standard normal interval ($/pm 2/sigma_E$) correctly estimates the experimental 95% confidence interval: 96% of all the observed errors fall within $/pm 2/sigma_E = /pm 1.04/% Cu$. If the predicted theoretical variance, $/sigma_E^2 = 0.261$ had been used, the result would have made little difference: $/pm 2/sigma_E = /pm 1.02/% Cu$.


next up previous contents
Next: Dispersion Variance Up: Variances and Regularization Previous: Variances and Regularization   Contents
Rudolf Dutter 2003-03-13