Validity range of the discrete Gaussian change-of-support model and its variant

 

 

J.-P. Chilès

Center of Geoscience and Geoengineering, MINES ParisTech, France

 

 

---

SYNOPSIS

The discrete Gaussian model is a very popular change-of-support model for estimating block grade distributions. It is designed for a stationary random function Z(x) that is regarded as the transform of a stationary random function Y(x) with Gaussian marginal. It is based on an assumption that is not fully true and thus constitutes an approximation to the exact solution. We examine the effect of this approximation on the modelling of the marginal distribution of block values of a lognormal random function. The initial discrete Gaussian model and a variant of that model are considered.

Keywords: discrete Gaussian model, DGM, change-of-support, grade distribution.

---

 

 

Introduction

In many situations decision-makers are faced with selectivity: extract ore above some cutoff grade, remediate soils whose pollution exceeds some threshold, restrict traffic speed when the ozone concentration exceeds some limiting value. In each case, data relates to a small support (e.g., a core, an air pollution sensor) that we regard as a point, whereas the decision is taken for a much larger support - selective mining units, decontamination units, average ozone concentration across a city in a given time interval (typically with a low threshold applied to the daily average concentration and a higher threshold for an hourly average concentration). To predict the effect of selectivity it is necessary to take into account that the grade or concentration distribution becomes less dispersed as the support becomes larger. In the framework of random function theory, the result depends on the whole spatial distribution. It can be obtained by means of Monte Carlo simulations (non-conditional simulations in the global case, conditional simulations in the local case). Even if Monte Carlo simulations are more accessible now than in the past, there are always situations where they require excessively intensive computation. It is therefore useful to have access to the approximate solutions provided by change-of-support models. Many such models have been developed by Georges Matheron (a synthesis is given in chapter 6 of Chiles and Delfiner, 2012). The most popular model is the discrete Gaussian model. It should not be used in just any situation because it has been developed for random functions such as transforms of stationary Gaussian random functions. The initial model was proposed in 1976 by Matheron (1976a). A simpler variant was proposed more recently by Emery (2007). We investigate here the accuracy of this change-of-support model for the modelling of the marginal distribution of block values. This investigation is carried out in the special case of lognormal random functions, which constitute a large and important class of random functions, studied notably by Krige (1978).

We first recall the assumptions of the discrete Gaussian model (DGM) and its variant, and explain why they lead to approximations to the true solution. We then describe the principle of the validation method. Finally, we explore the validity range of both models depending on space dimension, variance and covariance function of the logarithmic variable, and support size. In the conclusion we give some indications on the modelling of the local block distribution.

 

The discrete Gaussian model

Let us consider a stationary random function (SRF) Z(x) that can be expressed as the transform of an SRF Y(x) with standard normal marginal distribution. It is therefore of the form Z(x) = φ(Y(x)) with the transformation function φ = F-¹°G, where F is the marginal cumulative distribution function (c.d.f.) of Ζ(·) and G the standard normal c.d.f. Similarly, we can consider that the mean grade Z(v) of the block (v) is of the form Z(v) =

φ_ν(Υ_ν) where Y_v is a standard normal random variable and φ_v, the block transformation function that we want to determine (because we are here interested only in determining φ_ν, we consider a single block v but assume strict stationarity of Y to ensure that φ_ν does not depend on the block location).

Model DGM1

Consider a uniform random point x within the block v. The random variable Z(x) has for c.d.f. the marginal distribution F of the SRF Z(•) and can be expressed as the transform φ(Υ(χ)) of the random variable Y(x). The basic assumption of the discrete Gaussian model proposed by Matheron (1976a), hereafter referenced as model DGM1, is that the bivariate distribution of the (Y(x), Y_v) pair is Gaussian, with a positive correlation coefficient r. Matheron deduces from Cartier's relationship that the block transformation function φ_ν is given by

This defines the distribution of Z( v).

In practice φ is expressed through its Hermite polynomial expansion

where the χ_η are the normalized Hermite polynomials (see, e.g. Chiles and Delfiner, 2012, Appendix A.5), and the coefficients φ_n are given by

