How to handle extremely rare events in a meaningful way
Scott Ferson
Institute for Risk and Uncertainty, University of Liverpool
Nuclear Energy - GREEN CDT, 26 July 2023
Photo credit: Tatyana Ginzburg, used under license
This presentation has animations. Select Slide Show / From Beginning to see them.
Turn on recording and use microphone
Today
�NAFEMS is the International Association for the Engineering Modelling, Analysis and Simulation Community
https://en.wikipedia.org/wiki/Simulation_governance
Simulation governance is a managerial function concerned with assurance of reliability of information generated by numerical simulation. The term was introduced in 2011 [1] and specific technical requirements were addressed from the perspective of mechanical design in 2012.[2] Its strategic importance was addressed in 2015.[3][4] At the 2017 NAFEMS World Congress in Stockholm simulation governance was identified as the first of eight “big issues” in numerical simulation.
Edimu et al. (2010)
Risk analysis in engineering 26
So how are we doing? 7
Dependence 9
Kinds of uncertainty 8
Getting the inputs 18
Risk communication 3
Bayesian networks 4
Backcalculation 12
Conclusions 12
NAFEMS challenge
Fire alarm at 11 am
Look for the falconer on the roof next door
Open question
Suppose
A is in [2, 4]
B is in [3, 5]
What can be said about the sum A+B?
Deepwater Horizon Spill (2010)
Columbia Disaster (2003)
2007-2008 Financial Crisis
Fukushima Daiichi Nuclear Disaster (2011)
Air France Flight 447 Rio−Paris (2009)
Minneapolis Bridge Collapse (2007)
Credits (left to right):
Unknown photographer - US Coast Guard - 100421-G-XXXXL- Deepwater Horizon fire (Direct link), public domain.
Pawel Kierzkowski, https://en.wikipedia.org/wiki/Air_France_Flight_447, CC BY-SA 3.0.
Dominic Alves from Brighton, England - Northern Rock Queue, CC BY 2.0.
The Ponte Morundi collapsed 45 meters killing 43 people, Image credit: The Daily Sabah
Kevin Rofidal, United States Coast Guard, https://en.wikipedia.org/wiki/I-35W_Mississippi_River_bridge, public domain.
https://www.youtube.com/watch?v=aXiZ3RHR3bg
RAAC (reinforced autoclaved aerated concrete), collapse of a school roof in Kent (2018) image credit: no credit
Morry Gash, NPR, https://www.npr.org/2017/08/01/540669701/10-years-after-bridge-collapse-america-is-still-crumbling?t=1614450713959
The British bank Northern Rock failed as a result of a bank run precipitated by the global financial crisis. It was the first bank in 150 years to fail due to a bank run.
Don’t confuse the Fukushima Daiichi nuclear plant with a more dramatic picture of natural gas storage tanks burn at a facility in Chiba Prefecture, near Tokyo, which was another consequence of the same earthquake and tsunami that affected the Fukushima Daiichi nuclear power plant.
Hurricane Katrina
2011 natural gas storage tanks burn at a facility in Chiba Prefecture, near Tokyo, after !! It is not a nuclear power plant! But is another consequence of the same earthquake and tsunami that effected the Fukushima Daiichi nuclear power plant
Tōhoku earthquake and tsunami
2010 Failure of the blowout preventer led to eleven deaths, the largest-ever marine oil-spill, environmental catastrophe and costs to BP of more than $60B.
2009 crash with 228 fatalities of Air France flight 447 en route from Rio de Janeiro to Paris due to sensor failure and flawed human-software interaction
2011 explosive release of radioactive material in nuclear power plant meltdowns when cooling pumps lost battery power and backup generators were disabled by the tsunami resulting from the Tōhoku earthquake
https://programm.ard.de/TV/wdrfernsehen/kyrill---ein-orkan-fegt-durchs-land/eid_28111450974388
Ponte Morundi Collapse (2018)
Takoma Narrows Bridge failure (1940)
Gare Montparnasse (1895)
Boiler explosion (1880s—1910s)
Takoma Narrows (1940)
Challenger (1986)
Bhopal disaster (1984, 4000)
E. coli O157 (South Wales outbreak 2005; US 1993)
Fukushima Daiichi nuclear disaster 2011
Chernobyl 1986
Titanic (1912)
Minneapolis bridge (2004)
Mad cow
A radiation therapy machine caused 6 accidents in the late 1980s in which patients were given massive overdoses of radiation, because of a race condition program error. Now a standard case study in software engineering about the overconfidence of the engineers and lack of proper due diligence to resolve reported software bugs.
Therac-25 (late 1980s)
Bhopal disaster (1984)
Mad cow and E. coli O157 (late 1980s—today)
Hyatt Walkway Collapse (1981)
Challenger (1986)
Aloha Airlines Flight 243 (1988)
Chernobyl (1986)
Sampoong Department Store collapse (1995)
1988 explosive decompression at altitude of old aircraft with metal fatigue and possibly improper maintenance
Credits (left to right):
Seoul Metropolitan Fire & Disaster Headquarters, https://en.wikipedia.org/wiki/Sampoong_Department_Store_collapse, CC BY-SA 4.0, http://eungdapso.seoul.go.kr/Community/viewLongText.jsp?key=CD852A52B4760F9E09DBF27ABCB137A3AB33F5A7F1FF35C5B135421B870AD652.
Dennis Oda, AP, https://flightsafety.org/wp-content/uploads/2016/10/asw_mar11_p37-41.pdf, https://www.reddit.com/r/CatastrophicFailure/comments/euwgvz/on_28_april_1988_aloha_airlines_flight_243/, https://twitter.com/djsnm/status/996433183435509760.
Sally Banchbach, used under license.
Carl Montgomery, Flickr, https://en.wikipedia.org/wiki/Chernobyl_disaster#/media/File:Chernobylreactor_1.jpg, CC BY 2.0
FSIS, https://www.fsis.usda.gov/wps/wcm/connect/87caa3f9-0c76-45c7-be4e-84d73151ed9e/RD-2009-0011-072114.pdf?MOD=AJPERES, public domain
Catalina Márquez, https://en.wikipedia.org/wiki/Therac-25, CC BY-SA 4.0
Luca Frediani, https://en.wikipedia.org/wiki/Bhopal_disaster#/media/File:Bhopal-Union_Carbide_2.jpg, CC BY-SA 2.0
Dr. Lee Lowery, Jr., P.E., https://en.wikipedia.org/wiki/Hyatt_Regency_walkway_collapse#/media/File:Hyatt_Regency_collapse_end_view.PNG, public domain
https://en.wikipedia.org/wiki/Tacoma_Narrows_Bridge_(1940), public domain
Levy & fils, https://en.wikipedia.org/wiki/Gare_Montparnasse#/media/File:Train_wreck_at_Montparnasse_1895.jpg, public domain
Jon Bodsworth, Egypt Archive, https://en.wikipedia.org/wiki/Meidum#/media/File:02_meidum_pyramid.jpg, copyright free
Pyramid of Snefru at Maidum (2600 BCE)
Asteriod
Icelandic eruptions affecting European airtravel (2011)
Chilean eruption
California wildfires (2012)
Dought, pestilence
Landslide
Haiti earthquake (2010, 160,000)
Indian Ocean tsunami (2004, 280,000)
Yellow River floods (1931, ~4 million)
Yellow River floods (1931, ~4 million)
| |
Risk assessment applications
heavy metals, pesticides, PMx, ozone, PCBs, EMF, RF, etc., thermal pollution
traffic safety, bridge design, airplanes, spacecraft, nuclear plants, reliability of electricity distribution
portfolio planning, consultation, instrument evaluation
manufacturing and factory workers, farm workers, hospital staff
benzene in Perrier, E. coli in beef, salmonella in tomato, children’s toys
endangered species, fisheries and reserve management, entrainment and impingement
$
More and more important
Risk analysis
Throughout the nuclear industry
Risk
Risk as product
Risk curve
Badness
-1
10
0
1
Probability
0
Mortality
Financial loss
Sometimes called a Farmer curve (after the British nuclear regulator)
Probability unitless
Percentage per demand or outcome
Density per unit outcome
Frequency per unit time or outcome
Qualitative per unit
(e.g., ‘rare’, ‘plausible’, etc.)
Reg said
https://web.archive.org/web/20121218194620/https://netfiles.uiuc.edu/mragheb/www/NPRE%20457%20CSE%20462%20Safety%20Analysis%20of%20Nuclear%20Reactor%20Systems/The%20Risk%20Assessment%20Methodology.pdf
© M. Ragheb1/8/2011
In 1986, IAEA estimated 4,000 projected deaths due to Chernobyl
Risk assessment
Risk management
Old school risk analysis
X
X
X
X
X
X
Carelessness about ensembles is perhaps their fourth mistake
Traditional practice
What’s needed
An assessment should tell us
Probability distribution display
Pr(X ≤ x)
x
x
Pr(x ≤ X)
x
p(x)
0
1
0
1
CDF is not core damage frequency!
Probabilistic models
wrinkles
Ensemble
Specification of ensembles
Combining distributions (convolution)
X
X+Y
=
+
Y
Convolution is getting the output distribution from the input distributions. Monte Carlo is one way to do this.
Integrated calculations
Z
Y
X
W
f(W,X,Y,Z)
This is unlike intrusive analyses that solve operations sequentially
Sample realizations from each input distribution, and combine them to compute a scalar output value. Repeat many times to accumulate a distribution of output values.
Typical risk analysis problems
Case study: dike reliability
The inputs
Relative density of the revetment blocks
Block thickness
D = uniform(0.68, 0.72) meters
Slope of the revetment
Analysts’ “model uncertainty” factor
M = uniform(3.0, 5.2)
Significant wave height (average of the highest third of waves)
H = weibull(1.35 meters, 11)
Offshore peak wave steepness
s = normal(0.04, 0.0055)
How do we specify these distributions?
Analysis
0.6
0.7
0.8
D
0.2
0.3
0.4
0
1
α
3
5
2
4
6
M
0.02
0.04
0.06
s
1.5
1.6
1.7
1
0
Δ
Cumulative probability
Z
Probability of Z being negative
(i.e., dike failure) is very small
1
0
2
H
0
1
0
1
(meters)
Cumulative probability
Z
0
1
0
P = 0.000015
Script for the R programming language
# case study for reliability of the dike
m = 1000000 # number of replications
delta = runif(m, 1.60, 1.65) # density of the revetment blocks
D = runif(m, 0.68, 0.72) # block thickness
alpha = atan(runif(m, 0.32, 0.34)) # slope of the revetment
M = runif(m, 3.0, 5.2) # “model uncertainty” factor
H = rweibull(m, 11,1.35) # significant wave height
s = rnorm(m, 0.04, 0.0055) # offshore peak wave steepness
Z = delta * D - (H* tan(alpha)) / (cos(alpha)*M*sqrt(s))
plot(ecdf(Z)); sum(Z<0)/m
Google “download R”
So how are we doing?
Dramatic technocratic failures
These are not noble failures
Famous ‘misunderestimations’
1998 National Weather Service, US Army Corps
(RAWG 2005; USCPSOTF 2006)
We are understating risks and uncertainties
What are we doing wrong?
Wishful thinking
These mistakes can compound
Northeast Blackout of 2003
I reviewed a paper for Risk Analysis on reliability theory by candle light.
Almost 200,000 km of lines, operated by 500 separate companies
Dependence
Dependence
Dependence can be complex
X
Y
Dependence is subtler than correlation
Knowing the correlation doesn’t tell you the dependence. Each plot has (Pearson) correlation 0.816
op <- par(las=1, mfrow=c(2,3), mar=1.5+c(4,4,1,1), oma=c(0,0,0,0),
lab=c(6,6,7),
#cex.lab=2.0, cex.axis=1.3,
mgp=c(3,1,0))
ff <- y ~ x
for(i in 1:4) {
ff[[2]] <- as.name(paste("y", i, sep=""))
ff[[3]] <- as.name(paste("x", i, sep=""))
lmi <- lm(ff, data= anscombe)
xl <- substitute(expression(x[i]), list(i=i))
yl <- substitute(expression(y[i]), list(i=i))
plot(ff, data=anscombe, col="red", pch=21, cex=2.4, bg = "orange",
xlim=c(3,19), ylim=c(3,13)
# , xlab=eval(xl), ylab=yl # for version 3
, xlab='', ylab='' # for version 3
)
abline(lmi, col="blue")
}
x = c(4,8,13.269,17.5); y = c(4,10,8,12.5); plot(x,y
,col="red", pch=21, cex=2.4, bg = "orange",
xlim=c(3,19), ylim=c(3,13)
, xlab='', ylab='' # for version 3
); cor(x,y)
abline(lmi, col="blue")
x = c(6,4,8,10,13,17.5,16,14,15,17.3); y = c(5.468,4,6,8,13,9.241,11,11.5,12,10); plot(x,y
,col="red", pch=21, cex=2.4, bg = "orange",
xlim=c(3,19), ylim=c(3,13)
, xlab='', ylab='' # for version 3
); cor(x,y)
abline(lmi, col="blue")
par(op)
http://en.wikipedia.org/wiki/Correlation_and_dependence
http://en.wikipedia.org/wiki/File:Anscombe%27s_quartet_3.svg (four of these six graphs)
This is a featured picture, which means that members of the community have identified it as one of the best images on the English Wikipedia, adding significantly to its accompanying article. If you have a different image of similar quality, be sure to upload it using the proper free license tag, add it to a relevant article, and nominate it.
Case study: Vesely’s tank
What’s the chance the tank ruptures under pumping?
Vesely et al. 1981
Fault tree and its Boolean expression
E1
T
K2
E2
E3
S
E4
S1
E5
K1
or
and
or
or
or
R
E1 = T ∨ (K2 ∨ (S & (S1 ∨ (K1 ∨ R))))
Bill Vesely’s subjective probabilities
E1 = T ∨ (K2 ∨ (S & (S1 ∨ (K1 ∨ R))))
Component
Pressure tank T
Relay K2
Pressure switch S
Relay K1
Timer relay R
On-switch S1
Probability
5 × 10−6
3 × 10−5
1 × 10−4
3 × 10−5
1 × 10−4
3 × 10−5
Different dependency models
Vesely et al. (all variables independent)
E1 = T |∨| (K2 |∨| (S |&| (S1 |∨| (K1 |∨| R))))
Mixed dependencies
E1 = T |∨| (K2 ∨ (S &r (S1 |∨| (K1 /∨/ R))))
Correlated to no assumption about dependence
E1 = T |∨| (K2 ∨ (S & (S1 |∨| (K1 /∨/ R))))
No assumptions about any dependencies
E1 = T ∨ (K2 ∨ (S & (S1 ∨ (K1 ∨ R))))
Results
10−5
10−4
10−3
Probability of tank rupturing under pumping
Mixed dependencies
Relaxing S–(S1,K1,R) dependence
No dependence assumptions
All variables independent
3.5×10−5
[3.499×10−5, 3.504×10−5]
[3.50×10−5, 1.35×10−4]
[3×10−5, 1.4×10−4]
Can’t always assume independence
Kinds of uncertainty
Sources of uncertainties
Two kinds of uncertainty
Monte Carlo methods
and probability theory
Interval arithmetic and
constraint analysis
They can change in different contexts. Distinguishing them may be subtle, sometimes like dividing ice from snow. But it’s essential to do so, as getting hit with a snowball is different from getting hit with a block of ice.
Photo credit: Jason Hollinger, https://commons.wikimedia.org/wiki/File:Feathery_Snow_Crystals_(2217830221).jpg, behind translucent filter
Variability = aleatory uncertainty
Incertitude = epistemic uncertainty
Propagating incertitude
Suppose
A is in [2, 4]
B is in [3, 5]
What can be said about the sum A+B?
4
6
8
10
The right answer for
risk analysis is [5,9]
Must be treated differently
Propagation methods
��Relaxing all assumptions�� Laplace Dependence Shape
4
5
6
7
8
9
10
0
0.5
1
Cumulative distribution of Laplace’s unreasonable triangular density function
Any CDF inside this region is possible…same as interval
This contains any CDF arising from dependent uniforms
Aleatory or epistemic?
“These concepts only make unambiguous sense if they are defined with the confines of a model of analysis. In one model an addressed uncertainty may be aleatory, in another model it may be epistemic. So, the characterization of uncertainty becomes a pragmatic choice dependent on the purpose of the application.” (Der Kiureghian and Ditlevsen 2009)
“The probability of a component failure is 0.01”
(Jon Helton 2009)
Getting the inputs
Expert elicitation
Maximum likelihood
Bayesian inference
PERT
Method of moments
Maximum entropy
Estimation
methods
the great frontier of making things up
A two-front assault
How to specify input distributions
No data
Lots of data
Default distributions
Cooper’s quip
Arguments (frequently abused)
Engineers cannot always get data
Rare events
Random sample data
Just kidding
0
1
0
1
2
Probability density
0
0.002
0.004
0.006
0.008
0.01
0
1
p
Cumulative probability
Zellner
Bayes-Laplace
Jeffreys
Haldane
Means
0.001 Haldane
0.0015 Jeffreys
0.001986 Zellner
0.001996 Bayes-Laplace
Data is “once out of 1000”
Data is “zero out of 1000”
Means
0 Haldane
0.0005 Jeffreys
0.000942 Zellner
0.000998 Bayes-Laplace
0.000
0.002
0.004
0.006
0.008
0.010
0.0
0.2
0.4
0.6
0.8
1.0
p
Cumulative Probability
Random sample data
Confidence structure (c-box)
Parameter value
Probability
0
1
α
β
(α−β)100% confidence interval
θ
Confidence boxes
Estimators
Probability of rare event
p ~ [beta(k, n−k+1), beta(k+1, n−k)]
notation extends the use of tilda
Data
k = 2
n = 10
beta(2, 9)
0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
0.6
0.8
1
beta(3, 8)
Binomial rate
0
.4
.8
0
.2
.4
.6
.8
1
1
.2
.6
p
0/1
1/1
0/0
0/2
1/2
2/2
0/3
1/3
2/3
3/3
0/20
15/20
20/20
0/200
150/200
200/200
0/2000
1500/2000
2000/2000
…
…
…
…
…
…
…
…
…
…
…
…
Notice that the c-boxes in every row partition the vacuous ‘dunno’ interval
Zero out of 10k trials
0 8×10−k
“zero corner”
…
…
…
…
One out of 10k trials
0 101−k
“once corner”
…
…
…
…
Computing with c-boxes
0
0.2
0.4
0.6
0.8
1
0
1
0
0.2
0.4
0.6
0.8
1
0
1
0
0.2
0.4
0.6
0.8
1
0
1
Plan A
249 out of 1,014 failed
Plan B
39 out of 60 failed
Plan C
17 out of 20 failed
What if we tried all three plans independently?
Conjunction (AND)
a = balchbox(100,25)
b = balchbox(60,39)
c = balchbox(20,17)
ABC = a %&% b %&% c
abc = a %|&|% b %|&|% c
0
0.2
0.4
0.6
0.8
1
0
1
Probability that all three plans fail
Has the confidence performance!
Probability
Vesely’s pressurized tank
on-off switch S1 relay K1 timer relay R
0
0
0
0
0.06
0.03
0.006
0
0
0
0
0.015
0.002
0.03
0/2000 3/500 1/150
0/460 7/10000 0/380
tank T relay K2 pressure switch S
0.000
0.005
0.010
0.015
0.020
0.025
0.030
0
1
Top event, tank rupture under pressurization, assuming independence
Has the performance interpretation
Tight bounds or fat bounds
R script
# Vesely's pressurized tank system from sparse sample data assuming independence�many = 10000�constant <- function(b) if (length(b)==1) TRUE else FALSE�precise <- function(b) if (length(b)==many) TRUE else FALSE�leftside <- function(b) if (precise(b)) return(b) else return(b[1:many])�rightside <- function(b) if (precise(b)) return(b) else return(b[(many+1):(2*many)])�sampleindices = function() round(runif(many)*(many-1) + 1)�pairsides <- function(b) {i=sampleindices(); return(env(leftside(b)[i],rightside(b)[i]))}�env <- function(x,y) if ((precise(x) && precise(y))) c(x,y) else stop('env error’)�beta <- function(v,w) if ((v==0) && (w==0)) env(rep(0,many),rep(1,many)) else
if (v==0) rep(0,many) else if (w==0) rep(1,many) else sort(rbeta(many, v, w))�kn <- function(k,n) return(pairsides(env(beta(k, n-k+1), beta(k+1, n-k)))) ��orI <- function(x,y) return(1-(1-x)*(1-y))�andI <- function(x,y) return(x*y)��t = kn(0, 2000)�k2 = kn(3, 500)�s = kn(1, 150)�s1 = kn(0, 460)�k1 = kn(7, 10000)�r = kn(0, 380)��e1 = orI(t, orI(k2, andI(s, orI(s1, orI(k1, r)))))�
�plotbox <- function(b,new=TRUE,col='blue',lwd=2,xlim=range(b[is.finite(b)]),ylim=c(0,1),xlab='',ylab='Probability',...) {� edf <- function (x, col, lwd, ...) {� n <- length(x)� s <- sort(x)� if (2000<n) {s = c(s[1],s[1:(n/5)*5],s[n]); n <- length(s);} # subsetting for the plot to every 5th value� lines(c(s[[1]],s[[1]]),c(0,1/n),lwd=lwd,col=col,...)� for (i in 2:n) lines(c(s[[i-1]],s[[i]],s[[i]]),c(i-1,i-1,i)/n,col=col,lwd=lwd,...)� }� b = ifelse(b==-Inf, xlim[1] - 10, b)� b = ifelse(b==Inf, xlim[2] + 10, b)� if (new) plot(NULL, xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)� if (length(b) < many) edf(b,col,lwd) else� edf(c(min(b),max(b),b[1:min(length(b),many)]),col,lwd)� if (many < length(b)) edf(c(min(b),max(b),b[(many+1):length(b)]),col,lwd)� }��par(mfrow=c(2,4))�plotbox(t); title('tank T')�plotbox(k2); title('relay K2')�plotbox(s); title('pressure switch S')�plotbox(s1); title('on-off switch S1')�plotbox(k1); title('relay K1')�plotbox(r); title('timer relay R')�plotbox(e1); title('top event E1')
# Vesely's pressurized tank system from sparse sample data assuming independence
many = 10000
constant = function(b) if (length(b)==1) TRUE else FALSE
precise = function(b) if (length(b)==many) TRUE else FALSE
leftside = function(b) if (precise(b)) return(b) else return(b[1:many])
rightside = function(b) if (precise(b)) return(b) else return(b[(many+1):(2*many)])
sampleindices = function() round(runif(many)*(many-1) + 1)
pairsides = function(b) {i=sampleindices(); return(env(leftside(b)[i],rightside(b)[i]))}
env = function(x,y) if ((precise(x) && precise(y))) c(x,y) else stop('env error')
plotbox = function(x, xlab='', ylab='Cumulative probability', xlim=range(x), ylim=c(0,1), ...) {
plot(NULL,NULL,xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim, ...);
lines(ecdf(x[(many+1):(2*many)]),...); lines(ecdf(x[1:many]),...) }
beta = function(v,w) if ((v==0) && (w==0)) env(rep(0,many),rep(1,many)) else
if (v==0) rep(0,many) else if (w==0) rep(1,many) else sort(rbeta(many, v, w))
kn = function(k,n) return(pairsides(env(beta(k, n-k+1), beta(k+1, n-k))))
orI = function(x,y) return(1-(1-x)*(1-y))
andI = function(x,y) return(x*y)
t = kn(0, 2000); k2 = kn(3, 500)
s = kn(1, 150); s1 = kn(0, 460)
k1 = kn(7, 10000); r = kn(0, 380)
e1 = orI(t, orI(k2, andI(s, orI(s1, orI(k1, r)))))
plotbox(e1)
Experts opinions on rare event probabilities
In my experience
People using expert elicitation often take the pronouncements made by their informants exactly as they are specified
(If your experience is different, I’d like to chat)
Is there no uncertainty?
Risk communication
See also the Berlin presentation
0
.4
.8
0
.2
.4
.6
.8
1
1
.2
.6
p
0/1
1/1
0/0
0/2
1/2
2/2
0/3
1/3
2/3
3/3
0/20
15/20
20/20
0/200
150/200
200/200
0/2000
1500/2000
2000/2000
…
…
…
…
…
…
…
…
…
…
…
…
Confidence
0
1
Probability
0
1
Express risk as “k out of n”
Risk
k-n c-box
Match a calculated risk to a c-box
Equivalent binomial count
Confidence
0
1
Confidence
0
1
Probability
0
1
Probability
0
1
Amazon Mechanical Turk
Pair A was rated excellent by 10 out of 100 customers.
Pair B was rated excellent by 80 out of 100 customers.
Pair A was rated excellent by 66 out of 198 customers.
Pair B was rated excellent by 2 out of 6 customers.
Pair A was rated excellent by 66 to 88 out of 198 customers.
Pair B was rated excellent by 33 out of 100 customers.
50% of customers rated Pair A as excellent.
Pair B was rated excellent by
50 out of 100 customers.
50% of customers rated Pair B as excellent.
Pair A was rated excellent by 2 out of 4 customers.
Findings
80% 50% rational
10/100 80/100 same sample size
66/198 2/6 same magnitude
[66,88]/198 33/100 even with ambiguity
50% 50/100 prefer natural frequency
2/4 50% unless very uncertain
Sample sizes ~300 “master turkers”
Bayesian networks
Bayesian networks (Bayes nets)
Slide thanks to Diego Hector Estrada Lugo
Bayes nets
Slide thanks to Diego Hector Estrada Lugo
Many applications
Slide thanks to Diego Hector Estrada Lugo
Limited applicability
Probability
Slide thanks to Diego Hector Estrada Lugo
Backcalculation
Engineering design under uncertainty
101
A typical problem
Backcalculation
(⊆ means “is inside of” or “is a subset of”)
Backcalculation finds B such that A+B ⊆ C, from estimates for A and constraints on C
Planning a remediation program
Can’t just solve for concentration
When Concentration is put back into the forward equation, the resulting Dose is wider than planned
prescribed
known
unknown
forward equation
inverted equation
Dose = Concentration × Intake
Concentration = Dose / Intake
Can’t just solve for concentration
When Radiation is put back into the forward equation, the resulting Dose is wider than planned
prescribed
known
unknown
forward equation
inverted equation
Dose = Radiation × Exposure
Radiation = Dose / Exposure
Example
Suppose dose must be less than 3 milligrams
dose = [0, 3] mg
intake = [0.5, 4] liter
conc = dose / intake
[ 0, 6] mg liter-1
dose = conc * intake
[ 0, 24] mg
eight times larger than planned tolerable levels
Untangling backcalculations
B = backcalc(A, C) = [C1 − A1, C2 − A2]
B = factor(A, C) = [C1 / A1, C2 / A2]
Correct solution
Dose must be less than 3 milligrams
dose = [0, 3] mg
intake = [0.5, 4] liter
conc = factor(intake, dose) = [0/0.5, 3/4] mg liter-1
[ 0, 0.75] mg liter-1
Check:
dose = conc * intake
[ 0, 3] mg
Doesn’t work for distributions either
Dose = Concentration × Intake / Bodymass
Concentration = Dose × Bodymass / Intake
When concentration is put back into the forward equation, the resulting dose is wider than planned
prescribed
known
unknown
Planned dose
Body mass
Intake
Concentration
0
2
4
6
8
0.06
0.04
0.02
0.00
0
2
4
6
8
0.4
0.2
0.0
0
2
4
6
8
0.06
0.04
0.02
0.00
0
2
4
6
8
0.03
0.02
0.01
0.00
Concentration = Dose × Bodymass / Intake
Large doses in the realized distribution (arising from the allowed distribution of concentrations) are more common than were originally planned
0
2
4
6
8
0.06
0.05
0.04
0.03
0.02
0.01
0
Dose
planned
realized
Normal approximation
To simulate normal deviates with correlation r, compute
Y1 = Z1 × σ1 + μ1
Y2 = (rZ1 + Z2 √(1−r2)) × σ2 + μ2
where Z1 and Z2 are independent standard normal deviates, and μ and σ denote the respective desired means and standard deviations
Normal approximation with non-normal distributions
Iterative (trial & error) approach
By trial and error, you may be able to find a distribution of concentrations that yields something close to the planned distribution of doses. This is not always possible however.
0.06
0
2
4
6
8
0.05
0.04
0.03
0.02
0.01
0
Dose
Backcalculation with p-boxes
Suppose A + B = C, where
A = normal(5, 1)
C = {0 ≤ C, median ≤ 1.5, 90th %ile ≤ 35, max ≤ 50}
2
3
4
5
6
7
8
0
1
A
0
10
20
30
40
50
60
0
1
C
Getting the answer
-10
0
10
20
30
40
50
0
1
B
Algorithm for additive backcalculation of the kernel for B from the equation C = A + B, where the p-boxes for A and the constraint C are known
N = number of percentiles used in the representation; {e.g. 100}
{algorithm below finds the left side of the p-box for the distribution of B}
left limit on B at lowest percentile = left limit on C at lowest percentile – left limit on A at lowest percentile;
for i = 1 to N do begin
set flag to "not done";
for j = 0 to i–1 do begin
if (left limit on C at ith percentile) £ (left limit on A at [i–j]th percentile) + (left limit on B at the jth percentile)
then set flag to "done";
end; {if}
if flag is "done" then
left limit on B at ith percentile = left limit on B at [i–1]th percentile {the one right below it}
else left limit on B at ith percentile = left limit on C at ith percentile – left limit on A at 0th percentile;
end; {if}
end; {for j}
end; {for i, left bound}
{algorithm below finds the right side of the p-box for the distribution of B}
right limit on B at highest percentile = right limit on C at highest percentile – right limit on A at highest percentile;
for i = N–1 down to 0 do begin
set flag to "not done";
for j = N down to i+1 do begin
if (right limit on C at ith percentile) ³ (right limit on A at [i–j+N]th percentile) + (right limit on B at the jth percentile)
then set flag to "done";
end; {if}
if flag is "done" then
right limit on B at ith percentile = right limit on B at [i+1]th percentile {the one right above it}
else right limit on B at ith percentile = right limit on C at ith percentile – right limit on A at Nth percentile;
end; {if}
end; {for j}
end; {for i, right bound} end. {algorithm}
Check it by plugging it back in
A + B = C* ⊆ C
-10
0
10
20
30
40
50
60
0
1
C*
C
P-boxes
NASA shielding example
How are these constraints expressed with p-boxes?
P-boxes express these constraints
D [Gy]
R
[
s
-
1
Gy
]
0
300000
600000
T [s]
0
1000
2000
3000
0
1
0.0001
0.0011
What are the ordinate axes here?
D
R
T
Exceedance risk
Can’t just invert the equation
Dose = Radiation ×Time / Shielding
Shielding = Radiation × Time / Dose
When Shielding is put into the forward equation, the resulting Dose is wider than planned
prescribed constraint
unknown
known
forward equation
inverted equation
Wrong calculation
S = R × T / D
Tail doses more than 50 times larger than planned!
x
D = R × T / S
0
100000
200000
0
1
D
[
Gy
]
0
100000
200000
0
1
D
[
Gy
]
Exceedance risk
0
0.1
8
0
1
S
0
0.1
8
0
1
S
Correct way uses backcalculation
D = R × T / S
Wait, how could shielding vary?
The all-but-95% range of this p-box answer for S is about [0.95, infinity] (which is a lot thicker than the interval answer [0.72, infinity]) because the all-but-95% range of the allowed dose in this example was only 750 Gy, rather than 1000 Gy.
0
1000
2000
3000
0
1
D
[
Gy
]
0
1000
2000
3000
0
1
D
[
Gy
]
S = 1 / factor(R × T, D)
Exceedance risk
0
1
0
1
8
S
0
1
0
1
8
S
P-boxes better than worst case
Precise distributions don’t work
Backcalculation algebra
e.g., if A × B = C, then B = factor(A, C) = exp(backcalc(ln A, ln C))
When you need for | And you have estimates for | Use this formula to find the unknown |
A + B ⊆ C | A, C B ,C | B = backcalc(A,C) A = backcalc(B,C) |
A – B ⊆ C | A, C B ,C | B = –backcalc(A,C) A = backcalc(–B,C) |
A × B ⊆ C | A, C B ,C | B = factor(A,C) A = factor(B,C) |
A / B ⊆ C | A, C B ,C | B = 1/factor(A,C) A = factor(1/B,C) |
A ^ B ⊆ C | A, C B ,C | B = factor(log A, log C) A = exp(factor(B, log C)) |
to be inside
Conclusion
Conclusions
Decision making under risk and uncertainty
Don’t drink the Kool-Aid 1978
Don’t take the blue pill 1999
Be skeptical
Be sceptical
There’s a lot more
Complain if they shortchange you
Kinds of calculations
Answers can be very different with variability/uncertainty
Better safety assessments and reliability analyses
component failure rates, restoration times, performance indices
Base case Uncertainty
(no uncertainty) considered
Probability of failure 0.009 0.0072
Frequency of failure 3.284 1.6660
Duration of failure 76.799 63.465
Edimu, M., C.T. Gaunt and R. Herman (2011). Using probability distribution functions in reliability analyses. Electric Power Systems Research 81 (4): 915−921. https://doi.org/10.1016/j.epsr.2010.11.022 [Table 6]
Steps in a probabilistic risk analysis
Limitations of Monte Carlo
Common pitfalls of Monte Carlo
Monte Carlo simulation
How?
Why?
Why not?
Probability bounds analysis
How?
Why?
Why not?
Some of the available software tools
More information from Risk Institute
Questions?
NAFEMS challenge
�NAFEMS is the International Association for the Engineering Modelling, Analysis and Simulation Community
https://en.wikipedia.org/wiki/Simulation_governance
Simulation governance is a managerial function concerned with assurance of reliability of information generated by numerical simulation. The term was introduced in 2011 [1] and specific technical requirements were addressed from the perspective of mechanical design in 2012.[2] Its strategic importance was addressed in 2015.[3][4] At the 2017 NAFEMS World Congress in Stockholm simulation governance was identified as the first of eight “big issues” in numerical simulation.
Woodpile
Traditional practice
Traditional uncertainty analyses
Untenable assumptions
Need ways to relax assumptions
Incertitude in data measurements
When did the fish in my aquarium die during the night?
Coarse measurements, measurements from digital readouts
Chemical detection limits, studies prematurely terminated
Epidemiological or medical information, census data
Concentrations, solubilities, probabilities, survival rates
Presumed or hypothetical limits in what-if calculations
Simmer down…
We didn’t stop using Euclidean geometry when they invented non-Euclidean geometry
What I am saying is…
Bounding probability is an old idea
Several closely related ideas
Bounding
approaches
Euclid
Relax one axiom
Relax one axiom
Richer mathematics, wider array of applications
2
4
6
8
10
0
1
-10
0
10
20
0
1
0
10
20
30
40
50
0
1
2
4
6
8
10
0
1
2
4
6
8
10
0
1
2
4
6
8
10
0
1
0
10
20
30
40
50
0
1
2
4
6
8
10
0
1
2
4
6
8
10
0
1
min,
max,
mode
mean, std
range, quantile
min,
max
min,
max, mean
min, mean
Comparing p-boxes with maximum entropy distributions
sample data, range
min,max, mean,std
p-box
Maximum entropy solutions
When you know Use this shape
{minimum, maximum} uniform
{mean, standard deviation} normal
{minimum, maximum, mode} beta
{minimum, mean} exponential
{minimum, maximum, mean} beta
{min, max, mean, stddev} beta
{min, max, some quantiles} piecewise uniform
{mean, geometric mean} gamma
Comparing IP to Bayesian approach
Bayes fails incertitude
A few grossly uncertain inputs
Uniform
Uniform
Uniform
Uniform
Probability
A lot of grossly uncertain inputs...
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Uniform
Where does this surety come from?
What justifies it?
Probability
“Smoke and mirrors” certainty
Wishful thinking
Risk analysts often make assumptions that are convenient but are not really justified:
You don’t have to do this
Removing wishful thinking creates a p-box
Uncertainties expressible with p-boxes
Can uncertainty swamp the answer?