Inference of spatial indicator covariance parameters by maximum likelihood using MLREML

PII:

Computers & Geosciences Vol. 24, No. 5, pp. 453±464, 1998 # 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain S0098-3004(98)00015-6 0098-3004/98 $19.00 + 0.00

INFERENCE OF SPATIAL INDICATOR COVARIANCE PARAMETERS BY MAXIMUM LIKELIHOOD USING MLREML EULOGIO PARDO-IGUÂZQUIZA* TAMSAT, Department of Meteorology, University of Reading, Reading RG6 6BB, U.K. (Received 12 June 1997; revised 28 December 1997) AbstractÐThis paper discusses an application of a previously published code to show how maximum likelihood estimation may be used for the inference of spatial indicator covariance parameters. Although the maximum likelihood equations are derived on the assumption of a multivariate normal distribution for the experimental data, when the data are neither Gaussian nor can be transformed to normal scores (as is clearly the case with indicator data) maximum likelihood still provides a weighted least squares criterion of ®tting spatial covariance models. The estimator is still consistent, asymptotically unbiased and the estimates are asymptotically normally distributed. Maximum likelihood also provides the uncertainty of the estimates by assessing the standard errors of the estimates. The method is recommended for cases in which the number of experimental data is relatively small (several dozens) and irregularly distributed. In such situations, the experimental indicator variogram, calculated by classical non-parametric methods, is so noisy that it is extremely dicult to ®t it to a model. This diculty increases as the threshold moves further away from the median of the data. A simulated case study is used to demonstrate the application of the method and illustrate the results. The program MLREML (previously published in Computers and Geosciences) has been used to perform maximum likelihood estimation of spatial indicator covariances. # 1998 Elsevier Science Ltd. All rights reserved Code available at http://www.iamg.org/CGEditor/index.htm Key Words: Indicator variogram, Indicator kriging, Cumulative distribution function, Maximum likelihood estimation, Spatial statistics.

INTRODUCTION

Regionalized variables are distributed in space and are usually only known at a number of ®nite experimental points where measurements are taken. The value of the variable may be estimated at the nonsampled locations by, for example, ordinary kriging, where second-order stationarity is assumed for the random function Z(x) that models the regionalized variable. The result is that, at each point, both the expected value or kriging estimate, zÃ(x) and the variance of the estimation error, s^ 2 (x), are obtained. If the entire distribution function (rather than just mean and variance) is of interest, a model must be assumed for the distribution. An alternative, more ecient method is to estimate the distribution directly from the experimental information using the nonparametric method known as indicator kriging (Journel, 1983). With the indicator approach, the (conditional) cumulative distribution function (cdf) is estimated at every spatial location using the surrounding experimental data (thus the estimated cdf is conditional on these neighbourhood data, but this is understood and will not be used in *E-mail: [email protected]. 453

the subsequent notation). Some points to note about the indicator approach are: * The cdf is estimated only at K selected threshold values. The value of the cdf at values other than the thresholds is obtained by interpolation. * The core of the indicator approach is that the cdf appears as the expectation of a binary indicator transformation of the initial random function Z(x): F…zk † ˆ PfZ…x†Rzk g ˆ EfI…x; zk †g ˆ mI with

 I…x; zk † ˆ

1, 0,

if Z…x†Rzk otherwise:

…1†

…2†

The expectation is conditional on the experimental data. * The cdf may be estimated from the indicator data using indicator kriging (IK): ^ k† ˆ F…z

m X li …zk †
I…x i ; zk †,

k ˆ 1, . . . , K

…3†

iˆ1

This implies that K indicator covariances (or variograms) need to be inferred from the indicator data.

454

E. Pardo-IguÂzquiza

