1. Introduction to the MultiGaussian Model
1.1. Historic perspective
The multiGaussian model was firstly introduced in Geostatistics by Matheron in the 70’s [14] under the name of Disjunctive Kriging, as a way of directly estimating nonlinear functions of stationary grade values. Under this multivariate Gaussian framework, the derivation of conditional distributions, and hence of estimates for recovery functions (grade and tonnage above a cutoff), is particularly straightforward.
Several publications, since then, can be found in the literature which present the theory and applications of disjunctive kriging [10, 6, 16]. The problem considered is the estimation of the attribute in a volume (spatially discretized in an regular way by points ), from point samples with values , located at points , inside or outside the volume . If
is the estimator for the unknown value of the volume, it is usually desirable that this estimator presents properties such as unbiasedness, conditional unbiasedness and minimum error variance. Different estimators possess these properties, such as kriging, suitable if the sample values are normally distributed, and logarithmic kriging, when lognormality prevails. However, estimators derived from the multiGaussian model have the advantage of not only having these properties, but also of being distribution free, i.e., the model can be used whatever the multivariate probability distribution of the sample values
[18].1.2. Notion
The idea of the multivariate Gaussian approach is to transform the initial random function , into a random function with a standard Gaussian univariate distribution. Then, under the working hypothesis that is spatially multiGaussian, the conditional distributions of are determined. Derivation of conditional expectation in a nonlinear way, and estimates of recovery functions, are then obtained through inverse transformation from these conditional distributions.
In practice, the nonlinear estimator provided by this procedure differs from the true conditional expectation. However, and under the condition that the univariate distribution of is well known, practice has shown that the mean estimator of this nonlinear technique, , is generally better than the estimators provided by the direct linear kriging processes applied to the initial data (see [19, 7, 8]). Some caution, however, should be taken with very continuous variables, which may result in the overestimation of areas around highvalued samples and, consequently, in a dangerous bias if used for selecting mining blocks [20].
1.3. Range of validity
The multiGaussian approach is very convenient: the inference of the conditional distribution function (cdf) reduces to solving a simple kriging system at location u
. The tradeoff cost is the assumption that data follow a multiGaussian distribution, which implies first that the onepoint distribution of data (sample histogram) is normal.
Many variables in the earth sciences show an asymmetric distribution with a few very large values (positive skewness). Thus, the multiGaussian approach starts with an identification of the standard distribution and involves the “normalizing” of the skewed sample histogram. The common approach for this normalization is to apply a normal score transformation. This procedure will be described in detail in the following sections.
However, the normality of the onepoint cdf is a necessary but not sufficient condition to ensure that the random function model is multivariate Gaussian. Hence one should check whether the normalized data are also reasonably bivariate Gaussian. There are several ways to check that the twopoint distribution of data , is bivariate normal. For a detailed exposition, and many explanations, the reader can check [4, 3].
1.4. Conceptual problems
The normality of onepoint and twopoint cdfs are both necessary but not sufficient conditions to ensure that a multivariate Gaussian random function model is appropriate for modeling the spatial distribution of normal score data. One should also check the normality of the threepoint until the point cdfs. One type of check consist of comparing experimental multiplepoint frequencies with their theoretical Gaussian values through the comparison with analytical expressions. Although such an expression has been established, the main difficulty resides in the inference of such experimental frequencies. For example, the inference of a threepoint frequency such that
requires the availability of a series of triplet values with the same geometric configuration as the triplet .
Nonregular drilling grids and data sparsity prevent us from computing sample statistics involving more than two locations at a time. Therefore, in practice, because most of testing becomes tedious, if scatter plots do not invalidate the biGaussian assumption, the multiGaussian formalism is adopted.
2. Main result
Our principal concern, when conceiving a tool capable of the modeling spatial uncertainty at any support, is to take into consideration, on one hand, the spatial location of the samples in a way that the uncertainty is reduced when we approach to the samples and, on the other hand, to reproduce the heteroscedasticity that we find in real deposits, i.e., to model the zones of higher values with higher variance. This heteroscedastic behavior is commonly referred to as the proportional effect.
The main development, presented in what follows, has a rather simply but powerful origin, and is based in the following variance calculation for an average of values within a volume:
(2.1) 
with the amounts of points to average in . This result allows us to compute the variance of a volume in a simple way, by taking an “average” of the covariance between pairs of points within the volume.
The key step, that follows to this result, is to consider this covariance between the pairs of points , without the assumption of stationarity in the original values, but in their respective Gaussian values . Then, the covariance in the original values is retrieved as a nonlinear function of the covariance and also a function of the conditioning sample values. This nonlinear function will be derived in a straightforward way once the multiGaussian model is presented.
3. The MultiGaussian Model
In the multiGaussian model, the attribute under study is viewed as a realization of a random function that can be transformed into a Gaussian random field with zero mean and unit variance, in a certain position u
. This transformation is a quantile transformation (
u will be implicit most of the time hereafter):(3.1) 
where and , commonly written as (Figure 1). Hence,
or just . Lets define the cumulative density function (cdf) of a Gaussian random function with mean and variance , as , except for the standard case , written just as . Then , in shortened notation, is the anamorphosis function.
We notice that the quantile transformation (3.1), written in terms of the probability density and the anamorphosis function, is equivalent to:
leading to:
where and
is the probability density function (pdf) of a standard Gaussian distribution, commonly written as
. We will define, in general, as the pdf of a Gaussian random function with mean and variance (except for the standard case):(3.2) 
Knowing the anamorphosis of a variable is equivalent to knowing its distribution.
The multiGaussian assumption states that the distribution of any value conditioned by the sampled values is still Gaussian, with a mean and a variance equal to its simple kriging estimate and simple kriging variance respectively. Hence, the conditional cumulative probability function (ccdf) is
where , and where “data” represents the conditioning data . The conditional cdf of the original variable is then retrieved as:
where . This can be rewritten in terms of the probability densities:
=  
= 
from where follows:
(3.3) 
From here we can obtain the conditional density function without incurring in any quantile sampling of the posterior Gaussian distribution to construct the local distribution in the raw values:
(3.4) 
We will say that . Let us define hereafter.
Multiplying both sides by in (3.3) and integrating with respect to lead us to:
(3.5) 
By taking , we obtain :
(3.6) 
Let us insist that , however, from (3.6), it is clear that is not just equal to , which usually can lead to confusion and wrong results.
In the same manner, it is easy to obtain any moment (if they exist) under the multiGaussian model:
leading to an expression to compute the variance in a general way for the raw variable:
(3.7)  
We can extend the model into a variate model in the following way: the joint cumulative function of random variables (related to ) at locations , respectively, conditioned by the sampled values, follows a multivariate Gaussian distribution
with mean vector
and covariance matrix equal to , , such that and , .In the bivariate case the model looks as follows:
=  
=  
=  
= 
with:
=  
=  
=  
= 
This can be rewritten in terms of probability densities:
=  
= 
which leads to the following equality:
(3.8) 
Hence, we can obtain the probability distribution by:
(3.9) 
Multiplying both sides by in 3.8 and integrating with respect to and we obtain:
By taking , , we obtain :
Combining the previous results, we obtain a general expression for the covariance under the multiGaussian model. We can see explicitly the covariance as a measure of stochastic dependence between the locations and :
(3.10)  
Then, by plugging in this result in Equation 2.1, we obtain a general answer for the assessment of the variance in the volume under the multiGaussian framework.
In the nvariate case the mean vector and covariance matrix of the model looks as follows:
=  
=  

