Generating confidence intervals for spatial simulations –

determining the number of replications for spatial

terminating simulations
1
Robert M. Itami, 2Darrel Zell, 3Frank Grigel, 4Randy Gimblett
1
GeoDimensions Pty Ltd, E-Mail: bob.itami@geodimensions.com.au, 2Parks Canada, E-Mail:
Darrell.Zell@pc.gc.ca, 3Parks Canada, E-Mail: Frank.Grigel@pc.gc.ca, 4University of Arizona, E-
Mail: gimblett@ag.arizona.edu

Keywords: Spatial Simulation, Simulation Statistics, Simulation Output Analysis

EXTENDED ABSTRACT

A standard problem in the analysis of outputs from

terminating simulations is the need to determine
the number of replications needed to construct
confidence intervals for performance indicators
from the simulation (Law and Kelton, 2000). In
traditional industrial applications of simulation
such as manufacturing and queuing simulations a
single mean for each performance indicator is all
that is needed. In spatial simulations however, the
problem is more complex as performance
indicators can vary spatially as in the case of travel
simulations where performance indicators for each
destination must be analysed. This paper presents
three alternative methods recommended in the
simulation literature for determining the number of
replications required to obtain confidence intervals
based for a given alpha level and user defined
confidence interval half width or relative
preceision. The problem of measuring multiple
performance indicators is addressed with a short
discussion of the Bonferroni Correction. These
methods are then adapted to spatial simulations
using a travel simulation for Banff, Yoho,
Kootenay and Jasper National Parks as an
example. Outputs for daily link Use and daily link
encounters are examined applying different values
for absolute accuracy and relative precision.
Conclusions are then drawn on the relationship
between the sensitivity of performance indicators
to random variables in the simulation model and
the specification of absolute accuracy and relative
precision for spatial dynamic simulation models.

141
1. INTRODUCTION
S 2 ( n)
Spatial simulation models use random numbers to In equation 1 the expression: t1−α / 2 is
generate input variables such as arrival times,
n
durations at destinations, and trip selection. referred to as the confidence interval half width. If
Because of this, it is not recommended to draw any this value is less than the user specified accuracy
conclusions from the output of a single replication after the initial n replications for the “short run”
of a simulation model since it represents only one then there is no need for further replications.
realization of a stochastic process. This is a However if this value is larger then the user
common characteristic of all simulation models, specified accuracy, then n can be estimated using
however spatial simulations have added Equation 2.
complexity because output measures such as Equation 2.
number of visits or average duration at nodes are
measured across many different links and nodes. ⎡ ⎛ h ⎞2 ⎤
To solve the problem of variability in simulation n = Round ⎢ n × ⎜ * ⎟ ⎥
*

outputs from random inputs, the method of

independent replications (Goldman 1992) is used.
⎣⎢ ⎝ h ⎠ ⎦⎥
In this method one observation per replication is *
generated (usually the mean value of the WHERE n is the estimated number of
performance indicator). In the case of spatial replications needed
simulations you would have one observation per h is the half width from the sample run
link or node of interest. The simulation is
replicated a number of times to generate h* is the desired half width or absolute accuracy
confidence intervals that meet a certain reliability specified by the user.
(usually 0.90 or 0.95) and a given accuracy Law and Kelton (2000, p512) suggest a
(measured in the same units as the performance modification of the above estimate in a method
indicator). The problem is how to determine how they call the iterative method. In this case the
many replications are required to meet the user number of replications is increased by 1 each time
specified reliability and accuracy. and the confidence interval is recomputed after
each iteration until the desired accuracy is
2. STATISTICAL METHODS FOR achieved. This method assumes that the
DETERMINING THE NUMBER OF population variance will not change (appreciably)
REPLICATIONS FOR TERMINATING as the number of replications increase.
SIMULATIONS
Equation 3.
The method of independent replications requires
the model is run for a “small” number of ⎧⎪ S 2 (n) ⎫⎪
replications. In the case of probabilistic n ( β ) = min ⎨i ≥ n : ti −1,1−α / 2
*
a ≤ β⎬
i
⎩⎪ ⎭⎪
simulations, “small” may mean 10 to 15
replications.
The next step is to calculate the (1-α) confidence
interval using equation 1. (See Centeno and Reyes,
1998 and Law and Kelton, 2000 pp 253-259) a n* ( β )
WHERE is the estimated number of
Equation 1. replications needed with absolute accuracy
β
S 2 (n)
X ( n ) ± t1− α / 2 n is the number of replications from the “short
n run”

