A general approach to the assessment of uncertainty in volumes by using the multi-Gaussian model

TL;DR

This paper introduces a general, non-simulation-based method for assessing uncertainty in volumes conditioned on sampling data, using an extended multi-Gaussian model and Kriging to derive local distributions.

Contribution

It develops a numerical tool that explicitly relates volume uncertainty to data grades, spatial correlation, and conditioning values, applicable to any probabilistic distribution.

Findings

Provides explicit formulas for local uncertainty measures.

Enables straightforward computation of local distributions.

Applicable to arbitrary probabilistic data distributions.

Abstract

The goal of this research is to derive an approach to assess uncertainty in an arbitrary volume conditioned by sampling data, without using geostatistical simulation. We have accomplished this goal by deriving an numerical tool suitable for any probabilistic distribution of the sample data. For this, we have worked with an extension of the traditional multi-Gaussian model, allowing us to obtain a formulation that makes explicit the dependence of the uncertainty in the arbitrary volume from the grades within the volume, the spatial correlation of the data and the conditioning values. A Kriging of the Gaussian values is the only requirement to obtain not only conditional local means and variances but also the complete local distributions at any support, in an easy and straightforward way.

Figures

Equations

$$P(Y(\text{u})\le y_p,Y({\text{u}}^{\prime })\le y_p,Y({\text{u}}^{\prime \prime })\le y_p)$$

$$\mathbb{V}ar(Z_{{\mathbf{V}}_{\mathbf{k}}}|Z({\text{u}}_{\alpha })=z({\text{u}}_{\alpha }),\alpha \in \{1,\dots ,n(\text{k})\})=\frac{1}{|{\mathbf{V}}_{\mathbf{k}}{|}^2}\sum _{\begin{array}{c}i:{\text{u}}_i\in {\mathbf{V}}_{\mathbf{k}}\\ j:{\text{u}}_j\in {\mathbf{V}}_{\mathbf{k}}\end{array}}\mathbb{C}ov(Z^{\ast }({\text{u}}_i),Z^{\ast }({\text{u}}_j)),$$

$$F_Z(z)=F_Y(y),$$

$$Z(\text{u})=\varphi [Y(\text{u})]$$

$${\int }_{-\infty }^z{f}_Z(\text{z})d\text{z}={\int }_{-\infty }^{\varphi (y)}{f}_Z(\text{z})d\text{z}={\int }_{-\infty }^y{f}_Y(\text{y})d\text{y},$$

$$f_Z(z)=\frac{f_Y(y)}{{\varphi }^{{}^{\prime }}(y)}=\frac{g(y)}{{\varphi }^{{}^{\prime }}(y)},$$

$$g(y)=\frac{1}{\sqrt{2\pi }}\cdot e^{-\frac{y^2}{2}},g_{\mu }^{{\sigma }^2}(y)=\frac{1}{\sigma \sqrt{2\pi }}\cdot e^{-\frac{(y-\mu {)}^2}{2{\sigma }^2}}.$$

$$F_{Y|data}(y)=G\Bigg{(}\frac{y-y^{*}_{SK}}{\sigma_{SK}}\Bigg{)},$$

$$F_{Z|data}(z)$$

$$\displaystyle=\displaystyle G\Bigg{(}\frac{y-y^{*}_{SK}}{\sigma_{SK}}\Bigg{)}=G_{y^{*}_{SK}}^{\sigma_{SK}^{2}}(y),$$

$$f_{Z|data}(z)\cdot {\varphi }^{{}^{\prime }}(y)=g_{y_{SK}^{\ast }}^{{\sigma }_{SK}^2}(y).$$

$$f_{Z|data}(z)=\frac{1}{{\varphi }^{{}^{\prime }}(y)}\cdot g_{y_{SK}^{\ast }}^{{\sigma }_{SK}^2}(y)=\frac{1}{{\varphi }^{{}^{\prime }}({\varphi }^{-1}(z))}\cdot g_{y_{SK}^{\ast }}^{{\sigma }_{SK}^2}({\varphi }^{-1}(z)).$$

$${\int }_{-\infty }^{\varphi (y)}\text{z}\cdot f_{Z|data}(\text{z})d\text{z}$$

$$\displaystyle=\displaystyle\frac{1}{\sigma_{SK}}\cdot\int_{-\infty}^{y}\phi(\textbf{y})\cdot g\Bigg{(}\frac{\textbf{y}-y^{*}_{SK}}{\sigma_{SK}}\Bigg{)}d\textbf{y}$$