Although an estimate close to Equation (3) may be obtained using indicator cokriging (Journel and Alabert, 1989), the procedure is extremely demanding, requiring the modelling of K2 indicator covariances and cross-covariances. In practice, even IK may be considered cumbersome and replaced by median indicator kriging, but this is related to the diculty of inference of indicator variograms which is easily solved by maximum likelihood estimation as presented in this paper. * From the estimated cdf, di€erent estimates may be obtained for the unknown value z(x): point estimates such as the conditional expectation or the conditional median (Journel, 1983) and interval estimates, which are read directly from the distribution function. It is also possible to estimate the probability that z(x) will be greater than a given threshold zk by reading values directly from the cdf [actually as 1ÿFÃ(zk)]. The estimation of the cdf by IK as in Equation (3) at any experimental location x with n experimental data {z(xi); i = 1, . . . , n} requires the following sequence of steps: (i) Selection of the K thresholds {zk; k = 1, . . ., K} that discretize the cdf. (ii) For each threshold zk the indicator data {I(xi; zk); i = 1, . . . , n} are obtained from the original experimental data {z(xi); i = 1, . . ., n} using Equation (2). (iii) Inference of K indicator variograms from the indicator data {I(xi; zk); i = 1, . . ., n}. (iv) Resolution of K kriging systems, which gives the K estimates of the cdf by Equation (3) using the indicator data and the kriging weights. (v) The cdf is obtained from the K estimated values, subject to ensuring that the properties of cumulative distribution functions are veri®ed: FÃ(zk) $ [0, 1] and zkRzk'cFÃ(zk)RFÃ(zk'). (vi) Point estimates, interval estimates and recovery values are estimated from the estimated cdf. (vii) Steps (iv) to (vi) are repeated for every point x where the cdf is to be estimated. This paper focuses on step (iii) of the indicator approach, i.e. the inference of K indicator variograms from the indicator data. This step, known as indicator structural analysis, is considered next.

INDICATOR STRUCTURAL ANALYSIS

Given a second order stationary random function Z(x) and any two random variables Z(x) and Z(x + h) taken from the random function Z(x), the bivariate cumulative distribution function of {Z(x), Z(x + h)} is given by the noncentred indicator cross-covariance (Journel, 1983): KI …h; zk , zk 0 † ˆ EfI…x ‡ h; zk †
I…x; zk 0 †g 2

…4†

Thus the inference of K indicator covariances and cross-covariances amounts to the inference of the

bivariate cumulative distribution function, Equation (4), at discrete points given by the K thresholds {zk; k = 1, . . ., K}. For estimation of the cdf (Eq. (1)) by IK, the inference of K indicator (auto)covariances are sucient, the noncentred indicator covariances being de®ned by: KI …h; zk † ˆ EfI…x; zk †
I…x ‡ h; zk †g ˆ ProbfZ…x†Rzk , Z…x ‡ h†Rzk g:

…5†

The centred covariance is obtained by CI …h; zk † ˆ KI …h; zk † ÿ F 2 …zk †

…6†

and the indicator variogram is given by gI …h; zk † ˆ CI …0; zk † ÿ CI …h; zk † ˆ F…z† ÿ KI …h; zk †: …7† The inference of any one of these functions (noncentred covariance, centred covariance or variogram) amounts to the inference of the other two by Equation (7). In practice, the indicator variogram de®ned by gI …h; zk † ˆ 12 Ef‰I…x ‡ h; zk † ÿ I…x; zk †Š2 g

…8†

is estimated from the experimental indicator data using the method of moments through the wellknown equation: g^ I …h; zk † ˆ

1 X ‰I…x i ‡ h; zk † ÿ I…x i ; zk †Š2 …9† 2N…h† iˆ1 N…h†

where N(h) is the number of pairs of experimental points separated by a distance h. The procedure consists of the estimation of g^ (h; zk) for di€erent values of the distance h and the ®tting of a theoretical model to the estimated points, which can be done by eye or by some automatic procedure. Several problems appear with this approach: * In the frequent situation of a small number of irregularly distributed data, an angle tolerance and a lag tolerance need to be speci®ed in order to obtain, for a given distance h, a sucient number of couples N(h) in Equation (9). These approximations are required to ensure sucient pairs for g^ (h; zk) to be statistically reliable. The selection of the lag tolerance has a signi®cant e€ect on the smoothness of the experimental variogram, so that di€erent models may be considered for di€erent lag tolerances. * A second problem appears speci®cally in the estimation of indicator variograms. As the threshold zk moves further away from the median value of the experimental data, the binary indicator data become less balanced; i.e. for the median the same number of 1s and 0s appear in the experimental indicator data. If the threshold decreases from the median the number of 1s decreases and the reverse occurs if the threshold increases from the