and the ccdf is extended by the link with a variate Gaussian cdf:
=  
=  
=  
= 
3.1. Examples
3.1.1. The LogNormal Case
If is a normal random variable with mean and variance , then the random variable is said to be lognormal with parameters and , or just . Thus, a random variable is lognormal if is normal.
It follows that the anamorphosis function is obtained by using 3.1:
Hence:
(3.11) 
From here we obtain:
(notice that ) and that the conditioned local distribution is given by:
=  
=  
= 
This is a lognormal distribution with parameters
and for the mean and variance, respectively. Therefore, the local distribution at a certain point u conditioned by the data, which presents a lognormal prior distribution, preserves the lognormality.Following the same procedure and extending (3.9) to the dimesional case, lead us to compute that the local conditional distribution of the vector of random variables , each one with prior distribution, has a conditional distribution such that has an dimensional normal distribution with mean vector and covariance matrix , , such that and , , given the data.
In Figure 2, some possible posterior distributions and their bivariate behavior are presented.
3.1.2. The Exponential Case
is said to have an exponential distribution with parameter , , or just , if its probability density function is given by , or, equivalently, if its cdf is given by
It follows that the anamorphosis function is:
(3.12) 
From here we obtain:
We will not attempt to find out what kind of probability distribution results of the backtransformation, as we did with the lognormal case. However, numerical results of how some posterior distributions looks like are presented in Figure 3
. It is very interesting how the exponential distribution is not preserved any more in the posterior distributions (only when the posterior distribution follows a
).We present a synthetic case where the simulation procedure is compared with the results of the direct application of the formulations presented. The synthetic reality for the analysis is generated by LU nonconditional simulation with a known omni directional variogram in a 100x100 grid of 5x5 units spacing between nodes, for the Gaussian values. Then, the synthetic scenario is randomly sampled, obtaining 100 samples. An exponential distribution with parameter is used.
4. A Note Towards the Change of Support
In mining, the actual selection is made on panels, and not on core samples with size . Different panels may be blended into a larger volume (in downstream processes). Therefore, for the estimation of reserves, it is essential to take the support (size and geometry) of the selection unit into account (resulting in several hundred or thousand tones). Otherwise there is a risk of bias with serious economic consequences.
In general, this support will be considerably different from the support of the exploration unit (core samples for example). This is one of the most important aspects on the economic and technical context of a mining project, with consequences in the amount of recoverable reserves. The entire resources of a deposit are rarely mineable because of the variability of mining and treatment costs added to the spatial variability of ore quality. For a detailed explanation of the issue, the reader can consult the survey by Huijbregts [5]. In a nutshell, consider the histogram of dispersion of the grades as on Figure 5. The histogram of grades of all the mining units of larger support will have: (i) a mean equal to the mean of the core sample grades; (ii) an experimental dispersion variance less than the core support; and (iii) for high cutoff grades above the mean, , the hatched area of the histogram of core grades may seriously overestimate the possible mined recovery in both tonnage and mean grade and, correspondingly, underestimate the tonnage of metal left in the waste, i.e., underestimate the area corresponding to the panels with true grades .
In what follows, we present an extension of the multiGaussian model to address the distribution when a change of support is considered, locally. We are interested in finding the probability distribution of the main variable within a given volume of the spatial domain. For this, let’s say the volume is approximated by a discretization of points, located at positions, , each one with mass density , , we are interested in finding the probability distribution of the random variable:
(4.1) 
recalling that the pointsupport random variables are correlated, and they have to be conditioned by the sampled data, which, present themselves in different configurations and with different conditioning values, leading to different pointsupport distributions.
In this discretized framework, we can consider constant within the geological volume, and we omit, without loss of generality, this constant for the next steps. The problem can be addressed as follows.
Let’s begin by expressing the cumulative distribution of the sum:
=  
=  
=  
= 
This can be rewritten in term of the variable and the anamorphosis function, for a given value such that :
=  

