DeepAI

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

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.

• 1 publication
• 1 publication
09/20/2019

### The Locally Gaussian Partial Correlation

It is well known that the dependence structure for jointly Gaussian vari...
06/10/2020

### On the Conditional Expectation of Mean Shifted Gaussian Distributions

In this paper, we consider a property of univariate Gaussian distributio...
06/16/2021

### Estimating timber volume loss due to storm damage in Carinthia, Austria, using ALS/TLS and spatial regression models

A spatial regression model framework is presented to predict growing sto...
08/31/2020

### Direct Volume Rendering with Nonparametric Models of Uncertainty

We present a nonparametric statistical framework for the quantification,...
06/08/2022

### TreeFlow: Going beyond Tree-based Gaussian Probabilistic Regression

The tree-based ensembles are known for their outstanding performance for...
06/10/2019

### Multimodal Data Fusion of Non-Gaussian Spatial Fields in Sensor Networks

We develop a robust data fusion algorithm for field reconstruction of mu...
08/21/2020

### Testing for equality between conditional copulas given discretized conditioning events

Several procedures have been recently proposed to test the simplifying a...

## 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 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 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 , into a random function with a standard Gaussian univariate distribution. Then, under the working hypothesis that is spatially multi-Gaussian, the conditional distributions of are determined. Derivation of conditional expectation in a non-linear way, and estimates of recovery functions, are then obtained through inverse transformation from these 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 is well known, practice has shown that the mean estimator of this non-linear 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 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 , 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 -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

 P(Y(u)≤yp,Y(u′)≤yp,Y(u′′)≤yp)

requires the availability of a series of triplet values with the same geometric configuration as the triplet .

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 -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:

 (2.1) Var(ZVk|Z(uα)=z(uα),α∈{1,…,n(k)})=1|Vk|2∑i:ui∈Vkj:uj∈VkCov(Z∗(%ui),Z∗(uj)),

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 non-linear function of the covariance 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 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) FZ(z)=FY(y),

where and , commonly written as (Figure 1). Hence,

 Z(u)=ϕ[Y(u)]

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:

 ∫z−∞fZ(z)dz=∫ϕ(y)−∞fZ(z)dz=∫y−∞fY(y)dy,

leading to:

 fZ(z)=fY(y)ϕ′(y)=g(y)ϕ′(y),

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) g(y)=1√2π⋅e−y22,gσ2μ(y)=1σ√2π⋅e−(y−μ)22σ2.

Knowing the anamorphosis of a variable is equivalent to knowing its distribution.

The multi-Gaussian 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

 FY|data(y)=G(y−y∗SKσSK),

where , and where “data” represents the conditioning data . The conditional cdf of the original variable is then retrieved as:

 FZ|data(z) =FY|data(y) =G(y−y∗SKσSK)=Gσ2SKy∗SK(y),

where . This can be rewritten in terms of the probability densities:

 ∫ϕ(y)−∞fZ|data(z)dz = 1σSK⋅∫y−∞g(y−y∗SKσSK)dy = ∫y−∞gσ2SKy∗SK(y)dy,

from where follows:

 (3.3) fZ|data(z)⋅ϕ′(y)=gσ2SKy∗SK(y).

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) fZ|data(z)=1ϕ′(y)⋅gσ2SKy∗SK(y)=1ϕ′(ϕ−1(z))⋅gσ2SKy∗SK(ϕ−1(z)).

We will say that . Let us define hereafter.

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

 ∫ϕ(y)−∞z⋅fZ|data(z)dz =∫y−∞ϕ(y)⋅gσ2SKy∗SK(y)dy =1σSK⋅∫y−∞ϕ(y)⋅g(y−y∗SKσSK)dy (3.5) =∫y−∞ϕ(σSK⋅y+y∗SK)⋅g(y)dy.

By taking , we obtain :

 (3.6) [Z∗(u)]=∫∞−∞ϕ(y)⋅gσ2SKy∗SK(y)dy.

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 multi-Gaussian model:

 E[Zn(u)|data]=∫∞−∞ϕn(y)⋅gσ2SKy∗SK(y)dy,

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

 (3.7) Var[Z(u)|data] =E[Z2(u)|data]−E[Z(u)|data]2 =∫∞−∞ϕ2(y)⋅gσ2SKy∗SK(y)dy−[Z∗(u)]2.