Inference of spatial indicator parameters

median value. This means that as the threshold moves away from the median, the range of the variogram will decrease because it measures the statistically averaged length between boundaries of 1s and

455

0s. This may be seen as a destructuring e€ect (asymptotic tendency toward a pure nugget e€ect) of the covariance of the indicator variable. The dif®culty is that the appearance of the experimental

Figure 1. (A) Histogram of complete (2500) data set and (B) grey-scale map of complete data set, where all values equal or bigger than 10.0 (98.8 percentile) are represented as black.

456

E. Pardo-IguÂzquiza Table 1. Basic statistics of complete (2500) data set

Number of data Mean Variance Standard deviation Coecient of variation Range Skewness Kurtosis Minimum 10 percentile 25 percentile Median 75 percentile 90 percentile Maximum

2500 1.712 4.977 2.231 130.32% 31.446 4.506 36.023 0.025 0.241 0.493 1.053 2.033 3.698 31.472

variogram is so noisy that it is dicult to perceive any structure from the raw variogram, resulting in signi®cant uncertainty in the ®tted model. MAXIMUM LIKELIHOOD ESTIMATION

The method of maximum likelihood estimation has been widely studied in geostatistics (Kitanidis, 1983, 1987; Mardia and Marshall, 1984; Hoeksema and Kitanidis, 1985; Kitanidis and Lane, 1985;

Mardia and Watkins, 1989; Zimmerman, 1989; Dietrich and Osborne, 1991 among others). In all these cases maximum likelihood has been applied to continuous variables that are assumed to be a realization of a multivariate normal random function. Usually if the histogram of the data does not resemble a Gaussian model, a transformation such as normal scoring may be applied so that the marginal univariate distribution (estimated by the histogram) looks more like a Gaussian model. Of course this transformation does not ensure that the transformed data follow a multivariate Gaussian distribution, but in practice with only one realization of the random function it is not possible to estimate the multivariate distribution of the data. In this situation the choice of the multivariate Gaussian distribution is the more convenient for the reasons explained in Pardo-IguÂzquiza (1998). In the previous reference it is also shown how maximum likelihood may be seen as a weighted least squares procedure which still maintains the appropriate statistical properties of maximum likelihood, irrespective of the true multivariate distribution of the data (in practice unknown). In the present paper it is shown how the method may be applied to the inference of indicator covariance parameters where

Fig. 2(A±D). Continued.

Inference of spatial indicator parameters

457

Figure 2. ÐÐÐ: ``true'' indicator variogram using complete data set; - - -: model ®tted to experimental indicator variogram calculated using complete data set for thresholds (A) 10 percentile, (B) 25 percentile, (C) 50 percentile, (D) 75 percentile, (E) 90 percentile, (F) 95.7 percentile and (G) 98.8 percentile

the data are 0s and 1s and do not follow obviously a Gaussian distribution. A compilation of the equations used in maximum likelihood estimation may be found in Pardo-IguÂzquiza (1997), where there is a description of the program MLREML which performs maximum likelihood and restricted maximum likelihood estimation and evaluates the uncertainty of the estimated parameters. SAMPLING STATISTICS FROM A SIMULATION CASE STUDY

A realization of a regionalized random function on a 50  50 square grid of side length equal to

unity is perfectly known. The histogram of the 2500 data is shown in Figure 1(A); Figure 1(B) represents the grey scale map of the realization [where all the values equal or bigger than 10.0 (the 98.8 percentile) are represented in black]. The basic statistics are shown in Table 1. The indicator variograms calculated from the 2500 data set will be considered as the ``true'' indicator variograms. Seven thresholds have been selected equal to the percentiles 10, 25, 50, 75, 95, 95.7 and 98.8 of the complete data set. The z-thresholds are 0.241, 0.493, 1.053, 3.698, 6.0 and 10.0, respectively. The indicator variograms calculated with the complete data set and the previous thresholds are shown in

Table 2. Models ®tted to experimental indicator variograms calculated with complete (2500) data set Threshold 0.241 (percentile 10.0) 0.493 (percentile 25.0) 1.053 (percentile 50.0) 2.033 (percentile 75.0) 3.698 (percentile 90.0) 6.000 (percentile 95.7) 10.000 (percentile 98.8)

