Chapter Contents Previous Next
 The SIM2D Procedure

Example 58.1: Simulation

Continuing with the coal seam thickness example from the "Getting Started" section, this example asks a more complicated question. This question is economic in nature, and the (approximate) answer requires the use of simulation.

Simulating a Subregion for Economic Feasibility

The coal seam must be of a minimum thickness, called a cutoff value, for a mining operation to be profitable. Suppose that, for a subregion of the measured area, the cost of mining is higher than the remaining areas due to the geology of the overburden. This higher cost results in a higher thickness cutoff value for the subregion. Suppose also that it is determined from a detailed cost analysis that at least 60 percent of the subregion must exceed a seam thickness of 39.7 feet for profitability.

How can you use the SRF model ( and Cz(s)) and the measured seam thickness values Z(si),i = 1, ... ,75 to determine, in some approximate way, if at least 60 percent of the subregion exceeds this minimum?

Spatial prediction does not appear to be helpful in answering this question. While it is easy to determine if a predicted value at a location in the subregion is above the 39.7 feet cutoff value, it is not clear how to incorporate the standard error associated with the predicted value. The standard error is what characterizes the stochastic nature of the prediction (and the underlying SRF). It is clear that it must be included in any realistic approach to the problem.

A conditional simulation, on the other hand, seems to be a natural way of obtaining an approximate answer. By simulating the SRF on a sufficiently fine grid in the subregion, you can determine the proportion of grid points in which the mean value over realizations exceeds the 39.7 feet cutoff and compare it with the 60 percent value needed for profitability.

It is desirable in any simulation study that the quantity being estimated (in this case, the proportion exceeding the 39.7 feet cutoff) not depend on the number of simulations performed. For example, suppose that the maximum seam thickness is simulated. It is likely that the maximum value increases as the number of simulations performed increases. Hence, a simulation is not useful for such an estimate. A simulation is useful for determining the distribution of the maximum, but there are general theoretical results for such distributions, making such a simulation unnecessary. Refer to Leadbetter, Lindgren, and Rootzen (1983) for details.

In the case of simulating the proportion exceeding the 39.7 feet cutoff, it is expected that this quantity will settle down to a fixed value as the number of realizations increases. At a fixed grid point, the quantity being compared with the cutoff value is the mean over all simulated realizations; this mean value settles down to a fixed number as the number of realizations increases. In the same manner, the proportion of the grid where the mean values exceed the cutoff also becomes constant. This can be tested using PROC SIM2D.

A crucial, nonprovable assumption in applying SRF theory to the coal seam thickness data is that the values Z(si), i = 1, ... ,75 represent a single realization from the set of all possible realizations consistent with the SRF model ( and Cz(h)). A conditional simulation repeatedly produces other possible simulated realizations consistent with the model and data. However, the only concern of the mining company is with this single unique realization. It is not concerned with similar coal fields to be mined sometime in the future; it may never see another coal field remotely similar to this one, or it may not be in business in the future.

Hence the proportion found by generating repeated simulated realizations must somehow relate back to the unique realization that is the coal field (seam thickness). This is done by interpreting the proportion found from a simulation to the spatial mean proportion for the unique realization. The term "spatial mean" is simply an appropriate integral over the fixed (but unknown) spatial function z(s). (The SRF is denoted Z(s); a particular realization, a deterministic function of the spatial coordinates, is denoted z(s).)

This interpretation requires an ergodic assumption, which is also needed in the original estimation of Cz(s). Refer to Cressie (1993, pp. 53 -58) for a discussion of ergodicity and Gaussian SRFs.

Implementation Using PROC SIM2D

The subregion to be considered is the southeast corner of the field, which is a square region with length 40 distance units (in thousands of feet). PROC SIM2D is run on the entire data set for conditioning, while the simulation grid covers only this subregion. It is convenient to be able to vary the seed, the grid increment, and the number of simulations performed. The following macro implements the computation of the percent area exceeding the cutoff value by using the seed, the grid increment, and the number of simulated realizations as macro arguments.

The data set produced by PROC SIM2D is transposed so that each grid location is a separate variable. The MEANS procedure is then used to average the simulated value at each grid point over all realizations. It is this average that is compared to the cutoff value. The last DATA step does the comparison and determines the percent of the grid locations that exceed this cutoff value and writes the results to the listing file in the form of a report.

The macro is first invoked with a relatively coarse grid (grid increment of 10 distance units) and a small number of realizations (5). The next invocation uses a finer grid and 50 realizations, and the final invocation uses the same grid increment and 500 realizations. Each time, the macro is invoked with a different seed. The simulations indicate that around 87 percent of the subregion exceeds the cutoff value.

The number of grid points in the simulation increases with the square of the decrease in the grid increment, leading to long CPU processing times. Increasing the number of realizations results in a linear increase in processing times. Hence, using as coarse a grid as possible allows more realizations and experimentation with different seeds.

   /*- Set the covariance model parameters and cutoff value -*/
%let cc0=7.5;
%let aa0=30.0;
%let form=gauss;
%let cut=39.7;

