Alyson Wilson, Ph.D.
Technical Staff Member, Statistical Sciences Group
Los Alamos National Laboratory, Los Alamos, New Mexico
Stephen Pollock, Ph.D.
Professor Emeritus, Department of Industrial and Operations Engineering
University of Michigan, Ann Arbor, Michigan
This appendix illustrates two suggestions from Chapter 3 with illustrative R code. In particular, we consider:
the addition of an 18th stage to represent distributions of alternate consequences; and
replacing distributions of arc probabilities by expected values of the probabilities.
We work from the event tree in Figure C.1. For simplicity, we assume a single initiating event. For concreteness, we assign uncertainty distributions to each of the arc probabilities:
P_{A1} ~ Beta(2,2);
P_{T1} ~ Beta(4,1); and
P_{T2} ~ Beta(3,2).
In addition, we know the distributional form of each consequence distribution. Using the notation c(x | s_{1}) to denote the consequence distribution associated with the first arc, we assign the following distributions to consequences:
c(x | s_{1}) ~ Gamma(8000,2);
c(x | s_{2)} ~ Gamma(4500,1);
c(x | s_{3}) ~ Gamma(10000,2); and
c(x | s_{4}) ~ Gamma(5500,1).
We would like to know the form of the risk distribution. Summary statistics from this distribution (5th percentile, mean, 95th percentile) are used to summarize risk and present analyses in the Biological Threat Risk Assessment (BTRA) of 2006.
A simple way to simulate from the risk distribution is as follows:
Repeat n times;
Sample from each arc probability;
Calculate the probabilities for each scenario;
Choose a scenario using the calculated probabilities;
Sample from the consequence distribution for that scenario;
The n samples constitute a sample from the risk distribution; and
Summarize these samples using a histogram, empirical quantiles, and sample mean.
Below are the first 10 and last 10 pages of uncorrected machine-read text (when available) of this chapter, followed by the top 30 algorithmically extracted key phrases from the chapter as a whole.
Intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text on the opening pages of each chapter.
Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.
Do not use for reproduction, copying, pasting, or reading; exclusively for search engines.
OCR for page 80
Appendix C
Computational Example Illustrating the Replacement of
a Joint Distribution of Arc Probabilities with Marginal
Expected Values of Individual Arc Probabilities
Alyson Wilson, Ph.D.
Technical Staff Member, Statistical Sciences Group
Los Alamos National Laboratory, Los Alamos, New Mexico
Stephen Pollock, Ph.D.
Professor Emeritus, Department of Industrial and Operations Engineering
University of Michigan, Ann Arbor, Michigan
This appendix illustrates two suggestions from Chapter 3 with illustrative R code. In particular, we consider:
• the addition of an 18th stage to represent distributions of alternate consequences; and
• replacing distributions of arc probabilities by expected values of the probabilities.
We work from the event tree in Figure C.1. For simplicity, we assume a single initiating event. For concreteness,
we assign uncertainty distributions to each of the arc probabilities:
PA1 ~ Beta(2,2);
PT1 ~ Beta(4,1); and
PT2 ~ Beta(3,2).
In addition, we know the distributional form of each consequence distribution. Using the notation c(x | s1) to denote
the consequence distribution associated with the first arc, we assign the following distributions to consequences:
c(x | s1) ~ Gamma(8000,2);
c(x | s2) ~ Gamma(4500,1);
c(x | s3) ~ Gamma(10000,2); and
c(x | s4) ~ Gamma(5500,1).
We would like to know the form of the risk distribution. Summary statistics from this distribution (5th percentile,
mean, 95th percentile) are used to summarize risk and present analyses in the Biological Threat Risk Assessment
(BTRA) of 2006.
A simple way to simulate from the risk distribution is as follows:
• Repeat n times;
• Sample from each arc probability;
• Calculate the probabilities for each scenario;
• Choose a scenario using the calculated probabilities;
• Sample from the consequence distribution for that scenario;
• The n samples constitute a sample from the risk distribution; and
• Summarize these samples using a histogram, empirical quantiles, and sample mean.
0
OCR for page 80
APPENDIX C
PT1
c(x|s1)
PA1
1-PT1
c(x|s2)
PT2
c(x|s3)
1-PA1
1-PT2
c(x|s4)
FIGURE C.1 A simple event tree for two successive stages (events), each with two outcomes. For this example, each path through
the tree represents a unique scenario with its own consequence distribution.
R code implementing this algorithm follows.
n <- 1000000
consq <- rep(0,n)
for (i in 1:n) {
pa1 <- rbeta(1,2,2)
pt1 <- rbeta(1,4,1)
pt2 <- rbeta(1,3,2)
s1p <- pa1*pt1
s2p <- pa1*(1-pt1)
s3p <- (1-pa1)*pt2
s4p <- (1-pa1)*(1-pt2)
scen <- rmultinom(1,1,c(s1p,s2p,s3p,s4p))
if (scen[1] == 1) consq[i] <- rgamma(1,8000,2)
if (scen[2] == 1) consq[i] <- rgamma(1,4500,1)
if (scen[3] == 1) consq[i] <- rgamma(1,10000,2)
if (scen[4] == 1) consq[i] <- rgamma(1,5500,1)
}
hist(consq,freq=F,main=””,xlim=c(3500,6000),
xlab=”Consequence Distribution”,ylim=c(0,0.0035))
lines(density(consq))
quantile(consq,c(0.05,0.95))
mean(consq)
OCR for page 80
2 DEPARTMENT OF HOMELAND SECURITY BIOTERRORISM RISK ASSESSMENT
The histogram summarizing the risk distribution from this approach is given in Figure C.2, with an overlay of a
kernel density estimator of the risk distribution as the solid line.
The histogram and solid black line result from brute force sampling from the arc probability distributions and the
consequence distributions. The line with circles is the estimate from the methodology employed in the BTRA of 2006,
which can also produce risk curves. The line with triangles is the estimate from a greatly simplified algorithm that uses
only the marginal expected values of individual arc probabilities and simulations from the consequence distributions.
The line with crosses is calculated assuming a parametric (or tabular) form is known for the consequence distributions
and requires no simulation. Notice the good agreement between the four estimates.
For an event tree as complex as the one presented in the BTRA, this approach is infeasible. As we understand it,
the approach implemented in the BTRA is as follows:
• Draw 500 samples from each arc probability;
• Calculate 500 sets of scenario probabilities;
• Draw 1000 samples from each consequence distribution;
• Represent each consequence distribution as a histogram;
• For each of the 500 sets of scenario probabilities, calculate a weighted average of the mass in each bin of the
histogram, and call this one “sampled risk curve”;
• Calculate the average over all 500 risk curves. Use this as an approximation to the risk distribution and calculate
the mean, 5th percentile, and 95th percentile; and
• Also calculate the 5th and 95th percentiles for the entire set of risk curves.
FIGURE C.2 This plot illustrates estimates of the risk distribution for the simple event tree using three different algorithms.
R01268, Figure C-2
Fixed image, not changeable
OCR for page 80
APPENDIX C
R code implementing this algorithm follows.
nsampbr <- 500
pa1 <- rbeta(nsampbr,2,2)
pt1 <- rbeta(nsampbr,4,1)
pt2 <- rbeta(nsampbr,3,2)
s1p <- pa1*pt1
s2p <- pa1*(1-pt1)
s3p <- (1-pa1)*pt2
s4p <- (1-pa1)*(1-pt2)
nsampc <- 1000
cs1 <- rgamma(nsampc,8000,2)
cs2 <- rgamma(nsampc,4500,1)
cs3 <- rgamma(nsampc,10000,2)
cs4 <- rgamma(nsampc,5500,1)
bh1 <- hist(cs1,breaks=seq(3500,6000,length=101),plot=F)$density
bh2 <- hist(cs2,breaks=seq(3500,6000,length=101),plot=F)$density
bh3 <- hist(cs3,breaks=seq(3500,6000,length=101),plot=F)$density
bh4 <- hist(cs4,breaks=seq(3500,6000,length=101),plot=F)$density
qdm <- matrix(0,nsampbr,100)
for (i in 1:nsampbr) {
qdm[i,] <- s1p[i]*bh1 + s2p[i]*bh2 + s3p[i]*bh3 + s4p[i]*bh4
}
qdmean <- apply(qdm,2,mean)
qd5 <- apply(qdm,2,quantile,c(0.05))
qd95 <- apply(qdm,2,quantile,c(0.95))
x <- seq(3512.5,5987.5,by=25)
points(x,qdmean,type=”b”,pch=1)
The estimated risk distribution from this approach is given as the line with circles in Figure C.2.
As shown in Chapter 3, the risk distribution can be calculated without sampling from the arc probability distribu-
tions. For an event tree the size of the one used in the BTRA of 2006, this represents a significant computational
simplification. What is lost in the simplification is the family of risk curves—i.e., one curve for each possible outcome.
However, no analysis in the BTRA of 2006 and no improvement in analysis recommended by the committee can make
meaningful use of the information available in the family of risk curves, beyond that provided by their expectation.
Further, given the improvements proposed for the BTRA to incorporate additional consequence measures and utility
functions, the committee does not see upcoming analyses that require the family of risk curves.
Consider the following simplified algorithm:
• Draw 1000 samples from each consequence distribution;
• Represent each consequence distribution as a histogram; and
• Calculate a weighted average of the mass in each bin of the histogram using the expected arc probabilities and
use this as the estimated risk distribution.
OCR for page 80
DEPARTMENT OF HOMELAND SECURITY BIOTERRORISM RISK ASSESSMENT
R code implementing this algorithm follows.
ms1p <- (0.5)*(0.8)
ms2p <- (0.5)*(0.2)
ms3p <- (0.5)*(0.6)
ms4p <- (0.5)*(0.4)
nsampc <- 1000
cs1 <- rgamma(nsampc,8000,2)
cs2 <- rgamma(nsampc,4500,1)
cs3 <- rgamma(nsampc,10000,2)
cs4 <- rgamma(nsampc,5500,1)
bh1 <- hist(cs1,breaks=seq(3500,6000,length=101),plot=F)$density
bh2 <- hist(cs2,breaks=seq(3500,6000,length=101),plot=F)$density
bh3 <- hist(cs3,breaks=seq(3500,6000,length=101),plot=F)$density
bh4 <- hist(cs4,breaks=seq(3500,6000,length=101),plot=F)$density
erd <- ms1p*bh1 + ms2p*bh2 + ms3p*bh3 + ms4p*bh4
x <- seq(3512.5,5987.5,by=25)
points(x,erd,type=”b”,pch=2)
The estimated risk distribution from this approach is given as the line with triangles in Figure C.2.
If the conditional consequence distributions are given in parametric form, or in numerical look-up tables, calcula-
tion of the risk distribution can be done exactly, without resorting to estimating these distributions from the outputs
of Monte Carlo simulations. This method is simply:
• Calculate the expected arc probabilities; and
• Calculate the weighted average of the consequence distributions.
ms1p <- (0.5)*(0.8)
ms2p <- (0.5)*(0.2)
ms3p <- (0.5)*(0.6)
ms4p <- (0.5)*(0.4)
x <- seq(3512.5,5987.5,by=25)
points(x,ms1p*dgamma(x,8000,2) + ms2p*dgamma(x,4500,1) + ms3p*
dgamma(x,10000,2) + ms4p*dgamma(x,5500,1),type=”b”,pch=4)
The risk distribution (exact, and not an estimate) obtained using this approach is given as the line with crosses in
Figure C.2. This computation is both trivial and fast.