Cover Image


View/Hide Left Panel

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)

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, calculation 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.

The National Academies | 500 Fifth St. N.W. | Washington, D.C. 20001
Copyright © National Academy of Sciences. All rights reserved.
Terms of Use and Privacy Statement