data thick;
input east north thick @@;
datalines;
0.7  59.6  34.1   2.1  82.7  42.2   4.7  75.1  39.5
4.8  52.8  34.3   5.9  67.1  37.0   6.0  35.7  35.9
6.4  33.7  36.4   7.0  46.7  34.6   8.2  40.1  35.4
13.3   0.6  44.7  13.3  68.2  37.8  13.4  31.3  37.8
17.8   6.9  43.9  20.1  66.3  37.7  22.7  87.6  42.8
23.0  93.9  43.6  24.3  73.0  39.3  24.8  15.1  42.3
24.8  26.3  39.7  26.4  58.0  36.9  26.9  65.0  37.8
27.7  83.3  41.8  27.9  90.8  43.3  29.1  47.9  36.7
29.5  89.4  43.0  30.1   6.1  43.6  30.8  12.1  42.8
32.7  40.2  37.5  34.8   8.1  43.3  35.3  32.0  38.8
37.0  70.3  39.2  38.2  77.9  40.7  38.9  23.3  40.5
39.4  82.5  41.4  43.0   4.7  43.3  43.7   7.6  43.1
46.4  84.1  41.5  46.7  10.6  42.6  49.9  22.1  40.7
51.0  88.8  42.0  52.8  68.9  39.3  52.9  32.7  39.2
55.5  92.9  42.2  56.0   1.6  42.7  60.6  75.2  40.1
62.1  26.6  40.1  63.0  12.7  41.8  69.0  75.6  40.1
70.5  83.7  40.9  70.9  11.0  41.7  71.5  29.5  39.8
78.1  45.5  38.7  78.2   9.1  41.7  78.4  20.0  40.8
80.5  55.9  38.7  81.1  51.0  38.6  83.8   7.9  41.6
84.5  11.0  41.5  85.2  67.3  39.4  85.5  73.0  39.8
86.7  70.4  39.6  87.2  55.7  38.8  88.1   0.0  41.6
88.4  12.1  41.3  88.4  99.6  41.2  88.8  82.9  40.5
88.9   6.2  41.5  90.6   7.0  41.5  90.7  49.6  38.9
91.5  55.4  39.0  92.9  46.8  39.1  93.4  70.9  39.7
94.8  71.5  39.7  96.2  84.3  40.3  98.2  58.2  39.5
;

%macro area_sim(seed=,nr=,ginc=);

%let ngrid=%eval(40/&ginc+1);
%let tgrid=%eval(&ngrid*&ngrid);

proc sim2d data=thick outsim=sim1;
simulate var=thick numreal=&nr seed=&seed
scale=&cc0 range=&aa0 form=&form;
mean 40.14;
coordinates xc=east yc=north;
grid x=60 to 100 by &ginc
y=0 to 40 by &ginc;
run;

proc transpose data=sim1 out=sim2 prefix=sims;
by _iter_;
var svalue;
run;

proc means data=sim2 noprint n mean;
var sims1-sims&tgrid;
output out=msim n=numsim mean=ms1-ms&tgrid;
run;

/*- Determine the percentage of sites exceeding cutoff -*/
data _null_;
file print;
array simss ms1-ms&tgrid;
set msim;

/*- Loop over the grid sites to test cutoff  -*/
cflag=0;
do ss=1 to &tgrid;
tempv=simss[ss];
if simss[ss] > &cut then do;
cflag + 1;
end;
end;

area_per=100*(cflag/&tgrid);
put // +5 'Conditional Simulation of Coal Seam'
' Thickness for Subregion';
put / +5 'Subregion is South-East Corner 40 by 40'
' distance units';
put / +5 "Seed:&seed" +2 "Grid Increment:&ginc";
put / +5 "Total Number of Grid Points:&tgrid" +2
"Number of Simulations:&nr";
put / +5 "Percent of Subregion Exceeding Cutoff of
&cut ft.:"
+2 area_per 5.2;
run;
%mend area_sim;

%area_sim(seed=12345,nr=5,ginc=10);
%area_sim(seed=54321,nr=50,ginc=1);
%area_sim(seed=655311,nr=500,ginc=1);


Output 58.1.1: Conditional Simulation of Coal Seam Thickness

  Conditional Simulation of Coal Seam Thickness for Subregion Subregion is South-East Corner 40 by 40 distance units Seed:12345 Grid Increment:10 Total Number of Grid Points:25 Number of Simulations:5 Percent of Subregion Exceeding Cutoff of 39.7 ft.: 80.00 

  Conditional Simulation of Coal Seam Thickness for Subregion Subregion is South-East Corner 40 by 40 distance units Seed:54321 Grid Increment:1 Total Number of Grid Points:1681 Number of Simulations:50 Percent of Subregion Exceeding Cutoff of 39.7 ft.: 87.33 

  Conditional Simulation of Coal Seam Thickness for Subregion Subregion is South-East Corner 40 by 40 distance units Seed:655311 Grid Increment:1 Total Number of Grid Points:1681 Number of Simulations:500 Percent of Subregion Exceeding Cutoff of 39.7 ft.: 87.57 

 Chapter Contents Previous Next Top