U.S. Department of the Interior 

U.S. Geological Survey 

A Practical Primer on Geostatistics

By Ricardo A. Olea 

Open-File Report 2009–1103 

Version 1.4, December 2018
U.S. Department of the Interior 

RYAN K. ZINKE, Secretary 

U.S. Geological Survey 

James F. Reilly II, Director 

U.S. Geological Survey, Reston, Virginia 

First release: July 6, 2009 (ver. 1.0), online 

Revised: July 10, 2009 (ver. 1.0), online 

Revised: January 2010 (ver. 1.1), online 

Revised: July 2017 (ver. 1.2), online 

Revised: November 2017 (ver. 1.3), online 

Revised: December 2018 (ver. 1.4), online

For more information on the USGS—the Federal source for science about the Earth, 

its natural and living resources, natural hazards, and the environment—visit 

https://www.usgs.gov/ or call 1–888–ASK–USGS (1–888–275–8747). 

For an overview of USGS information products, including maps, imagery, and publications, 

visit https://store.usgs.gov/. 

Any use of trade, firm, or product names is for descriptive purposes only and does not imply 

endorsement by the U.S. Government. 

Although this information product, for the most part, is in the public domain, it also may 

contain copyrighted materials as noted in the text. Permission to reproduce copyrighted items 

must be secured from the copyright owner. 

Suggested citation: 

Olea, R.A., 2018, A practical primer on geostatistics (ver. 1.4, December 2018): U.S. Geological 

Survey Open-File Report 2009–1103, 346 p., https://doi.org/10.3133/ofr20091103. 

ISSN 2331-1257 (online) 
1 

A practical primer 

on geostatistics

By Ricardo A. Olea

U.S. Department of the Interior

U.S. Geological Survey

Open-File Report 2009–1103 

Version 1.4, December 2018
2 

CONTENTS

1. Introduction 3 

2. Univariate Statistics 14

3. Bivariate Statistics 53

4. SGeMS Software 67

5. Spatial Statistics 85

6. Transformations 98

7. Semivariogram 108

8. Simple Kriging 127

9. Ordinary Kriging 155

12.Cokriging 225

13.Crossvalidation 240

14.Critical Review 254

15.Sequential Gaussian

Simulation 266

16.Simulated Annealing 282

17.Filter Simulation 298

18.Reliability 318

19.Navigation Charts 333

10.Universal Kriging 192

11.Block Kriging 211

Selected Bibliography 340

Index 342
3 

1. INTRODUCTION
4 

THE CHALLENGE

• Most geological phenomena are extraordinarily

complex in their interrelationships and vast in their

geographical extension.

• Ordinarily, engineers and geoscientists are faced with

corporate or scientific requirements to properly

prepare geological models with measurements

involving a small fraction of the entire area or volume

of interest.

• Exact description of a system such as an oil reservoir

is neither feasible nor economically possible.

• The results are necessarily uncertain.

Note that the uncertainty is not an intrinsic property

of the systems; it is the result of incomplete knowledge 

by the observer.
5 

THE AIM OF GEOSTATISTICS

• The main objective of geostatistics

is the characterization of spatial

systems that are incompletely

known, systems that are common

in geology.

• A key difference from classical

statistics is that geostatistics uses

the sampling location of every

measurement. 

 Unless the measurements show spatial correlation, the 

application of geostatistics is pointless.

 Ordinarily the need for additional knowledge goes 

beyond a few points, which explains the display of results 

graphically as fishnet plots, block diagrams, and maps.

Two-dimensional example

Easting

Northing

Data

How much?

0.86

0.77

0.45

0.63

0.49

0.51

? 0.32
6 

GEOSTATISTICAL METHODS

Geostatistics is a collection of numerical techniques 

for the characterization of spatial attributes using 

primarily two tools: 

• probabilistic models, which are used for spatial data

in a manner similar to the way in which time-series

analysis characterizes temporal data, or

• pattern recognition techniques.

The probabilistic models are used as a way to

handle uncertainty in results away from sampling 

locations, making a radical departure from alternative 

approaches like inverse distance estimation methods.
7 

DIFFERENCES WITH TIME SERIES

On dealing with time-series 

analysis, users frequently 

concentrate their attention 

on extrapolations for 

making forecasts.

 Although users of 

geostatistics may be 

interested in extrapolation, 

the methods work at their 

best interpolating.

Time

Attribute

0 

Distance

Attribute

? 

 This simple difference has significant methodological 

implications.
8 

HISTORICAL REMARKS

• As a discipline, geostatistics was firmly established

in the 1960s by the French engineer Georges

Matheron, who was interested in the appraisal of ore

reserves in mining.

• Geostatistics did not develop overnight. Like other

disciplines, it has built on previous results, many of

which were formulated with different objectives in

various fields.
9 

PIONEERS 

Seminal ideas conceptually related to what today we 

call geostatistics or spatial statistics are found in the 

work of several pioneers, including:

• 1940s: A.N. Kolmogorov in turbulent flow and N.

Wiener in stochastic processing

• 1950s: D. Krige in mining

• 1960s: B. Mathern in forestry and L.S. Gandin in

meteorology
10

CALCULATIONS

Serious applications of geostatistics require the use of 

digital computers.

 Although for most geostatistical techniques 

rudimentary implementation from scratch is fairly 

straightforward, coding programs from scratch is 

recommended only as part of a practice that may help 

users to gain a better grasp of the formulations.
11

SOFTWARE 

For professional work, the reader should employ 

software packages that have been thoroughly tested to 

handle any sampling scheme, that run as efficiently as 

possible, and that offer graphic capabilities for the 

analysis and display of results.

 This primer employs primarily the package Stanford 

Geomodeling Software (SGeMS)—recently developed at 

the Energy Resources Engineering Department at 

Stanford University—as a way to show how to obtain 

results practically.

 This applied side of the primer should not be 

interpreted as the notes being a manual for the use of 

SGeMS. The main objective of the primer is to help the 

reader gain an understanding of the fundamental 

concepts and tools in geostatistics.
12

ORGANIZATION OF THE PRIMER

The chapters of greatest importance are those 

covering kriging and simulation. All other materials 

are peripheral and are included for better 

comprehension of these main geostatistical 

modeling tools.

 The choice of kriging versus simulation is often a 

big puzzle to the uninitiated, let alone the different 

variants of both of them. Chapters 14, 18, and 19 

are intended to shed light on those subjects.

 The critical aspect of assessing and modeling 

spatial correlation is covered in chapter 7.

 Chapters 2 and 3 review relevant concepts in 

classical statistics.
13

COURSE OBJECTIVES

This course offers stochastic solutions to common 

problems in the characterization of complex geological 

systems. 

At the end of the course, participants should have:

• an understanding of the theoretical foundations of

geostatistics;

• a good grasp of its possibilities and limitations; and

• reasonable familiarity with the SGeMS software,

thus opening the possibility of practically applying

geostatistics.
14

2. UNIVARIATE

STATISTICS
15

EVERYTHING AND A PIECE

In statistics, a population is the collection of all possible 

outcomes or individuals composing the complete system 

of interest; for example, all people in the United States.

 Populations may be hard or impossible to analyze 

exhaustively. In statistics, a limited collection of 

measurements is called a sample; for example, a Gallup 

Poll.

 Unfortunately, the term “sample” is employed with 

different meanings in geology and statistics.

 The statistical usage of the term “sample” is observed in 

what follows.

Geology Statistics

collection sample

sample observation
16

RANDOM VARIABLE

Variate

Frequency

A random variable or 

variate is a quantity that 

may take any of the values 

within a given set with 

specified relative 

frequencies.

 The concept is used in geostatistics extensively to 

characterize a population or convey the unknown value 

that an attribute may take at any spatiotemporal 

location.
17

DESCRIPTIVE ANALYSIS

• A sample of an attribute ordinarily comprises several

measurements, which are best understood when

organized in some way. This is an important aspect

of statistics.

• The number of measurements in a sample is the

sample size.

• There are multiple options to make the data more

intelligible. Some of these are more convenient than

others, depending on factors such as the sample size

and the ultimate objectives of the study.
18

SAMPLE VISUALIZATION
19

FREQUENCY TABLE

Given some numerical information, if the interval 

of variation of the data is divided into class 

intervals—customarily of the same lengths—and 

all observations are assigned to their 

corresponding classes, the result is a count of 

relative frequency of the classes.
20

UNCF EXAMPLE FREQUENCY TABLE

Class Count Frequency, %

1 7,680-7,710 1 1.43

2 7,710-7,740 1 1.43

3 7,740-7,770 1 1.43

4 7,770-7,800 2 2.86

5 7,800-7,830 5 7.14

6 7,830-7,860 6 8.57

7 7,860-7,890 10 14.29

8 7,890-7,920 11 15.71

9 7,920-7,950 13 18.57

10 7,950-7,980 7 10.00

11 7,980-8,010 7 10.00

12 8,010-8,040 4 5.71

13 8,040-8,070 2 2.86

Total 70 100.00

This example from a major oil company relates to depth in 

feet to an unconformity (UNCF) in an undisclosed area.

 It will be used as a common reference to graphically 

illustrate other definitions in this chapter.
21

HISTOGRAM

0

4

8

12

16

20

7695 7755 7815 7875 7935 7995 8055

Depth, ft

Frequency, percent

A histogram is a graphical representation of a 

frequency table.
22

CUMULATIVE FREQUENCY

Summaries based on frequency tables depend on the 

selection of the class interval and origin.

 Given a sample of size n, this 

drawback is eliminated by 

displaying each observation zi

versus the proportion of the 

sample that is not larger than zi. 

Each proportion is a multiple of 

100/n. The vertical axis is 

divided in n intervals and the 

data are displayed at the center 

of the corresponding interval.

 Customarily, the vertical axis is scaled so that data from 

a normal distribution (page 49) display as a straight line.
23

SUMMARY STATISTICS
24

SUMMARY STATISTICS

Summary statistics are complementary or alternative 

descriptors to histograms and cumulative 

distributions. 

 A statistic is a synoptic value that is calculated from 

a sample of observations, which is usually, but not 

necessarily, an estimator of some population 

parameter.

 Generally, summary statistics are subdivided into 

three categories:

• Measures of location or centrality

• Measures of spread or dispersion

• Measures of shape
25

MEASURES OF LOCATION

Measures of location give an idea about the central 

tendency of the data. They are: 

• mean

• median

• mode
26

MEAN 

The arithmetic mean or 

simply the mean, , of a 

sample of size n is the 

additive average of all 

the observations, zi

: 

 The mean of the UNCF 

sample is 7,912.2 ft. 

mˆ

. 1 ˆ

1

∑=

=

n

i

zi nm Mean
27

MEDIAN 

The median, Q2, of a 

sample is the value that 

evenly splits the number 

of observations zi into a 

lower half of smaller 

observations and an 

upper half of larger 

measurements. 

 If zi is sorted by 

increasing values, then

( )

( )  ( ) 







+⋅

=

+

+

5.0 , isif even.

, isif odd,

12/2/

2/1

2 nzz

nz

Q

n n

n

The median of the UNCF sample is 7,918 ft.

Median
28

MODE

The mode of a sample is 

the most probable or 

frequent value, or, 

equivalently, the center 

point of the class containing 

the most observations.

 For the UNCF sample, 

the center of the class with 

the most observations is 

7,935 ft.

0

4

8

12

16

20

7695 7755 7815 7875 7935 7995 8055

Depth, ft

Frequency, percent

Mode
29

ROBUSTNESS 

Robustness denotes the ability of statistical methods 

to work well not only under ideal conditions but also in 

the presence of data problems, mild to moderate 

departures from assumptions, or both.

 For example, in the presence of large errors, the 

median is a more robust statistic than the mean.
30

MEASURES OF SPREAD

Measures of spread provide an idea of the 

dispersion of the data. The most common measures 

are:

• variance

• standard deviation

• extreme values

• quantiles

• interquartile range
31

VARIANCE 

The variance, , is the average squared dispersion 

around the mean:

expressions that are commonly restricted to estimate 

variances of finite populations.

 When dealing with samples, the denominator is 

often changed to n - 1.

 Because this is a quadratic measure, it is less 

robust than most other measures of spread. 

The variance of the UNCF sample is 5,474 sq ft.

2 σˆ

( ) ,

1 1 2

1

2

1

2 2





 



 ∑ =−= ∑ ⋅− = =

mnz

n

mz

n

n

i

i

n

i

σ i


32

STANDARD DEVIATION 

The standard deviation is 

the positive square root 

of the variance.

 It has the advantage of 

being in the same units 

as the attribute.

 The standard deviation 

of the UNCF sample is 

74.5 ft.

σˆ σˆ σˆ σˆ

theorem, for any sample and , the proportion of data 

that deviates from the mean at least is at most : 

t > 1

( ) 2

1 Prop ˆ ˆ t

tmX σ ≤⋅≥−

mˆ t ⋅σˆ

According to Chebyshev’s

−2 t
33

EXTREME VALUES

The extreme values 

are the minimum and 

the maximum. 

 For the UNCF 

sample, the minimum 

value is 7,696 ft and 

the maximum value is 

8,059 ft.

This measure is not 

particularly robust, 

especially for small samples.
34

QUANTILES 

The idea of the median splitting the ranked sample 

into two equal-size halves can be generalized to any 

number of partitions with equal numbers of 

observations. The partition boundaries are called 

quantiles or fractiles. The names for the most 

common quantiles are:

• Median, for 2 partitions

• Quartiles, for 4 partitions

• Deciles, for 10 partitions

• Percentiles, for 100 partitions

The number of boundaries is always one less than 

the number of partitions.
35

UNCF QUARTILES

Q1 = 7,871.75 ft

Q2 = 7,918 ft

Q3 = 7,965.75 ft Q1

Q2

Q3

Q2 coincides with the median.
36

INTERQUARTILE RANGE

The interquartile range, 

iqr, is the difference 

between the upper and 

the lower quartiles

thus measuring the 

central spread of the 

data.

 For the UNCF sample, 

the interquartile range is 

94 ft.

, QQ 13 iqr −=

 The interquantile range is more robust than the variance 

but insensitive to values in the lower and upper tails.

iqr

Q1 Q3
37

OUTLIER 

 A practical rule of thumb is to 

regard any value deviating 

more than 1.5 times the 

interquartile range, iqr, from the 

median as a mild outlier and a 

value departing more than 3 

times iqr as an extreme outlier.

Outliers are values so markedly different from the rest of the 

sample that they raise the suspicion that they may be from a 

different population or that they may be in error, doubts that 

frequently are hard to clarify. In any sample, outliers are 

always few, if any.

 For the UNCF sample, all mild outliers seem to be legitimate 

values, while the extreme outlier of 8,240 ft is an error.

0

4

8

12

16

20

76957755781578757935

79958055

81158175

8235

Depth, ft

Frequency, percent

mild

3iqr 1.5iqr

mild
38

BOX-AND-WHISKER PLOT

The box-and-whisker plot is a 

simple graphical way to summarize 

several of the statistics:

• Minimum

• Quartiles

• Maximum

• Mean

Variations in this presentation

style abound. Extremes may

exclude outliers, in which case

the outliers are individually plotted 

as open circles. Extremes 

sometimes are replaced by the 

5th and 95th percentiles. 

7600

7700

7800

7900

8000

8100

Depth, ft

Max.

Min.

Q1

Q2

Q3

UNCF 

Mean +
39

MEASURES OF SHAPE

The most commonly used measures of shape in 

the distribution of values are:

• Coefficient of skewness

• Quartile skew coefficient

• Coefficient of kurtosis
40

COEFFICIENT OF SKEWNESS

If : , the left tail is 

longer;

 , the distribution is 

symmetric;

 , the right tail is 

longer.

( )

3

1

3

1

1

σ

∑=

−

=

n

i

i mz

n

B

B1 = 0

B1 > 0

B1 < 0

0

4

8

12

16

20

7695 7755 7815 7875 7935 7995 8055

Depth, ft

Frequency, percent

B1 = -0.38

The coefficient of skewness is a measure of asymmetry 

of the histogram. It is given by:

The UNCF coefficient of skewness is -0.38.
41

QUARTILE SKEW COEFFICIENT

The quartile skew coefficient serves the same purpose as 

the coefficient of skewness, but it is more robust, yet only 

sensitive to the central part of the distribution. Its 

definition is:

( ) ( )

iqr

QQQQ

qs 23 −−− 12

=

0

4

8

12

16

20

7695 7755 7815 7875 7935 7995 8055

Depth, ft

Frequency, percent

qs=0.02

If : , the left tail is longer;

 , the distribution is 

symmetric;

 , the right tail is longer.

The UNCF quartile skew coefficient is 0.02.

qs > 0

qs = 0

qs < 0
42

COEFFICIENT OF KURTOSIS

This statistic measures the concentration of values around 

the mean. Its definition is:

( )

4

1

4

2

1

σ

∑=

−

=

n

i

i mz

n B

If : , the distribution is more 

 peaked than the 

 Gaussian distribution;

 , it is as peaked as 

 the Gaussian;

 , it is less peaked 

 than the Gaussian.

B2 < 3

 The UNCF coefficient of kurtosis 

is 3.08.

0

0.1

0.2

0.3

0.4

0.5

-5 -2.5 0 2.5 5

x

f(x)

Gaussian distribution

B2 = 3

0

4

8

12

16

20

7695 7755 7815 7875 7935 7995 8055

Depth, ft

Frequency, percent

B2=3.08

B2 = 3

B2 > 3
43

MODELS
44

PROBABILITY 

• The probability of an impossible event is zero, Pr[A] = 0.

It is the lowest possible probability.

• The maximum probability is Pr[A] = 1, which denotes

certainty.

• When two events A and B cannot take place

simultaneously, Pr[A or B] = Pr[A] + Pr[B].

• Frequentists claim that Pr[A] = NA/N, where N is total

number of outcomes and NA the number of outcomes of

A. The outcomes can be counted theoretically or

experimentally.

• For others, a probability is a degree of belief in A, even

if no random process is involved nor a count is possible.

Probability is a measure of the likelihood that an event, A, 

may occur. It is commonly denoted by Pr[A].
45

JOINT AND CONDITIONAL PROBABILITIES

Let A and B be two random events. The probability of their 

joint (simultaneous) occurrence, , is

where is the conditional probability of A given that 

B has occurred.

If then the result of one event affects the 

probability of the other. Such types of outcomes are said to 

be statistically dependent. Otherwise, they are statistically 