ti −1,1−α / 2
WHERE X is the mean of the performance is the is the (1-α) percentile of the t-
indicator for the current replication student distribution with i-1 degrees of
freedom
t1−α / 2 is the is the (1-α) percentile of the t-student
distribution with n-1 degrees of freedom S 2 (n) is the sample variance from n replications

S 2 (n) is the sample variance i is the iteration (greater than n)

A third method described by Law and Kelton
(2000, p. 513) uses a measure called “relative
accuracy”. Relative accuracy is the Confidence

142
Interval Half Width from Equation 1 divided by accept that some indicators may not be used
X . is an estimate of the actual relative error. In reliably in drawing conclusions from the
this method the user specifies a desired relative simulation.
error λ' . 4. DETERMINING THE NUMBER OF
Equation 4. REPLICATIONS FOR SPATIAL
SIMULATIONS
⎧ ti −1,1−α / 2 S 2 ( n ) / i ⎫
⎪ ⎪ In traditional industrial applications of simulation
nr* ( λ ) = min ⎨i ≥ n : ≤ λ' ⎬ such as manufacturing and queuing simulations a
⎪⎩ X (n) ⎪⎭ single mean for each performance indicator is all
that is needed. In spatial simulations however, the
WHERE nr ( λ ) is the estimated number of
*
problem is more complex as performance
replications needed with relative accuracy indicators can vary spatially as in the case of travel
λ simulations where performance indicators for each
destination must be analysed simultaneously.
N is the number of replications from the “short Essentially, the approach required is to apply the
run” same statistical methods described in section 2 to
each and every location in the spatial simulation
λ' is the user-specified relative accuracy where performance indicators are to be measured.
For instance, in travel simulations we may be
t1−α / 2 is the is the (1-α) percentile of the t-student
interested in:
distribution with n-1 degrees of freedom
• the total visits per destination,
S 2 (n) is the sample variance from n replications
• average visit duration per destination, and
i is the iteration (greater than n) • average queuing times at parking
facilities.
3. MULTIPLE PERFORMANCE
INDICATORS AND THE BONFERRONI In this case, we have three performance indicators.
CORRECTION Suppose our network has 10 destinations we wish
to evaluate. We must first determine the alpha
In a typical spatial simulation, there is normally level we wish to test for. If the overall alpha is
more than one performance indicator being 0.10 then, according to the Bonferroni Correction
measured. If simulations are viewed as we must use an alpha level of 0.10/3 or 0.03.
experiments where we are testing hypotheses about
the system under study, then the alpha levels of the Next we decide which method will be used to
statistical tests applied in the previous section must determine the number of replications either by
be adjusted using the Bonferroni Correction specifying the confidence interval half width for
(MathWorld, 2005). The Bonferroni Correction is equations 2 and 3 or by specifying the relative
used when several tests are being performed accuracy for equation 4.
simultaneously. Where a given alpha level for a The simulation is then replicated for a “short run”
single performance indictor may be appropriate, it of say, 10 replications and the outputs from each of
is not for the set of all comparisons. The simplest the 10 destinations for the 3 performance
form of the Bonferroni Correction is to take the indicators are gathered. Using this output we then
desired alpha level and divide by the number of apply the corresponding method (equations 2, 3 or
performance indicators being tested. Thus, if the 4) for each of the 10 destinations for the 3
desired alpha level is 0.10 and there are five performance indicators using an alpha level of
performance indicators, the adjusted alpha level 0.03.
would be 0.10/5 or 0.02 for each test.
Law and Kelton (2000) suggest another approach 5. EXAMPLE: PATTERN OF USE
in which the sum of the alpha levels for each test SIMULATION FOR CANADAS 4
equals the desired alpha level, suggesting that each MOUNTAIN PARKS
performance indicator can have a unique alpha A probabilistic spatial simulation of tourist flows
level in the statistical test. They also suggest that in Banff, Jasper, Kootenay and Yoho National
more than 10 performance indicators is impractical Parks in the Canadian Rocky Mountains was
and that, given the stochastic nature of simulation developed using RBSim software. Data from
it may be impractical to meet statistical independent travelers was collected in 2003 in a
requirements for all performance indicators survey conducted by Parks Canada over the entire
simultaneously and that one must just have to year. The 2003 Mountain Parks Study provided