Variogram model

Variance

0.02 + 0.07Sph(7.0) 0.035 + 0.155Sph(9.0) 0.06 + 0.19Sph(10.0) 0.05 + 0.15Sph(10.0) 0.02 + 0.08Sph(8.0) 0.012 + 0.033Sph(6.0) 0.0045 + 0.008Sph(5.0)

0.09 0.19 0.25 0.20 0.10 0.045 0.0125

Where the variogram model: A + BSph(C) represents a spherical model with sill B, range C and nugget e€ect A.

458

E. Pardo-IguÂzquiza

Figure 2(A)±(G), respectively. In Figure 2(A)±(G) the ®tted models have been represented as a dashed line. These models have been summarised in Table 2. The ®rst experiment in inference consists in obtaining the sampling statistics for the maximum likelihood estimator, using di€erent sample sizes and with random location of the data. The question is, what estimates are obtained if we only know the value of the variable for a sample of size n (where n has been chosen equal to 15, 30 and 60 in this experiment) with random location? Given the grid of 2500 data, there are …2500 15 † di€erent ways of selecting 15 data. Some of the con®gurations will give better estimates than others due both to the geometry of the con®guration and the actual values of the variable at those locations. Because an exhaustive analysis of all the con®gurations is impractical, m con®gurations at random have been selected (in this experiment m is a set equal to 100). For each sample the estimates of range, nugget e€ect and variance are obtained by maximum likelihood, giving 100 estimates of each parameter. The histogram of the estimates is an approximation of the sampling distribution of the estimator. Some statistics are calculated to summarise the sampling distribution: mean, bias, variance and mean square error (mse). The mean and variance are estimated by: m 1X  y^ i yˆ m iˆ1

^ ˆ var…y†

…10†

m 1X …y^ i ÿ  y†2 , m iˆ1

respectively.

…11†

The bias b is the di€erence between the mathematical expectation of the sampling distribution and the true value of the parameter: ^ ˆ Efy^ i g ÿ y ˆ  b…y† yÿy

…12†

where y is the true value of the parameter. The mse is the dispersion of the estimates with respect to the true value of the parameter: ^ ˆ mse…y†

m 1X …y^ i ÿ y†2 m iˆ1

…13†

that can be expressed in terms of bias and variance as ^ ‡ var…y†: ^ ^ ˆ b2 …y† mse…y†

…14†

In Equations (10)±(14), y represents a single parameter to be estimated. Mean, variance, bias and mean square error for the range, variance and nugget e€ect parameters have been calculated for the di€erent thresholds and the results are shown in Tables 3±5 for sample sizes of 15, 30 and 60 data, respectively. As the sample size increases, bias decreases (estimator asymptotically unbiased), variance decreases (estimator consistent) and the sampling histogram tends towards a Gaussian bell shape (estimates asymptotically normally distributed). In Tables 3±5, the last column with the label ``data'' refers to the number of the 100 random locations that gives estimates (i.e. the indicator data are not all 1s or all 0s, in which case the indicator variance is 0.0). In the same column and for the range parameter, the number of times that the estimated range reaches an upper limit of 19.1 appears in brackets, i.e. all the estimates of the range larger than this number are assigned to this class. With this we avoid allowing a small number of large estimates to in¯uence the sampling statistics. The important fact is the

Table 3. Sampling statistics with sample size of 15 Threshold

Parameter

Mean

Variance

Bias

mse

Data

0.241 0.493 1.053 2.033 3.698 6.0 10.0 0.241 0.493 1.053 2.033 3.698 6.0 10.0 0.241 0.493 1.053 2.033 3.698 6.0 10.0

range range range range range range range variance variance variance variance variance variance variance nugget e. nugget e. nugget e. nugget e. nugget e. nugget e. nugget e.

9.64 9.88 9.99 10.76 8.81 8.90 8.05 0.1016 0.1689 0.2265 0.1576 0.0915 0.0742 0.0632 0.0047 0.0147 0.0267 0.0162 0.0156 0.0013 0.0000

41.59 41.56 48.86 38.27 35.47 35.02 25.57 0.0019 0.0021 0.0006 0.0026 0.0012 0.0004 0.0000 0.0005 0.0018 0.0033 0.0017 0.0001 0.0001 0.0000