We can extend the model into a -variate model in the following way: the joint cumulative function of (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 bi-variate case the model looks as follows:

 FZ(u0),Z(u1)|data(z0,z1) = P(Z(u0)≤z0,Z(u1)≤z1|data) = P(Y(u0)≤y0,Y(u1)≤y1|data) = FY(u0),Y(u1)|data(y0,y1) = GΣSKμSK(y0,y1),

with:

 μSK = (y∗SK(u0)y∗SK(u1)) ΣSK = (E{(Y∗(u0)−Y(u0))2}E{(Y∗(u0)−Y(u0))⋅(Y∗(u1)−Y(u1))}E{(Y∗(u0)−Y(u0))⋅(Y∗(u1)−Y(u1))}E{(Y∗(u1)−Y(u1))2}) = ⎛⎝σ2−kT0⋅C−1⋅k% 0σu0,u1−kT0⋅C−1⋅k1σu1,u0−kT1⋅C−1⋅k0σ2−kT1⋅C−1⋅k1⎞⎠ = (σ2SK(u0)σSK(u0,u1)σSK(u1,u0)σ2SK(u1)).

This can be rewritten in terms of probability densities:

 ∫z0−∞∫z1−∞fZ0,Z1|data(z0,z1)dz0dz1 = ∫ϕ(y0)−∞∫ϕ(y1)−∞fZ0,Z1|data(z0,z1)dz0dz1 = ∫y0−∞∫y1−∞gΣSKμSK(y0,y1)dy0dy1,

which leads to the following equality:

 (3.8) fZ0,Z1|data(ϕ(y0),ϕ(y1))⋅ϕ′(y0)⋅ϕ′(y1)=gΣSKμSK(y0,y1).

Hence, we can obtain the probability distribution by:

 (3.9) fZ0,Z1|data(z0,z1)=1ϕ′(y0)⋅1ϕ′(y1)⋅gΣSKμSK(y0,y1).

Multiplying both sides by in 3.8 and integrating with respect to and we obtain:

 ∫ϕ(y0)−∞∫ϕ(y1)−∞z0z1fZ0,Z1|data(z0,z1)dz0dz1=∫y0−∞∫y1−∞ϕ(y0)ϕ(y1)gΣSKμSK(y0,y1)dy0dy1.

By taking , , we obtain :

 E[Z(u0)∗Z(u1)∗]=∫∞−∞∫∞−∞ϕ(y0)ϕ(y1)gΣSKμSK(y0,y1)dy0dy1.

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 and :

 (3.10) Cov(Z∗(u0),Z∗(u1)) =E[Z(u0)∗Z(u1)∗]−[Z∗(u0)][Z∗(u1)] =∫∞−∞∫∞−∞ϕ(y0)ϕ(y1)gΣSKμSK(y0,y1)dy0dy1 −∫∞−∞ϕ(y)⋅gσ2SK(u0)y∗SK(u0)(y)d%y⋅∫∞−∞ϕ(y)⋅gσ2SK(u1)y∗SK(u1)(y)dy.

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 multi-Gaussian framework.

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

 μSK = ⎛⎜ ⎜⎝y∗SK(u1)⋮y∗SK(un)⎞⎟ ⎟⎠ ΣSK = ⎛⎜ ⎜ ⎜⎝σ2SK(u1)⋯σSK(u1,un)⋮⋱⋮σSK(un,u1)⋯σ2SK(un)⎞⎟ ⎟ ⎟⎠,

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

 FZ(u1),…,Z(un)|data(z1,…,zn) = P(Z(u1)≤z1,…,Z(un)≤zn|data) = P(Y(u1)≤y1,…,Y(un)≤yn|data) = FY(u1),…,Y(un)|data(y1,…,yn) = GΣSKμSK(y1,…,yn).

### 3.1. Examples

#### 3.1.1. The Log-Normal Case

If is a normal random variable with mean and variance , then the random variable is said to be log-normal with parameters and , or just . Thus, a random variable is log-normal if is normal.

The cumulative distribution function is given by:

 FZ(z)=G(ln(z)−μσ).

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

 G(ln(z)−μσ)=G(y),

Hence:

 (3.11) z=ϕ(y)=eμ⋅eσy,y=ϕ−1(z)=ln(z)−μσ.

From here we obtain:

 ϕ′(y)=σ⋅eμ⋅eσy=σϕ(y),

(notice that ) and that the conditioned local distribution is given by:

 fZ|data(z) = 1ϕ′(y)⋅gσ2SKy∗SK(y) = 1σ⋅eμ+σy⋅1σSK√2π⋅e−(y−y∗SK)22σ2SK = 1zσσSK√2π⋅e−(ln(z)−μ−σy∗SK)22σ2σ2SK.

This is a log-normal 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 log-normal prior distribution, preserves the log-normality.

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 bi-variate 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

 FZ(z)=1−e−λz,z≥0.

It follows that the anamorphosis function is:

 (3.12) z=ϕ(y)=−1λ⋅ln(1−G(y)).

From here we obtain:

 ϕ′(y)=1λ⋅11−G(y)⋅g(y),

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

).

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 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 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 points, located at positions, , each one with mass density , , we are interested in finding the probability distribution of the random variable:

 (4.1) Z∗V=ρ(u1)Z∗(u1)+ρ(% u1)Z∗(u2)+⋯+ρ(un)Z∗(un),

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 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:

 FZ∗V(z) = P(Z∗V≤z) = P(ZV≤z|data) = P(Z(u1)+⋯+Z(un)≤z|data) = ∫⋯∫\displaylimitsz1+⋯+zn≤zfZ1,…,Zn|data(z1,…,zn)dz1…dzn