$$={\int }_{-\infty }^y\varphi ({\sigma }_{SK}\cdot \text{y}+y_{SK}^{\ast })\cdot g(\text{y})d\text{y}.$$

$$[Z^{\ast }(\text{u})]={\int }_{-\infty }^{\infty }\varphi (\text{y})\cdot g_{y_{SK}^{\ast }}^{{\sigma }_{SK}^2}(\text{y})d\text{y}.$$

$$\mathbb{E}[Z^n(\text{u})|data]={\int }_{-\infty }^{\infty }{\varphi }^n(\text{y})\cdot g_{y_{SK}^{\ast }}^{{\sigma }_{SK}^2}(\text{y})d\text{y},$$

$$\mathbb{V}ar[Z(\text{u})|data]$$

$$={\int }_{-\infty }^{\infty }{\varphi }^2(\text{y})\cdot g_{y_{SK}^{\ast }}^{{\sigma }_{SK}^2}(\text{y})d\text{y}-[Z^{\ast }(\text{u}){]}^2.$$

$$f_{Z_0,Z_1|data}(\varphi (y_0),\varphi (y_1))\cdot {\varphi }^{\prime }(y_0)\cdot {\varphi }^{\prime }(y_1)=g_{{\mathbf{\mu }}_{SK}}^{{\mathbf{\Sigma }}_{SK}}(y_0,y_1).$$

$$f_{Z_0,Z_1|data}(z_0,z_1)=\frac{1}{{\varphi }^{\prime }(y_0)}\cdot \frac{1}{{\varphi }^{\prime }(y_1)}\cdot g_{{\mathbf{\mu }}_{SK}}^{{\mathbf{\Sigma }}_{SK}}(y_0,y_1).$$

$${\int }_{-\infty }^{\varphi (y_0)}{\int }_{-\infty }^{\varphi (y_1)}{\mathbf{z}}_0{\mathbf{z}}_1{f}_{Z_0,Z_1|data}({\mathbf{z}}_0,{\mathbf{z}}_1)d{\mathbf{z}}_{0}d{\mathbf{z}}_1={\int }_{-\infty }^{y_0}{\int }_{-\infty }^{y_1}\varphi ({\mathbf{y}}_0)\varphi ({\mathbf{y}}_1)g_{{\mathbf{\mu }}_{SK}}^{{\mathbf{\Sigma }}_{SK}}({\mathbf{y}}_0,{\mathbf{y}}_1)d{\mathbf{y}}_{0}d{\mathbf{y}}_1.$$