2.64 0.88 ÿ0.01 ÿ0.76 0.81 2.90 3.05 0.0116 ÿ0.0210 ÿ0.0234 ÿ0.0423 ÿ0.0084 0.0292 0.0507 ÿ0.0153 ÿ0.0203 ÿ0.0332 ÿ0.0337 ÿ0.0184 ÿ0.0106 ÿ0.0045

48.57 42.33 48.86 38.84 36.13 43.43 34.90 0.0019 0.0026 0.0011 0.0043 0.0013 0.0013 0.0025 0.0008 0.0022 0.0043 0.0029 0.0004 0.0002 0.00002

74 (11) 99 (19) 100 (23) 99 (19) 77 (10) 50 (5) 22 (1) 74 99 100 99 77 50 22 74 99 100 99 77 50 22

Inference of spatial indicator parameters

459

Table 4. Sampling statistics with sample size of 30 Threshold

Parameter

Mean

Variance

Bias

mse

Data

0.241 0.493 1.053 2.033 3.698 6.0 10.0 0.241 0.493 1.053 2.033 3.698 6.0 10.0 0.241 0.493 1.053 2.033 3.698 6.0 10.0

range range range range range range range variance variance variance variance variance variance variance nugget e. nugget e. nugget e. nugget e. nugget e. nugget e. nugget e.

8.90 10.24 11.23 8.64 9.19 8.14 6.30 0.0901 0.1755 0.2367 0.1758 0.0874 0.0531 0.0362 0.0075 0.0249 0.0444 0.0198 0.0118 0.0058 0.0000

33.72 22.68 24.86 24.65 21.24 25.78 21.33 0.0012 0.0015 0.0003 0.0015 0.0019 0.0006 0.0001 0.0005 0.0022 0.0043 0.0017 0.0007 0.0004 0.0000

1.90 1.24 1.23 ÿ1.36 1.19 2.14 1.30 0.0001 ÿ0.0144 ÿ0.0132 ÿ0.0241 ÿ0.0125 0.0081 0.0235 ÿ0.0124 ÿ0.0100 ÿ0.0155 ÿ0.0301 ÿ0.0815 ÿ0.0061 ÿ0.0045

37.33 24.22 26.37 26.50 22.67 30.37 23.03 0.0013 0.0017 0.0005 0.0017 0.0016 0.0007 0.0006 0.0007 0.0023 0.0046 0.0026 0.0008 0.0004 0.00002

96 (9) 100 (11) 100 (12) 100 (4) 93 (2) 69 (5) 25 (1) 96 100 100 100 93 69 25 96 100 100 100 93 69 25

way in which this modal class decreases as the sample size increases (consistency of the estimator). In Figure 3(A)±(C) there is a representation of the sampling histogram of the range for threshold 6.0 (95.7 percentile) for sample sizes of 15 (Fig. 3A), 30 (Fig. 3B) and 60 (Fig. 3C). As the sample size increases, the dispersion of the estimates decreases and are concentrated around the true value equal to 6.0. Taking the mean values from the sampling distribution with sample size 60 for range, variance and nugget e€ect as the parameters of a model, it is possible to compare these models with the experimental indicator variograms obtained from the complete (2500) data set; they have been represented in Figure 4(A)±(G) for the di€erent

thresholds. Apart from the last threshold (98.8 percentile) there is a good agreement between both variograms for the di€erent thresholds. The di€erences between the ``true'' variograms and the mean variograms from the sampling distribution gives the bias in the di€erent parameters that will decrease as the sample size increases. In a second experiment a sample of 60 data with random location is selected (in fact this is just the ®rst of the 100 random selection used in the previous experiment). The basic statistics of this sample data set are shown in Table 6. The data ®le is given in Table 7, so that the results may be replicated using MLREML (Pardo-IguÂzquiza, 1997). The estimated parameters for the di€erent thresholds are given in Table 8. For the same data

Table 5. Sampling statistics with sample size of 60 Threshold

Parameter

Mean

Variance

Bias

mse

Data

