1994 INTERNATIONAL SYSTEM DYNAMICS CONFERENCE
Analytic Uncertainty Modeling
Henry Neimeier
The MITRE Corporation
7525 Colshire Drive
McLean, Virginia, 22102, USA
Abstract
The analytic uncertainty is useful analysis is imp It
provides the entire resulting probability distribution instead of a single uncertain point estimate of
the mean. Both analytic costs, and ion costs are far less than in discrete
event simulation. The price paid is some lack in modeling flexibility.
Discrete simulation requires multiple long simulation runs to obtain a statistically significant point
estimate. The different result values from multiple runs with identical parameter values but different
random number seeds, are averaged to obtain the point estimate of the mean result value.
Conversely, the analytic solution gives the entire resulting probability distribution with minimal
calculation. The analytic solution also considerably simplifies sensitivity analysis. A single analytic
tun is done for each input parameter setting. Discrete event simulation requires multiple runs for each
input to obtain a statistically i mean result.
In functional economic analysis we are interested in the relative future costs of alternative systems.
There are uncertainties in process performance, resource cost esti
required, workload, interest and inflation rates. There is also uncertainty in the future projection of
these el Analytic 'y modeling provides a simple way of calculating output measure
uncertainty from model input parameter uncertainties.
System Dynamics : Methodological and Technical Issues, page 195
1994 INTERNATIONAL SYSTEM DYNAMICS CONFERENCE
Analytic Uncertainty Modeling
Overview
The beta distribution analytic technique outlined in this paper is useful whenever sensitivity analysis
is important. It provides the entire resulting probability distribution vice a single uncertain point
estimate of the mean. Both analytic costs, and ‘ion costs are far less than
in discrete event simulation. The price paid is some lack in modeling flexibility.
With the appropriate choice of parameter values, the beta distribution closely fits all the
classical probability distributions. The sums and products of beta variates are also approximately beta
distributed. This paper documents the error in beta approximation. The beta distribution can be fit
based on the mini mean, i and standard deviation statistics. In a complex results
calculation all that is required is to keep track of these statistics as the calculation proceeds. At any
point in a calculation, the probability distribution of the result, can be derived by fitting a beta
distribution based on the four statistics.
Discrete Event Simulation Times
Discrete simulation requires multiple long simulation runs to obtain a statistically significant point
estimate. The different result values from multiple runs with identical parameter values but different
random number seeds, are averaged to obtain the point estimate of the mean result value.
Conversely, the analytic solution gives the entire resulting probability distribution with minimal
calculation. The analytic solution also considerably simplifies sensitivity analysis. A single analytic
run is done for each parameter. setting vice multiple runs for a statistically significant result. The
simulation time (T) required to be 95 percent confident in a relative error (e) is approximated by the
following equation for open GI/G/1 (general independent inter arrival times, general service times, 1
server) queuing networks:
T= B(C,24C2) 274 (1-12) e2)
Where:
T =simulation time for a specified relative error
t =service time
oP = square coefficient of variation in inter arrival time
(variance in inter arrival time divided by the mean
inter arrival time squared)
g? = square coefficient of variation in service times
Z = =unit normal deviate (Z=2 for 95 percent confidence
r =utilization (service time divided by inter arrival time)
e = tolerated relative error
Figure | is a semi-log plot of simulation time required for 95 percent confidence in a specified relative
error as a function of utilization. It represents the exponential inter arrival and service time case
(Ca=Cs=1). Note that at high utilization and low relative errors extremely long simulation times are
required. To achieve 5 percent relative error in the mean on an 80 percent utilized queuing network
requires one million service times.
In functional economic analysis we are interested in the relative future costs of alternative
systems. There are uncertainties in process performance, resource requirements, cost estimates,
investment required, workload, interest and inflation rates. There is also uncertainty in the future
projection of these elements. Thus there is uncertainty in the discounted present value cost
distribution for each alternative system. A plot of cumulative probability versus cost, aids the
decision process. The entire range in cost distribution is of interest. Figure 2 shows the expected
System Dynamics : Methodological and Technical Issues, page 196
1994 INTERNATIONAL SYSTEM DYNAMICS CONFERENCE
number of simulation events required to obtain an event in the tail of the result distribution when
using discrete event simulation. The equation plotted is:
E=1/p© Where:
E= expected number of simulation events
P = distribution tail probability
C = uncertain model components
The lower the tail probability, and the more components in the model, the more events are required.
For example, an average of one million simulation events-are Tequired i in a six ‘Component model to
simultaneously be in the 10 percent tail of all distrib To sii ly be in the 1
percent tail requires an average of one trillion simulation events. In the limit it requires an infinite
number of simulation events to capture the entire range of results. Thus discrete event simulation is
not practical if one is interested in the entire result distribution in other than very small models with
few P If the mini and istribution values are not needed then discrete
event simulation is practical. However, even in this case the model development, execution, and
sensitivity analysis costs are higher.
Relative Error
no
- -
: a
3 --"
0.1
Hier Millions
|
---=-0.01
seeee ees 0.005
2
Es
é
8
z
3
0.75 +
Utilization
Figure 1. Simulation Time For Specified Relative Error
System Dynamics : Methodological and Technical Issues, page 197
1994 INTERNATIONAL SYSTEM DYNAMICS CONFERENCE
to frail y (P)
10 3 —— 0.010
7
1o 2 0.025
—- —-— 0.100
P > —-- 0.150
0.200
Mean Simulation Events
Figure 2. Si ion Events For A ified Tail Probability Versus Model Components ares
Beta distribution
The beta density distribution is bounded by high (max) and low (min) values.
P=C {(x-min)/(max-min)} ael {1-(x-min)/(max-min)}
ere:
P is the probability density
C normalizes the area under the distribution to unity
(C= VU," x Pl ayy
min is the minimum variable value
max is the maximum variable value
If a random variable is bounded then its distribution function is uniquely determined by its moments
(reference 3, p 126). The beta distribution can be fit based on the minimum, maximum, mean and
variance. The fit parameters a and b are based on a range variable r and a skewness variable s, and are
defined in standard terms as:
r= variance / (max - min)”
s = (mean - min) / (max- min)
a=s"(1-s)/r-s
b=s(l-s)/r-l-a
If (a-1)(b-1)>0 then the beta distribution has a mode at
min + (max - min) (a- 1)/(a+b-2)
Figure 3 plots the beta distribution a and b parameter values as a function of the r and s
parameters. The ratio of standard deviation to range (square root of r) selects the appropriate curve.
System Dynamics : Methodological and Technical Issues, page 198
1994 INTERNATIONAL SYSTEM DYNAMICS CONFERENCE
Separate curves are presented for a and b parameters at a selected set of "r" values. The symbol
legend is at the right. The symbol for the a curve is a darkened version of the b
symbol. The "s” sk (am init ) maxis ini )) is plotted along the
abscissa. Select the appropriate value and read the a and b values off the ordinate. Note that the
smaller the standard deviation the larger the a and b parameters. For the mean at the midpoint
between minimum and maximum values (s=.5) we have a symmetric distribution. The a and b
parameters are equal in this symmetric case.
V r=Std.Dev/Range
é —h—a.t6
hats
5 a2
PA —'—2.22
4 i ——a.24
H Y, ——a.26
3 fi —+a.30
/. A —_—~a.38
2 LA ——b.16
7; 7 ——p.18
1 4 “7 —i—b.2
iG CAg =e.
Agee: NS ——b22
Oey — 4 —t—b.24
0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95|—0—b.26
s=(Mean-Minimum) / beso
(Maximum-Minimum) ——b.38
Figure 3. Beta Distribution a and b Parameters
Figure 4 shows examples of different beta distribution shapes. For the mean midway between
the minimum and maximum values, the distribution is symmetric with equal a and b parameter values.
Unity a and b parameters yield a uniform distribution (1,1). Values less than 1 lead to a "U" shaped
distribution (.5,.5). With a and b equal 2 (2,2) the distribution has a parabolic shape. At higher a and
b values (6,6; 12,12) the distribution has a shape similar to a normal distribution. The lower the
standard deviation relative to the range, the higher the a and b parameter values, and the more peaked
the distribution shape. If a is greater than b then the distribution has a negative skew (9,2).
Conversely, if a is less than b the distribution has a positive skew. Triangular distributions (1,3) and
left (.5,3) and right "J" shaped distribution shapes are also possible.
System Dynamics : Methodological and Technical Issues, page 199
1994 INTERNATIONAL SYSTEM DYNAMICS CONFERENCE
a,b Parameters
art 8S
> =r
ss
Ey sane
&
2 re
a ‘cnr 7 |
s
8
= 722
oe
aia
Figure 4. Beta Distribution Shapes
Spanning the classical distribution space
Figure 5 is a skew (+b1) kurtosis (b2) plot of the popular classical continuous statistical distributions.
The beta distribution covers most of the classical distribution area except the log normal line. In the
case of the log normal line a log transformation is performed before the beta distribution is fit. The
beta distribution shapes are shown on the chart regions. Beta distributions can fit uniform, triangular,
J shaped, U shaped, exp ial, Erlang, hyp p ial, gamma, and normal distributions.
Exponential
Point
yy
© | ve cont
ong vie oan 6 os
@ Fi
8 eg
iS I~ gan
x 1
ah |
4
car:
pon
aor”
Skew
Figure 5. Skew Kurtosis Plot of Classical Distributions
Beta approximation errors
The beta distribution closely fits all the classical distributions. The) sums and products of beta variates
are also closely app’ d by a beta distribution. The g beta distribution can easily be
calculated. One does not have to resort to putati i ‘ion for the sums of
System Dynamics : Methodological and Technical Issues, page 200
1994 INTERNATIONAL SYSTEM DYNAMICS CONFERENCE
random variates. All that is required is to calculate the mini i mean and standard
deviation statistics. The following paragraphs calculate the error in using the beta distribution to fit
sums or products of uniform random variables. Then the error in approximating the normal
distribution is calculated. Finally the mean and variance formulas for different types of calculations
are presented.
Sums and products of uniform random variables
The uniform distribution is fit exactly by the beta distribution (a=b=1). The sum of two uniform
variates is a triangular distribution with a discontinuity at the mode. In the case of discontinuities, the
error in beta approximation is minimal. The cumulative probability distribution function of the sum
of n independent uniform random variables between 0 and 1 is :
Pr{Sn£x]=S (-1) (n!/G! (nj)! )) &f)" /n!
Where:
x isa random variate between 0 and n
nis the number of uniform (0-1) random variates that are summed
Pr[Sn] is the cumulative probability of the sum of n independent
uniform random variates
The summation (S) is over all integers j <x
The minimum, maximum, mean, and variance of the distribution of the sum distribution were used to
fit a beta distribution. The difference between the fit beta distribution and the actual sum distribution
was calculated over the entire range of distribution values. Figure 6 gives the results for sums and
products of 2,3, or 4 uniform distributions.
The maximal error for sums ranges from 1.2 percent for 2 variates down to less than 1/2
percent for the sum of four uniform variates. The error rates in the tails of the distribution are far less.
This is the area we are most interested in when making confidence statements.
Since the product of variates tend to be log normally distributed (outside the skew kurtosis
beta region) the beta approximation for the product is not as good as the sum. Log transformation of
the data before fitting would give a closer fit, but this would increase the calculation complexity.
Using the uniform distribution as an example, the cumulative distribution function of the
product of two uniform distributions is:
Prix £ p] =p - p In(p)
System Dynamics : Methodological and Technical Issues, page 201
1994 INTERNATIONAL SYSTEM DYNAMICS CONFERENCE
Beta Error Sum Of 2 Uniform Variables Beta Error Product Of 2 Uifora Variables
0.04 0.01
0.005+ 0.0076;
[\ 0,005;
5 as 2
0.0025;
0.005
“0.04 02 o4 0.6 0.8 1
0.0025:
Beta Error Sum Of 3 Uniform Variables Beta Efrox Product Of 3 Unitora Yarisbles
0. 006
0.045:
0. 004
0.002 [\ oe
05/ 1 16 2/25 3 om
0.002
~0.004 0.2 4 06 08 1
0. 006 -0.005;
Beta Error Sum Of 4 Uniform Variables Bete Efror Product Of 3 Unitora Variables
0.004 /\ 0.015;
0.002
aN 0.08
3, 4 0.006:
0,002
oz fa os 08 1
0.004 -0. 005;
Figure 6. Beta Error Functions
The cumulative distribution function of the product of three random variates is:
Pr{x £ p]= p-p In(p) +p In(p)* /2
The cumulative distribution function of the product of four random variates is:
Pr[x £ p] = p-pin(p) +p In(p)* /2- p in(p)? / 6
The error increases as the number of terms in the product increases. However the error in the upper
tail is very reasonable (<1/2 percent). Kotlarski and t than have ii igated general di
under which products of independent variables have a beta distribution. Springer uses Mellin
transforms to calculate the products and quotients of beta random variables.
System Dynamics ; Methodological and Technical Issues, page 202
1994 INTERNATIONAL SYSTEM DYNAMICS CONFERENCE
Beta approximation to normal distribution
The sum of independent variates hes a normal distribution as more variates are summed
(central limit theorem) The symmetric normal distribution tails extend from minus infinity to plus
infinity. However the area in the tails beyond 3 or 4 standard deviations from the mean is minimal
(0.27 percent for 3 standard deviations, 0.006 percent for 4 standard deviations). Figure 7 shows the
error in approximating the normal distribution with a Beta distribution.
Beta Error Beta Error
088 0.0075:
0.04 0.005
ue /\ 0.0025
a4 o4f 06 Oe oa Vo fos lot
0.005 -0,.0025
0.04 0.008
-0.0075
0.045
Range = 6 Standard Deviations Range = 8 Standard Deviations
(a=b=4) (a=b=7.5)
Figure 7. Beta Errors In Approximating Normal Distributions
The maximal error is 1.5 percent when the beta range is set to 6 standard deviations (+-3) or 0.75
percent when the beta range is set to 8 standard deviations (+-4). Note that again the error in the high
interest tails of the distribution is far less.
Mean and variance statistics
In practice any calculation with beta variates usually yields beta variates. Thus if one encodes
parameter uncertainty with a beta distribution, the results of a calculation with many uncertain
parameters will also be beta distributed. This allows analytic solution for the resulting distribution
without resorting to discrete simulation. This greatly reduces the calculation requirement and
simplifies parametric sensitivity analysis.
In a process cost or performance calculation, operations must be performed on uncertain
parameter values to determine expected results and uncertainty in results. Operations include:
addition, sub i iplication, power, pol ial function, general function. Using the Beta
distribution we must keep track of the calculation minimum, maximum, mean (u) and variance (s2)
values. The variance of the sum or diffe of two indep is just the sum of their
component variances. If the are lated the following equation is used:
Fy4g7s? +82 425, oa =52 + 22 S19
Where s}2 is the covariance of parameters | and 2. The variance of the product of two uncertain
independent parameters | and 2 is given by the following formula:
s?yg=s2] 829 +uy 2829 +uy 2s?)
If parameter a is a constant with value a (minimum = mean = maximum = a, and variance = 0) then
the variance of the product of parameters 1 and 2 is:
24902 52,
‘System Dynamics : Methodological and Technical Issues, page 203
1994 INTERNATIONAL SYSTEM DYNAMICS CONFERENCE
In the case of division, the maximum result is the maximum quotient divided by the minimum divisor.
Conversely the minimum result is the minimum quotient divided by the maximum divisor. The
variance of the result and the corrected mean are obtained from:
s wigh is an 2 2 au 2)
2 102 101
uyy = (uy 0g) (H8/19)2)
Ify is a linear function of n input variables (y=agta, x) tajxot..+a,X,) then the variance of y is
given by:
n
s 3 (atioxy2 si+28 S$ (flax; otax; )sjj2
i=l i=] jritl
If the input variables are independent the second term is ignored. If y is not a linear function of the
input variables the above equation is only an approximation. In the case of a general single variable
function (y=f(x)) with known derivatives, the variance of the result is approximated by the following:
(@ylax)* s,2+ 1/20? y/ax2)2 5,443 yax3) @y/Ax) 5x44.
Usually only the first term provides a bi imati Additional terms can be found in the
Tukey report. Approximations for functions of several variables (z=f(x,y,...)) are derived based on
the multidimensional Taylor series. | Note that the mean y (y=f(x)) is not equal to the result of
substituting the mean x value into the function, if the function is non-linear. The shift in mean value
is given by:
172 @yldx2) 5,2+1/8 (4 ylax4) 54+
For a general function, a first-cut approximation is to the output for
mean, and mean plus standard deviation inputs. The squared difference of these outputs is an
approximation of the output variance.
Conclusions
In functional economic analysis we are interested in the relative future costs of alternative
systems. There are uncertainties in process performance, cost estit and
workload. Thus there is uncertainty in the discounted Present value cost distribution for each
alternative system. This paper p dan analytic t to calculate the result distribution for
alternatives from estimates of p ies. The technique has wide applicati It
considerably simplifies sensitivity analysis. It should be considered when the probability distribution
of a result is desired rather than a single point estimate of the mean. Both analytic development
costs, and computer execution costs are far less than in discrete event simulation. The price paid is
some lack in modeling flexibility.
References
Whitt, W. 1989.Planning Queueing Simulations, Management Science,
Vol 35, No.11, November 1989.
Johnson, N., Kotz,S. 1970. Continuous Univariate Distributions-2,
John Wiley & Sons, New York.
Wilks, S. 1962 Mathematical Statistics, New York, John Wiley & Sons,
System Dynamics : Methodological and Technical Issues, page 204
1994 INTERNATIONAL SYSTEM DYNAMICS CONFERENCE
New York. ;
Tukey,J. The Prop ion of Errors, Fli ions, and Tole Basic Ge lized Formulas,
AD155082, Princeton University, Princeton,
New Jersey.
Koltarski,I.1962 On Groups of N Independent Random Variables Whose Product Follows The Beta
Distribution, Warsaw, COLLOQUIUM MATHEMATICUM, vol IX, FASC.2,pp325-332, Warsaw,
Poland.
Jambunathan, M.1954.Some Properties Of Beta And Gamma Distributions, The Annals of
Mathematical Statistics, Vol.25,pp 401-405,1954.
Springer,M.1979. The Algebra of Random Variables, John Wiley & Sons, New York.
Seiler,F. 1987. Error Propagation for Large Errors, Risk Analysis, Vol.7, No.4,pp 509-518.
System Dynamics : Methodological and Technical Issues, page 205