$$\mathbb{E}[Z({\text{u}}_0{)}^{\ast }Z({\text{u}}_1{)}^{\ast }]={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\varphi ({\mathbf{y}}_0)\varphi ({\mathbf{y}}_1)g_{{\mathbf{\mu }}_{SK}}^{{\mathbf{\Sigma }}_{SK}}({\mathbf{y}}_0,{\mathbf{y}}_1)d{\mathbf{y}}_{0}d{\mathbf{y}}_1.$$

$$\mathbb{C}ov(Z^{\ast }({\text{u}}_0),Z^{\ast }({\text{u}}_1))$$

$$={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\varphi ({\mathbf{y}}_0)\varphi ({\mathbf{y}}_1)g_{{\mathbf{\mu }}_{SK}}^{{\mathbf{\Sigma }}_{SK}}({\mathbf{y}}_0,{\mathbf{y}}_1)d{\mathbf{y}}_{0}d{\mathbf{y}}_1$$

$$-{\int }_{-\infty }^{\infty }\varphi (\text{y})\cdot g_{y_{SK}^{\ast }({\text{u}}_0)}^{{\sigma }_{SK}^2({\text{u}}_0)}(\text{y})d\text{y}\cdot {\int }_{-\infty }^{\infty }\varphi (\text{y})\cdot g_{y_{SK}^{\ast }({\text{u}}_1)}^{{\sigma }_{SK}^2({\text{u}}_1)}(\text{y})d\text{y}.$$

$$F_{Z}(z)=G\Bigg{(}\frac{ln(z)-\mu}{\sigma}\Bigg{)}.$$

$$G\Bigg{(}\frac{ln(z)-\mu}{\sigma}\Bigg{)}=G(y),$$

$$z=\varphi (y)=e^{\mu }\cdot e^{\sigma y},y={\varphi }^{-1}(z)=\frac{ln(z)-\mu }{\sigma }.$$

$${\varphi }^{{}^{\prime }}(y)=\sigma \cdot e^{\mu }\cdot e^{\sigma y}=\sigma \varphi (y),$$

$$F_Z(z)=1-e^{-\lambda z},z\ge 0.$$

$$z=\varphi (y)=-\frac{1}{\lambda }\cdot ln(1-G(y)).$$

$${\varphi }^{{}^{\prime }}(y)=\frac{1}{\lambda }\cdot \frac{1}{1-G(y)}\cdot g(y),$$

Full text

1. 1 Introduction to the Multi-Gaussian Model
2. 2 Main result
3. 3 The Multi-Gaussian Model
4. 4 A Note Towards the Change of Support
5. 5 Hermite Polynomials
6. 6 Application of Hermite Polynomials Expansions in a General Case
7. 7 Conclusion

Keywords. Geostatistics, conditional expectation, conditional variance, simulations, multi-gaussian model, Hermite polynomials.

The Multi-Gaussian model is a very simple non linear technique used in the context of spatial estimation of variables, introduced by Matheron in the 70’s [12, 16, 14]. It solves many classic problems regarding spatial estimation, when the linear estimation does not provide a good solution.

One of the most successful applications of the model has relationship with the simulation of ore deposits, producing different plausible scenarios, conditioned to a set of samples (i.e. matching the data in their positions and values), reproducing their spatial variability. These scenarios are post-processed, using non-linear transfer functions (usually related with the application of cutoff grades in open-pit design, mine planning and exploration) in order to assess the impact of the geological variability on mineral resources and reserves.

The multi-Gaussian model also helps obtaining estimators of local conditional expectations and variances, at point and block support, including arbitrary volumes [16, 17, 11, 18]. The discrete Gaussian model [15, 2] can also be used to provide the change of support distribution locally, however it does not consider the proportional effect or the blending when an arbitrary number of different zones is considered.

The goal of this research is related with this last problem. Simulations allow to address straightforwardly the probability distribution of variables over any volume, just by changing the support of the outcomes over different scenarios. However, its result depend on the number of realizations and, therefore, provides an inexact result.

We present next the steps to close this gap by analyzing the inference of the mean estimator, then the variances and finally the complete conditional distributions, at point support and when change of support is considered. This is done considering both the conditioning effect of sample locations and the proportional effect, without resorting to the use of simulations.

The plan of the paper is as follows. In Section 1 we provide an introduction related with the history and the range of validity of the model. Then in Section 2, we comment on our procedure to obtain the variance for any support. In Section 3, we formally establish the multi-Gaussian model. After the statement, it is rather easy to derive point $n$-variate probability distributions. Some examples are presented at the end of that section, of the application of the model and its comparison with standard simulation results. In Section 4, we propose our methodology for change of support with an example, providing a general and simple approach. Finally, Section 5 provides an introduction and classic results of Hermite Polynomials. Hermite Polynomials are highly suitable to be incorporated into the Gaussian model. This section serves as a compilation of proofs to the reader and can be read independently of the rest. An elegant result of the application of Hermite Polynomials is presented in Section 6. Final comments of these results and future work are provided in Section 7.

1. Introduction to the Multi-Gaussian Model

1.1. Historic perspective

The multi-Gaussian 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 non-linear 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 $Z$ in a volume $\boldsymbol{V}$ (spatially discretized in an regular way by points $\textbf{u}_{i},i\in\{1,\dots,|\boldsymbol{V}|\}$), from $n$ point samples with values $z(\textbf{u}_{\alpha}),\alpha\in\{1,\dots,n\}$, located at points $\textbf{u}_{\alpha},\alpha\in\{1,\dots,n\}$, inside or outside the volume $\boldsymbol{V}$. If $Z^{*}_{\boldsymbol{V}}$ 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 multi-Gaussian 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 $Z(\textbf{u})$, into a random function $Y(\textbf{u})$ with a standard Gaussian univariate distribution. Then, under the working hypothesis that $Y(\textbf{u})$ is spatially multi-Gaussian, the conditional distributions of $Y(\textbf{u})$ are determined. Derivation of conditional expectation in a non-linear way, and estimates of recovery functions, are then obtained through inverse transformation from these $Y(\textbf{u})$ conditional distributions.

In practice, the non-linear estimator provided by this procedure differs from the true conditional expectation. However, and under the condition that the univariate distribution of $Z(\textbf{u})$ is well known, practice has shown that the mean estimator of this non-linear technique, $[Z^{*}(\textbf{u})]$, 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 high-valued samples and, consequently, in a dangerous bias if used for selecting mining blocks [20].

1.3. Range of validity

The multi-Gaussian approach is very convenient: the inference of the conditional distribution function (cdf) reduces to solving a simple kriging system at location u. The trade-off cost is the assumption that data follow a multi-Gaussian distribution, which implies first that the one-point 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 multi-Gaussian 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 one-point 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 two-point distribution of data $y(\textbf{u}_{\alpha}),\alpha\in\{1,\dots,n(\textbf{u})\}$, is bivariate normal. For a detailed exposition, and many explanations, the reader can check [4, 3].

1.4. Conceptual problems

The normality of one-point and two-point 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 three-point until the $N$-point cdfs. One type of check consist of comparing experimental multiple-point 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 three-point frequency such that

[TABLE]

requires the availability of a series of triplet values with the same geometric configuration as the triplet $(\textbf{u},\textbf{u},\textbf{u}^{\prime\prime})$.

Non-regular 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 $h$-scatter plots do not invalidate the bi-Gaussian assumption, the multi-Gaussian 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:

[TABLE]

with $|\boldsymbol{V_{k}}|$ the amounts of points to average in $\boldsymbol{V_{k}}$. 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 $\mathbb{C}ov(Z^{*}(\textbf{u}_{i}),Z^{*}(\textbf{u}_{j}))$, without the assumption of stationarity in the original values, but in their respective Gaussian values $\mathbb{C}ov(Y^{*}(\textbf{u}_{i}),Y^{*}(\textbf{u}_{j}))$. Then, the covariance in the original values is retrieved as a non-linear function of the covariance $\mathbb{C}ov(Y^{*}(\textbf{u}_{i}),Y^{*}(\textbf{u}_{j}))=\mathbb{C}ov(\textbf{u}_{i}-\textbf{u}_{j})$ and also a function of the conditioning sample values. This non-linear function will be derived in a straightforward way once the multi-Gaussian model is presented.

3. The Multi-Gaussian Model

In the multi-Gaussian model, the attribute under study is viewed as a realization of a random function $Z(\textbf{u})$ that can be transformed into a Gaussian random field $Y(\textbf{u})$ 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):

[TABLE]

where $F_{Z}(z)=P(Z\leq z)$ and $F_{Y}(y)=P(Y\leq y)$, commonly written as $G(y)$ (Figure 1). Hence,

[TABLE]

or just $z=\phi(y)$. Lets define the cumulative density function (cdf) of a Gaussian random function with mean $\mu$ and variance $\sigma^{2}$, as $G_{\mu}^{\sigma^{2}}$, except for the standard case $G_{0}^{1}$, written just as $G$. Then $\phi=F^{-1}_{Z}\circ G$, 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:

[TABLE]

leading to:

[TABLE]

where $\displaystyle\phi^{{}^{\prime}}(y)=\frac{d\phi(y)}{dy}$ and $f_{Y}(y)$ is the probability density function (pdf) of a standard Gaussian distribution, commonly written as $g(y)$. We will define, in general, $g_{\mu}^{\sigma^{2}}$ as the pdf of a Gaussian random function with mean $\mu$ and variance $\sigma^{2}$ (except for the standard case):

[TABLE]

Knowing the anamorphosis of a variable $Z(\textbf{u})$ is equivalent to knowing its distribution.

The multi-Gaussian assumption states that the distribution of any value $Y(\textbf{u})$ conditioned by the sampled values is still Gaussian, with a mean and a variance equal to its simple kriging estimate $y^{*}_{SK}(\textbf{u})$ and simple kriging variance $\sigma_{SK}^{2}(\textbf{u})$ respectively. Hence, the conditional cumulative probability function (ccdf) is

[TABLE]

where $F_{Y|data}(y)=P(Y\leq y|data)=G(y|data)$, and where “data” represents the conditioning data $Y(\textbf{u}_{\alpha})=y(\textbf{u}_{\alpha}),\alpha\in\{1,\dots,n(\textbf{k})\}$. The conditional cdf of the original variable is then retrieved as:

[TABLE]

where $F_{Z|data}(z)=P(Z\leq z|data)$. This can be rewritten in terms of the probability densities:

[TABLE]

from where follows:

[TABLE]

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:

[TABLE]

We will say that $Z|data\sim\boldsymbol{\phi}(y^{*}_{SK},\sigma_{SK}^{2})$. Let us define $Z^{*}:=Z|data$ hereafter.

Multiplying both sides by $\phi(y)$ in (3.3) and integrating with respect to $y$ lead us to:

[TABLE]

By taking $y\to\infty$, we obtain $[Z^{*}(\textbf{u})]:=\mathbb{E}[Z(\textbf{u})|data]$:

[TABLE]

Let us insist that $[Y^{*}(\textbf{u})]=y^{*}_{SK}$, however, from (3.6), it is clear that $[Z^{*}(\textbf{u})]$ is not just equal to $\phi([Y^{*}(\textbf{u})])$, 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 multi-Gaussian model:

[TABLE]

leading to an expression to compute the variance in a general way for the raw variable:

[TABLE]

We can extend the model into a $n$-variate model in the following way: the joint cumulative function of $n$ random variables $Y(\textbf{u}_{1}),\dots,Y(\textbf{u}_{n})$ (related to $Z(\textbf{u}_{1}),\dots,Z(\textbf{u}_{n})$) at locations $\textbf{u}_{1},\dots,\textbf{u}_{n}$, respectively, conditioned by the sampled values, follows a multivariate Gaussian distribution $\mathcal{N}(\boldsymbol{\mu}_{SK},\boldsymbol{\Sigma}_{SK})$ with mean vector $\boldsymbol{\mu}=(y^{*}_{SK}(\textbf{u}_{1}),\dots,y^{*}_{SK}(\textbf{u}_{n}))^{\prime}$ and covariance matrix equal to $\boldsymbol{\Sigma}=(\boldsymbol{\Sigma})_{ij}$, $i,j\in\{1,\dots,n\}$, such that $(\boldsymbol{\Sigma})_{ij}=(\boldsymbol{\Sigma})_{ji}=\sigma_{SK}(\textbf{u}_{i},\textbf{u}_{j})=\sigma_{SK}(\textbf{u}_{j},\textbf{u}_{i})$ and $(\boldsymbol{\Sigma})_{ii}=\sigma_{SK}^{2}(\textbf{u}_{i})$, $i\in\{1,\dots,n\}$.

In the bi-variate case the model looks as follows:

[TABLE]

with:

[TABLE]

This can be rewritten in terms of probability densities:

[TABLE]

which leads to the following equality:

[TABLE]

Hence, we can obtain the probability distribution by:

[TABLE]

Multiplying both sides by $\phi(y_{0})\cdot\phi(y_{1})$ in 3.8 and integrating with respect to $y_{0}$ and $y_{1}$ we obtain:

[TABLE]

By taking $y_{0}$, $y_{1}\to\infty$, we obtain $\mathbb{E}[Z(\textbf{u}_{0})^{*}Z(\textbf{u}_{1})^{*}]:=\mathbb{E}[Z(\textbf{u}_{0})Z(\textbf{u}_{1})|data]$:

[TABLE]

Combining the previous results, we obtain a general expression for the covariance under the multi-Gaussian model. We can see explicitly the covariance as a measure of stochastic dependence between the locations $\textbf{u}_{0}$ and $\textbf{u}_{1}$:

[TABLE]

Then, by plugging in this result in Equation 2.1, we obtain a general answer for the assessment of the variance in the volume $\boldsymbol{V}$ under the multi-Gaussian framework.

In the n-variate case the mean vector and covariance matrix of the model looks as follows:

[TABLE]

and the ccdf is extended by the link with a $n$-variate Gaussian cdf:

[TABLE]

3.1. Examples

3.1.1. The Log-Normal Case

If $X$ is a normal random variable with mean $\mu$ and variance $\sigma^{2}$, then the random variable $Z=e^{X}$ is said to be log-normal with parameters $\mu$ and $\sigma^{2}$, or just $Z\sim Logn(\mu,\sigma^{2})$. Thus, a random variable $Z$ is log-normal if $ln(Z)$ is normal.

The cumulative distribution function is given by:

[TABLE]

It follows that the anamorphosis function is obtained by using 3.1:

[TABLE]

Hence:

[TABLE]

From here we obtain:

[TABLE]

(notice that $\phi^{{}^{\prime}}(\phi^{-1}(z))=\sigma\phi(\phi^{-1}(z))=\sigma z$) and that the conditioned local distribution is given by:

[TABLE]

This is a log-normal distribution with parameters $\mu+\sigma y^{*}_{SK}$ and $\sigma^{2}\sigma_{SK}^{2}$ for the mean and variance, respectively. Therefore, the local distribution at a certain point u conditioned by the data, which presents a log-normal prior distribution, preserves the log-normality.

Following the same procedure and extending (3.9) to the $n$-dimesional case, lead us to compute that the local conditional distribution of the vector of random variables $(Z_{1},\dots,Z_{n})^{\prime}$, each one with $Logn(\mu,\sigma^{2})$ prior distribution, has a conditional distribution such that $(X_{1},\dots,X_{n})^{\prime}=(ln(Z_{1}),\dots,ln(Z_{n}))^{\prime}$ has an $n$-dimensional normal distribution with mean vector $\boldsymbol{\mu}=(\mu+\sigma y^{*}_{SK}(\textbf{u}_{1}),\dots,\mu+\sigma y^{*}_{SK}(\textbf{u}_{n}))^{\prime}$ and covariance matrix $\boldsymbol{\Sigma}=(\boldsymbol{\Sigma})_{ij}$, $i,j\in\{1,\dots,n\}$, such that $(\boldsymbol{\Sigma})_{ij}=(\boldsymbol{\Sigma})_{ji}=\sigma^{2}\sigma_{SK}(\textbf{u}_{i},\textbf{u}_{j})$ and $(\boldsymbol{\Sigma})_{ii}=\sigma^{2}\sigma_{SK}^{2}(\textbf{u}_{i})$, $i\in\{1,\dots,n\}$, given the data.

In Figure 2, some possible posterior distributions and their bi-variate behavior are presented.

3.1.2. The Exponential Case

A continuous random variable $Z$ is said to have an exponential distribution with parameter $\lambda$, $\lambda>0$, or just $Z\sim Exp(\lambda)$, if its probability density function is given by $f_{Z}(z)=\lambda e^{-\lambda z}$, $z\geq 0$ or, equivalently, if its cdf is given by

[TABLE]

It follows that the anamorphosis function is:

[TABLE]

From here we obtain:

[TABLE]

We will not attempt to find out what kind of probability distribution results of the backtransformation, as we did with the log-normal 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 $\mathcal{N}(0,1)$).

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 non-conditional simulation with a known omni directional variogram $\gamma(\mathbf{h})=0.1\cdot Nugg+0.9\cdot Sph(100)$ 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 $\lambda=1/0.8$ 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 $\sim\textbf{u}$. Different panels may be blended into a larger volume $\boldsymbol{V}$ (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 $Z(\textbf{u})$ 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, $z>m$, 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 $z<m$.

In what follows, we present an extension of the multi-Gaussian 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 $n$ points, located at $\textbf{u}_{i}$ positions, $i=1,\dots,n$, each one with mass density $\rho(\textbf{u}_{i})$, $\sum_{i=1}^{n}\rho(\textbf{u}_{i})=1$, we are interested in finding the probability distribution of the random variable:

[TABLE]

recalling that the point-support 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 point-support distributions.

In this discretized framework, we can consider $\rho(\textbf{u})$ 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:

[TABLE]

This can be rewritten in term of the $Y$ variable and the anamorphosis function, for a given value $a$ such that $z=\phi(a)$:

[TABLE]

Then, the probability density can be determined joining the following two developments. On one hand:

[TABLE]

and on the other hand:

[TABLE]

This shows that we are moving from the volume of integration $\mathcal{V}(a)$ parameterized by $\mathcal{V}(a)=\{z_{1},\dots,z_{n}|z_{1}+\dots+z_{n}\leq\phi(a)\}$, when we compute $F_{Z^{*}_{\boldsymbol{V}}}$, to the surface of this volume $\mathcal{S}(a)=\partial\mathcal{V}(a)$ when we compute $f_{Z^{*}_{\boldsymbol{V}}}$, which defines planes in term of the $Z$ variable but not necessarily in terms of the Gaussian variable $Y$.

Notice that $\mathbf{y_{n}}$ is now determined by the rest of the variables. Thus, we can just rewrite $\phi^{-1}\{\phi(a)-\phi(\mathbf{y_{1}})+\dots+\phi(\mathbf{y_{n-1}})\}=\mathbf{y_{n}}(a,\mathbf{y_{1}},\dots,\mathbf{y_{n-1}})={y_{n}}$. Therefore,

[TABLE]

The procedure described is summarized in Figure 6, where we want to compute the sum of two correlated lognormal random variables. The straight lines $z_{2}=z-z_{1}$ where integration needs to be done, in the bi-lognormal distribution space, are “bent”, together with the density itself, being transformed into curved lines $\phi(y_{2})=\phi(a)-\phi(y_{1})$ 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 multi-Gaussian density in $\mathcal{S}(a)\subset\mathbb{R}^{n-1}$. However, a highly suitable alternative could be to consider Monte-Carlo 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 cut-off value z, $P(Z^{*}_{\boldsymbol{V}}>\textit{{z}})$, by computing $f_{Z^{*}_{\boldsymbol{V}}}(z)$ in a discrete set of points in the interval $(0,\textit{{z}})$, 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 $n=4$, 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:

[TABLE]

with $\int_{\boldsymbol{V}}\rho(\textbf{u})d\textbf{u}=1$. Supposing we are considering an anamorphosis transformation given by an initial log-normal distribution with mean [math] and variance $1$, then 4.3 can be expressed as:

[TABLE]

where $Y^{*}(\textbf{u})$ 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 $P$), suitable for avoiding numerical integration in some of the formulations already presented. Firstly, we begin by introducing Hermite Polynomials, and secondly, we present a self-contained 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 $H_{n}(Y(\textbf{u}))$ are polynomials that have special properties related to the normal distribution. Hermite Polynomials are defined by Rodrigues’ formula:

[TABLE]

where $g^{(n)}(y)=\frac{d^{n}g(y)}{dy^{n}}$, $\displaystyle\frac{1}{\sqrt{n!}}$ is a normalization factor, $y$ is a Gaussian value, and $g(y)$ is the standard Gaussian probability distribution function (defined in Equation 3.2). The Hermite Polynomial $H_{n}(Y(\textbf{u}))$ is a polynomial of degree $n$. More specifically:

[TABLE]

Theorem 5.1**.**

$H_{n}$* is a polynomial of degree $n$ that satisfy the recurrence relations:*

[TABLE]

[TABLE]

Proof.

Calculating the successive derivatives of the Gaussian density:

[TABLE]

[TABLE]

[TABLE]

and, in general

[TABLE]

Noticing that definition (5.1) implies:

[TABLE]

and plugging this relation in (5.2) (and then dividing by $g(y)$) leads to:

[TABLE]

Dividing by $\sqrt{n!}$ leads to the first relation. From here, we obtain a way to generate the Hermite polynomials:

[TABLE]

On the other hand:

[TABLE]

which gives us a second recurrence relation:

[TABLE]

Comparing between these two recurrence relations, (5.3) and (5.4), follow. ∎

Theorem 5.2**.**

The Hermite polynomials are orthogonal by the Gaussian distribution:

[TABLE]

Proof.

Without loss of generality suppose $n\leq p$, and lets define

[TABLE]

By the substitution of $H_{p}(y)g(y)$ with $\displaystyle\frac{g^{(p)}(y)}{\sqrt{p!}}$ we get:

[TABLE]

Integrating by parts gives:

[TABLE]

By using the recurrence formula:

[TABLE]

Iterating in the same way $n$ times we get:

[TABLE]

We get two cases:

$\bullet$ if $n=p$, then $I_{n,p}=\displaystyle\frac{\sqrt{n!}}{\sqrt{n!}}\int_{-\infty}^{\infty}g(y)dy=1$

$\bullet$ if $n\neq p$, then I_{n,p}=\displaystyle\frac{\sqrt{n!}}{\sqrt{p!}}\Big{[}g^{(p-n-1)}(y)\Big{]}^{+\infty}_{-\infty}=\displaystyle\frac{\sqrt{n!}}{\sqrt{p!}}\Big{[}H_{(p-n-1)}(y)g(y)\Big{]}^{+\infty}_{-\infty}=0 ∎

If $Y\sim\mathcal{N}(0,1)$, the expressions above can be written in probabilistic terms:

[TABLE]

where the mean of $H_{n}(y)$ is:

[TABLE]

and the variance