143
the trip itinerary data for a trace simulation. Out of 1. The week of arrival is determined (1 through
the 13,373 first time arrivals to the Mountain 53)
parks, 9348 respondents agreed to complete the
2. The week day is determined (Monday through
questionnaire. Out of the 9348 respondents, 2383
Sunday).
questionnaires were returned, resulting in a 25.5
percent return rate. Out of those 2383 3. The total number of arrivals for the week for
questionnaires, 1982 respondents actually the current gate is determined from the
completed the trip diary component (necessary to weighted traffic count data for the week.
develop the trip itineraries for the trace
simulation), resulting in a 21.2 percent return rate. 4. The weekly arrival distribution for the current
week is selected
The survey trips represent a “single realization” of
the full variation of travel patterns in the four 5. The total number of arrivals for the week is
parks. To use a pattern of use simulation for multiplied by each day in the weekly arrival
management, a probabilistic model of travel distribution. This determines the number of
patterns must be developed from analysis of arrivals for each day of the week.
survey trips. Trip itineraries vary within the 6. The hourly distribution is then selected
constraints of seasonal patterns, which will
constrain certain activities and access to 7. The exact time of arrival for the current trip is
destinations because of weather. Trips must be selected from the hourly arrival distribution by
grouped or clustered in order to classify trip generating a random number and calculating
itineraries according to seasonal variations in the exact minute of arrival by piecewise linear
pattern of use. In order to do this, a number of interpolation from the hourly arrival
different cluster analysis techniques were applied distribution.
to the 1620 survey trips using a number of 8. Finally, the trip clusters that occur in this
different clusters and against different definitions week are looked up and then all trips that fall
of a “season”. The strongest clustering results in these clusters for the current entry gate are
came from a technique called “kmeans”. In this selected to create the pool of trips. Each trip
technique, the analyst must nominate the number in the pool has an equal chance of being
of clusters. 3 and 4 clusters were tested with 3 selected. A random number is generated and a
clusters differentiating the winter season trips most single trip is selected from the pool
distinctly. The other trips were clustered on the
basis of the type of activity: either active (hiking, 9. The trip selection is now complete and the
biking etc) or inactive. process is repeated for the next trip until all
trips for the current day are selected. If all
Based on analysis of the survey trips looking at trips are scheduled, the process progresses to
sample size during each week of the year and the the next day of the simulation and the process
results of the cluster analysis, two periods of the is repeated until the last day of the simulation.
year were selected to develop a probabilistic
simulation from survey trip itineraries and traffic The first two weeks of the simulation outputs must
counts. The winter period selected is the month of be discarded because the simulation starts with no
January 2003. The summer period selected is 4 visitors so any outputs during an initial “warm up”
weeks beginning June 23, 2003 and ending July period are not representative of system behavior.
24, 2003. We report on the January results in this To correct this problem, the simulation is allowed
paper. to run until the system comes to full capacity. In
this case the period that was selected was two
The next step in the development of the weeks since the longest visits are 11 days. The
probabilistic simulation is to randomly assign, for warm up period gives the simulation time to
each day of the simulation, the correct number of populate the system to capacity before collecting
trip arrivals to each day of the simulation and then statistics from the simulation. For the winter
scheduling the exact minute of arrival based on the simulation 1 week of simulation was obtained for
hourly arrival distribution. A trip itinerary must be analysis.
randomly selected from a pool of trips specified
from the cluster analysis for the entry gate and Two performance indicator were measured for
time of year. This is done for each gate by links: Link Use and Link Encounters. Link use is
proceeding from the first day of the simulation to a frequency count for the number of parties
the last day of the simulation according to the visiting each link for each day of the simulation.
following procedure: Link Encounters is the number of direct contacts
between parties along a link.
For each day,