0.241 0.493 1.053 2.033 3.698 6.0 10.0 0.241 0.493 1.053 2.033 3.698 6.0 10.0 0.241 0.493 1.053 2.033 3.698 6.0 10.0

range range range range range range range variance variance variance variance variance variance variance nugget e. nugget e. nugget e. nugget e. nugget e. nugget e. nugget e.

9.34 9.35 10.39 9.01 8.59 6.65 5.07 0.0884 0.1808 0.2448 0.1799 0.0838 0.0411 0.0188 0.0215 0.0292 0.0467 0.0348 0.0162 0.0073 0.0000

27.19 15.04 15.32 11.46 16.64 16.13 8.45 0.0009 0.0009 0.0002 0.0008 0.0009 0.0004 0.00003 0.0009 0.0015 0.0026 0.0021 0.0285 0.0003 0.0000

2.34 0.35 0.39 ÿ0.99 0.59 0.65 0.7 ÿ0.0015 ÿ0.0091 ÿ0.0051 ÿ0.0200 ÿ0.0161 ÿ0.0038 0.0063 0.0015 ÿ0.0058 ÿ0.0123 ÿ0.0151 ÿ0.0037 ÿ0.0045 ÿ0.0045

32.67 15.16 15.47 12.44 16.98 16.57 8.46 0.0009 0.0010 0.0003 0.0012 0.0011 0.0004 0.00007 0.0009 0.0016 0.0027 0.0021 0.0008 0.0004 0.00002

100 (10) 100 (6) 100 (3) 100 (1) 100 (4) 90 (3) 58 (0) 100 100 100 100 100 90 58 100 100 100 100 100 90 58

460

E. Pardo-IguÂzquiza

Figure 3. Sampling histogram for range parameter when threshold is 6.0 (95.7 percentile) and with sample sizes of (A) 15, (B) 30 and (C) 60.

Inference of spatial indicator parameters

Figure 4. ÐÐÐ: ``true'' indicator variogram using complete data set; - - -: models obtained from mean of the sampling distribution for range, variance and nugget e€ect and with sample size equal to 60.

461

462

E. Pardo-IguÂzquiza Table 6. Basic statistics of sample data set

Number of data Mean Variance Standard deviation Coecient of variation Range Skewness Kurtosis Minimum 10 percentile 25 percentile Median 75 percentile 90 percentile Maximum

60 1.769 3.435 1.853 104.79% 9.559 2.217 8.718 0.097 0.197 0.570 1.082 2.249 3.742 9.657

set, the indicator variograms considering the di€erent thresholds have been calculated by the classical non-parametric approach using Equation (9). Figure 5(A)±(F) represent the experimental variograms and the models estimated by maximum likelihood. [The variogram for the 10.0 (98.8 percentile) threshold is not given because for this threshold all the sample indicator data are equal to 1, i.e. this particular data set is not informative of the indicator covariance for that threshold]. For most of the experimental variograms the model ®tted by eye or any automatic procedure (in the sense of curve ®tting) would provide di€erent models for those obtained by maximum likelihood. Some of the experimental variograms are particularly noisy and is not easy to guess a model. But more important than that is the fact that maximum likelihood assesses the uncertainty of the estimates that can be used to obtain interval estimates if a probability distribution function is assumed for the error of estimation (for example Gaussian). The standard errors and the estimates intervals (68% con®dence intervals assuming a Gaussian distribution for the estimation error) are given in Table 9 as obtained using MLREML. The relatively inferior results for the nugget e€ect parameter are easily explained. As can be seen in Figure 2(A)±(G), models have been ®tted by eye to the experimental variograms obtained from the complete (2500) data set. They have been considered as the ``true'' values of the parameters. Although the ranges and sills of the experimental variograms seem relatively clear, the nugget e€ect is not so clear and could have been modelled with lower values as is suggested by the sampling statistics obtained by maximum likelihood.