independent. For example:

• Heads and tails in successive flips of a fair coin are

independent events.

• Being dealt a king from a deck of cards and having two

kings on the table are dependent events.

BA ]|Pr[

≠ [ ] ABA ,Pr]|Pr[

Pr[ and ⋅=∩= BBABABA ]Pr[]|Pr[]Pr[]

Pr[ and BA ]
46

Number of balls

Box Blue Red Total

#1 20 5 25

#2 12 18 30

32 23 55

BAYES’S THEOREM

This is a widely used relationship for the calculation of 

conditional probabilities. Given outcomes A and B, because

. :

[ ] [ ]

[ ] [ ] A

B

AB BA Pr

Pr

|Pr |Pr =

• If one only knows that there are 2 boxes, Pr[A] = 1⁄2 = 0.5.

• Now, if the table is available, Pr[B|A] = 20/25 = 0.8.

• Pr[B] = 32/55 = 0.59.

[ ] .05.0 69

.0 59

8.0 |Pr BA = =

Example

Suppose there are two boxes. 

A blue ball is drawn (event B). 

What is the probability the ball 

came from box #1 (event A)?

Hence:

Pr[ ] [ ] Pr ∩=∩ ABBA
47

PROBABILITY FUNCTIONS

Analytical functions approximating experimental 

fluctuations are the alternative to numerical descriptors 

and measures. They provide approximations of 

general conditions. Their drawback is the loss of fine 

detail in favor of simpler models.

 Models approximating histograms are called 

probability density functions.

 Variations in the parameters of a probability density 

function allow the generation of a family of distributions, 

sometimes with radically different shapes.

 The main subdivision is into discrete and continuous 

distributions, of which the binomial and normal 

distribution, respectively, are typical and common 

examples.
48

BINOMIAL DISTRIBUTION

This is the discrete 

probability density function, 

. , of the number of

successes in a series of 

independent (Bernoulli) 

trials, such as heads or 

tails in coin flipping. If the 

probability of success at 

every trial is p, the 

probability of x successes 

in n independent trials is

( ) ( ) ( ) xpp n

xnx

n npxf x xn ,,2,1,0,1 !!

!

,; =−  − = −

( ) ,; npxf

f (x ; 0.5, 12)

0

0.1

0.2

0.3

0.4

Number of successes

Probability

0 1 2 3 4 5 6 7 8 9 10 11 12

f (x ; 0.1, 12)

0

0.1

0.2

0.3

0.4

Number of successes

Probability

0 1 2 3 4 5 6 7 8 9 10 11 12
49

NORMAL DISTRIBUTION

The most versatile of all 

continuous models is the normal 

distribution, also known as the 

Gaussian distribution. Its 

parameters are μ and σ, which 

coincide with the mean and the 

standard deviation.

( )

( )

= ∞<<∞−

− −

xf e x

x

, 2

1

,;

2

2

2σ

μ

πσ

σμ

0

0.1

0.2

0.3

0.4

0.5

-5 -2.5 0 2.5 5

x

f(x; 0,1) If X = log(Y) is normally distributed, Y is said to follow 

a lognormal distribution. Lognormal distributions are 

positively defined and positively skewed.
50

PROBABILITY FROM MODELS

[ ] ( ) xfxX dx

x

∫

∞−

=≤

1

Prob 1

[ ] ( )

.0 88

Prob .1 17 Normal 1,0;

.1 17

=

=≤ ∫

∞−

X x dx

[ ] ( ) Prob 1 1 =≤ xFxX

Examples:

[ ] ( )

.0 88

.1Prob 17 .1 17

=

X =≤ F
51

EXPECTED VALUE 

Let X be a random variable having a probability distribution 

. and let be a function of x. The expected value of 

. is denoted by the operator and it is the

probability weighted average value of

If X is continuous, such as temperature:

and if it is discrete, like in coin flipping:

In the latter case, for the trivial example of if all 

values are equally probable, the expected value turns into

which is exactly the definition of the mean.

( ) xu

( ) xf ( ) xu

E[ ] ( ) xu

E[ ] ( ) ( ) ( ) , ∫

∞

∞−

= xfxuxu dx

E[ ] ( ) = ∑ ( ) ( ). x

xfxuxu

( ) xu .

( ) = xxu ,

[ ] ,

1 E = ∑

x

x

n

x
52

MOMENT

Moment is the name given to the expected value 

when the function of the random variable, if it exists, 

takes the form , where k is an integer larger

than zero, called the order.

 If a is the mean, then the moment is a central 

moment.

 The central moment of order 2 is the variance. For 

an equally probable discrete case,

( )k − ax

∑( ) =

−==

n

i

i mx

n

M

1

2 2

2

1

σ


53

3. BIVARIATE STATISTICS
54

TOOLS

Frequently there is interest in comparing two or more 

measurements made for the same object or site. 

Among the most common classical alternatives, we 

have:

• Scatterplot

• Correlation coefficient

• Regression

• Quantile-quantile plot

• Probability-probability plot

Some of these concepts can be generalized to deal

with more than two variables.
55

SCATTERPLOT 

A bivariate scatterplot is 

a Cartesian posting in 

which the abscissa and 

the ordinate are any two 

variables consistently 

measured for a series of 

objects.

Scatterplots are 

prepared for exploring 

or revealing form, 

direction, and strength 

of association between 

two attributes.

Mecklenburg Bay sea floor, Germany
56

COVARIANCE 

The covariance is a measure of joint variation. 

 Given two random variables X and Y with means and 

. , their covariance is the expected value: 

 The covariance estimator when using a sample of point 

measurements is: 

m x

mY

   ,YX

ECov XYX  mm Y

     











n

i

i

n

i

i

n

i

YX ii

yx

nn

yx

n 1 11

,

1

1

1

1

Cˆvo
57

CORRELATION COEFFICIENT

This coefficient is the number most 

commonly used to summarize bivariate 

comparisons. If and are the standard 

deviations for two variables, their correlation 

coefficient, ρ, is given by:

ρ

YX

YX

σσ

ρ ⋅ =

Cov ,

• It only makes sense to employ ρ for

assessing linear associations.

• varies continuously from -1 to 1:

1, perfect direct linear correlation 

0, no linear correlation

-1, perfectly inverse correlation

σ X σ Y

-1

-0.5

0 

1 
58

REGRESSION 

Regression is a method for establishing analytical 

dependency of one or more variables on another 

mainly to determine the degree of their dependency, 

to estimate values not included in the sample, or to 

summarize the sample.

• The variable in the abscissa is called the

regressor, independent, or explanatory variable.

• The variable in the ordinate is the regressed,

dependent, or response variable.

• In many studies, which variable goes into which

axis is an arbitrary decision, but the result is

different. Causality or the physics of the process

may help in resolving the indetermination.
59

REGRESSION MODEL

• The model is:

• is any continuous

function of x that is judiciously 

selected by the user. are 

unknown parameters.

• Term is a random variable 

accounting for the error.

• Having selected ,

parameters are calculated by 

minimizing total error, for which 

there are several methods.

( ) xfy ii i ;θ += ε

( ) xf i

;θ

( ) xf i

;θ

Explanatory variable (x)

Response variable (

y )

εi 

( ) xf i

;θ

θ

θ

ε

 Avoid applying the model outside the extreme values 

of the explanatory variable.
60

LINEAR REGRESSION

The simplest regression 

case is one in which the 

parameters are estimated 

by minimizing the mean 

square error, , associated 

with a linear model, 

 In this special situation, 

accounts for the proportion 

of variation accounted for by 

the regression.

 In the example, ρ = 0.94. 

Hence, in this case, the

∑=

=

n

i

i n

r

1

2 1 2 ε

linear regression explains 88% of the variation.

Mecklenburg Bay seafloor, Germany

ρ = 0.94

( ) 2 100 ⋅ ρ

2 ρ

2

r
61

NONLINEAR REGRESSION

• In theory, the

higher the

polynomial degree,

the better the fit.

• In practice, the

higher the

polynomial, the less

robust the solution.

• Overfitting may

capture noise and

not systematic

variation.

Linear

25.47 2 r =

Quadratic polynomial

20.56 2 r =

Cubic polynomial

20.33 2 r =

Sixth degree polynomial

19.98 2 r =
62

IMPLICATIONS

• allows prediction of one

variable when only the other

is known, but the inference

may be inaccurate,

particularly if the correlation

coefficient is low;

• means the variables are

related, but the association

may be caused by a

common link to a third

lurking variable, making the

relationship meaningless;

• does not necessarily imply

cause and effect.

High to good correlation: Countries with a population of 

more than 20 million in 1990

50

55

60

65

70

75

80

1 10 100 1000

Persons per television set

Life expectancy, year

ρ = -0.81
63

QUANTILE-QUANTILE PLOT

• A quantile-quantile or Q-Q plot is a

scatterplot based on ranked data.

• The pairing is independent of the

object or site where the

observations were taken. The first

pair of coordinates has the minimum

value for each attribute, the second

pair is made of the second smallest

readings, and so on until finishing

with the two maximum values.

Interpolations are necessary for

different size samples.

• Q-Q plots are sensitive to a shift and

scaling of the distributions.
64

STANDARDIZED VARIATE

If X is a random variable 

with mean μ and 

standard deviation σ, 

the standardized 

variate, Z, is the 

transformation:

σ

− μ = X Z

 A standardized variate always has a mean of zero and 

a variance of one.

A standardized Gaussian distribution is a called a 

standard normal distribution. It is often denoted by N( ).1,0

σ

− μ

=

X Z
65

PROBABILITY-PROBABILITY (P-P) PLOT

A P-P plot is another scatterplot prepared by extracting 

information from the cumulative distributions of two variates.

• If the variates are in different units, preliminary

standardization is necessary.

• For given thresholds, the axes show the cumulative

probabilities for the two distributions being compared.

Density

Velocity
66

Q-Q AND P-P PLOTS

Main use of the Q-Q and P-P plots is 

as a quick way to decide the degree of 

similarity between two distributions.

• If the distributions are the same, the

points align along the main diagonal.

• There is no statistic or level of

significance for evaluating the results.

• P-P plots are insensitive to shifting

and scaling, and the vertical scale is in

probability units.

• The Q-Q plot in the example illustrates

its potential for calling the user’s

attention to the fact that normal

distribution with the same parameters

can take negative values.

P-P

Q-Q

West Lyons field, Kansas
67

4. SGeMS SOFTWARE
68

GENERAL DESCRIPTION

The Stanford Geostatistical Modeling Software (SGeMS)

is a general-purpose, user-friendly, state-of-the-art 

geostatistical software package. 

The software code is in the public domain, 

downloadable from http://sgems.sourceforge.net/

(accessed in October 2017).

The code is in C++ and runs interactively under 

Windows.

A manual published by Cambridge University Press is 

listed in the selected bibliography.

SGeMS is an outgrowth of the Stanford Geostatistical 

Library (GSLIB), with which it shares several features. 

For details about GSLIB, please consult the book listed 

in the bibliography.
69

MAIN MENU
70

OBJECTS 

A. Data, when each record specifies location.

B. Grid, when locations are in a regular rectangular

(usually square) or parallelepipedal (usually cubic)

pattern. Location is implicit and specified in a

separate file. A grid file contains attribute values only.

Typically they are the result of some calculation, but

they can be used as input for subsequent processing.

Objects are files with numerical information that can 

be of two types:

 Evidently, no process is possible without prior 

loading of an object, which is done by: 

• clicking on <Objects> in the upper left corner of

main menu, or

• dragging the file icon into the Object area of the

main menu.
71

A. DATA FILES

• The first record is intended to be a header to identify

the data. It is often automatically used to label the

output.

• The second record must be the number of fields per

data entry.

• Starting with the third record, there should be one

record per attribute containing a short description. A

logical choice is to provide attribute type and unit of

measurement, if any.

Data may be up to three-dimensional and have to be 

coded into plain ASCII files. No special formatting is 

necessary, yet they must observe the following 

structure:
72

DATA FILES

o The entries must be strictly numeric. No letters or

special symbols are allowed.

o There must be exactly as many columns as

indicated in the second record. A special code

chosen by the user, such as -999, should be used

if there are missing values.

o It is required to provide at least two Cartesian

coordinates specifying location.

• If the file has v attributes, the numerical values for the

sample start at record v + 3, one for each site, under

the following restrictions:

• The number of records is unlimited.
73

EXAMPLE OF DATA FILE

UNCF unconformity

4

Identification

Easting, m

Northing, m

Depth, ft

1 32000 87015 7888

2 38600 87030 8020

3 44000 86400 7949



6200 85050 87000 8003

6202 49400 84000 8032
74

DATA LOADING

• Click on <Objects>.

• Click on <Load Object> at new menu.

• Search for the file of interest.

• Highlight file name and click on <OK>.

Simply:

• Choosing <Point Type> at <Select object type>.

• Clicking on <Next>.

The top of the file will be displayed. Finish by:

• Drag the file icon into the object area of the main

menu

or:
75

1. COMPLETING DATA LOADING

• Type an object name at

<Pointset name>.

• Indicate columns for

easting and northing

and also depth in threedimensional studies.

• Specify missing data

code, if any.

• Click on <Finish>.

The loading has been

completed and the data 

file will show in the object 

panel of the main menu.
76

B. GRID LAYOUT

Results from modeling of 

spatial attributes come in 

the form of arrays of 

values at regular 

intervals.

Coordinates for the nodes 

are not printed; they are 

implicit in the regular 

arrangement. 1 2 

134 135
77

IMPLICIT ARRANGEMENT

• In three-dimensional modeling, layer values go from

the lowermost layer to the uppermost layer.

• Within a layer, the first value to go into the file is the

one in the lower left corner.

• Then come the rest of the values for the lowermost

row, from left to right.

• Within a layer, last value to enter the file is the one

in the extreme upper right corner.

• Specification of grid geometry goes into a separate

file.
78

GRID SPECIFICATION FILE

Remember that grids 

must be specified 

before attempting 

any processing that 

will result in the 

generation of a grid, 

such as estimation 

or simulation.

For specifying a 

grid, click <Object> 

and then select 

<New Cartesian 

Grid>. 
79

EXAMPLE OF GRID

To minimize storage space, SGeMS stores grid values in 

binary form. Grids can only be displayed in graphical form.

 Numerical display of a grid requires saving it in ASCII 

format and further viewing it with a text editor. To do such 

saving:

• Click on <Object>.

• Select <Save Object>.

• Name the file.

• Select place for storage

at <Look in>.

• Make sure to select

GSLIB as the format.

• Click on <Save>.
80

VIEW OF GRID WITH TEXT EDITOR

UNCFgrid (211x281x1)

2 

OK

OK_krig_var

7773.04 2075.6 

7766.52 1902.38 

7763.12 1745.86 

7759.76 1595.19 

7756.46 1450.97

 . 

. 

.

8017.93 2229.09 

8018.47 2313.04 

8018.98 2404.05 

8019.46 2501.91 

8019.92 2606.35 

SGeMS automatically 

creates the first record, 

which includes the 

number of columns, 

rows, and layers.
81

OBJECT SAVE

In SGeMS, all objects in a 

session make up a project.

 The user can save all 

objects in a project for 

archiving or further use. 

To do that:

• Click on <File> at the upper right corner of the

SGeMS screen;

• Click on <Save Project>. This is all that it is

necessary to do for an existing project.

• If the project is new, navigate directories to locate a

storage folder, select name for the project file, and

click on <Save>.
82

HISTOGRAM PREPARATION

• Click on <Data Analysis> at main menu.

• In the new menu, click on <Histogram>.

• Select the object.

• Choose the property.

• Set the number of classes.

• If there is a need to change the axis extreme values,

click on <Display Option> and enter the settings.

As a first practical application of SGeMS, let us prepare 

a histogram of the object previously loaded. For that:
83

HISTOGRAM DISPLAY
84

SAVING IMAGE

• Click on <Save as

Image>.

• Select a folder.

• Give a name to the

file.

• Select type.

• Decide on grid and

summary statistics.

• Click on <Save>.

The histogram graphical file is now at the specified folder.
85

5. SPATIAL STATISTICS
86

RANDOM FUNCTIONS

A random function is a collection of random variables, 

one per site of interest in the sampling space.

 A realization is the set of values that arises after 

obtaining one outcome for every distribution.

 Geostatistics relies heavily on random functions to 

model uncertainty.

s

t
87

MAIN DIFFERENCES BETWEEN CLASSICAL 

GEOSTATISTICS AND CLASSICAL STATISTICS 

• Geographical location is an integral part of any sample.

• Measurements have an associated volume called support.

For example, in the case of density, the support is the

specimen volume; a depth measurement has a point support.

• Geostatistics does not assume that the variables are

independent and identically distributed.

• A spatial sample is regarded as a single realization of a

random function, instead of multiple outcomes of a single

random variable.

• Most geostatistical formulations do not require any particular

probability distribution, although some work better under

normality conditions.

• The most efficient forms of spatial sampling are those

following a regular pattern—cubic or hexagonal, for example.
88

ERGODIC ASSUMPTION

We have seen that stochastic methods regard data as 

one of an infinite number of possible realizations.

 Ordinarily, one needs multiple realizations to infer 

properties of the complete ensemble.

 The ergodic assumption basically states that just 

one realization is indeed sufficient to make reliable 

assessments about ensemble properties.

It is an assumption because it is not possible to test it.
89

STATIONARITY 

This is another fundamental assumption, presuming 

invariance through space of properties of the random 

function.

 The most general case assumes that the joint 

probability distribution of any number of random 

variables is the same.

 Most geostatistical formulations require 

independence of location only for moments and just 

up to order two.
90

BIAS 

Bias in sampling denotes preference in taking the 

measurements.

 Bias in an estimator implies that the calculations are 

preferentially loaded relative to the true value, either 

systematically too high or too low.

• When undetected or not compensated, bias induces

erroneous results.

• The opposite of the quality of being biased is to be

unbiased.
91

PREFERENTIAL SAMPLING

Although estimation and simulation methods are 

robust to clustered preferential sampling, parameters 

that need to be inferred from the same sample prior 

to estimation or simulation can be seriously distorted, 

especially for small samples.

 The solution to preferential sampling is preparation 

of a compensated sample to eliminate the clustering, 

for which there are several methods.

 Declustering is important for the inference of global 

parameters, such as any of those associated with the 

histogram.
92

EXAMPLE OF SAMPLE WITH 

PREFERENTIAL CLUSTERING
93

DECLUSTERING (1) 

Detect presence of clusters by preparing a cumulative 

distribution of distance to nearest neighbor.
94

DECLUSTERING (2)

Decompose the clustering. 

S1 S2
95

DECLUSTERING (3) 

Preferential sampling shows as poor overlapping 

between the two histograms.

Histograms of the attribute by distance class
96

DECLUSTERING (4)

a few observations (S3) from the 

clustered subset (S2).

• Transfer observations from

(S2) by decreasing distance to

nearest neighbor in S4.

• Stop transferring points when

the distribution of distances for

the original subset S1 is about

the same as that for the

transferred observations (S3)

within S4.

S1 

S3 

One possibility is to obtain the declustered subset (S4) by 

expanding the subset without clusters (S1) by transferring
97

(S3) 

(S1) 

(S1) 

(S3) 

(S1) 

Best results are obtained by randomly repeating the 

drawing from S2 to generate multiple declustered sets. 

The final answer is obtained by taking median values 

from the cumulative distributions of all declustered sets.

POSTING OF A 

DECLUSTERED SAMPLE (S4) 
98

6. TRANSFORMATIONS
99

COMPARING DISTRIBUTIONS

distribution becomes 

mandatory. 

The simplest way to 

compare a sample, , 

to an hypothesized 

distribution, , is to 

superimpose their 

cumulative distributions.

 By transforming the 

data, one can force the 

data to follow almost any 

other distribution.

( ) xF

n ( ) xF

The distribution followed by some data and the one that may 

be required by a method may be two different distributions.

 If a method requires a particular distribution, such
100

KOLMOGOROV-SMIRNOV TEST

The maximum discrepancy, D, is the measure of similarity in 

the Kolmogorov-Smirnov test.

 The significance of the discrepancy depends on the sample 

size and on whether the parameters of the hypothesized 

distribution are known. A statistic called p-value, commonly 

read from tables, is used to decide the significance of the 

discrepancy. Customarily, if the p-value is above 0.1 or 0.05, 

the equality of the distributions is accepted.

 Standard tables available to run the significance of the 

discrepancy are not valid in geostatistics because the tables 

do not take into account spatial dependence.

 It has been found, however, that the standard tables may 

be a limiting case for the situation considering spatial 

correlation. In that case, each time two distributions pass the 

Kolmogorov-Smirnov test without considering spatial 

correlation, they also pass it considering that correlation. 

The same is not true in case of rejection.
101

WEST LYONS FIELD, KANSAS

Remember that the cumulative 

probability of any normal 

distribution follows a straight 

line in normal probability scale.

 At 10 ft, the thickness 

distribution of Lyons field 

reaches a maximum 

discrepancy of 0.09 against a 

normal distribution with the

Because the pu-value is already above the threshold of

0.05–0.1, it is not necessary to go into the complications of 

considering spatial correlation. Without further ado, it can 

be accepted that the thickness is normally distributed.

same mean and standard deviation as the sample. The 

p-value without considering spatial correlation, pu, is 0.14.

D = 0.090

pu = 0.14
102

NORMAL SCORES

In geostatistics, as in classical statistics, the most 

versatile distribution is the normal distribution.

 If the sample strongly deviates from normality, it may 

be convenient or required to transform the data to follow 

a normal distribution, so as to properly apply a specific 

method.

 The transformation is fairly straightforward. It is 

performed by assigning to the cumulative probability of 

every observation the value of the standard normal 

distribution for the same cumulative probability.
103

NORMAL SCORES

For the value zi of 

the attribute, the 

normal score is 1.45.

 The transformation 

also can be used 

backwards to bring 

calculated values in 

the normal score 

space to the original 

attribute space.
104

LOGNORMAL TRANSFORMATION

• This is a special case of normal score transformation.

The transformed value is simply the logarithm of the

original observation.

• If the original distribution is exactly lognormal, the

distribution of the transformed values will be exactly

normal. Because this is never the case, the

distribution for a transformed variable is only

approximately normal.

• Most minerals and chemical elements have

distributions of concentration close to lognormal,

which explains the great attention paid to this

distribution since the early days of geostatistics.
105

INDICATORS

Given a continuous random function , where s

denotes geographical location, its indicator is

the binary transformation:

The potential of indicators lies in the possibility to 

handle imprecise data and the fact that the expected 

value of an indicator is equal to the cumulative 

probability for the threshold value:

Z( ) s

( ) s,uI

( ) ( )

 ( ) 



≤

>

=

uZ

uZ uI

s

s

s

,1 if

,0 if

,

[ ] ( )( ) [ ( ) ( )]

( )( ) nuF

nuIE nuZ

|,

|, Prob |

s

s s

=

= ≤
106

INDICATOR CODING FOR DIFFERENT 

TYPES OF CONTINUOUS DATACourtesy of A. Journel, Stanford University
107

DISCRETE INDICATORS

The indicator of a categorical attribute is:

 In this case, the indicator is now equal to the 

probability of belonging to class k. 

 In this case, the indicator transformation offers the 

possibility of digital processing for information that is 

commonly qualitative.

( ) u,kI

( ) 



 ∈

=

,0 otherwise

,1 if

,

ku

u kI

[ ] u |, ( )( ) nkIE [ ] ∈= |Prob ( ) nku
108

7. SEMIVARIOGRAM
109

THE OBJECTIVE

The objective of structural 

analysis is to capture the 

style of the fluctuations 

through an analytical 

function.
110

CONTINUITY

Spatial continuity involves the concept that small 

values of an attribute are in geographical proximity 

to other small values, while high values are close to 

other high values. Transitions are gradual.

 Assessment of covariance or its close equivalent, 

the semivariogram, has been the classical way in 

geostatistics to measure spatial correlation.
111

RELEVANCE 

• The semivariogram or the covariance are

necessary in the formulation and application of

most estimation and simulation methods.

• Semivariogram analysis leads to some standalone applications, such as sampling design.
112

DEFINITIONS

Given a spatial distance h involving magnitude and 

direction, the semivariogram is:

and the spatial covariance is:

both with a constant mean 

 Different from the classical covariance (page 56), the 

spatial covariance is grouped by distance and applies to 

the same attribute. Therefore, a more descriptive yet 

seldom used term for the covariance is autocovariance.

( ) h = [ ( ) ( ) ZZVar +− hss ] 2

1

γ

γ ( ) h

[ ] ( ) [ ] ZZm ( ) +== hss .EE

Cov( ) ( )( ) ( )( E[ ZmZ hssh −+−= m)],

m
113

EQUIVALENCE 

If the mean is constant and 

the covariance is 

independent of location, 

then always

which makes it immaterial 

which one to use. 

 In general, what is 

estimated is the 

semivariogram because its 

estimation does not require 

knowledge of the mean. 

γ ( ) ( ) ( ) h Cov 0 −= Cov h ,

 Lag a is the range and the semivariogram asymptote 

is called the sill, which is equal to the variance Cov(0).

Lag

Second order moment

Semivariogram

Covariance

Cov(0)

a 

0 
114

SEMIVARIOGRAM ESTIMATOR

 : estimated or experimental semivariogram; 

 : lag;

 : number of pairs units apart;

 : observation at site . 

Always remember that for the estimator to be valid, the 

mean must be constant; hence, there must be no trend.

For any direction, it is recommended that and 

the estimation be restricted to h no more than half the 

extension of the sampling domain in such direction.

( ) ( ) [ ( ) ( )] ( )

∑=

= +−

h

hss

h

h

n

i

i i zz

n 1

2

2

1

γˆ

h

γˆ( ) h

n( ) h h

( ) z si si

n( ) h ≥ 30
115

TRIVIAL EXAMPLE

0

1

2

3

4

5

0 1 2 3 4 5

Lag

Semivariogram

∆ 5.0 ∆ 0.1 ∆ 5.1

si ( )i z s 2 ∆0

2 ∆ 5.02 ∆ 0.12 ∆ 5.12 ∆ 0.22 ∆ 5.2

2 ∆ 0.32 ∆ 5.32 ∆ 0.4

2 ∆ 5.4

7.0 3.2 0.0 1.21 3.24 10.89 22.09 24.01 18.49 16.81 12.25 3.76

7.5 4.3 0.0 0.49 4.84 12.96 14.44 10.24 9.00 5.76 2.25

8.0 5.0 0.0 2.25 8.41 9.61 6.25 5.29 2.89 0.64

8.5 6.5 0.0 1.96 2.36 1.0 0.64 0.04 0.49

9.0 7.9 0.0 0.04 0.16 0.36 1.44 4.41

9.5 8.1 0.0 0.36 0.64 1.96 5.29

10.0 7.5 0.0 0.04 0.64 2.89

10.5 7.3 0.0 0.36 2.25

11.0 6.7 0.0 0.81

11.5 5.8 0.0

∑∆2

h 0.0 7.52 22.74 38.97 50.15 43.99 30.87 23.21 14.50 6.76

n( ) h 10 9 8 7 6 5 4 3 2 1

λ( ) h 0.0 0.42 1.42 2.78 4.18 4.40 3.85 3.86 3.62 3.38
116

CASE OF IRREGULAR PATTERNS

If the observations, are not regularly spaced, they are 

grouped into distance classes of equal radial thickness, 

which customarily is set equal to twice 

z( ) s ,

si

sj

Easting

Northing

Each class contributes one value to the experimental 

semivariogram,

 Commonly used 

parameters are:

 : lag tolerance 

 : angular tolerance

 : bandwidth

 only operates when the 

cone is wider than twice 

Otherwise, prevails. 

ht

δ

bt

bt

. bt

δ

γ ( ) h . 

. ht
117

MOST COMMON PERMISSIBLE MODELS

Spherical

Exponential

Gaussian

 































 

 ah

ah

a

h

a

h

h

, if

, if

2

1

2

3

3

C

C



  ]1[

2

3 















a

h

 h eC

  ]1[

3 















a

h

 h eC

A sum of permissible models is also permissible. 

Permissible semivariogram models are a class of models 

guaranteeing a unique solution to an estimation system of 

equations and a nonnegative estimation variance.

a
118

ASSESSMENT OF SPATIAL 

CORRELATION

Spatially correlated attributes have 

gradually increasing semivariograms that 

reach the range at a distance clearly 

greater than zero. 

If the semivariogram does not seem to 

be zero for a lag close to zero, the value 

of convergence is called nugget effect. 

Yet, 

Lack of spatial correlation is revealed by 

a semivariogram that remains constant, 

following what is called a pure nuggeteffect model, which is permissible.

Mere assessment of spatial correlation does not require 

semivariogram modeling.

Semivariogram

Lag

Proof of spatial correlation

Nugget

effect {

Lag

Lack of spatial correlation

Semivariogram

 ( ) = .00
119

ESTIMATION BASIC STEPS

A. Make sure that the data are error free.

B. Correct preferential clustering, if any.

C. Perform transformations if required or desired.

D. Estimate semivariograms along at least 3 different

directions.

E. Detect trend. In the presence of trend, the

experimental semivariogram is an artifact that

increases monotonically, reaching values orders of

magnitude higher than the sample variance.

F. If necessary, model semivariogram along the trendfree direction.

 When the attribute is isotropic, parameters do not change 

with h. A semivariogram model prepared disregarding 

azimuthal variation is said to be omnidirectional.
120

A. DATA LOADING

If data have not yet been loaded into SGeMS, load 

them as explained in chapter 4 (page 74).
121

B. MODELING WITH SGeMS

First click on <Data Analysis> and then on <Variogram>.

See page 187 for more information about this dataset.
122

C. SELECT FOLDER AND ATTRIBUTE

• Select <Head Property> equal to <Tail Property>.

• After making selections, click on <Next>.
123

D. DIRECTIONS AND LAGS

For grids, the azimuth is specified by directional cosines. 

For example, north in a horizontal place is (0, 1, 0).

 For datasets, azimuth is specified as clockwise degrees 

from the north; so, for example, NE is 45o. 

 Upon 

specifying 

lags and 

directions, 

click on 

<Next>.
124

E – F. TREND DETECTION AND MODELING 

• Hide the main

menu screen;

• expand sideways

semivariogram

menu;

• click on

<Window>, and

then click on

<Tile>;

• pick best fitting

model parameters

by trial and error.

Here there is no trend because the semivariogram tends

to stabilize at a value close to the sample variance. See 

pages 202–203 for modeling of a nonstationary dataset.
125

SAVING SEMIVARIOGRAMS

• Click on the <File> option of semivariogram modeling.

• Pick an option, for example, <Export Plot as Images>.

• Check <Show Grid> if it is desired to display it.

• Click on the three dots <...>.

• Select storage folder.

• Give generic name to semivariograms.

• Click on <Save> in both menus.
126

OTHER DIRECT APPLICATIONS

Direct applications of the semivariogram, not necessarily 

requiring to go through modeling, include:

• Detection of anisotropy by

finding out if the

semivariogram is different

along different directions.

• Comparative studies of

origins in the processes

behind attributes based on

the premise that different

geneses lead to different

semivariograms.

• Sampling design. In any direction, sampling interval

should not exceed half the range.

Lag

Semivariogram

Direction 1

Direction 2

Direction 3

Direction 4
127

8. SIMPLE KRIGING
128

KRIGING

Kriging is a form of generalized 

linear regression for the 

formulation of an optimal spatial 

estimator in a minimum meansquare-error sense.

Given some observations at 

locations si, kriging provides an 

estimate and a standard error at 

locations, such as s0, typically 

not considered in the sampling.

 In contrast to classical linear regression, kriging takes 

into account observation volume and stochastic 

dependence among data.

 The method works best inside the convex hull determined 

by the peripheral data.

s1

s3

s4

s2

s6

s5

s0

s7

s8
129

VARIANTS 

There are several forms of kriging, each making 

different assumptions to adjust to different styles of 

attribute fluctuations.

• Simple kriging

• Ordinary kriging

• Universal kriging

• Lognormal kriging

• Indicator kriging

The last two types are simply any of the first three 

forms of kriging applied to transformed data as 

explained in chapter 6.
130

SIMPLE KRIGING

Simple kriging is the estimation of

( ) ( )( ),

1

0

*

∑=

−⋅+=

n

i

SK s mz λ sii mz

where:

 is the estimator at geographical location

 denotes the location of measurement ;

 is the number of observations to consider;

 is the mean of 

is a real number, commonly referred to as weight. 

n

λi

si

( ) 0

* zSK s

m

; s0

Z( ) s ;

i
131

NORMAL EQUATIONS

The originality of kriging is the use of weights, such 

that they minimize the estimation error variance. Such 

weights turn out to be the solution of the following system 

of equations, called the normal system of equations: 

( ) ( )

( ) ( )

( ) ( )



















=

=

=

∑

∑

∑

=

=

=

k

k

i

i ki

k

i

i i

k

i

i i

ss ss

ss ss

ss ss

Cov , Cov ,

Cov , Cov ,

Cov , Cov ,

0

1

20

1

2

10

1

1

λ

λ

λ



, λi
132

DERIVATION OF NORMAL EQUATIONS

The steps are:

A. Make necessary assumptions.

B. Calculate error variance.

C. Express error variance in terms of covariances.

D. Find the weights that minimize the error variance.

E. Calculate the minimum error variance.
133

A. ASSUMPTIONS

1. The sample is a partial realization of a random

function , where denotes spatial location.

2. is second order stationary, which implies:

 where denotes expected value, is a scalar 

constant, and is a vectorial distance.

3. Unique to simple kriging is the assumption that the

mean is not only constant but known.

Z( ) s s

E[ ] ( ) s = mZ

E[ ( )( ) ( )(ZmZ hss m)] =−+− Cov( ) h ,

E[ ]⋅

h

m

Z( ) s

 These assumptions make simple kriging the most 

restricted form of kriging.

 The weakness of this simplicity is the same as with any 

simple model: limitation in its applicability or suboptimal 

results if its use is forced.
134

B. ESTIMATION VARIANCE

By definition, the estimation variance is:

 

We know neither the random function nor the 

estimator , so we cannot directly calculate the 

variance of the difference.

 The way to solve this conundrum is to transform the 

expression for the variance to have it in terms of 

quantities that we know.

( ) ( ) ( ) [ ] 0 0

*

0

2 σ s = Var SK − ZZ ss

( ) Z s0

( ) 0

* ZSK s

. 
135

C. ALGEBRAIC MANIPULATIONS (1)

The residual of a random function is the difference 

between the random function and its expected value. In 

this case:

Expressing this estimation variance in terms of the 

residuals, we have:

if one defines .

      E    ssss  mZZZY

       























     

k

i

ii

k

i

ii

YY Y

0

0

1

0

2

 s Var  ss Var  s

Y  s

0

 1
136

ALGEBRAIC MANIPULATIONS (2)

From the properties of the variance (pages 31 and 52):

Expanding the squares and introducing the expectations 

inside the summations:

( ) E ( ) E ( ) .

2

0

2

0

0

2

















 



 −





























= ∑ ∑ = =

k

i

ii

k

i

σ s λY sii λY s

( ) [ ] ( ) ( )

[ ] ( ) [ ] ( ).EE

E

0 0

0 0

0

2

∑∑

∑∑

= =

= =

−

=

k

i

k

j

ji i j

k

i

k

j

ji ji

YY

YY

ss

s ss

λλ

σ λλ
137

ALGEBRAIC MANIPULATIONS (3)

Grouping the terms by coefficients:

Backtransforming to the random function and expanding: 

Because of the cancelation of terms,

( ) {E[ ] ( ) ( ) [ ] ( ) [ ] ( )}.EE

0 0

0

2

∑∑= =

= −

k

i

k

j

σ s λλ ji ji i YYYY ssss j

( ) {E[ ] ( ) ( ) [ ] ( ) [ ] ( )},EE

0 0

0

2

i j

k

i

k

j

s = ∑∑ ji ji − ZZZZ ssss = =

σ λλ

( ) { [ ( ) ( ) ( ) ( )]

E[ ] ( ) [ ] ( ) [ ] ( ) [ ] ( )}.EEE

E

2

0 0

2

0

2

i j i j

k

i

k

j

ji i j ji

ZZZmZmm

mZm ZZmZ

s sss

s ssss

+− + −

= ∑∑ +−− = =

σ λλ
138

ALGEBRAIC MANIPULATIONS (4)

which is the same as:

Finally, remembering that 

which is an expression in terms of covariance, a 

measurable property, hence a useful expression.

   .,Cov

0 0

0

2

ji

k

i

k

j

s  jiss

 

  

0

 :1

     

  



 





 

k

i

k

j

ji ji

k

i

i i

1 1

1

0 00 0

2

Cov ,,

Cov 2, Cov ,

ss

sss ss



 
139

D. OPTIMAL WEIGHTS (1)

The expression for is valid for any weights. Among 

the infinite weight combinations, we are interested in those 

that minimize the estimation variance associated with the 

estimator

The expression for is quadratic in the unknowns, 

the weights. Under these circumstances, according to an 

operations research theorem, the necessary and sufficient 

condition to have a unique global nonnegative minimum is 

that the quadratic term must be larger than or equal to zero. 

This implies that the covariance function must be such that 

any combination, like the one in the general expression for 

. cannot be negative. In mathematics, such is the 

property of being positive definite. This is the condition a 

covariance model has to satisfy to be permissible. 

( ) 0

2 σ s

( ) 0

2 σ s

( ). 0

* ZSK s

( ), 0

2 σ s
140

OPTIMAL WEIGHTS (2)

The optimal weights are those that make zero all partial 

derivatives of with respect to the weights:

for all locations .

Cancellation of factor 2 and rearrangement give:

( ) 2Cov( ) 2, Cov( ) 0,

1

0

0

2

−= + =

∂

∂ ∑=

k

i

j i ji

j

ss ss

s λ

λ

σ

( ) 0

2 σ s

= ,,2,1, kjj

( ) ( )

( ) ( )

( ) ( )



















=

=

=

∑

∑

∑

=

=

=

k

k

i

i ki

k

i

i i

k

i

i i

ss ss

ss ss

ss ss

Cov , Cov ,

Cov , Cov ,

Cov , Cov ,

0

1

20

1

2

10

1

1

λ

λ

λ


141

E. OPTIMAL ESTIMATION VARIANCE

Multiplying each optimality condition by and adding 

them together, we have a new optimality relationship:

Introducing under the second summation and 

replacing the resulting double summation in the general 

expression for :

Cov( ) , Cov( )., 0

1 1 1

i

k

i

k

i

i

k

j

∑ ∑ ji ss ji ∑ ss = = =

λλ = λ

λi

λi

( ) 0

2 σ s

( ) ( ) ( ) Cov , Cov .,

1

0 00 0

2

∑=

= −

k

i

σ SK s ss λi ss i
142

SUMMARY

Simple kriging normal equations for optimal weights : 

Simple kriging estimation:

Simple kriging mean-square-error variance:

( ) ∑ ( )( ) =

+= −

n

i

SK ii mzmz

1

0

* s λ s

( ) ( )

( ) ( )

( ) ( )



















=

=

=

∑

∑

∑

=

=

=

k

k

i

i ki

k

i

i i

k

i

i i

ss ss

ss ss

ss ss

Cov , Cov ,

Cov , Cov ,

Cov , Cov ,

0

1

20

1

2

10

1

1

λ

λ

λ



( ) ( ) ( ) ∑=

= −

k

i

SK i i

1

0 00 0

2 σ s Cov ,ss λ Cov ,ss

λi
143

MATRICES 

A matrix is a rectangular 

array of numbers, such as A. 

When n = m, A is a square 

matrix of order n.

Transposing all rows and 

columns in A is denoted as

 Matrices are a convenient 

notation heavily used in 

dealing with large systems of 

linear equations, which 

notably reduce in size to just 

A X = B or X = A-1B, where

























=

nn nm

m

m

aaa

aaa

aaa









21

2221 2

1211 1

A

The main diagonal of a square matrix is the sequence

of elements a11, a22, ..., ann from upper left to lower right.

A-1 denotes the inverse of matrix A.

. T A

[ ]

T B = 21  bbb m

[ ]

T X = 21  xxx m
144

SIMPLE KRIGING MATRIX FORMULATION 

( ) ( )

( ) ( )

( )

( )



















=









































k kkk K

k

ss

ss

ss ss

ss ss

Cov ,

Cov ,

Cov , Cov ,

Cov , Cov ,

0

1 10

1

11 1

 







λ

λ

( )

( )

( ) 







































−

−

+=

k

T

k

SK

mZ

mZ

Z m

λ

λ



1 1

0

*

s

s

s

( ) ( )

( )

( ) 







































= −

k

T

k

SK

λ

λ

σ 

1

0

10

0 00

2

Cov ,

Cov ,

Cov ,

ss

ss

s ss
145

EXERCISE 

0

50

100

150

200

250

300

0 100 200 300 400

Easting

Northing

3 

2 

1 

4 ? 

m=110

SAMPLE

Index Easting Northing Attribute

1 10 20 40

2 30 280 130

3 250 130 90

4 360 120 160

? 180 120

( ) / 250 Cov 2000 h h − = e
146

EXERCISE 

For the previous configuration, using simple kriging:

A. Discuss the weight values.

B. Find the estimate.

C. Compute the estimation variance.

D. Calculate the estimation weights.
147

A. CALCULATION OF WEIGHTS (1)

Distance matrix among observations:

Distance from estimation location to each observation:

























364 0. 366 7. 110 04.

264 0. 266 03.

260 08.

0

[ ]

T 197.2 219.3 70.7 180.0
148

CALCULATION OF WEIGHTS (2)

Weights:

























−

=















































 −

.0 001

.0 646

.0 128

.0 185

973 6.

1507 2.

831 8.

908 7.

466 4. 461 2. 1285 8. 2000

695 6. 689 4. 2000

704 8. 2000

2000 1
149

B. REMARKS

• Weight 4 is almost zero

despite being closer to

the estimation location

than weight 1. This

reduction known as

screen effect is caused

by the intermediate

location of weight 3.

• Weight 4 is negative.

• Negative weights open

the possibility to have

estimates outside the

extreme values of the

data (nonconvexity).

0

50

100

150

200

250

300

0 100 200 300 400

Easting

Northing

1 

0.185

2 

0.128

 3 

0.646

 4 

-0.001
150

C. ESTIMATE

( )

( ) 180, 86120 7.

.0 001

.0 646

.0 128

.0 185

110160

11090

110130

11040

180, 110120

*

*

=

























− 























−

−

−

−

+=

SK

T

SK

z

z

( )

( )

( ) 







































−

−

+=

k

T

k

SK

mZ

mZ

mz

λ

λ



1 1

0

*

s

s

s
151

D. ESTIMATION VARIANCE

( )

( ) 180, 752120 9.

.0 001

.0 646

.0 128

.0 185

973 6.

1507 2.

831 8.

908 6.

180,120 2000

2

2

=

























− 























−=

SK

T

SK

σ

σ

( ) ( )

( )

( ) 







































= −

k

T

k

SK

λ

λ

σ 

1

0

10

0 00

2

Cov ,

Cov ,

Cov ,

ss

ss

s ss
152

PROPERTIES OF SIMPLE KRIGING

• Nonconvexity

• Screen effect

• Declustering

• Global unbiasedness

• Minimum mean-square-error optimality

• Exact interpolator

• When interpolation is exact, the error variance is 0

• Intolerance to duplicated sites

• The estimate is orthogonal to its error

• Independence to translation of Cartesian system
153

PROPERTIES OF SIMPLE KRIGING

• Independence of estimator to covariance scaling

• Scaling of covariance changes the estimation

variance by same factor

• Estimation variance depends on data configuration

• Estimation variance does not depend directly on

individual sample values

• If the weights sum up to 1, the estimator is

independent of the value of the mean
154

MULTIGAUSSIAN KRIGING

• This is kriging of normal scores.

• The univariate marginal distribution of a multivariate

normal distribution is normal, yet the converse is not true.

• Under multivariate normality, simple kriging is the best of

all possible estimators, linear or nonlinear, biased or

unbiased.

• Multivariate normality cannot be verified; at most, it

cannot be rejected.

• Considering that the mean of normal scores is zero, if

they are second order stationary, they provide a rare

opportunity to apply simple kriging.

• Backtransformation of estimated normal scores outside

the interval of variation for the normal scores of the data

is highly uncertain.
155

9. ORDINARY KRIGING
156

BRILLIANT TRICK

One can rewrite the simple kriging estimator as

 If the weights sum up to 1, then the estimator is 

independent of m, thus resulting in a more general 

formulation applicable to samples with constant, but 

unknown mean (no trend). 

 It can be proved than the constraint has the additional 

benefit of making the estimator unbiased.

 Ordinary kriging was a first formulation of improved 

simple kriging. Simple kriging roots predate 

geostatistics. We will see that one can still find weights 

minimizing the mean square estimation error.

( ) 1 ( ). 1 1

0

*

∑∑= =

+





 







−=

n

i

ii

n

i

ZSK s m i λλ Z s
157

ORDINARY KRIGING

The ordinary kriging estimator is:

where:

 denotes estimation at geographic location

 is the location of measurement ;

 denotes the number of observations to consider;

 is the mean of and

is a real weight.

k

λi

si

( ) 0

* zOK s

m

; s0

Z( ) s ;

i

( )

( )















=

=

∑

∑

=

=

k

i

i

k

i

ii

OK

Z

Z

1

1

0

*

λ ,1

λ s

s
158

DERIVATION OF EQUATIONS

A. Make necessary assumptions;

B. Find expression for estimation error variance;

C. Modify estimation error variance equation to have

it in terms of covariances;

D. Out of all possible weights, find those that

minimize the error variance.

E. Replace the optimal weights in the expression for

error variance. This is the kriging variance.

The whole purpose of ordinary kriging is to find 

optimal weights that minimize the mean square 

estimation error. The steps are the same as those 

required to formulate simple kriging:
159

A. STRONGER ASSUMPTIONS

1. The sample is a partial realization of a random

function , where denotes spatial location.

2. is second order stationary, that is:

 where denotes expected value, is a scalar 

constant, and is a vectorial distance.

 

Z( ) s s

E[ ] ( ) s = mZ

E[ ( )( ) ( )(ZmZ hss m)] =−+− Cov( ) h ,

E[ ]⋅

h

m

Z( ) s

 As in simple kriging, the mean is still constant, but 

now it can be unknown.
160

WEAKER ASSUMPTIONS

1. The sample is a partial realization of a random

function , where denotes spatial location.

2. honors the intrinsic hypothesis, meaning that:

Z( ) s s

E[ ] ( ) s = mZ

Var[ ( ) ( ) ZZ ] =+− γ ( ) hhss ,2

Z( ) s

The following is a set of slightly weaker assumptions, 

in the sense that they do not require that the variance 

be finite. 

where is the semivariogram. γ ( ) h

 The mean must be constant, but because it will 

not enter into the calculations, it can be unknown.
161

B. ESTIMATION VARIANCE

The estimation variance is the error variance:

where is the random function at the estimation 

location and is the ordinary kriging estimator.

 We know neither of these; thus, it is not possible to 

calculate . The way out is to transform the

expression for the variance to have it in terms of 

quantities that we can calculate. This will require even 

more laborious algebraic manipulations than for simple 

kriging. Nothing comes for free!

( ) ( ) ( ) [ ] 0 0

*

0

2 σ s = Var OK − ZZ ss

( ) Z s0

( ) 0

* ZOK s s0

( ) 0

2 σ s

, 
162

LEMMA 1 (1)

Let be a random variable of a continuous random 

function. Then for any coefficient :

Proof:

One can regard as a new constant and 

 as a new variable.

( ) Z si

E ( ) E[ ] ( ) ( ). 1 1

2

1

∑ ∑∑ = = =

=



























 k

i

k

j

i ji ji

k

i

λiZ s λλ ZZ ss

λi

λλ ji

( ) ( ) ZZ ss ji
163

LEMMA 1 (2)

Then, by the distributive and commutative properties of 

the expectation:

and the proof follows by condensing the expression 

using a second summation index.

( ) [ ] ( ) ( )

∑ [ ] ( ) ( ) ∑ [ ] ( ) ( )

∑ ∑

= =

= =

+ ++

=





























k

i

ik ik

k

i

i i

k

i

i i i

k

i

i

ZZ ZZ

Z ZZ

1 1

2 2

1

1 1

2

1

E E

E E

ss ss

s ss

λλ λλ

λ λλ


164

LEMMA 2 (1)

Let be a random variable of a continuous random 

function. Then for any coefficient :

Proof:

By the definition of variance (pages 31 and 52):

( ) Z si

Var ( ) Cov[ ].,

1 1

2

1

∑ ∑∑ = = =

=



























 k

i

k

j

i ji ji

k

i

λi

Z s λλ ss

λi

Var ( ) E ( ) E ( ) .

2

1

2

1 















 



 −





























=



 





∑ ∑ ∑ = = =

k

i

ii

k

i

ii

k

ii

λ Z sii λ Z s λ Z s
165

LEMMA 2 (2)

In the right side, by Lemma 1 and expanding second the term:

Factoring,

The fact that the difference is equal to the covariance of 

 and proves the lemma.

( ) [ ] ( ) ( )

[ ] ( ) [ ] ( ).EE

Var E

1 1

1 1 1

∑∑

∑ ∑∑

= =

= = =

−

=



 





k

i

k

j

ji i j

k

i

k

j

i ji ji

k

i

i

ZZ

Z ZZ

ss

s ss

λλ

λ λλ

( ) Z si ( ) Z sj

( ) { [ ] ( ) ( )

[ ] ( ) [ ] ( )}.EE

Var E

1 1 1

i j

k

i

k

j

i ji ji

k

i

i

ZZ

Z ZZ

ss

s ss

−

=



 





∑ ∑∑ = = =

λ λλ
166

LEMMA 3 (1)

Let be an intrinsic random function. Then:

Proof:

By the stationarity assumption:

which does not change by adding and subtracting a third 

variate

( ) ( ) ( )

Cov( ssss ).,

ssssss

−−−

−+−=−

ji

ji i j γ γγ

Z( ) s

( ) E[{ ( ) ( )} ], 2

1 2 ss ji =− i − ZZ ss j γ

( ) E[{{ } ( ) ( ) { } ( ) ( ) } ]. 2

1 2

γ ss ji =− i j −−− ZZZZ ssss
167

LEMMA 3 (2)

Expanding:

and again, by the stationarity assumption:

The proof follows from the definition of the covariance 

because of the special fact that the mean of the differences 

is zero because the mean is constant.

( ) [{ } ( ) ( ) ] [{ } ( ) ( ) ]

E[ ( ) ( ), ( ) ( )],

E

2

1 E

2

1 2 2

ssss

ss ss ss

ZZZZ

ZZ ZZ

i j

ji i i

−− −

 =− +− −

( ) ( ) ( )

E[ ( ) ( ), ( ) ( ) ssss ].

ssssss

i j ZZZZ

ji i j

−− −

  −+−=−
168

C. NEW ERROR VARIANCE EXPRESSION (1)

Replacing the estimator by its definition in the expression 

for the error variance:

because From Lemma 2,

which, by Lemma 3, is:

( ) ( ) ( ) [ ] ( ( ) ( ))



 





= =− ∑ − =

k

i

OK ZZ ii ZZ

1

0 0 0

*

0

2 σ s Var ss Var λ ss

( ) Cov( ),,

1 1

0 0 0

2

∑∑= =

= −−

k

i

k

j

σ s λλ ji i j ssss

( ) { ( ) ( ) ( )}. 1 1

0 0 0

2

∑∑= =

= −−−+−

k

i

k

j

s iji j ssssss ji σ γλλ γ γ

.1

1

∑ = =

k

i

λi
169

NEW ERROR VARIANCE EXPRESSION (2)

Expanding, because of the constraint :

The first two terms in the right-hand side of the equation 

are the same, thus:

1

1

∑ = =

k

i

λi

( ) ( ) ( )

( ). 1 1

1

0

1

0 0

2

∑∑

∑ ∑

= =

= =

− −

= +− −

k

i

k

j

jiji

k

j

jj

k

i

ii

ss

s ss ss

γλλ

σ γλ γλ

( ) ( ) 2 ( ). 1 1 1

0 0

2

∑ ∑∑ = = =

= −− −

k

i

k

i

k

j

s ii ss ss jiji σ γλ γλλ
170

D. OPTIMAL WEIGHTS

The expression for is quadratic in the unknowns, 

the weights. In such a case, the necessary and sufficient

condition to have a unique global minimum is that the 

quadratic term must be equal to or larger than zero, which 

implies that the semivariogram must be negative definite. 

Thus, we need to employ only permissible models. 

The expression for is valid for any weights. 

Among the infinite weight combinations, we are interested 

in those minimizing the estimation variance associated 

with the estimator . This is a constrained

operations research problem, best solved by the Lagrange 

method of multipliers.

( ) 0

* ZOK s

( ) 0

2 σ s

( ) 0

2 σ s
171

LAGRANGE METHOD 

The Lagrange method of multipliers is a procedure for the 

optimization of a function of several variables

 constrained by .

 The solution is found by equating to 0 the partial 

derivatives of an unconstrained auxiliary function

with respect to all , regarding the Lagrange multipliers,

, as constants.

The resulting equations plus the original constraints 

provide a system of equations from which the unknown 

and the auxiliary variables, , can be calculated.

( ) yyyf n ,,, 21  φφφ m ,,, 21 

n

( ) ∑=

+

m

i

yyyf n ii

1

21 ,,, φμ

yi

μi

n

i y

μi

m
172

TRIVIAL LAGRANGE EXAMPLE 

To find the maximum area of a rectangle whose perimeter 

is the given value , it is necessary to find the maximum 

value of the product of the sides , subject to the 

constraint .

The auxiliary function is in this case

which after differentiation with respect to the unknowns 

and gives

whose solution is 

p

pba =−+ 022

⋅ ba

u ab μ( ) −++= pba ,22











=+

=+

=+

pba

b

a

22

02

02

μ

μ

a

b

== .0 25pba μ −= .0, 125.
173

NORMAL EQUATIONS

The Lagrange auxiliary function in this case is:

The derivatives with respect to the unknowns are:

Thus

2 ( ) ( ) ,02 for .,,1

1

0 kj L k

i

j jii

j

−−= ==+−  ∂

∂ ∑=

γ γλ μ λ

ss ss

( ) ,2 ( ) 2, .1

1 1 1 1

0 















= ∑ ∑∑ − + ∑ − = = = =

k

j

j

k

j

k

i

k

j

L γλ jj ss γλλ ss jiji λμ

( ) ( )

( ) ( )





















=

=−

=−

∑

∑

∑

=

=

=

1

, ,

, ,

1

0

1

01

1

1

k

i

i

k

k

i

iki

k

i

i i

λ

γμγλ

γμγλ

ss ss

ss ss


174

E. OPTIMAL ESTIMATION VARIANCE

Multiplying each optimality condition by and adding them 

together, we have the new optimality relationship:

Because the weights sum up to 1, introducing under the 

second summation and replacing the resulting double 

summation in the general expression for the estimation 

variance, one obtains:

or

λi

( ) , ( ) , . 1 1

0

1 1

∑ ∑ ∑ ∑ = = = =

= +

k

i

i

k

i

ii

k

i

k

j

γλλ ss jiji ss λμγλ

λi

σ ( ) ( ) ( ) = ∑ ∑ γλ − − μγλ = =

k

i

k

i

OK ii ii

1 1

0 0 0

2 s ,2 ss ,ss

( ) ( ) .,

1

0 0

2 σ = ∑ − μγλ=

k

i

OK s ii ss
175

SUMMARY

Ordinary kriging normal equations for optimal weights : 

Ordinary kriging estimate:

Ordinary kriging mean-square-error variance:

( ) ( ) ∑=

=

n

i

zOK z ii

1

0

* s λ s

λi

σ ( ) ( ) = ∑ − μγλ=

k

i

OK ii

1

0 0

2 s ,ss

( ) ( )

( ) ( )





















=

=−

=−

∑

∑

∑

=

=

=

1

, ,

, ,

1

0

1

10

1

1

k

i

i

k

k

i

iki

k

i

i i

λ

γμγλ

γμγλ

ss ss

ss ss


176

INTRINSIC MATRIX FORMULATION

( ) ( )

( ) ( )

( )

( )























=

























− 























1

,

,

0111

, 1,

1

, 1,

0

1 10

1

11 1

k kk k K

k

ss

ss

ss ss

ss ss

γ

γ

μ

λ

λ

γ γ

γ γ









( )

( )

( ) 























− 























=

μ

λ

λ

k

T

k

OK Z

Z

Z 

1 1

0

*

0

s

s

s

( )

( )

( ) 























− 























=

μ

λ

λ

γ

γ

σ

k

T

k

OK s



1

0

10

0

2

1

,

,

ss

ss
177

STATIONARY MATRIX FORMULATION

( ) ( )

( ) ( )

( )

( )























=

















































1

Cov ,

Cov ,

0111

Cov , Cov 1,

1

Cov , Cov 1,

0

1 10

1

11 1

k kk k K

k

ss

ss

ss ss

ss ss

 







μ

λ

λ

( )

( )

( ) 















































=

μ

λ

λ

k

T

k

OK Z

Z

Z 

1 1

0

*

0

s

s

s

( ) ( )

( )

( ) 















































= −

μ

λ

λ

σ

k

T

k

OK



1

0

10

0 00

2

1

Cov ,

Cov ,

Cov ,

ss

ss

s ss
178

EXERCISE 

0

50

100

150

200

250

300

0 100 200 300 400

Easting

Northing

3 

2 

1 

4 ? 

SAMPLE

Index Easting Northing Attribute

1 10 20 40

2 30 280 130

3 250 130 90

4 360 120 160

? 180 120

( ) / 250 Cov 2000 h h − = e
179

EXERCISE 

A. Find the estimation weights.

B. Note any interesting peculiarities in the weight

values.

C. Use the data and weights to evaluate the

estimate.

D. Calculate the estimation variance.

Apply ordinary kriging to the previous sample to do 

the following:

Compare to the results obtained using simple

kriging.
180

A. ORDINARY KRIGING WEIGHTS

Distance matrices:

Weights and Lagrange multiplier:

[ ]

T 197.2 219.3 70.7 180.0

364 0. 366 7. 110 04.

264 0. 266 03.

260 08.

0

























































−

=































































 −

42.714

.0 011

.0 650

.0 141

.0 198

1

973 6.

1507 2.

831 8.

908 7.

01111

466 4. 461 2. 1285 8. 2000

695 6. 689 4. 2000

704 8. 2000

2000 1
181

B. WEIGHT COMPARISON

Weight SK OK

1 0.185 0.198

2 0.128 0.141

3 0.646 0.650

4 -0.001 0.011

Sum 0.958 1.000

0

50

100

150

200

250

300

0 100 200 300 400

Easting

Northing

1 

0.198

2 

0.141

3 

0.650

4

0.011
182

C. ESTIMATES

Using the subscripts OK to denote ordinary kriging and 

SK for simple kriging,

( )

( )

( ) 180,120 was 86 7.

180, 86120 6.

.0 011

.0 650

.0 141

.0 198

160

90

130

40

180,120

*

*

*

SK

OK

T

OK

z

z

z

=

















































=

( )

( )

( ) 























− 























=

μ

λ

λ

k

T

k

OK Z

Z

z



1 1

0

*

0

s

s

s
183

D. ESTIMATION VARIANCE

( )

( ) 180, 754120 7. ( ) 180,120 was 752 9.

42.714

.0 011

.0 650

.0 141

.0 198

1

973 6.

1507 2.

831 8.

908 6.

180,120 2000

2 2

2

OK SK

T

OK

 



=

































 −































−=

 (−42.714) provides the variation in per unit of change in 

the constraint. The change in the sum of weights from OK to SK 

is 0.042. Hence the reduction of −1.8 (752.9 − 754.7) agrees 

with the fluctuation in (0.042·(−42.724)). 

2 

2 
184

PROPERTIES 

Ordinary kriging shares the same properties with simple 

kriging. In addition,

• . 

• For the trivial case of considering only one observation

in the calculation of the estimate,

 where is the distance from the observation, , 

to the estimation location, .

( ) ( ) 0

2

0

2

OK ≥ σσ SK ss

h01

( ) ( ) 0 01 ,2 2 σ OK ⋅= γ hs

s0

( ) s1 z
185

SEARCH NEIGHBORHOOD

• Because of the screen

effect and numerical

instabilities, it is

recommended that only

the closest observations

to the estimation location

be used.

• Three observations are a

reasonable bare minimum

and 25 are more than

adequate.

• Use octant search to

further ensure good radial

distribution.
186

SUBSAMPLE SIZE SENSITIVITY FOR THE 

CASE OF THE UNCF DATASET

5 0.351 0.311 0.043 0.157 0.137 --- 0.766 0.159

7 0.358 0.328 0.047 0.157 0.137 -0.047 0.750 0.159

8 0.356 0.325 0.050 0.178 0.165 -0.074 0.767 0.158

10 0.355 0.324 0.050 0.178 0.170 -0.077 0.756 0.158

15 0.355 0.324 0.051 0.178 0.169 -0.077 0.753 0.158

20 0.355 0.324 0.052 0.178 0.169 -0.078 0.754 0.158

25 0.355 0.324 0.052 0.178 0.169 -0.078 0.758 0.158

λ1 λ 2

*

OK ∑ z =

k

i

i

6

λ

4 λ 5 λ λ 3

* σ OK k
187

 The sampling is perfectly regular on a 286 by 286 grid, 

making 81,796 observations with a Gaussian anisotropic 

semivariogram (p. 124): 

 The exhaustive sample will be used solely for comparisons. 

All calculations will be done employing the same subset of 

size 40. Based partly on crossvalidation, its semivariogram is

THE ELY WEST SAMPLE

γ ( ) h .0 1 13 ⋅+= Gaussian(N15E,80,35).

This is a semi-artificial, exhaustive 

dataset that will be used extensively to 

illustrate and compare several methods.

 The data are transformed digital 

elevations in an area in Nevada with 

extreme cases of continuity: flat dry 

lakes surrounded by rugged sierras 

running roughly north-south.

( ) h 1.0 += 11⋅Gaussian(N10E,80,34). E γ
188

ORDINARY KRIGING WITH SGeMS

A. LOAD SAMPLE

Load data according to instructions in 

chapter 4.

 If data have been previously loaded 

and saved:

• click on <File>;

• click on <Open Project>;

• select folder and click <OK>.

Loading is complete.

 To post data, operate on 

<Object> and <Preferences>.
189

B. CREATE GRID FOR STORING RESULTS

• Click on <Objects>;

• click on <New

Cartesian Grid>;

• fill in options;

• click on <Create

Grid>.
190

C. SET PARAMETERS

• Click on

<Estimation> and

then on <kriging>.

• Complete the form.

To be able to do

that, it is required to

previously have

modeled the

semivariogram.

• Click on <Run

Algorithm> when

ready.

You may want to click on <Save> to store all parameters 

as *.par file for further reference. 
191

D. RESULTS

*

zOK

2 σ OK

Click on camera icon ( ) at lower right corner of screen 

to save maps as electronic files.
192

10. UNIVERSAL KRIGING
193

INTRODUCTORY REMARKS

There are several instances of attributes—water depth 

near the shore, temperature in the upper part of the 

Earth’s crust—that have a clear, systematic variation. 

Hence, models that presume constancy of the mean are 

inadequate for the characterization of such attributes.

 Universal kriging is a full generalization of simple 

kriging that takes to the limit the ordinary kriging 

improvement. Universal kriging removes both the 

requirement that the attribute must have a constant mean 

and that the user must know this constant mean.

 As the model complexity goes up, unfortunately so does 

uncertainty. Thus, universal kriging should not be used 

indiscriminately.

Universal kriging is also known as kriging with trend.
194

ESTIMATOR 

The universal kriging estimator has the same form as that 

for ordinary kriging:

where

 denotes spatial location,

 are real weights,

 is the number of observations to consider, and

 is a variate at site .

( )

( )















=

=

∑

∑

=

=

k

i

i

k

i

ii

UK

Z

Z

1

1

0

*

λ ,1

λ s

s

s

λi

k

( ) Z si si
195

ASSUMPTIONS 

1. The residuals are a partial realization of a 

random function , where is 

the trend or drift of the random function .

2. As in ordinary kriging, one can select second order

stationarity or intrinsic hypothesis, depending on

whether the derivations will be done on terms of

covariances or semivariograms of the residuals.

3. The mean is now neither known nor constant.

4. The trend in the mean is considered a deterministic

component amenable to analytical modeling. The

most common practice, yet not the only alternative, is

to employ polynomials of the geographical

coordinates. The coefficients become additional

unknowns in the trend expression.

( ) ( ) ( ) −= mZY Z sss

Y( ) s

mZ ( ) s

Z( ) s
196

TREND (DRIFT)

( ) ∑=

+=

p

i

Z iifaam

1

s 0
197

NORMAL EQUATIONS

( ) ( ) ( )

( ) ( ) ( )

( ) ( ) ( )

( ) ( )

( ) ( )

( ) ( ) 













































=

=

=

=

++ =

++ =

++ =

∑

∑

∑

∑

∑ ∑

∑ ∑

∑ ∑

=

=

=

=

= =

= =

= =

0

1

02

1

2

01

1

1

1

0

1

0

1

20

1

0 2

1

2

10

1

0 1

1

1

1

Cov , Cov ,

Cov , Cov ,

Cov , Cov ,

ss

ss

ss

ss s ss

ss s ss

ss s ss

p

k

i

ipi

k

i

ii

k

i

ii

k

i

i

Y k

p

l

kll

k

i

kiYi

Y

p

l

ll

k

i

iYi

Y

p

l

ll

k

i

iYi

ff

ff

ff

f

f

f

λ

λ

λ

λ

λ μμ

λ μμ

λ μμ





Note that the covariance is that of the residuals.
198

ESTIMATION ERROR VARIANCE

( ) ( ) ∑ ( ) ∑ ( ) = =

= − −−

p

j

jj

k

i

UK Y iYi f

1

0 0

1

0 0

2 σ s Cov 0 λ Cov ,ss μμ s
199

SOME MATRIX DEFINITIONS

( ) ( ) ( ) ( )

( ) ( ) ( ) ( )

( ) ( )

( ) ( ) 











































=

00 0

00 0

1 001 0

Cov , Cov 1,

Cov , Cov 1,

1

11 1

1 1

11 1 11 1

 

   

 

 

 

   

 

p kp

k

kY kkY k kp

Y kY p

f f

f f

f f

f f

s s

s s

ss ss s s

ss ss s s

K

[ ( ) ( ) ( ) ( )]

T

Y Y k p f f 10 Cov 0 1, 01 0 k = Cov ,ss  ss s  s

[ ]

T

x = 1  k −− μμλλ 10  − μ p

[ ( ) ( ) ]

T

= Z sZ 1  Z sk  000
200

MATRIX FORMULATIONS

Kx = k

( ) xZs T ZUK 0 =

*

( ) xkT

UK CovY 0 −=

2 σ

Normal equations

Estimator

Estimation variance
201

NUMERICAL EXERCISE

Elevation of water table in northwestern Kansas is a 

clear example of an attribute with a linear trend. Use 

the sample to:

A. Apply the semivariogram estimator on page 114 along

different directions.

B. Obtain the semivariogram for the residuals by fitting a

model to the experimental semivariogram along the

trend-free direction—in this case, the only

experimental semivariogram that is not an artifact.

C.Use the closest observations to estimate the weights

and Lagrange multipliers at location (60, 193).

D.Find the estimate and its kriging variance.

E. Produce a map of water table elevations using

SGeMS.
202

A. EXPERIMENTAL SEMIVARIOGRAM

FOR THE ATTRIBUTE
203

B. SEMIVARIOGRAM MODEL ALONG

TREND-FREE DIRECTION (N20W)

 ( ) h 119 += Gaussian(sill = 3948, = 29range 8. )

A low angular tolerance of 10° was used to properly 

filter out the trend. 
204

C. WEIGHTS

Well

Easting

mi

Northing

mi

Elevation

ft

993 61.56 197.85 3065.0

1002 62.94 194.81 3099.4

1003 55.68 193.56 3200.0

1502 64.96 189.77 3114.9

1504 54.80 190.60 3217.1

1505 59.12 189.47 3189.7

60.0 193.0 ?

187

189

191

193

195

197

199

54 56 58 60 62 64 66

Easting, mi

Northing, mi

993

1002

1003

1502

1504

1505

Data
205

DISTANCES 









































.8 72 .6 57 .5 34 .5 85 .4 47 0

.9 91 .9 17 .3 1009 .20 0

.8 77 .5 1044 .03 0

.7 27 .7 37 0

.3 34 0

0

[ ]

T .5 09 .3 45 .4 36 .5 92 .5 73 .3 64

Observation–observation:

Estimate–observation:
206

COVARIANCES AND UNKNOWNS

























































=

187. 19485 . 19381 . 18956 . 19077 . 18960 .47 000

61. 6256 . 5594 . 6468 . 5496 . 5980 .12 00

0111111

2716 0. 3084 4. 3262 7. 3192 0. 3372 7. 3710 5.

2494 0. 2633 9. 3511 0. 2438 4. 3710 5.

2707 1. 3250 3. 2471 3. 3710 5.

2970 1. 2954 3. 3710 5.

3489 4. 3710 5.

3710 5.

K

[ ]

T k = 3296 1. 3478 8. 3385 5. 3181 3. 3209 4. 3461 14. 60 0. 193 0.

[ ]

T

x = .0 073 .0 348 .0 274 .0 012 − .0 056 .0 33350 .95 .0 27 .0 53
207

D. ESTIMATE AND VARIANCE

( )

( ) 60,193 3149 4.

.0 53

.0 27

33.95

.0 350

.0 056

.0 012

.0 274

.0 348

.0 073

0

0

0

3189 7.

3217 1.

3114 9.

3200 0.

3099 4.

3065 0.

60,193

*

*

=

























































−

























































=

UK

T

UK

z

z ( )

( ) 60, 111193 .71

.0 53

.0 27

33.95

.0 350

.0 056

.0 012

.0 274

.0 348

.0 073

193.00

60 0.

1

3461 4.

3209 4.

3181 3.

3385 5.

3478 8.

3296 1.

60,193

2

2

=

























































−

























































=

UK

T

UK

σ

σ
208

WEIGHTS, ESTIMATES, AND VARIANCES

187

189

191

193

195

197

199

54 56 58 60 62 64 66

Easting, mi

Northing, mi

3149.6

111.7

3065

0.073

3099

0.348

3115

0.012

3200

0.274

3217

-0.056 3189

0.350
209

E. WATER TABLE ELEVATION WITH SGeMS

• Click on

<Estimation and

<kriging>;

• Fill in form and

modify default

values. Trend

model is:

where:

 trend

a, b, c: coefficients

X: easting

Y: northing

• Click on <Run

Algorithm>.

( ) s ⋅+⋅+= YcXbam

m( ) s :
210

MAPS 

2 σ UK

* ZUK

ft

sq ft
211

11. BLOCK KRIGING
212

MOTIVATION 

• Block kriging is the generic name given to any form of

kriging in which one is interested in estimates of the

linear average of an attribute inside supports that are

intermediate in size between the support of the

observations and the sampling domain.

• This formulation was originally the most widely used

form of kriging in mining applications.

• Any form of point kriging is a limiting form of an

equivalent formulation of block kriging.

• Ordinary block kriging under second order stationarity

is the most common form of block kriging.
213

LINEAR AVERAGE

Let be a point support random function in and 

let be a support of finite size centered around 

. . The true block average at location is:

( ) ( )

( )

. 1

0

s0 ss

s

dZ

S

Z

S

S = ∫ S

Z( ) s D

( ) S s0 S

s0 ( ) ZS s0 s0
214

ESTIMATOR 

The purpose of block kriging is to produce direct block 

averages from point or quasi-point measurements, not 

by averaging point estimates.

Let be variates of a random function and let 

. be the support for the estimator at 

site . Then the ordinary block kriging estimator is:

( ) Z si k

( ) s0 ⊂ DS ( ) 0

* ZS s

s0

( ) ( )















=

=

∑

∑

=

=

1

1

1

0

*

k

i

i

k

i

ZOBK sZ ii

λ

s λ
215

POINT TO BLOCK COVARIANCE

The main peculiarity of block kriging is the use of 

supports of two different sizes for the observations and 

the estimator, and hence covariances between two 

different volume sizes.

Let be covariances between the variate 

and all the variates inside the block . Then the 

covariance between the site and the block is:

Cov( ) 0,ss ( ) Z s0

S

( ) Cov( ) .,

1 Cv , 0 = ∫ 0

S

d

S

s S sss

Z( ) s

Cv( ) s0,S
216

BLOCK TO BLOCK COVARIANCE

Let be covariances between a variate 

in block and variates inside the block . Then 

the block covariance between two blocks is:

Cov( ) ,us

Si Z( ) u

( ) Cov( ) .,

1 CV , ∫ ∫ =

S Si j ji

ji dd

SS

SS usus

Z( ) s

( ) SS ji CV ,

Sj
217

ASSUMPTIONS 

The minimum set of assumptions required to formulate 

ordinary block kriging are the same as those for ordinary 

kriging:

1. The sample is a partial realization of a random 

function, Z (s ) , where s denotes spatial location.

2. Z (s ) is second order stationary, which implies:

where E [ ⋅] denotes expected value; m is a scalar,

unknown constant; and h is a vectorial distance.

3. The mean is constant, but unknown.

The SGeMS implementation requires that the block size

be of the same dimensions as the grid cell size.

E[ ] ( ) s = mZ

E[ ( )( ) ( )(ZmZ hss m)] =−+− Cov( ) h ,
218

NORMAL EQUATIONS

In an entirely similar way to the derivation of the normal 

equations for point ordinary kriging, one can obtain the 

following system of equations that provides the ordinary 

block kriging weights :

( ) ( )

( ) ( )

( ) ( )



























=

=+

=+

=+

∑

∑

∑

∑

=

=

=

=

1

Cov , Cv ,

Cov , Cv ,

Cov , Cv ,

1

1

2

1

2

1

1

1

k

i

i

k

k

i

i ki

k

i

i i

k

i

i i

S

S

S

λ

λ μ

λ μ

λ μ

ss s

ss s

ss s



λi
219

ESTIMATION VARIANCE

The minimum mean square error for ordinary block 

kriging, is:

The form of the expression remains the same as the one 

for point ordinary kriging, but all point covariances have 

been replaced by block-to-block or point-block 

covariances.

σ ( ) ( ) ( ) = −∑ − μλ =

k

i

OBK SS ii S

1

0

2 s CV , Cv s ,

,

2 σ OBK
220

POINT COVARIANCE MODELING

As in other forms of kriging, performance of block kriging 

presumes knowledge of some covariance terms, which in 

this case, however, are of different types.

 The solution to this requirement is not different from what 

we have seen before: model the covariances and ignore 

the effects that the use of models instead of the true 

covariances have in the normal equations. The common 

practice is to first model the semivariogram and then use

to obtain the point-point covariance model. The blockblock and point-block covariances are derived from the 

point-point covariance.

Cov( ) ( ) ( ) h Cov 0 −= γ h
221

NON-POINT COVARIANCE MODELING

In early days of geostatistics there was great effort 

devoted to solve exactly the integrals for point-block and 

block-block covariances for multiples shapes of blocks.

 Today the simplest and prevailing practice is to calculate 

the integrals by numerical approximation. The block is 

tessellated into smaller units that are treated as points and 

then the block covariance is calculated as the average 

between all possible point-point covariances. The 

following table gives the recommended minimum number 

of subdivisions.

Dimension Subdivisions

1 10

2 6 by 6

3 4 by 4 by 4
222

EXAMPLE OF NUMERICAL 

CALCULATION

When the point-point covariance is the exponential 

model e-h/0.8, the values show the point covariance 

between the point on the left and the center of squares 

0.05 miles on the side. Averaging the 25 numbers 

gives a value of 0.682 for Cv(0.3).
223

AVERAGE POROSITY

WEST LYONS FIELD, KANSAS

Ordinary kriging

Cell size: 0.05 mi

Ordinary block kriging

Cell size: 0.25 mi
224

STANDARD ERROR OF POROSITY 

Ordinary kriging

Cell size: 0.05 mi

Ordinary block kriging

Cell size: 0.25 mi
225

12. COKRIGING
226

INTRODUCTION

• Simple, ordinary, and universal kriging are not

multivariate models in the usual statistical sense of the

term. Despite the fact that they employ random function

models comprising infinite numbers of random variables,

they are employed in the modeling of one attribute at a

time.

• The main purpose of cokriging is to work as a true

multivariate estimation method able to deal

simultaneously with two or more attributes defined over

the same domain, which is called a coregionalization.

• Not all attributes must be measured at all locations, but a

bare minimum of collocated measurements per pair of

attributes is necessary for structural analysis.
227

VECTORIAL NOTATION

Like in spatial univariate geostatistics, geographical 

location is the vectorial term:

 Cokriging makes heavy use of vectorial and matrix 

notation for conciseness in the notation and similarity with 

the univariate formulations. Let us start with the definition 

of a vectorial random function.

If is a set of random 

functions defined over the same domain, then the vectorial 

random function is the vectorial matrix: 

( ) ( ) ( ) ( ) 21  k ,,,,, ZZZZ p ssss

[( ,, )] .

T si = eastingi northingi elevationi

p

( ) sZ

( ) ( ) ( ) ( ) ( ) [ ] . 21

T

= k  ZZZZ p sssssZ
228

COEFFICIENT MATRIX

The coefficient matrix, , is the square matrix:

where are real numbers and is the number of 

attributes to model.

Λi





























=Λ

i

pp

i

p

i

p

i

p

ii

i

p

ii

i

λλλ

λλλ

λλλ









21

2221 2

1211 1

i λ jk p
229

ORDINARY COKRIGING ESTIMATOR

The objective of cokriging is to find the minimum meansquare-error weights for the linear estimator, :

or

both subject to

( ) 0

* sZ

( ) ( ) ∑=

Λ=

n

i

ii

1

0

* sZ sZ

. 1

∑ =Λ I =

n

i

i

( )

( )

( )

( )















































Λ

Λ

Λ

=

n

T

n sZ

sZ

sZ

sZ 

2

1

2

1

0

*
230

VECTORIAL MOMENTS

For the purpose of listing the assumptions, we need to 

define the vector of means, 

and the matrix covariance

m:

[ ]

T m = 21  mmm p

( ( ) ( ))

( ( ) ( )) ( ( ) ( ))

( ( ) ( )) ( ( ) ( ))

.

Cov , Cov ,

Cov , Cov ,

,

1 1

11 11 11





















+ +

+ +

=+

hsZsZ hsZsZ

hsZsZ hsZsZ

Cov hsZsZ

pp pp p p

p p



  


231

ORDINARY COKRIGING ASSUMPTIONS

1. The sample is a partial realization of a vectorial

random function, .

2. is second order stationary, which implies:

3. Each term in m is unknown, but constant.

4. None of the variables is a linear combination of the

others:

E[ ] ( ) = msZ

E[ ( )( ) ( )( mhsZmsZ )] =−+− Cov( ) h .

( ) ( ) .,,2,1, 0 ZaaZ for k p

kj

k +≠ ∑ jj =  ≠

s s

( ) sZ

( ) sZ

The minimum set of assumptions required to formulate 

ordinary cokriging under second order stationarity are:
232

NORMAL EQUATIONS

where is a square matrix of Lagrange multipliers of the 

same order p as the number of attributes in the 

coregionalization.

( ) ( )

( ) ( )

( ) ( )



























=Λ

Λ =+

Λ =+

Λ =+

∑

∑

∑

∑

=

=

=

=

I

Cov μss Cov ss

Cov μss Cov ss

Cov μss Cov ss

n

i

i

p

k

i

p pi

k

i

i

k

i

i

1

0

1

'

02

1

2

'

2

01

1

1

'

1

, ,

, ,

, ,



μ
233

ESTIMATION VARIANCE

If 

then the estimation variance for ordinary cokriging 

is:

[ ]

T

W n μ ' '

2

'

1  ΛΛΛ=

( ) diagonal( ( ) ). 0

2 s Cov Wk0 T σ OC = −

[ ( ) ( ) , ] ,, 01 0

T

n k = Cov ss  Cov Iss

( ) 0

2 σ OC s
234

STRUCTURAL ANALYSIS

• Besides the covariances, the user must model

additional cross covariances.

- must be positive definite to properly solve the 

normal equations. Individual selection of permissible 

models for all covariances and cross covariances is 

not sufficient to produce a positive definite 

p

( ) pp − 2/1

Cov( )⋅

Cov( )⋅ .

The structural analysis required by cokriging is much 

more demanding than the one for kriging because:

 The common way to conduct structural analysis is 

through the application of the so-called linear 

coregionalization model.

•
235

LINEAR COREGIONALIZATION MODELING

Linear coregionalization modeling is a solution based on 

the decomposition of all structural functions in terms of a 

basic set of permissible models, first done in terms of 

semivariograms and cross semivariograms and then 

converted to covariances and cross covariances.

For all structural functions to be permissible, all coefficient 

matrices of order must be positive semidefinite. If a 

matrix is positive semidefinite, its determinant, all its 

principal minor determinants, and all its eigenvalues are 

nonnegative. A cross semivariogram, , is the moment:

( ) ( ) ∑=

=Γ

L

l

ll

1

γ hBh

Bl p

( ) E[{ ( ) ( )} ( ) ( ) { }]. 2

1

γ ij h = i i +− j ZZZZ j +− hsshss

γ ij( ) h
236

EXAMPLE OF LINEAR 

COREGIONALIZATION MODEL

The following would be an example of isotropic linear 

coregionalization for the case of two attributes considering 

nugget effect, one exponential model with a range of 200, 

and one spherical model with a range of 150:

This notation implies, for example, that the cross 

semivariogram between attributes 1 and 2 is:

  .1,

2 150

1

2 150

3

21

11

1

32

24

31

12

3

200

3





























 





























 

























 hh

e

h

h

  .1,

2 150

1

2 150

3

121

3

200

3

12 















































 hh

e

h

 h
237

EFFECTIVE APPLICATION OF COKRIGING

• The secondary variables must be significantly

correlated to the primary variable.

• A significant proportion of secondary measurements

must be at locations other than those of the primary

variable.

 A corollary is that when all attributes are measured at 

the same sites without missing values, cokriging provides 

similar results to kriging of individual attributes. 

Abundant information on additional variables is not a 

guarantee of superior results of multivariate cokriging 

over just kriging for the variable of interest. The 

following are two prerequisites:
238

POROSITY OF SEA FLOOR,

MECKLENBURG BAY, GERMANY

Sample Kriging estimation
239

COKRIGING FOR SEA-FLOOR POROSITY 

MECKLENBURG BAY, GERMANY

Cokriging porosity map using porosity 

Sample for grain size median and grain size median data 
240

13. CROSSVALIDATION
241

ACCOUNTABILITY 

• Upon selecting any estimator plus its parameters and

producing estimations for a given spatial attribute,

crossvalidation addresses the curiosity and need to know

about the quality of the job done.

• The ultimate approach would be to run an exhaustive

sampling and then assess the quality of the estimation by

comparing the new data with the collocated estimates

produced with the original data. Unquestionably, the best

estimator would be the one with the best statistics on the

basis of a previously agreed criterion.

• As conclusive as this approach may be, its

implementation ordinarily is either impossible to execute

or prohibitively expensive and defeats the purpose of the

estimation in practical situations.
242

ALTERNATIVE EVALUATION METHOD

• Given a sample of size , crossvalidation is based on the

dropping of each observation with replacement. For 

each ignored observation, an estimate is computed at 

the location of the discarded observation by using at 

most the remaining measurements.

• By pretending that a drawn observation was never taken,

one can genuinely attempt to produce an estimate, and

by comparing it to the deleted measurement, one can

have an instant and cost free error measurement.

• The discarded observation is replaced in the sample

before continuing with the analysis.

n

n −1

Crossvalidation is an ingenious alternative to exhaustive 

sampling.
243

WORD OF CAUTION

• Crossvalidation does not indicate whether an

observation, estimate, parameter, or assumption is

inadequate.

• Complex interdependencies among the errors have

precluded finding their distribution, thus hindering a

rigorous analysis or testing of the results.

• Nonetheless, crossvalidation remains a useful tool

primarily to dispel blunders and assist in drawing

honest conclusions by comparison, instead of making

arbitrary decisions or assumptions.

Crossvalidation is a fast and inexpensive way for 

indirectly testing anything about kriging or the data. 

However, one must apply crossvalidation with caution.

 Crossvalidation is not part of SGeMS, but it can be 

found in GSLIB.
244

UNFAVORABLE CIRCUMSTANCES

Special circumstances deserving extra caution are:

• In the case of kriging, the estimate is insensitive to

scaling of the semivariogram, thus making the

estimates insensitive to proportional changes in both

sill and nugget.

• Irregular sampling patterns may lead to

unrepresentative errors. A typical situation is the

sampling in clusters, which results in abnormally low

errors at all locations making the clusters. For better

results, decluster first, then crossvalidate.
245

DIAGNOSTIC STATISTICS

Comparisons or remedial actions based on crossvalidation 

are more straightforward and conclusive if the collection of 

errors is reduced to some key statistics such as:

all statistics do not always point in the same direction.

• Some quantiles of the error distribution, such as the 5th

percentile, the median, and the 95th percentile.

• The mean and standard deviation for the square of the

errors.

• The slope and intercept of the regression line between the 

measured and estimated values.

• The correlation coefficient between the measured and

estimated values.

There is plenty of redundancy in the statistics. In addition, 
246

CRITERIA

Guidelines to follow based on the statistics include:

• The ideal estimation method should result in no errors

and a correlation coefficient between measurements and

estimates of 1.

• Large individual errors may be the effect of blunders or

outliers. Data preparation and coding, including

geographical location, should be checked carefully.

• A slope different from 45o in the regression line between

measurement and estimates denotes conditional bias.

• Ideally, errors should be spatially uncorrelated, resulting

in a pure nugget effect semivariogram.

• If the z-scores—the ratio of the error over the standard

error—follow a standard normal distribution, one cannot

discard the possibility that the errors are multinormally

distributed.
247

PRACTICAL EXAMPLE: UNCF

A. The sample size is too

small for a directional study

of the semivariogram.

Employ crossvalidation for

investigating anisotropy.

B. Utilize optimal parameters

to run ordinary kriging.

C. For an 8-km search radius,

employ crossvalidation for

finding the optimal number

of estimation points.

D. Analyze crossvalidation

errors.
248

A. OPTIMAL NUMBER ESTIMATION

POINTS 

For an 8-km search radius, the mean square error of 

ordinary kriging is minimum when the maximum 

number of estimation points is 11.

Maximum number of 

observations

Mean square error

sq ft

3 589.6

4 632.1

5 658.7

6 657.4

7 593.9

8 545.2

9 539.5

10 451.0

11 361.2

12 376.0

13 385.1

300

350

400

450

500

550

600

650

700

2 4 6 8 10 12 14

Maximum number of observations

Mean square error, ft
249

B. INVESTIGATION OF ANISOTROPY

For anisotropic semivariograms, the practice is to vary 

the range as the length of a cord from the center to the 

perimeter of an ellipse. Below are the results of an 

analysis keeping the average of the axes equal to 

isotropic range.

Axes size

m

Major axis 

orientation

Mean square 

error, sq ft

 8107 8107 N 361.2

 8500 7714 N 302.8

 9000 7214 N 257.7

 9300 6914 N 247.0

 9400 6810 N 246.0

 9500 6714 N 246.2

 9600 6614 N 247.5

 9400 6814 N20E 253.8

 9400 6814 N10E 243.6 ( ) 















−+=





 



 −

2

.8 11

3

10 .5 64 1

h

γ h e

Isotropic semivariogram
250

C. OPTIMAL ORDINARY KRIGING

PARAMETERS 

• Minimum number of observations: 3

• Maximum number of observations: 11

• Search neighborhood radius: 8,000 m

• Best semivariogram model:

 Type: Gaussian

 Nugget: 10 sq ft

 Sill: 5,640 sq ft

 Major range: 9,400 m

 Direction: N10E

 Minor range: 6,814 m
251

ORDINARY KRIGING MAPS

Estimated depth Kriging standard deviation

m 

m m 

m 
252

D. CROSSVALIDATION ERRORS (1)

• Errors are slightly positive.

• Distributions for estimated

and true values are

extraordinarily similar.

• Crossvalidation errors are

spatially correlated.
253

CROSSVALIDATION ERRORS (2)

• The Kolmogorov-Smirnov pu-value is 0.59, which is

significantly larger than 0.1. Thus, according to the criteria

in p. 100, one cannot reject the normality of errors

independently of the influence of the spatial correlation.

• The larger the standard error, the larger the error

dispersion, but without being correlated.

D = 0.09

pu = 0.59
254

14. CRITICAL REVIEW
255

SAMPLING 

• Minimum mean square error methods—kriging

included—produce realistic estimations when the

average sampling interval is one order of magnitude

smaller than the minor axis of the semivariogram

range.

• As the sampling becomes sparser, maps based on

kriging show increasing smoothing, which can be

appreciated in various ways.

We explore here the weak sides of kriging as a

motivation to go into stochastic simulation.
256

EVALUATION OF ORDINARY KRIGING 

THROUGH THE EXHAUSTIVE SAMPLE 

AT ELY WEST

Exhaustive sample

81,796 measurements

Ordinary kriging map

40 measurements
257

CROSSVALIDATION 

Estimated values 

are close to being 

pure noise, with 

estimates for low 

measurements 

above the main 

diagonal and 

conversely for the 

highest data values.

 This a form of 

smoothing.
258

SEMIVARIOGRAMS OF 

ORDINARY KRIGING MAP

There is wide disagreement between the actual 

semivariogram and that of the ordinary kriging map. 

This is a consequence of smoothing.
259

HISTOGRAMS

The ordinary kriging histogram is not bimodal like the 

histogram of the exhaustive sample, which is also a 

consequence of smoothing.
260

CUMULATIVE DISTRIBUTIONS

Discrepancy in the proportion of zeros results in 

clearly different distributions.
261

ERRORS ACCORDING TO 

EXHAUSTIVE SAMPLE

In agreement with the global unbiasedness of kriging, the 

errors are symmetric around a value approximately equal 

to zero.
262

SEMIVARIOGRAM OF ERRORS 

ACCORDING TO EXHAUSTIVE SAMPLE

Errors are spatially correlated, with a range remarkably 

similar to that of the sample, but with a significantly 

smaller sill.
263

A VIEW AT ERRORS AND 

KRIGING STANDARD DEVIATIONS 

Errors Kriging standard deviation

4 

3 2 2 3 

4 

Most kriging standard deviations are no 

more than twice the error, but there is 

no correlation, nor the type of assurance 

provided by the Chebyshev’s theorem.
264

REMEDIAL ALTERNATIVES

There have been efforts along two different lines to 

eliminate smoothing:

• Compensation procedures

• Simulation

Success of compensation procedures has been

limited.

 By far the bulk of the effort and results have been in 

the generation of equiprobable realizations using 

simulation.

When simulation is used, uncertainty is modeled by 

generating several realizations. The wider the 

discrepancies, the larger the uncertainty.
265

SIMULATION METHODS

The most common methods are:

• Turning bands method

• Lower-upper triangular method

• Truncated plurigaussian methods

• Object oriented methods

• Sequential Gaussian simulation

• Simulated annealing

• Filter simulation

When a realization by any method reproduces the

value and location of all observations, the procedure is 

said to be a conditional simulation. Otherwise, it is an 

unconditional simulation. 
266

15. SEQUENTIAL GAUSSIAN

SIMULATION
267

BACKGROUND

Sequential Gaussian simulation is a stochastic method that 

is basically multigaussian kriging with feedback.

Let be a subset of variates of a random function 

and let be a sample of size . By the multiplication 

rule in p. 45, one can break a bivariate cumulative 

frequency distribution function into the product of two 

univariate cumulative frequency distribution functions:

where is the probability 

conditional to the original sample plus a value drawn 

from the distribution 

Z( ) s N

z( ) s n

( ) ( ) ( ) 1|,|,|,,,ss 2121 = s 11 ⋅ s 22 ntFntFnttF +

( ) 1|, s 22 ntF + Pr ( )( ), 22 s ≤ tZ

( ).|, s 11 ntF

( ) z s1
268

EXTENSION OF THE DECOMPOSITION 

By extension, the decomposition of a bivariate distribution 

can be generalized to any number of random variables:

 Sequential Gaussian simulation is a numerical 

application of this recursive decomposition for the special 

case when the joint distribution . 

is multivariate normal, in which case the marginal univariate 

distribution and all the conditional univariate 

distributions are also normal.

In practice, the locations are the N nodes of a regular grid.

F( nttt ) N N |,,,,,,, 21  sss 21 

( ) s1 N |, ntF

( intF −+ 1|, ) s ii

( )

( ) ( ) ( 1|,|, ).1|,

|,,,,,,,

11 22

21 21

⋅ ⋅⋅+ −+

=

NntFntFntF

F nttt

NN

N N

s s s

sss



 
269

TEN-STEP APPROACH (1)

1. In case the sample is not univariate normal, transform

the data into normal scores. SGeMS can do it, going to

<Algorithms> then <Utilities> then <transf>.

2. Using standard semivariogram modeling techniques,

find the best semivariogram model for the transformed

data.

3. Define a grid, pick nodes at random, and stack them all

in a visitation queue.

4. For faster execution, move the data to the closest node.

In case two or more observations land in the same

node, implement a criterion to end up with one value

per node, such as averaging or discarding.

To generate a partial realization of a random function at the 

nodes of a regular grid, do the following:
270

TEN-STEP APPROACH (2)

5. From the node sequence queue, draw the next location

. due for the calculation of a simulated value. Employ 

multigaussian kriging to find the estimate and the 

corresponding estimation variance using an 

expanded sample comprising all original data plus all 

values that may have been already simulated.

si

( )i z s*

( ) si

2 σ

n

6. Use and to 

define a normal distribution with mean 

and variance , 

 . 

( )i z s

* ( ) si

2 σ

( ) z si

*

( ) si

2 σ

( ( ) ( )) i i zN ss* 2

,σ

x

z(

s,x)

* z

σ
271

TEN-STEP APPROACH (3)

9. If is not the last node in the queue, then go back to

Step 5.

10. In case Step 1 was necessary, backtransform all the

values in the partial realization to the original sampling

space.

7. Draw at random a value, 

. , from 

The number is the 

simulated value at location 

. 

8. Add to the expanded 

sample.

. si

( ) z sis

( ( ) ( )).,

* 2

i i z ( ) sis zN σ ss

si

x

z(

s,x)

* z

σ

( )is z s
272

EXERCISE 

A. Generate 100 realizations by sequential Gaussian

simulation; compare the first one to the exhaustive sample.

B. Examine the semivariograms for the first realization along

the north-south and east-west directions and compare

them to the exhaustive semivariogram.

C. Compare the histograms for the exhaustive sample and

the first realization.

D. Prepare a Q-Q plot to compare the same distributions.

E. Use the exhaustive sample to obtain the errors for the first

realization. Map them and study the semivariogram.

F. Investigate errors in the first realization and the standard

deviation that results from all realizations at every node.

Use SGeMS and the Ely West sample of size 40 (page 188) 

to do the following:
273

A. SGeMS CODE FOR GENERATING 100

SEQUENTIAL GAUSSIAN SIMULATIONS

• Note that the

attribute is the

data, not the

normal scores.

• However, the

semivariogram

is for the

normal scores,

not for the

attribute.

• No backtransformation is

necessary.
274

TWO REALIZATIONS FOR ELY WEST

These are conditional realizations honoring the same 40 

data. Thus, although hard to note, the values of both 

realizations at those 40 nodes are the same.
275

BACK TO THE EXHAUSTIVE CASE

Exhaustive sample

81,796 measurements

First realization

40 measurements
276

B. FIRST REALIZATION

SEMIVARIOGRAMS

Agreement is not perfect but is significantly better 

than in the case of kriging.
277

C. HISTOGRAMS

The realization follows the sample histogram, not that of 

the exhaustive sample.
278

D. CUMULATIVE DISTRIBUTIONS

According to the Q-Q plot, agreement between distributions 

is also better than in the case of ordinary kriging. 
279

E. ERRORS ACCORDING TO

EXHAUSTIVE SAMPLE

Errors are unbiased but are larger than those for 

kriging.
280

SEMIVARIOGRAM FOR ERROR IN 

FIRST REALIZATION

The error semivariogram resembles that of the 

attribute.
281

F. A LOOK AT FIRST REALIZATION ERRORS

AND NODAL STANDARD DEVIATIONS

Errors Nodal standard deviation

4

3 2 2 3

4

Given a node, the nodal standard deviation 

is the standard deviation of the values from 

all realizations at that node.

 Larger standard deviations are associated 

with large errors, but there is no correlation.
282

16. SIMULATED ANNEALING
283

THE GENERAL IDEA

Simulated annealing was one of the first methods for the 

modeling of spatially correlated attributes that moved away 

from generating values for uninformed nodes by combining 

the sample values.

 Simulated annealing started as an approach to model the 

cooling of a piece of hot metal, but it is now more widely 

used as an optimization tool instead.

 In geostatistics, simulated annealing is used to move 

around node values in order to have a realization with 

some property as close as possible to a target function, 

which commonly is the semivariogram.

 Its main feature is the lack of any distributional 

assumptions.
284

“COOLING” A GRID

A simulation grid is populated by placing each sample 

value at the collocated or closest node. The unsampled 

nodes are given values randomly taken from a probability 

distribution, typically that of the sample.

 By so doing, the realizations honor the data and the 

selected histogram, but in a typical situation, unless the 

sample is large relative to the grid size, the grid is totally 

chaotic and lacking spatial correlation.

 The aim of simulated annealing is to transform the grid 

to gain additional properties, customarily to be spatially 

correlated by following a semivariogram model.
285

OBJECTIVE FUNCTION

When the target function is a semivariogram model , the 

objective function G is:

where is the semivariogram of the realization.

 Reduction in the value of G is achieved through swapping 

pairs of values at nodes which are not sampling locations.

 Key to a successful process is not to reject 100 percent of 

the exchanges of nodal values that do not result in an 

improvement of the objective function G. 

 The cooling schedule controls the proportion of 

unfavorable swaps as a function of total swaps attempted.

[ ] ( ) ( )

( ) ∑ − =

h h

hh

2

2

γ

γγ G G

γ ( ) h

γ G ( ) h
286

COMPONENTS OF AN 

ANNEALING SCHEDULE

t is the temperature, starting at .

Other critical components are (p. 289):

 temperature reduction factor;

 is the maximum number of all swaps for a given 

temperature, after which it is reduced by ;

 is the maximum number of accepted swaps for a given 

temperature, after which it is reduced by ; 

 is the maximum number of times for reaching ;

 is the target value for the objective function.

[ ] 





 ≤

= −

, otherwise.

i,1

Pr

t

GG

new old

old new

e

GGf

accept swap

Ka

λ

0t

λ

λ

min OG min )(

Ka

Ka

e

K
287

ANNEALING IN ACTION

G=1

Swaps=0

G=0.0001

Swaps=10,895,372
288

FEATURES 

• No computational power at all is used to calculate node

values.

• As all nodal values are drawn from the sample

distribution, its representativity is more critical than for

the other simulation methods.

• All efforts are devoted to moving around node values.

• Given a sample, processing time is in general the

longest among simulation methods, which increases

more steeply than linearly with the number of grid nodes.

• Despite its novel computational approach, simulated

annealing continues to rely on a sample and global

properties that can be derived from the sample, most

commonly the semivariogram.

• The most salient property is independence from any

distributional assumptions.
289

REPETITION OF PAGE 272 EXERCISE 

NOW FOR SIMULATED ANNEALING

The method is not part of SGeMS, but it is in GSLIB.

Parameters for SASIM (1992)

***************************

START OF PARAMETERS:

s40.dat \conditioning data (if any)

1 2 0 4 \columns: x,y,z,vr

-1.0e21 1.0e21 \data trimming limits

0 \0=non parametric; 1=Gaussian

s40.dat \non parametric distribution

4 0 \columns: vr,wt

0 14.0 \minimum and maximum data values

1 1.0 \lower tail option and parameter

4 2.0 \upper tail option and parameter

sasim.out \output File for simulation

sasim.var \output File for variogram

3 500000 \debug level, reporting interval

sasim.dbg \output file for debugging

1 \annealing schedule? (0=auto)

1.0 0.12 4000000 400000 5 0.0001 \manual schedule: t0,lambda,ka,k,e,Omin

2 \1 or 2 part objective function

112063 \random number seed

100 \number of simulations

286 0 1 \nx,xmn,xsiz

286 0 1 \ny,ymn,ysiz

1 0.0 1.0 \nz,zmn,zsiz

150 \max lags for conditioning

1 0.1 0 \nst, nugget, (1=renormalize)

3 45.0 13.0 \it,aa,cc: STRUCTURE 1

15.0 0.0 0.0 0.45 0.0 \ang1,ang2,ang3,anis1,anis2:
290

A. TWO REALIZATIONS FOR ELY WEST
291

BACK TO THE EXHAUSTIVE CASE

Exhaustive sample

81,796 measurements

First realization

40 measurements
292

B. SEMIVARIOGRAMS FOR

THE FIRST REALIZATION

Semivariogram reproduction is good.
293

C. HISTOGRAMS

Exhaustive sample First realization

There is also good reproduction of the sample histogram.
294

D. CUMULATIVE DISTRIBUTIONS

Agreement between the distributions of realizations 

and the population is limited only by the unbiasness 

of the empirical sample.
295

E. ERRORS IN FIRST REALIZATION

ACCORDING TO EXHAUSTIVE SAMPLE

Errors are comparable to those of sequential Gaussian 

simulation.
296

ERROR SEMIVARIOGRAMS 

FOR FIRST REALIZATION

Errors have a semivariogram resembling the 

semivariogram of the attribute.
297

F. TRYING TO RELATE ERRORS AND

NODAL STANDARD DEVIATIONS

Nodal standard deviation

Errors in first 

realization

4 

3 2 2 3 

4 

The tendency of errors to stay within 

four nodal standard deviations is still 

true, but it is less pronounced than 

for kriging and sequential Gaussian 

simulation.
298

17. FILTER SIMULATION
299

CLASSICAL GEOSTATISTICS

• Until recently, geostatistics had evolved to provide

modeling tools relying on methods typical of

stochastic processing.

• Central to past developments is the heavy use of

moments up to order 2, particularly the

semivariogram or its equivalent the covariance.
300

MULTIPLE POINT STATISTICS

• The first attempts to expand classical geostatistics

aimed to calculate statistics involving three or more

points.

• These efforts have had limited success because of

the high demands in terms of sample size required

for the estimation of higher order moments. 
301

MULTIPLE POINT RELATIONSHIPS

• However, the idea of considering relations among

several points has stayed.

• The concept of moment has received less

attention than techniques related more to image

analysis and pattern recognition, such as filtering

and cluster analysis.

Filter simulation is a member of the multipoint

statistics family of simulation methods. It is based on 

pattern recognition and does not require the 

semivariogram at all, even though one can still 

estimate it for purpose of evaluation or comparison.
302

MULTIPOINT SOURCES

Sources of multiple point 

relations are:

• Either images borrowed

from sites considered

analogous to the one of

interest or

• Synthetic models

generated on the basis of

fundamental geological

principles.

These pictures are called 

training images. They are 

absolutely critical to the 

success of the simulation.

Courtesy of A. Journel, Stanford University
303

FILTERING 

Filters are used to extract essential pattern information, 

such as average, gradient, and curvature. SGeMS 

employs six filters for two-dimensional patterns and nine 

in three-dimensional cases.

applied

to

applied

to

kth filter weights pattern kth score value

results

in

Filter weights a pattern a score

generate

Stanford UniversityCourtesy of A. Journel, 
304

CLUSTER ANALYSIS

The general idea is to group objects 

in attribute space into clusters as 

internally homogeneous as possible 

and as different from the other 

clusters as possible.

 Different types of approaches, 

distances, and proximity criteria have 

resulted in several methods. Filter 

simulation offers two:

• Cross partition, where classes are defined by the

quantiles of each filter score. Two patterns belong to the

same cluster if they are bounded by the same quantiles

for all filters. Results depends on the number of quantiles.

• K-means, which operates by trial and error. Cluster

membership evidently varies with the number of clusters.

0 2 4 6 8 10

Attribute 1

0

2

4

6

8

10

Attribute 2

1 

2 

3 

4 

5 

6 

A 

B 

C 

Score 1

Score 2
305

EXAMPLE OF CATEGORICAL CLUSTER

Prototype

• The prototype here

is the average of 56

template patterns

from a training

image involving two

lithologies.

• By the nature of

cluster analysis, all

the patterns are

quite similar.

Courtesy of A. Journel, Stanford University
306

FILTER SIMULATION

Filter simulation is closer to putting together a puzzle than 

to estimation. The steps to generate a realization are:

1. Define a grid and select a training image at least as large

as the study area.

2. Tessellete the training image into small templates.

3. Use different filters to generate scores to characterize

template patterns.

4. Use the filtering multivariate information to classify the

templates into clusters that group pattern styles.

5. Define a random visitation sequence.

6. Select at random a template from the cluster best

matching the conditioning data or select a template at

random if there are no conditioning data.

7. Paste the central portion of the template (inner patch) in

the realization.

8. Go back to Step 5 until the modeling is completed.
307

DIAGRAMMATIC GENERATION OF AN 

UNCONDITIONAL REALIZATION

Training image

Note that the progression is by patches, not by pixels.

Courtesy of A. Journel, Stanford University
308

BACK TO THE ELY WEST EXAMPLE

Training image

The training image must at least extend over the area 

of interest.
309

PAGE 272 EXERCISE ONCE MORE, NOW 

EMPLOYING FILTER SIMULATION
310

A. TWO REALIZATIONS FOR ELY WEST
311

COMPARISON TO EXHAUSTIVE SAMPLE

Exhaustive sample

81,796 measurements

First realization, based on

40 measurements plus 

the training image
312

B. SEMIVARIOGRAMS

FOR FIRST REALIZATION

Reproduction of exhaustive semivariogram is not as 

good as that of simulated annealing, but comparable 

to the overall approximation by sequential Gaussian 

simulation.
313

C. HISTOGRAMS

Filter simulation does a superior job approximating the 

histogram of the exhaustive sample because of the 

positive influence of the training image.
314

D. Q-Q PLOT

The distribution for the first realization is in close 

agreement with that of the exhaustive sample. 
315

E. ERRORS IN FIRST REALIZATION

ACCORDING TO EXHAUSTIVE SAMPLE

Errors are unbiased and the lowest among the simulation 

methods considered in this primer.
316

ERROR SEMIVARIOGRAMS

In common with all other methods, errors here are 

correlated. The semivariogram resembles the one for 

the attribute being modeled.
317

F. A VIEW OF ERRORS AND NODAL

STANDARD DEVIATIONS

Error Nodal standard deviation

5 5 

3 2 2 3 

Above a nodal standard deviation of 

1.5, the dispersion of errors tends to 

be the same regardless of the value of 

the nodal standard deviation.
318

18. RELIABILITY
319

UNCERTAINTY

From the outset, the 

quantitative 

assessment of the 

uncertainty associated 

with the resulting 

models has been an 

important objective of 

geostatistics.

 Basic to the analysis is having information about 

the probability distribution of the errors.

 This chapter explores the different ways in which 

geostatistics uses results from previous methods for 

uncertainty assessment.

x

z(

s,x)
320

THE RANDOM FUNCTION MODEL

By collecting the values of multiple realizations at the 

same node, it is possible in practice to model the 

concept of multiple variates over a sampling domain.

s

t

Diagrammatic idea

Actual implementation:

central node, Ely West
321

MAPS OF MOMENTS

Although nodal distributions have all the information for 

assessing uncertainty, customarily there are so many 

nodes that their individual inspection is impractical, 

creating, once more, a need to prepare summaries.

 The easiest solution is to map some consistent 

summary statistic for the nodal distributions. The most 

common options are:

A. Nodal mean, also known as E-type

B. Nodal standard deviation

Both options share the advantage of reducing all

information to one map, no matter how many the 

realizations.
322

A. NODAL MEAN MAP

• E-type maps tend to remind us of kriging maps, sharing

their smooth appearance.

• As in kriging, smoothing disqualifies E-type maps from

being realistic models.

• In the event that an E-type map is a good model, it is

always less demanding to prepare it by using kriging.

Ordinary kriging estimation Sequential Gaussian simulation

The link between E-type maps and kriging is an interesting

fact to remember. 
323

B. NODAL STANDARD DEVIATION MAPS

Ordinary kriging estimation Sequential Gaussian simulation

• The resemblance between standard deviations from

simulation and kriging can be good, but never as good as

the one for the means.

• Although nodal standard deviation provides a semiquantitative feeling about uncertainty, other summary

statistics are more efficient for summarizing the

information contained in the realizations.
324

NORMALITY

Kriging can be used in reliability 

estimation, but it is necessary to 

assume the form of the nodal 

distributions. The common 

practice is to assume normality, 

with a mean equal to the 

estimate, and a standard x z(

s,x)

* z

* σ

 When the errors indeed follow a multivariate normal 

distribution reasonably closely, results can be quite 

reliable, as we will see for the case of the UNCF sample.

 Whereas kriging may not be great at characterizing the 

true attribute, it can do a much better job, comparable to 

that of simulation, when it comes to assessing confidence 

intervals.

deviation equal to the kriging standard error, .

*

σ

,

* z
325

BEST DISPLAYS FOR SUMMARIZING 

UNCERTAINTY

A. The probability distribution for some global attribute

after reducing each realization to one single number,

which can be of the same nature or different from

the sampled attribute.

B. The kth percentile map displaying the threshold for

which there is a probability that the

true value is less than or equal to the threshold. The

units of this type of map are the same as those of

the attribute.

C. A probability map showing the chances that the true

value be equal to or less than a given threshold.

D. Interval probability map displaying fluctuations in the

chances that the true value be between a pair of

fixed values.

k /0 100 ≤≤ 1
326

A. SUMMARY DISTRIBUTIONS

ASSESSING GAS POTENTIAL IN 

GREATER NATURAL BUTTES, UTAH
327

B. BACK TO THE UNCF EXAMPLE.

10th PERCENTILE MAP

Ordinary kriging estimation

Prob[Z( ) s ≤ node _value( ) s ] = .0 10

Sequential Gaussian simulation
328

C. PROBABILITY MAP FOR TWO

THRESHOLDS USING SEQUENTIAL 

GAUSSIAN SIMULATION

Prob[Z( ) s ≤ 7872] Prob[Z( ) s ≤ 7966]
329

PROBABILITY MAP FOR TWO 

THRESHOLDS UNDER NORMALITY

Prob[Z( ) s ≤ 7872] Prob[Z( ) s ≤ 7966]
330

D. INTERVAL PROBABILITY

Ordinary kriging estimation

Prob[7872 Z( ) s ≤≤ 7966]

Sequential Gaussian simulation
331

BACK TO ELY: WHERE ARE TRUE VALUES 

LANDING AT THE NODAL DISTRIBUTIONS?

Each node shows 

the cumulative 

probability below 

the true value.

 Blanks denote 

that the true value 

is outside the 

nodal distribution, 

which is certainly 

not a good 

situation.
332

ACCURACY OF CONFIDENCE INTERVALS

Using the nodal 

distributions, for a 

given probability, the 

scatterplots show how 

frequently the true 

values are indeed 

below the quantile 

associated with the 

given probability.

 Note that these plots 

are global assessments 

comprising all nodes 

for all realizations.
333

19. NAVIGATION CHARTS
334

KRIGING OR SIMULATION?

• Perform kriging if the paramount criterion is individual

minimization of prediction error, if the data spacing is low

relative to the semivariogram range, or if the main interest is

an exploratory mapping of a poorly known attribute through

a deliberately smooth, thus simplified, version of reality.

• Go for simulation if it is not acceptable to assume normality

of the errors for reliability assessment, if correct modeling of

spatial continuity is more important than local accuracy, or if

there is interest in a stochastic assessment of some global

attribute, such as the examples in page 326.

Although there are no clear excluding rules, common criteria 

to choose one approach over the other are:

 Always remember that, unless you are in a terrible hurry, 

nothing prevents trying both kriging and simulation, as done 

with the Ely data, and compare the results.
335

RANKING OF FOUR METHODS

IN THE SPECIAL CASE OF ELY WEST

1 2

3 4

Exhaustive sample

Method

Reproduction

of

distribution

Reproduction

of

semivariog. Texture

Issues 

with

errors

Mean

square

error

Confidence

interval

Process.

time Total

1 Filter simulation 1 2 1 1 2 2 3 12

2 Seq. Gaussian simul. 2 2 2 1 3 4 2 16

3 Simulated annealing 2 1 3 4 3 1 4 18

4 Ordinary kriging 4 4 4 3 1 3 1 20
336

WHAT TYPE OF KRIGING?

• Use simple kriging if there is no trend and the population

mean is known, a rare situation.

• Run ordinary kriging if the mean is unknown and the

trend is approximately constant within the search

neighborhood (page 185). This is the predominant case.

Ordinary kriging is actually a special name for universal

kriging for a polynomial of degree 0.

• Otherwise, employ universal kriging proper of degree 1 or

2, whose main feature is extrapolation, which is still poor.

Selection of the type of point kriging is ruled by the 

knowledge of and type of mean (trend).

 Use the same criteria to decide the type of block kriging 

and cokriging.

 As the complexity of the trend model goes up from simple 

to universal kriging, so does the estimation variance.

 Use of semivariogram or covariance is immaterial 

because, for all practical purposes, they are equivalent.
337

WHAT KIND OF SIMULATION?

• Apply filter simulation if you already have a good idea

about the spatial variation of your attribute, good

enough to be able to provide a training image. The

quality of the results is crucially dependent on the

quality of the training image.

• In terms of usage, the popularity of the methods is

divided among the other two types covered here. Go

for simulated annealing if you want to stay away from

any form of normality assumption and the grid size is

less than 100,000 nodes.

• Otherwise, try sequential Gaussian simulation.
338

ACKNOWLEDGMENTS
339

Teaching is the best way to learn. Many thanks to all those participants of my courses who helped 

shape the present and earlier versions of these notes with their criticism. 

 I am grateful to Alexander Boucher (Stanford University) for productive discussions during my 

evaluation of SGeMS and for the care he took two years later to scrutinize this primer. Michael Friedel 

(U.S. Geological Survey), and Marlene Noble (U.S. Geological Survey) contributed to the polishing of 

this report with their thorough reviews. Marlene deserves extra credit for an earlier round of remarks 

while taking the latest offering of the course.

 U.S. Geological Survey editor Peter Stauffer did a careful job scrutinizing final details. Eric Morrissey 

(U.S. Geological Survey) inserted all the tags necessary to facilitate the use of the report by visually 

impaired readers.
340

SELECTED BIBLIOGRAPHY
341

Caers, J., 2011, Modeling uncertainty in the Earth Sciences. Chichester, UK, John Wiley & Sons 

Ltd, 229 pp.

Chilès, J.-P. and Delfiner, P., 1999, Geostatistics: modeling spatial uncertainty. New York, WileyInterscience, 695 pp.; 2nd ed. 2012, 734 pp.

Christakos, G., Olea, R.A., Serre, M.L., Yu, H.-L., and Wang, L.-L., 2005, Interdisciplinary public 

health reasoning and epidemiologic modelling. Berlin, Springer, 319 pp.

Deutsch, C.V. and Journel, A.G., 1992, GSLIB: geostatistical software library and user’s guide. 

New York, Oxford University Press, 340 pp. and 2 diskettes; 2nd ed., 1998, 369 pp. and CDROM.

Goovaerts, P., 1997, Geostatistics for natural resources evaluation. New York, Oxford University 

Press, 483 pp.

Isaaks, E.H. and Srivastava, R.M., 1989, An introduction to applied geostatistics. New York, Oxford 

University Press, 561 pp.

Mariethoz, G. and Caers, J., 2015, Multiple-point geostatistics. Chichester, UK, John Wiley & Sons 

Ltd, 264 pp.

Moore, D.S. and Notz, W.I., 2006, Statistics: concepts and controversies (6th ed.). New York, 

Freeman, 561 pp.

Olea, R.A., ed., 1991, Geostatistical glossary and multilingual dictionary. New York, Oxford 

University Press, 177 pp.

Olea, R.A., 1999, Geostatistics for engineers and earth scientists. Norwell, Mass., Kluwer Academic 

Publishers, 313 pp.

Remy, N., Boucher, A., and Wu, J., 2009, Applied geostatistics with SGeMS: a user’s guide: New 

York, Cambridge University Press, 284 pp. and CD-ROM.
342

INDEX
343

anisotropy 126, 187, 247, 249

annealing schedule 285–286

assumption 29, 243

block kriging 217

cokriging 231

ergodic 88

multivariate normality 154

ordinary kriging 159–160

simple kriging 133

simulated annealing 283, 288

stationarity 89

universal kriging 195

autocovariance 112

backtransformation 103, 154, 271, 

273

Bayes’s theorem 46

bias 90, 154, 246

binomial distribution 47–48

block kriging 212–224

box-and-whisker plot 38

causality 58, 62

Chebyshev’s theorem 32, 263

cluster 91–97, 119, 244, 301, 

304–306

coefficient

kurtosis 42

quartile skew 41

skewness 40, 49

coin flipping 45, 48, 51

cokriging 226–239

comparison

of block and ordinary kriging 

223–224

of cokriging and kriging 226, 

238–239

of covariance and 

semivariogram 112–113

of geostatistics and time-series 

analysis 6–7

of kriging and simulation 

322–323, 327–332, 334–335

of kriging methods 336

of ordinary and simple kriging 

156, 181–184

of simulation methods 335, 337

of statistics and geostatistics 5, 

87

conditional probability 45–46

conditional realization 274

conditional simulation 265

coregionalization 226, 232

model 234–236

correlation 246, 263, 281

coefficient 57, 60, 62, 245

spatial 118, 284

covariance 56, 110–111, 133, 

138–139, 153, 159, 197, 

215–216, 220–222, 299, 336

block-to-block 216, 219–221

point-to-block 215, 219–222

cross partition 304

cross semivariogram 235–236

crossvalidation 241–253

cumulative distribution 22, 65, 93, 

97, 99, 101, 103, 105, 267

dataset

Ely West 187

UNCF 20, 247

decile 34

declustering 91–97, 119, 152, 244

dependent variable 58

diagonal of matrix 143, 233

distribution 99, 283, 321

binomial 47–48

continuous 49, 51, 105

discrete 48, 51, 107

Gaussian 22, 42, 47, 49–50, 

64, 87, 154, 270, 324, 334

lognormal 49

normal 22, 42, 47, 49–50, 64, 

87, 154, 268, 270, 324, 334

standard normal 64, 102, 246

drift 114, 195–196, 201

ergodicity 88

error 59, 119, 152, 242, 244–249, 

252–253, 279–281, 295-297, 

315–317, 319, 324, 334

mean square 60, 128, 152, 

158, 219, 229

estimation variance 117, 134, 270, 

336

block kriging 219, 224

cokriging 233

simple kriging 134, 138–142, 

146, 151, 153, 183–184
344

ordinary kriging 161–170, 174, 

184

universal kriging 198, 200, 

207–208

estimator 24, 90, 128, 241, 270, 324

block kriging 214, 223

cokriging 229

covariance 56

ordinary kriging 157, 183

semivariogram 114

simple kriging 130, 142, 150, 

153–154, 156, 183

universal kriging 194, 207–208

E-type map 321–322

expected value 51–52, 56, 105,

133, 135, 159, 217

experimental semivariogram 

114–116, 119, 201–202

explanatory variable 58–59

exponential semivariogram 117, 236

extrapolation 7

extreme value 33, 37–38, 59, 82

filter 301, 303, 306

filter simulation 299–317, 337

fractile 34

frequency table 19–20, 22

Gandin, L.S. 9

Gaussian distribution 22, 42, 47, 49, 

50, 64, 87, 154, 270, 324, 334

Gaussian semivariogram 117, 187, 

249

genesis 126

grid 70, 76–80, 189, 268–269, 284, 

306

GSLIB 68, 79, 243, 268, 289

simulated annealing 289

histogram 21, 24, 40, 47, 82–84, 

91, 272, 284

historical remark 8–9

independence 45, 48, 87, 89, 

152–153, 156, 288

independent variable 58–59

indicator 105–107

interpolation 7, 63, 152

interquartile range 36–37, 41

intrinsic hypothesis 160, 166, 195

inverse matrix 143

isotropy 119, 236

joint probability 45–46, 89, 268

k-means 304

Kolmogorov, A.N. 9

Kolmogorov-Smirnov test 100–101,

253

Krige, D. 9

kriging 128–129, 237–238, 

243–244, 255, 322–324, 336

block 212–224

multigaussian 154, 267, 270

ordinary 156–191, 193, 226, 

247–251, 256–263, 336

simple 130–154, 156, 

181–184, 193, 226, 336

universal 193–210, 226, 336

with trend 193

kurtosis 42

Lagrange method 171–173, 201, 

232

lognormal distribution 49

lognormal transformation 104

lurking variable 62

main diagonal of matrix 143

Mathern, B. 9

Matheron, G. 8

matrix 143–144, 176–177, 

199–200, 227–232, 235–236

inverse 143

main diagonal 143, 233

transposed 143

maximum 33, 38, 63

mean 26, 29, 31–32, 38, 42, 49, 

51–52, 64, 101, 112–114, 133, 

153–154, 156–157, 159–160, 

193, 195, 230, 245, 270, 

321–324,336

square error 60, 128

measure

centrality 24–25

dispersion 24, 30

shape 24, 39

median 27, 29, 34–35, 37, 97, 245

minimum 33, 38, 63, 139, 170, 185

mode 28, 259

moment 52, 89, 235, 299–301, 321

multigaussian kriging 154, 267, 270

node 76, 269–272, 281, 283–285, 

288, 320–321, 323, 332
345

noise 61, 257

normal distribution 22, 42, 47, 49, 

50, 64, 87, 154, 268, 270, 324, 

334

normal equation

block kriging 218

cokriging 232, 234

ordinary kriging 173

simple kriging 131, 142

universal kriging 197

normal score 102–103, 154, 269

nugget effect 118, 236, 244, 246

object 70–81, 188–189

omnidirectional semivariogram 119

ordinary kriging 156–191, 193, 226, 

247–251, 256–263, 336

orthogonality 152

outlier 37–38, 246

P-P plot 65–66

pattern recognition 301

percentile 34, 38, 245, 325

permissible model 117–118, 139,

170, 234–235

plot

P-P 65–66

Q-Q 63, 66, 272

population 15–16, 24, 31, 294, 336

positive definite 139, 234

probability 44–46, 50, 325

conditional 45–46

cumulative 22, 64, 99, 102, 105

density function 47–50

distribution 47–50, 87, 284, 

319–321, 324–326

joint 45–46

probability-probability plot 65–66

project 81, 188

property

ordinary kriging 184

simulated annealing 288

simple kriging 152–153

Q-Q plot 63, 66, 272

quantile 34, 245, 304

quantile-quantile plot 63, 66, 272

quartile 34–36, 38

random function 86–87, 89, 105,

133–135, 159–160, 195, 213–214,

217, 226–227, 267, 269, 320

random variable 16, 59, 64, 86–87, 

89, 226

range 113, 118, 126, 236, 249–250, 

255, 334

realization 86–88, 133, 159–160, 

195, 217, 231, 264–265, 

268–269, 271–272, 274, 

283–285, 306–307, 320–323, 332

regressed variable 58

regression 58–61, 128, 245–246

regressor 58–59

reliability 319–332

residual 135, 195–197, 201

response variable 58

robustness 29, 31, 33, 36, 41, 61, 

91

sample 15, 17, 24, 26–28, 31–34, 

56, 58–63, 72, 87, 91, 99–102, 

119, 124, 133, 153, 159–160, 

217, 231, 242, 267, 269–272, 

283–284, 287–288, 325

Ely West 187

size 17, 63, 242, 247, 267, 272, 

300

UNCF 20, 247

saving

grid display 191

grid values 79

histogram display 84

object 81, 188

parameter file 190

semivariogram 125

scaling 63, 66, 153, 244

scatterplot 55, 332

screen effect 149, 152, 185

search 185–186, 247–248

semivariogram 109–126, 170, 190, 

195, 201–203, 220, 235–236, 

244, 246–247, 249, 255, 269, 

272, 283–285, 288, 299, 336

experimental 114–116, 119, 

201–202

exponential model 117, 236

Gaussian model 117, 187, 249

omnidirectional 119

pure nugget effect 118, 236, 

244, 246

spherical model 117, 236
346

sequential Gaussian simulation 

267–281, 322–323, 330–332, 337

SGeMS 11, 13, 68–84

availability 68

data 70–75

filter simulation 309

grid 70, 76–80

histogram 82–84

image save 84, 191

loading 74–75, 188

main menu 69

object 70–82, 188–189

ordinary kriging 188–191

project 81, 188

save 79, 81, 84, 125, 188, 

190–191

semivariogram 120–125

sequential Gaussian simulation 

273

universal kriging 209–210

shift 63, 66

sill 113, 244, 250

simple kriging 130–154, 156, 

181–184, 193, 226, 336

simplicity 133

simulated annealing 283–297, 

331–332, 335, 337

simulation 78, 255, 264–265, 337

conditional 265

filter 299–317, 337

sequential Gaussian 267–281, 

322–323, 330–332, 337

simulated annealing 283–297, 

331–332, 335, 337

unconditional 265

skewness 40–41, 49

smoothing 255, 257–259, 264, 334

spatial continuity 110, 187

spatial correlation 5, 100, 110, 118, 

128, 284

spherical semivariogram 117, 236

standard deviation 32, 49, 57, 64, 

101, 245, 263, 272, 281, 321, 

297, 317, 323

standardization 64–65

standard normal distribution 64, 

102, 246

stationarity 89, 113–114, 133, 154, 

159, 166–167, 195, 212, 217, 231

statistic 24, 66, 84, 100, 241, 

245–246, 300, 321, 323

statistical dependence 45, 58, 243

statistical independence 45, 48, 87

structural analysis 109, 226, 

234–236

support 87, 212–216

theorem

Bayes’s 46

Chebyshev’s theorem 32, 263

time series analysis 6–7

training image 302, 305–308, 337

transformation 64, 99–107, 119, 

269, 271

Indicator 105–107

lognormal 104

normal 102–103, 269

transposed matrix 143

trend 114, 119, 195–196, 201, 209

unbiasedness 90, 152, 154, 156, 

261, 279, 294, 315

uncertainty 4, 6, 86, 193, 264, 319, 

321, 323, 325

unconditional simulation 265

universal kriging 193–210, 226, 336

variable

dependent 58

explanatory 58–59

independent 58–59

lurking 62

random 16, 59, 64, 86–87, 89, 

226

regressed 58

response 58

variance 31–32, 36, 52, 64, 113, 

117, 119, 124, 136, 164, 270

variate 16, 267–268, 320

visitation schedule 269, 306

weight 130–132, 139–140, 142, 146,

148–149, 153, 156–158, 170, 174,

180, 208, 218, 229, 303

Wiener, N. 9

z-score 246, 253