[TABLE]

Theorem 5.3**.**

All function $\phi$ square integrable by $g$, say that $\int_{-\infty}^{\infty}\phi^{2}(y)g(y)dy$, can be expanded into Hermite polynomials:

[TABLE]

where the coefficients $\phi_{n}$ are given by

[TABLE]

Also we have

[TABLE]

We will omit this proof of the expansion into Hermite polynomials. The next two formulas can be obtained by using the Theorem (5.2).

Let $\lambda$ be any value. Let’s consider the Taylor expansion of the function $g(y-\lambda)$ in a point $y$:

[TABLE]

Also:

[TABLE]

Equating both sides gives:

[TABLE]

Proposition 5.4**.**

[TABLE]

Proof.

Consider the Fourier transform of the standard Gaussian density:

[TABLE]

The integral can be repeatedly differentiated with respect to $y$. Then we have:

[TABLE]

Using the definition (5.1) of the Hermite polynomials:

[TABLE]

and (5.10), Hermite polynomials can be rewritten as:

[TABLE]

Observe that the term $t^{n}$ in the integrand suggests that we consider:

[TABLE]

By (5.9), this leads to

[TABLE]

Notice that the exponential can be rewritten as:

[TABLE]

By equating both expression and noticing that $n=2q+p$, it leads to