Relationship [1] then implies that φ_ν can be expressed in the form

The variance of Z(v) can be derived from the coefficients φ_nrη of the expansion of φ_ν or from the covariance C(h) of the SRF Z. For consistency, r is obtained by equating these two expressions, that is

Considered as an equation in r, Equation [5] has a unique solution between 0 and 1. The correlation coefficient r is called the change-of-support coefficient.

Extensions of the model (not considered here) enable the local estimation of a block by disjunctive kriging or in a multivariate Gaussian framework. A step further makes it possible to estimate the ore tonnage and metal quantity that will be obtained when selecting blocks on the basis of future grade estimates (once blast-holes become available). This is the so-called information effect. On these points, see Matheron (1976a) and Chiles and Delfiner (2012, pp. 455-466).

Model DGM2

The variant DGM2 proposed by Emery (2007) is simpler but requires the additional assumption that the bivariate distri- bution of Y(x) and Y(x') for two independent random points within the same block v is Gaussian. In that case, Emery shows that r² is the variance of the average Y(v) of Y(•) in the block v:

where ρ(h) is the covariance (here a correlogram) of the SRF Y(•). Moreover, Y_v is simply the average Y(v) rescaled to a unit variance by the change-of-support coefficient r:

This induces large simplifications in the extensions of the model to local estimation (with and without information effect), notably in the framework of a multivariate Gaussian assumption (Emery, 2005, 2008; Chiles and Delfiner, 2012, pp. 455-466).

Discussion of the assumptions

Do there exist random functions satisfying the assumptions of models DGM1 and/or DGM2? It is easy to simulate Gaussian samples with pair-wise correlation r, thus corresponding to samples of Y with independent random locations in the block v (see Emery, 2007). However, the author is not aware of any random function model for fixed locations leading to such correlations for random locations in v. Should such a model exist, it would be specific to that support v. The above assumptions should therefore be considered as approximate only.

Let us first consider the situation where Z is a Gaussian SRF, that is, φ amounts to an affine transformation (in that case we have no need of the discrete Gaussian model, but it is interesting to see what it would mean). Models DGM1 and DGM2 yield the same value for r because C(h) is proportional to ρ(h), and relationship [7] is exact. Since x is random in v, the bivariate distribution of Y(x) and Y_v is a mixture of standard bivariate normal distributions with correlations

for x having any possible location in v. The average value of ω(x) when x scans v is r, so that the average correlation between Y(x) and Y_v is r but, depending on r, it can be an average of very different values. A mixture of such bivariate Gaussian distributions is not a bivariate Gaussian distribution: it is a Hermitian bivariate distribution (Matheron, 1976b; see also Chiles and Delfiner, 2012, pp. 418 and 423). In the one-dimensional case, that is, when v is a segment of length L, and for an exponential covariance with unit sill and with scale parameter a, the change-of-support coefficient r and the function ω(x) are respectively given by

ω(x) is minimal at the extremes (x = 0 or L) and maximum for x = L / 2. The contrast between the largest correlation and the smallest one is

This ratio is equal to 1 for L/ a = 0, remains close to 1 when L / a remains moderate (with a value of about 1.245 for L / a = 1), and increases to 2 when L / a increases to infinity. Figure 1 shows the graph of ω (x) for some values of L / a. The approximation of a constant correlation ω (x) thus seems to be acceptable when L / a is not too large.

 

 

The situation is more critical for the bivariate distribution of Y(x) and Y(x') because it is a mixture of standard bivariate Gaussian distributions with correlations ρ(x' - x), thus ranging from 1 when x = x' to ρ(L) when x and x' are at opposite corners of the block v and L denotes their separation vector (we assume here isotropy and monotonicity of the correlogram). In the above one-dimensional case, L is the length of the segment and this minimum correlation is equal to exp(-L / a): it can be very small if L is large with respect to a.

When Z is not Gaussian, the transformation function φ is not linear. The variogram of Z is therefore not proportional to the variogram of Y. Note that the change-of-support coefficient of model DGM2 depends only on the variogram of Y and does not depend on the transformation φ. This is not the case for model DGM1, because the covariance C of Z depends also on the transformation φ:

The correlograms ρ(h)ⁿ are less and less structured when n increases. For example, if ρ(h) is an exponential covariance with scale parameter a, ρ(h)ⁿ remains exponential, but with scale parameter a / n, and therefore tends to a pure nugget effect when n tends to infinity. As a consequence, the change-of-support coefficient r of model DGM1 is larger than that of model DGM2, with a difference that increases as the block v becomes larger and the terms φ_η with large n dominate in the development of the transformation function φ. The efficiency of model DGM2 will be similar to that of model DGM1 when the block size is small and the φ_n decrease rapidly with n.

Note that the change-of-support coefficient r cannot simultaneously satisfy relationships [5] and [6] (except if Z is bivariate Gaussian). Unlike model DGM1, model DGM2 therefore does not preserve the variance of Z(v) expressed in the right-hand side of relationship [5], which is an important parameter of the block distribution. This may be critical if the ratio r₁/r₂ of the coefficients r provided by models DGM1 and DGM2 is significantly larger than 1.

 

Validation of the discrete Gaussian model

We will check the validity of the assumptions of the discrete Gaussian models for modelling the marginal block distribution. This corresponds to the so-called global change-of-support, in contrast with the local change-of-support which consists of predicting the distribution of block values conditional on data available in the block and its neighbourhood.

The result depends on the spatial distribution of the random function Z and we consider the ideal situation where this spatial distribution is known. In practice, we will consider the case where Z is a lognormal SRF (that is, its logarithm is a Gaussian SRF). We can assume without lack of generality that it has a unit mean. It is then characterized by its logarithmic standard deviation σ and by the correlogram ρ(h) of its logarithm. The random function Z(x) is thus of the form

It can also be expressed as which corresponds to an expansion of the transformation function φ with

The φ_η decrease rapidly when n increases if σ is small, but less rapidly if σ> 1. For σ = 3, for example, φ_η increases up to n = 8 and then decreases slowly. The logarithmic assumption thus includes very contrasted behaviours. This can also be seen on the variance of Z: it is given by

Since Z(x) has unit mean, its coefficient of variation is close to σ when σ is small (0.53 for σ = 0.5) but takes large values when σ is large: 1.31 for σ = 1, 7.32 for σ = 2, 90.01 for σ = 3. A logarithmic standard deviation as large as 3 is far beyond what is seen in mining applications (such a value has been reported for permeability in hydrogeology (e.g. Zimmerman et al., 1998) but the main variable in that case is log-permeability rather than permeability). When dealing with ore grade or pollution concentration, we are of course interested in Z(x) itself rather than its logarithm, and such a logarithmic standard deviation would be extreme. Would it correspond to an actual situation, we would be in a very uncomfortable position because a simple parameter such as the arithmetic mean of Z would require a large number of observations to be estimated reliably. Of course, it is possible to deduce the mean of Z from the mean and the variance of log Z, but then we rely heavily on the assumption of a lognormal distribution, which cannot be taken for granted.

A specificity of lognormal SRFs is that DGM models lead to a lognormal distribution for Z(v), with logarithmic variance r²σ²: This is the well-known permanence of lognormality. It can be shown easily by application of Equations [1] or [4]. Another specificity is that the correspondence between the covariance C(h) of Z and the correlogram ρ(h) of Y takes the simple form

,

which facilitates the computation of the right-hand side of Equation [5].

Validation method

Several authors have already checked the validity of model DGM1 or DGM2 in the lognormal case (Matheron, 1981; Cressie, 2006; Emery, 2007; Chiles and Delfiner, 2012). The principle is to compare the block transformation function given by the DGM model with the 'exact' block transformation function obtained with a large number of simulations.

The first exercise was conducted by Matheron (1981). Owing to the computation capabilities available at that time, this check was limited to the one-dimensional case with an exponential covariance (the random function is then a Markov random process), and used a limited number of simulations. We extend it to two and three dimensions, with a finer discretization of the block and a much larger number of simulations, enabling relatively large values for σ to be considered. The principle of the checking is as follows:

► Consider a block v of the d-dimensional space, defined as the union of M discrete points forming a regular grid

► Build Nunconditional simulations Y_k(x): k = 1, N, of a point-support stationary Gaussian random function Y( x) with zero mean, unit variance, and covariance ρ(h)

► Build the corresponding simulations Z_k(x) = exp (σΥ_k(x)-σ²/2)

► Calculate the simulated block values Z_k(v), namely the average values of Z_k(x) among all x defining the block v

► Sort the Z_k(v) by increasing values; let W_k(v) denote the sorted series

► Use (Wk(v) + Wk₊₁(v)) / 2 as the k/ N th quantile of the distribution of Z(v), that is as the value of φv(y) for y = G^-1(k / N)

► Compare to the values predicted by models DGM1 and DGM2, that is, to exp(r σy - r² σ² / 2) with the corresponding change-of-support coefficient r.

Note that we do not consider the average value at M grid nodes as an approximation to the average in a block v of the d-dimensional space: we substitute the problem in the discrete space to the problem in the continuous case, so that there is no approximation in the approach. The simulations are generated with the discrete spectral method, which produces perfectly Gaussian simulations (up to the quality of the pseudo-random number generator). If it were not possible to exactly reproduce the desired covariance for Y, we would replace it by a covariance as close as possible to it and check this model (this is likely to occur with Gaussian covariances, for example). Note that an alternative would be to build simulations based on covariance matrix decomposition (e.g. Chiles and Delfiner, pp. 493-494); this method also produces perfectly Gaussian simulations, meets any covariance model, but has stronger limitations than the discrete spectral method in terms of size of the grid. For consistency, the integrals in Equations [5] and [6] are replaced by discrete sums. We use a large number of simulations (up to 100 000).

Note that this approach, fully similar to that of Matheron (1981), is slightly different from that used by Emery (2007) to check model DGM2. Indeed, Emery simulated standard Gaussian values with pairwise correlations all equal to r, thus corresponding to the values of Y at independent random points in v.

 

Results

We focus on a spherical correlogram for Y, and on three contrasted block sizes L with respect to the range a of the correlogram: L / a = 0.1, 1, and 10.

When L / a = 0.1, the change-of-support coefficient r₂ of model DGM2 is close to 1 (0.980 in 1D, 0.951 in 2D, 0.933 in 3D) and the coefficient r₁ of model DGM1 remains very close to r₂ even for a large σvalue (for σ = 3 we obtain 0.981 in 1D, 0.955 in 2D, 0.938 in 3D). Both models lead to similar block distributions, very close to the true one (see Chiles and Delfiner, 2012, p. 454 for the 2D case with σ = 1.5).

When L / a = 1, the coefficient r₂ is equal to 0.79 in 1D, 0.59 in 2D, and only 0.46 in 3D. Moreover, r₁ increases significantly with σ when σ exceeds 1, as can be seen in Figure 2.

 

 

Figures 3 and 4 show the results obtained with 100 000 simulations for σ = 1 and 2 respectively. Model DGM1 quite perfectly reproduces the true transformation function - and thus the block distribution. Model DGM2 gives good results as far as large values are not considered but comes with a slight bias for y > 2 when σ = 1, and a significant bias when σ = 2.

 

 

 

 

When L / a = 10, the coefficient r₂ is equal to 0.271 in 1D, 0.077 in 2D, and only 0.022 in 3D, and like in the preceding case r₁ can have much larger values when σ exceeds 1. Model DGM1 presents a slight bias at the extremes of the distribution, whereas DGM2 is biased everywhere but close to the median (see Chiles and Delfiner, 2012, p. 454 in the 2D case with σ = 1.5).

 

Conclusion

This validation exercise shows that the original DGM1 model of Matheron (1976a) gives a very good approximation to the true block distribution, except for extreme values when the logarithmic standard deviation is very large. The variant DGM2 of Emery (2007) is also a very good approximation to the true distribution provided that we are in one of the following situations: (i) the block size is small with respect to the range, (ii) the logarithmic standard deviation is not too large, or (iii) we are not interested in the distribution of high grades. In such a case, this variant can be applied safely. Application of the models, and especially of DGM2, should be done carefully with highly skewed grade distributions and/or very large blocks or panels.

Further checks should be carried out. The approximation of the DGM models has to be quantified also in 1D (time series) and in 3D (the approximation is less valid as the space dimension increases). It is also interesting to examine other covariance models for ρ(h). When σ is fixed, the solution provided by DGM models depends only on r, but two covariance models that would give the same r value do not necessarily give the same true block distribution. Finally, the presence of a nugget effect extends the validity of DGM models, but this has to be quantified.

We have addressed only the global change-of-support problem. The local change-of-support (prediction of the block distribution conditional on neighbouring data) has been examined by Cressie (2006), who considers unbiased lognormal estimators that are exponentials of the simple or ordinary kriging estimators of Y(v), which amounts to the assumptions of model DGM2. The experiment compares this ordinary lognormal kriging with the optimal solution provided by conditional expectation, obtained by a Monte Carlo method. The results indicate that ordinary lognormal kriging performs well in situations where the block size is small with respect to the range, the lognormal standard deviation is not too large, and the neighbourhood is sparse. The first two conditions are required for the global model to be efficient. The third expresses the fact that conditional expectation makes better use of numerous data than an estimator whose form is limited to the exponential of a linear combination of the logarithms of the data.

 

Acknowledgements

An initial version of this paper was presented at the Ninth International Geostatistics Congress, Oslo, June 2012.

 

Bibliography

Chilès, J.-P. and Delfiner, P. 2012. Geostatistics: Modeling Spatial Uncertainty. 2nd edn. John Wiley & Sons, New York, USA.         [ Links ]

Cressie, N. 2006. Block kriging for lognormal spatial processes. Mathematical Geology, vol. 38, no. 4. pp. 413-443.         [ Links ]

Emery, X. 2005. Simple and ordinary multigaussian kriging for estimating recoverable reserves. Mathematical Geology, vol. 37, no. 3. pp. 295-319.         [ Links ]

Emery, X. 2007. On some consistency conditions for geostatistical change-of-support models. Mathematical Geology, vol. 39, no. 2. pp. 205-223.         [ Links ]

Emery, X. 2008. Change of support for estimating local block grade distributions. Mathematical Geology, vol. 40, no. 6. pp. 671-688.         [ Links ]

Krige, D.G. 1978. Lognormal-de Wijsian Geostatistics for Ore Evaluation. Monograph Series, Geostatistics 1. South African Institute/or Mining and Metallurgy, Johannesburg.         [ Links ]

MATHERON, G. 1976a. Forecasting block grade distributions: the transfer functions. Advanced Geostatistics in the Mining Industry. Guarascio, M., David, M., and Huijbregts, C. (eds.). Reidel, Dordrecht, Holland. pp. 237-251.         [ Links ]

MATHERON, G. 1976b. A simple substitute for conditional expectation: the disjunctive kriging. Advanced Geostatistics in the Mining Industry. Guarascio, M., David, M., and Huijbregts, C. (eds.). Reidel, Dordrecht, Holland. pp. 221-236.         [ Links ]

Matheron, G. 1981. Remarques sur le changement de support. Technical Report N-690, Centre de Géostatistique, Fontainebleau, France.         [ Links ]

Zimmerman, D.A., De Marsily, G., Gotway, C.A., Marietta, M.G., Axness, C.L., Beauheim, R., Bras, R., Carrera, J., Dagan, G., Davies, P.B., Gallegos, D.P., Galli, A., Gómez-Hernández, J., Grindrod, P., Gutjahr, A.L., Kitanidis, P.K., Lavenue, A.M., McLaughlin, D., Neuman, S.P,. Ramarao, B.S,. Ravenne, C., and Rubin, Y. 1998. A comparison of seven geostatistically-based inverse approaches to estimate transmissivities for modeling advective transport by groundwater flow. Water Resources Research, vol. 34, no. 6. pp. 1373-1413.         [ Links ] ♦