Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.
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. 80
APPENDIX C 81 PT1 c(x|s1) PA1 1-PT1 c(x|s2) PT2 1-PA1 c(x|s3) 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)
82 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
APPENDIX C 83 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.
84 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.