Table 7. Sample data set. First, second and third columns are X coordinate, Y coordinate and value of variable. The di€erent sets of indicator data are obtained from this data set by applying corresponding threshold 4.000000 50.000000 27.000000 25.000000 35.000000 34.000000 2.000000 8.000000 6.000000 42.000000 30.000000 26.000000 12.000000 4.000000 6.000000 19.000000 16.000000 50.000000 24.000000 42.000000 12.000000 20.000000 9.000000 28.000000 18.000000 48.000000 41.000000 12.000000 38.000000 32.000000 27.000000 22.000000 41.000000 21.000000 13.000000 3.000000 28.000000 36.000000 15.000000 49.000000 9.000000 36.000000 36.000000 31.000000 19.000000 21.000000 40.000000 22.000000 43.000000 11.000000 20.000000 5.000000 49.000000 1.000000 32.000000 46.000000 23.000000 29.000000 36.000000 15.000000

37.000000 32.000000 29.000000 24.000000 8.000000 7.000000 34.000000 30.000000 10.000000 43.000000 36.000000 26.000000 42.000000 27.000000 20.000000 30.000000 44.000000 44.000000 14.000000 29.000000 2.000000 47.000000 47.000000 22.000000 1.000000 14.000000 21.000000 27.000000 34.000000 50.000000 34.000000 10.000000 16.000000 46.000000 13.000000 13.000000 9.000000 30.000000 23.000000 8.000000 38.000000 5.000000 35.000000 43.000000 8.000000 36.000000 2.000000 8.000000 8.000000 37.000000 18.000000 42.000000 33.000000 29.000000 42.000000 11.000000 48.000000 45.000000 46.000000 21.000000

9.554886E ÿ 01 2.851474E ÿ 01 4.308289 3.370415 1.504360 1.166946 2.424769 1.049588 1.078215 8.238431 1.851702E ÿ 01 3.443186 9.450347E ÿ 01 4.372724E ÿ 01 5.701655E ÿ 01 2.188935 1.851650 2.308521E ÿ 01 6.274740E ÿ 01 4.011951E ÿ 01 1.181852 2.715797 1.014867 1.065530 2.145306 9.748030E ÿ 02 2.249006 9.494542E ÿ 01 1.441131E ÿ 01 9.288399E ÿ 01 1.924591 2.505698 1.989438 4.435335 1.839142 1.167101 5.447964 4.429599E ÿ 01 1.086624 6.922982E ÿ 01 2.015690 7.094427E ÿ 01 1.523585E ÿ 01 1.970460E ÿ 01 9.656535 2.524806 2.698479 4.041601 6.055389E ÿ 01 5.207313 1.930894E ÿ 01 8.856252E ÿ 01 2.584733E ÿ 01 1.961285 2.488525E ÿ 01 2.386725E ÿ 01 1.241654 6.613083E ÿ 01 2.404311 9.302798E ÿ 01

CONCLUSIONS

Where the number of experimental data is relatively small (a few dozen) and the data are irregularly distributed, the estimation of the indicator variograms by the classical method of moments is dicult because of the noise present in the exper-

imental raw variogram, especially for thresholds di€erent from the median of the data. The ®tting of a theoretical model to these variograms becomes uncertain. The maximum likelihood estimation method for the inference of indicator variograms

Inference of spatial indicator parameters

463

Table 8. Parameters estimated by maximum likelihood with sample (60) data set Threshold 0.241 (percentile 10.0) 0.493 (percentile 25.0) 1.053 (percentile 50.0) 2.033 (percentile 75.0) 3.698 (percentile 90.0) 6.000 (percentile 95.7) 10.000 (percentile 98.8)

Variogram model

Variance

0.054 + 0.060Sph(11.0) 0.0 + 0.172Sph(12.44) 0.032 + 0.229Sph(12.69) 0.075 + 0.136Sph(13.20) 0.078 + 0.026Sph(8.49) 0.0 + 0.033Sph(4.44) n.a.

0.114 0.172 0.261 0.211 0.104 0.033 n.a.

n.a.: not applicable (with the threshold 10.0, all the 60 indicator data are equal to 1).

Figure 5. Experimental variogram and models estimated by maximum likelihood from sample data set for thresholds (A) 10 percentile, (B) 25 percentile, (C) 50 percentile, (D) 75 percentile, (E) 90 percentile and (F) 95.7 percentile.

464

E. Pardo-IguÂzquiza

Table 9. Parameters estimated by maximum likelihood with sample (60) data set, interval estimates with 68% con®dence level Threshold 0.241 0.493 1.053 2.033 3.698 6.0 10.0 0.241 0.493 1.053 2.033 3.698 6.0 10.0 0.241 0.493 1.053 2.033 3.698 6.0 10.0