=  = 
Then, the probability density can be determined joining the following two developments. On one hand:
and on the other hand:
This shows that we are moving from the volume of integration parameterized by , when we compute , to the surface of this volume when we compute , which defines planes in term of the variable but not necessarily in terms of the Gaussian variable .
Notice that is now determined by the rest of the variables. Thus, we can just rewrite . Therefore,
(4.2) 
The procedure described is summarized in Figure 6, where we want to compute the sum of two correlated lognormal random variables. The straight lines where integration needs to be done, in the bilognormal distribution space, are “bent”, together with the density itself, being transformed into curved lines in the Gaussian framework.
Remark 4.1.
Although this procedure gives exact results, it is unfeasible in practice because of its computational complexity, requiring to integrate the multiGaussian density in . However, a highly suitable alternative could be to consider MonteCarlo integration of the integral 4.2. One can argue that this would be similar to using geostatistical simulations. It is not, since we can obtain immediately answers such as the probability above a cutoff value z, , by computing in a discrete set of points in the interval , in contrast with geostatistical simulations, where there is no control in the outcomes range and values.
We show an example of reblocking, where the exact distribution of the points and the average of points are presented using expression 4.2, together with the outcomes of the simulation, at point and block support. We have been able to implement it with , in the exponential example of Figure 4. The results are summarized in Figure 7.
Remark 4.2.
If the spatial discretization is infinitesimal, expression 4.1 can be formally expressed as:
(4.3) 
with . Supposing we are considering an anamorphosis transformation given by an initial lognormal distribution with mean and variance , then 4.3 can be expressed as:
where is the conditioned Gaussian random variable. Some approximations to this form of the problem can be found in [13].
5. Hermite Polynomials
In this section Hermite Polynomials are introduced as an infinite sum alternative (truncated at some term ), suitable for avoiding numerical integration in some of the formulations already presented. Firstly, we begin by introducing Hermite Polynomials, and secondly, we present a selfcontained list of proofs for most of the important results that we will use in the following section. The proofs are taken from [9, 14, 1, 2].
Hermite Polynomials are polynomials that have special properties related to the normal distribution. Hermite Polynomials are defined by Rodrigues’ formula:
(5.1) 
where , is a normalization factor, is a Gaussian value, and is the standard Gaussian probability distribution function (defined in Equation 3.2). The Hermite Polynomial is a polynomial of degree . More specifically:
=  
=  
= 
Theorem 5.1.
is a polynomial of degree that satisfy the recurrence relations:
Proof.
Calculating the successive derivatives of the Gaussian density:
and, in general
(5.2) 