[TABLE]

∎

The Hermite polynomials are not orthonormal for the nonstandard Gaussian distribution. We will need to compute expressions for shifts and stretches of the polynomials, in terms of non-shifted/non-stretched polynomials.

Proposition 5.5**.**

[TABLE]

with $\displaystyle C_{p}^{n}=\frac{p!}{n!(p-n)!}$.

Proof.

By using (5.7), one can write:

[TABLE]

Equating the coefficients associated with $\lambda$ in (5.14) and (5.15), and noticing that $p=n+k+2q$, leads to:

[TABLE]

By re-arranging (5.8):

[TABLE]

and plugging $k=p-n$ follows:

[TABLE]

∎

6. Application of Hermite Polynomials Expansions in a General Case

In this section we show an example, taken from [2], of the application of the results previously shown, when Hermite Polynomials are introduced for the computation of the conditional local mean, appropriate for implementation on any probability distributions presented in earth sciences databases.

The idea, in a first step, is to decompose any transfer function into Hermite Polynomials, in which an expansion up to a degree $P$ is used to give a close approximation to the shape of the function:

[TABLE]

Procedures to calculate the coefficients are derived from the properties already summarized in the previous section, and can be found in [17, 21]. Once the coefficients are computed, $\phi^{\prime}[Y(\textbf{u})]$ follows directly by using the recurrence relationships:

[TABLE]

Then, if we want to calculate $[Z^{*}(\textbf{u})]$, by plugging the Hermite expansion of the anamorphosis in equation 3.5 leads to:

[TABLE]

As we see, in order to assess the mean value in the raw variable, Hermite polynomials should be computed in a shifted/stretched version. Thus, we can use equation 5.13, which leads to:

[TABLE]

Since the Hermite polynomials have a zero mean except the first one (equation 5.5), all the terms in the sum vanish except the one associated with $n=0$, reducing the computation of the mean to:

[TABLE]

7. Conclusion

The aim of the study was to develop a methodology that allow us the quantification of the uncertainty in an arbitrary volume conditioned by sampling data, without the use of the traditional geostatistical simulation, which is time consuming and hard to manage, specially for large grid sizes. We have accomplished this objective, concluding with a tool suitable for any probabilistic distribution of the sample data.

For this, we have proceeded with an extension of the traditional multi-Gaussian model, allowing us to obtain the formulations needed, and making explicit the dependence of the uncertainty from the grades, the spatial correlation and conditioning data. Only a Kriging of the Gaussian values is needed to obtain conditional local means and variances.

Also we were able to obtain complete local distributions at any support, in an easy and straightforward way, although with the trade-off cost of resorting to an approximate version of the distribution when several number of points are considered for the change of support.

Future work should consider this last point as a root of study, improving such approximations, if it is the case; the search of simplified versions of the results obtained, by the use of polynomials; and extending the framework to multiple variables, taking advantage of correlated heterotopic data.