144
There were 658 links visited and 493 links with confidence interval half width with the mean. This
encounters recorded in the simulation. Since daily shows the value of this measure to indicate the
link use and link encounters were generated, there amount of variation between simulation runs.
were 3799 link use days and 2083 link encounter
days recorded with confidence intervals calculated Date Link Use Link
for each. Encounters
Link Jan- CI Half Rel CI Half Rel
For purposes of understanding the implications of 03 Prec Prec
the three methods for estimating the number of 116 15 14.82 0.04 2.37 0.83
replications needed to obtain user specified
116 16 11.02 0.03 4.14 1.07
measures of confidence and reliability from the
simulation, the results of two links are shown in 116 17 18.76 0.05 2.10 0.57
the tables below. Confidence intervals for the 7 116 18 16.00 0.04 3.90 1.05
replications for daily link use and daily link 116 19 15.33 0.03 5.01 0.80
encounters were calculated using different alpha
values and different user-specified confidence 126 15 13.52 0.02 11.92 0.36
half-width values to investigate the impact these 126 16 14.49 0.02 23.09 0.47
values have on the three methods for estimating 126 17 20.99 0.04 26.41 0.82
the number of replications.
126 18 12.05 0.02 54.57 0.63
Date Link Use Link 126 19 13.02 0.02 21.10 0.84
Encounters
Link Jan 03 Mean StDev Mean StDev Table 2. Confidence Half Intervals and
Relative Precision for Link Use and Link
116 15 380.71 20.18 2.86 3.23
Encounters, Alpha = 0.10
116 16 387.00 15.00 3.86 5.64
116 17 401.57 25.55 3.71 2.86
Date Link Use Link Encounters
116 18 432.86 21.79 3.71 5.31
Link Jan Eq. Eq. Eq. Eq. Eq. Eq.
116 19 448.71 20.88 6.29 6.82 03 2 3 4 2 3 4
126 15 559.14 18.41 32.86 16.23 116 15 1583 1113 7 39 31 156
126 16 580.14 19.73 48.71 31.44 116 16 850 616 7 120 88 260
126 17 595.43 28.58 32.29 35.96 116 17 2465 1784 7 31 25 74
126 18 679.00 16.41 86.71 74.31 116 18 1792 1297 7 106 79 249
126 19 670.71 17.73 25.14 28.74 116 19 1646 1192 7 176 128 144
Table 1. Means and Standard Deviations for 126 15 1280 927 7 994 720 32
Link Use and Link Encounters 126 16 1470 1064 7 3731 2701 152
Table 1 shows the means and standard deviations 126 17 3084 2232 7 4883 3534 152
for link use and link encounters for two links. 126 18 1016 736 7 20848 15089 91
Link 116 and Link 26. Link 116 is typified by low
number of encounters as compared to Link 126. 126 19 1187 860 7 3118 2257 159
Note that standard deviations for link encounters
Table 3. Estimated replications using
are relatively high compared to link use. This is
Equations 2, 3 and 4 for Link Use and Link
because link encounters are much more sensitive
Encounters, Alpha = 0.10, User CI Half Width =
to random variations than link use since small
1, User Relative Precision = 0.15
changes in arrival times and volumes can change
the number of encounters on a link.
Table 3 shows the estimated number of
Table 2 shows the confidence interval half width
replications for link use using the three different
(CI Half) and relative precision (Rel prec) from
methods described in this paper. For equations 2
confidence intervals calculated with an alpha of
and 3 an absolute accuracy (the user specified
0.10 (90% confidence). Law and Kelton (2000)
confidence interval half width) of 1 is used. For
suggest that an alpha of 0.10 is a reasonable level
Equation 4 a user specified relative accuracy of
of confidence given that random number are used
0.15 is used. For link use equation 4 shows that
to generate the performance indicators. Note the
there is no need for further replications whereas
large values for relative precision for link
equations 2 and 3 show that many more
encounters as compared to link use. Remember
replications are needed to reduce the confidence
that the relative precision is the ratio of the
interval to a half width of 1. This is instructive,