Parameter (and true value)

e

s.e.

[e ÿ s.e., e + s.e.]

range (7.0) range (9.0) range (10.0) range (10.0) range (8.0) range (6.0) range (5.0) variance (0.09) variance (0.19) variance (0.25) variance (0.20) variance (0.10) variance (0.045) variance (0.125) nugget e€ect (0.02) nugget e€ect (0.0350 nugget e€ect (0.06) nugget e€ect (0.05) nugget e€ect (0.02) nugget e€ect (0.012) nugget e€ect (0.0045)

11.00 12.44 12.69 13.20 8.49 4.44 n.a. 0.114 0.172 0.261 0.211 0.104 0.033 n.a. 0.054 0.0 0.032 0.075 0.078 0.0 n.a.

4.60 2.14 2.72 4.13 8.48 2.07 n.a. 0.0208 0.0313 0.0475 0.0384 0.0189 0.0059 n.a. 0.0139 0.0084 0.0176 0.0202 0.0189 0.0054 n.a.

[6.4, 15.6] [10.3, 13.58] [9.97, 15.41] [9.07, 17.33] [0.01, 16.97] [2.37, 6.51] n.a [0.093, 0.135] [0.141, 0.203] [0.214, 0.309] [0.173, 0.294] [0.085, 0.123] [0.012, 0.039] n.a. [0.040, 0.068] [0.0, 0.008] [0.014, 0.050] [0.055, 0.095] [0.059, 0.097] [0.0, 0.005] n.a.

e: estimate; s.e.: standard error; n.a. not applicable.

o€ers an alternative that gives an estimation of the variogram parameters as well as an estimation of their uncertainty. This case where we are using the maximum likelihood equations derived for multivariate Gaussian data, with indicator data, may be seen as a weighted least square procedure that still maintains good statistical properties of the maximum likelihood method. The program MLREML may be used for an ecient inference of the covariance indicator parameters by maximum likelihood. AcknowledgmentsÐThe work of the author has been supported by a European Community Training Fellowship in the Department of Meteorology (Reading University, U.K.) as part of the Training and Mobility of Researchers (TMR) program. The anonymous reviewers are thanked for their constructive comments. REFERENCES Dietrich, C. R. and Osborne, M. R. (1991) Estimation of the covariance parameters in kriging via the restricted maximum likelihood. Mathematical Geology 23(1), 119±135. Hoeksema, R. J. and Kitanidis, P. K. (1985) Analysis of the spatial structure of properties of selected aquifers. Water Resources Research 21(4), 563±572.

Journel, A. G. (1983) Non-parametric estimation of spatial distributions. Mathematical Geology 15(3), 445±468. Journel, A. G. and Alabert, F. (1989) Non-Gaussian data expansion in the Earth sciences. Terra Nova 1 , 123± 134. Kitanidis, P. K. (1983) Statistical estimation of polynomial generalized covariance functions and hydrologic applications. Water Resources Research 19(2), 909±921. Kitanidis, P. K. (1987) Parametric estimation of covariances of regionalized variables. Water Resources Bulletin 23(4), 671±680. Kitanidis, P. K. and Lane, R. W. (1985) Maximum likelihood parameter estimation of hydrologic spatial process by the Gauss±Newton method. Journal of Hydrology 73, 53±71. Mardia, K. V. and Marshall, R. J. (1984) Maximum likelihood estimation of models for residual covariance in spatial regression. Biometrika 71(1), 135±146. Mardia, K. V. and Watkins, A. J. (1989) On multimodality of the likelihood in the spatial linear model. Biometrika 76(2), 289±295. Pardo-IguÂzquiza, E. (1997) MLREML: a computer program for the inference of spatial covariance parameters by maximum likelihood and restricted maximum likelihood. Computers & Geosciences 23(2), 153±162. Pardo-IguÂzquiza, E. (1998) Maximum likelihood estimation of spatial covariance parameters. Mathematical Geology 30(1), 95±108. Zimmerman, D. L. (1989) Computationally ecient restricted maximum likelihood estimation of generalized covariance functions. Mathematical Geology 21(7), 655±672.