This can be rewritten in term of the variable and the anamorphosis function, for a given value such that :

FZ∗V(z) = = ∫⋯∫\displaylimitsz1+⋯+zn≤zfZ1,…,Zn|data(z1,…,zn)dz1…dzn = FZ∗V(ϕ(a)).

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

 ddaFZ∗V(ϕ(a))=fZ∗V(ϕ(a))⋅ϕ′(a),

and on the other hand:

 ddaFZ∗V(ϕ(a)) =∫∞−∞…∫∞−∞gΣSKμSK(y1,…,yn−1,ϕ−1{ϕ(a)−ϕ(y1)+⋯+ϕ(yn−1)})ϕ′[ϕ−1{ϕ(a)−ϕ(y1)+⋯+ϕ(yn−1)}]⋅ϕ′(a)dyn−1…dy1.

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) fZ∗V(ϕ(a))=∫S(a)gΣSKμSK(y1,…,yn−1,yn)ϕ′(yn)dyn−1…dy1.

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 bi-lognormal 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 multi-Gaussian density in . 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, , 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) Z∗V=∫Vρ(u)Z∗(u)du,

with . Supposing we are considering an anamorphosis transformation given by an initial log-normal distribution with mean and variance , then 4.3 can be expressed as:

 Z∗V=∫Vρ(u)eY∗(u% )du,

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 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 are polynomials that have special properties related to the normal distribution. Hermite Polynomials are defined by Rodrigues’ formula:

 (5.1) Hn(y)=1√n!⋅g(n)(y)g(y)∀n≥0

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:

 H0(y) = 1 H1(y) = −y H2(y) = 1√2⋅(y2−1)
###### Theorem 5.1.

is a polynomial of degree that satisfy the recurrence relations:

 √(n+1)⋅Hn+1(y)+yHn(y)+√n⋅Hn−1(y)=0
 H′n(y)=−√n⋅Hn−1(y)
###### Proof.

Calculating the successive derivatives of the Gaussian density:

 g′(y)+yg(y)=0
 g′′(y)+yg′(y)+g(y)=0
 g(3)(y)+yg′′(y)+2g(y)=0

and, in general

 (5.2) g(n+1)(y)+yg(n)(y)+ng