145
because it may indicate that the criteria for Date Link UseLink
absolute accuracy may be unrealistic given the Encounters
high standard deviations for link use. Note that Link Jan- Eq. Eq. Eq. Eq. Eq. Eq.
equation 3 (the iterative method) shows a 03 2 3 4 2 3 4
consistently lower estimate for the number of 116 15 24 19 7 7 7 221
replications required to generate confidence
116 16 13 12 7 7 7 370
intervals with a half width of 1.
116 17 39 28 7 7 7 104
In Table 3 the number of replications required for
link encounters is greater for Equation 4 than 116 18 28 21 7 7 7 354
equations 2 and 3. This is because there are 116 19 26 20 7 7 7 204
relatively few encounters per link (see Table 1) for 126 15 20 16 7 16 13 44
link 116. However for relative accuracy, we see 126 16 23 18 7 59 41 74
that we would require anywhere from 74 to 260
replications in order to reach a relative precision of 126 17 49 34 7 77 53 215
0.15 for link 116. This reflects the high standard 126 18 16 13 7 331 200 128
deviations as compared to the means for link 126 19 19 15 7 49 35 226
encounters for link 116 and the resulting high
relative precision. For link 126 there were many Table 5. Estimated replications using
more encounters per day than link 116. The result Equations 2, 3 and 4 for Link Use and Link
is that we have a larger confidence interval half Encounters, Alpha = 0.05, User CI Half Width =
width, which requires many more replications to 10, User Relative Precision = 0.15
achieve our user specified absolute accuracy of 1.
This shows how important it is to carefully select
appropriate standards of accuracy for each In Table 5 we see a significant drop in the number
performance indicator. of replications required from equations 2 and 3 for
link use because we have now increased the user
To examine the impact of the Bonferroni specified half width from 1 to 10. Equation 4 for
Correction, we now assume that in order to link use shows no need for more replications
achieve an overall alpha of 0.10 for the simulation, because of the low relative precision values
we need to estimate replications using an alpha obtained as the result of the high use levels and
level of 0.05 for each of our two performance small confidence interval half widths. For link
indicators. encounters no further replications are required
from equations 2 and 3 because of the low number
Date Link Use Link
of encounters and the high user specified
Encounters
confidence interval half width of 10. Equation 4
Link Jan-03 CI Half Rel CI Half Rel for node 116 still shows a large number of
Prec Prec replications due to the high standard deviations as
116 15 18.67 0.05 2.98 1.04 compared to the small confidence interval half
116 16 13.88 0.04 5.22 1.35 widths (see Table 4). For Node 126 we see a large
116 17 23.63 0.06 2.65 0.71 reduction in the number of replications as
compared to table 3 primarily because of the
116 18 20.15 0.05 4.91 1.32
increase in the user specified confidence interval
116 19 19.31 0.04 6.31 1.00 half width of 10.
126 15 17.03 0.03 15.01 0.46
Now we look at the estimated replications for the
126 16 18.25 0.03 29.07 0.60 entire network. Tables 6 and 7 show the summary
126 17 26.43 0.04 33.26 1.03 of results for all links for link use and link
126 18 15.17 0.02 68.73 0.79 encounters for an alpha of 0.05, and a confidence
interval half width of 10. Note in table 6 it is
126 19 16.40 0.02 26.58 1.06 possible to achieve the user desired accuracy for
90% of the links with less than 60 replications for
Table 4. Confidence Half Intervals and
link use whereas in table 7 it takes over 180
Relative Precision for Link Use and Link
replications to achieve the same coverage for link
Encounters, Alpha = 0.05
encounters.

Table 4 shows, as we would expect, that the

confidence interval half widths and relative
precisions have increased because of the lower
alpha value.

146
Link Use Absolute Accuracy = 10
Reps Count %Links Accum% Link Use - Relative Precision = 0.15
1: 20 2926 77.02% 77.02% Reps Count % Accum%
21: 40 428 11.27% 88.29% 1: 40 3493 91.95% 91.95%
41: 60 111 2.92% 91.21% 41: 80 87 2.29% 94.24%
61: 80 12 0.32% 91.52% 81:120 44 1.16% 95.39%
81:100 16 0.42% 91.95% 121:160 23 0.61% 96.00%
101:120 6 0.16% 92.10% 201:240 14 0.37% 96.37%
121:140 17 0.45% 92.55% 281:320 10 0.26% 96.63%
141:160 16 0.42% 92.97% 401:440 30 0.79% 97.42%
181:200 267 7.03% 100.00% 481:500 23 0.61% 98.03%
Total 3799 100.00% 501: 75 1.97% 100.00%
Number of replications for link use using equation Total 3799 100.00%
3 for alpha = 0.05, CI Half Width = 10 for entire
network Table 7. Number of replications for link use
using equation 4 for alpha = 0.05 relative
Link Encounters Absolute Accuracy = 10
accuracy = 0.15 for entire network
Reps Count %Links Accum%
1: 20 1120 53.77% 53.77% Link Encounters - Relative Precision = 0.15
21: 40 169 8.11% 61.88% Reps Count % Accum%
41: 60 95 4.56% 66.44% 2.74% 2.74%
1: 40 57
61: 80 66 3.17% 69.61% 4.42% 7.15%
41: 80 92
81:100 61 2.93% 72.54% 6.10% 13.25%
81:120 127
101:120 45 2.16% 74.70% 6.67% 19.92%
121:160 139
121:140 36 1.73% 76.43% 5.62% 25.54%
161:200 117
141:160 33 1.58% 78.01% 5.95% 31.49%
201:240 124
161:180 26 1.25% 79.26% 5.62% 37.11%
241:280 117
181:200 432 20.74% 100.00% 4.66% 41.77%
281:320 97
Total 2083 100.00% 2.59% 44.36%
321:360 54
Table 6. Number of replications for link 361:400 28 1.34% 45.70%
encounters using equation 3 for alpha = 0.05, CI 401:440 147 7.06% 52.76%
Half Width = 10 for entire network
441:480 131 6.29% 59.05%
481:500 48 2.30% 61.35%
Tables 7 and 8 show the results of equation 4 with
a relative accuracy of 0.15. A similar pattern is 501: 805 38.65% 100.00%
seen here with less than 40 replications needed for Total 2083 100.00%
link use to achieve the relative precision for 0.15
for link use and over 500 replications required for Table 8. Number of replications for link
link encounters. Since both link use and link encounters using equation 4 for alpha = 0.05
encounters are generated in the same simulation relative accuracy = 0.15 for entire network
the obvious method for selecting the required
simulations is simply to take the maximum value
from all estimates, which is this case is over 500
replications. However this ignores the costs of 6. CONCLUSIONS
processing time. In this simulation, each These results show a number of important
replication takes around an hour to run. If we run characteristics of performance indicators and the
the simulation for 500 replications this means 500 impact on the number of replications required to
hours of computer time and with limitations on file meet user-specified reliability measures. First the
sizes for output databases, it is likely the more sensitive a performance measure is to
simulation will fail from reaching file size limits. random variation, the higher the variances and the

147
wider the confidence interval and therefore the
greater number of replications required. Second, if
using equations 2 and 3 for estimating the number
of replications, it is important to carefully select
the absolute accuracy (desired CI half width) for
each performance indicator – it may be helpful to
use the output analysis for the “short run” to help
determine reasonable values for each performance
indicator. Third, absolute accuracy and relative
precision are two very different measures yielding
very different results in terms of calculating the
number of replications. It may be useful to
evaluate both measures when estimating the
number of replications for terminating simulations
using the methods described in this paper.
When estimating the number of replications to
obtain reliable results for management there is a
tradeoff between the different performance
measures described here and the processing time
and limitations of output file sizes. The reliability
of results will therefore vary with each
performance indicator and with each location. If a
trade-offs need to be made between accuracy and
computing resources it is advisable to map out the
reliability of the results for each performance
indicator. The analyst and manager may simply
have to live with a degree of uncertainty or accept
that confidence interval widths will vary across
space. At least, the methods presented in this
paper give the analyst the tools needed to
accurately measure this variation spatially.

7. REFERENCES
Centeno, Martha. and Reyes, M. Florencia (1998),
"So you have your model: What to do next.
A tutorial on simulation output analysis" in
Proceedings of the 1998 winter Simulation
conference, D.J. Medeiros, E.F. Watson,
J.S. Carson and M.S. Manivannan, editors.
http://www.informs-cs.org
Goldsman, D. (1992), Simulation Output Analysis.
In Proceedings of the 1992 Winter
Simulation Conference, ed. J. J. Swain, D.
Goldsman, R. C. Crain, and J. R. Wilson,
97-103.
Law, Averill M. and Kelton, W. David (2000)
Simulation Modeling and Analysis, Third
Edition, McGraw Hill, Boston.

148