Embedded Files
Post date: Mar 23, 2014 9:15:40 PM
#################################################################################
# Compute reliability of Vesely's pressurized tank system from sparse sample data assuming independence among the components
# wipe everything from memory
rm(list=ls())
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)
# env <- function(x,y) if ((precise(x) && precise(y))) c(x,y) else stop('env error') # only works for precise distributional inputs
# this is a MUCH better env(); should check that the bad version wasn't screwing up calculations!
env <- function(x,y) c(pmin(sort(leftside(x)), sort(leftside(y)))[sampleindices()],pmax(sort(rightside(x)), sort(rightside(y)))[sampleindices()])
#x = runif(many)
#y = runif(many) * 3 - 1
#plotbox(y,lwd=4)
#plotbox(x,0,lwd=4)
#e = env(x,y)
#plotbox(e,0,col='red')
#plotbox(leftside(e),0,col='green')
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))
pairsides <- function(b) {i = sampleindices(); return(env(leftside(b)[i],rightside(b)[i]))}
# c-box for Bernoulli parameter: x[i] ~ Bernoulli(parameter), x[i] is either 0 or 1, n = length(x), k = sum(x)
kn <- function(k,n) return(pairsides(env(beta(k, n-k+1), beta(k+1, n-k))))
# logical operators
orI <- function(x,y) return(1-(1-x)*(1-y))
andI <- function(x,y) return(x*y)
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)
}
#################################################################################
#################################################################################
#################################################################################
# example of c-boxes with Vesely's risk of the tank rupturing under pressurization
################
# hypothetical data in the form of k failures observed out of n trials are used to create c-boxes using the formula kn(k,n)
t = kn(0, 2000)
k2 = kn(3, 500)
s = kn(1, 150)
s1 = kn(0, 460)
k1 = kn(7, 10000)
r = kn(0, 380)
################
# compute failure probability of top event (tank rupturing under pressurization)
# assuming mutual independence
#e2 = t %|||% (k2 %|||% (s %|&|% (s1 %|||% (k1 %|||% r))))
e1 = orI(t, orI(k2, andI(s, orI(s1, orI(k1, r)))))
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')
We can also check whether these calculations with confidence actually exhibit the promised confidence. And the simulations described below indicate that indeed they do.
We evaluated the statistical performance of using c-boxes to compute the risk of the tank rupturing under pressure. We found that the calculation singhs [lower, right plot], which means that the confidence intervals derived from the top event c-box [third of the lower plots] covers the hypothetically true probability at least at the nominal rate. The R code to show this is to the right of the displayed area of the slide. This analysis is based on fixed hypothetical values for the probabilities of each of the six inputs which are depicted as vertical lines in the plots. For reference, the analogous Bayesian posteriors based on Jeffreys priors for the size inputs are shown as gray beta distributions, and the resulting posterior predictive distribution for the probability of the top event (rupture under pressure) is also shown in gray in the top event plot.
We evaluated the statistical performance of using c-boxes to compute the risk of the tank rupturing under pressure. We found that the calculation singhs [lower, right plot], which means that the confidence intervals derived from the top event c-box [third of the lower plots] covers the hypothetically true probability at least at the nominal rate. The R code to show this is to the right of the displayed area of the slide. This analysis is based on fixed hypothetical values for the probabilities of each of the six inputs which are depicted as vertical lines in the plots. For reference, the analogous Bayesian posteriors based on Jeffreys priors for the size inputs are shown as gray beta distributions, and the resulting posterior predictive distribution for the probability of the top event (rupture under pressure) is also shown in gray in the top event plot.
The small violations are the result of using small numbers of replications in the Monte Carlo simulation used to make the singh calculations. If the integer variables ‘many’ and ‘some’ are increased, the violations disappear, and also the ‘teeth’ of the curves approach the red uniform distribution more closely.
The small violations are the result of using small numbers of replications in the Monte Carlo simulation used to make the singh calculations. If the integer variables ‘many’ and ‘some’ are increased, the violations disappear, and also the ‘teeth’ of the curves approach the red uniform distribution more closely.
Here is the R code used to make these plots:
###########################################################################################
# Bill Vesely's pressurized tank, which is a BOOLEAN EXPRESSION UNDER INDEPENDENCE
# wipe everything from memory
rm(list=ls())
# replications, trials, and sample sizes
some = 10000 # increase for more accuracy
many = 10000 # increase for more accuracy
# infrastruture
constant <- function(b) if (length(b)==1) TRUE else FALSE
precise <- function(b) if (constant(b) || (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) c(pmin(sort(leftside(x)), sort(leftside(y)))[sampleindices()],pmax(sort(rightside(x)), sort(rightside(y)))[sampleindices()])
# distribution constructors
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))
bernoulli <- function(p) runif(many) < p
betabinomial <- function(size,v,w) rbinom(many, size, beta(v,w))
negativebinomial <- function(size,prob) rnbinom(many, size, prob)
gamma <- function(shape,rate=1,scale=1/rate) {rate = 1/scale; rgamma(many, shape,rate)}
gammaexponential <- function(shape,rate=1,scale=1/rate) {rate = 1/scale; rexp(many,gamma(shape,rate))}
normal <- function(m,s) rnormal(many, m,s)
student <- function(v) rt(many, v)
chisquared <- function(v) rchisq(many, v)
lognormal <- function(m,s) rlnorm(many, m,s)
uniform <- function(a,b) runif(many, a,b)
histogram <- function(x) x[trunc(runif(many)*length(x))+1]
mixture <- function(x,w) {r=sort(runif(many)); x=c(x[1],x); w=cumsum(c(0,w))/sum(w); u=NULL; j=length(x); for (p in rev(r)) {repeat {if (w[[j]] <= p) break; j = j - 1}; u=c(x[[j+1]],u); }; u[order(runif(length(u)))] }
# confidence structures
bern1 <- function(x, n=length(x), k=sum(x)) return(pairsides(env(beta(k, n-k+1), beta(k+1, n-k)))) # x[i] ~ Bernoulli(probability), x[i] is either 0 or 1
bern2 <- function(x, n=length(x), k=sum(x)) return(env(beta(k, n-k+1), beta(k+1, n-k)))} # x[i] ~ Bernoulli(probability), x[i] is either 0 or 1
Bzero <- 1e-6; Bone <- 1-Bzero;
Bzero = 0; Bone = 1
km <- function(k,m) {
if ((k < 0) || (m < 0)) stop('Improper arguments to function km')
return(pmin(pmax(env(beta(k,m+1),beta(k+1,m)),Bzero),Bone))
}
KN = function(k,n) {
if ((k < 0) || (n < k)) stop('Improper arguments to function KN')
return(pmin(pmax(env(beta(k,n-k+1),beta(k+1,pmax(0,n-k))),Bzero), Bone))
}
# other estimators
freqkm <- function(k,m) return(1/(1+m/k))
freqKN <- function(k,n) return(k/n)
jeffkm = function(k,m) return(beta(k+0.5, m+0.5))
jeffKN = function(k,n) return(beta(k+0.5, n-k+0.5))
# logical operators
orI <- function(x,y) return(1-(1-x)*(1-y))
andI <- function(x,y) return(x*y)
# arithmetic operations
OPi = function(x,y,op) { # op can be '+', '-', '*', '/', '^', 'pmin', or 'pmax'
nx = length(x)
ny = length(y)
if (op=='-') return(OPi(x, c(-y[(many+1):ny], -y[1:many]),'+'))
if (op=='/') return(OPi(x, c(1/y[(many+1):ny], 1/y[1:many]),'*'))
if ((nx==1) && (ny==1)) return(do.call(op,list(x,y)))
if (nx==1) return(c(do.call(op,list(x,y[1:many])), do.call(op,list(x,y[(many+1):ny]))))
if (ny==1) return(c(do.call(op,list(x[1:many],y)), do.call(op,list(x[(many+1):nx],y))))
c(do.call(op,list(x[1:many],y[1:many])), do.call(op,list(x[(many+1):nx],y[(many+1):ny])))
}
opi = function(x,y,op) { # in case it is not obvious what OPi is doing
nx = length(x)
ny = length(y)
if ((nx==1) && (ny==1)) return(do.call(op,list(x,y)))
if (nx==1) return(opi(rep(x,2*many),y))
if (ny==1) return(opi(x,rep(y,2*many)))
if (op=='+') return(c(x[1:many]+y[1:many], x[(many+1):nx]+y[(many+1):ny]))
if (op=='-') return(opi(x, c(-y[(many+1):ny], -y[1:many]),'+'))
if (op=='*') return(c(x[1:many]*y[1:many], x[(many+1):nx]*y[(many+1):ny]))
if (op=='/') return(opi(x, c(1/y[(many+1):ny], 1/y[1:many]),'*'))
if (op=='^') return(c(x[1:many]^y[1:many], x[(many+1):nx]^y[(many+1):ny]))
if ((op=='min') || (op=='pmin')) return(c(pmin(x[1:many],y[1:many]), pmin(x[(many+1):nx],y[(many+1):ny])))
if ((op=='max') || (op=='pmax')) return(c(pmax(x[1:many],y[1:many]), pmax(x[(many+1):nx],y[(many+1):ny])))
stop('ERROR unknown operator in opi')
}
ci <- function(b, c=0.95, alpha=(1-c)/2, beta=1-(1-c)/2) {
left = sort(b[1:many])[round(alpha*many)]
if (many < length(b)) right = sort(b[(many+1):length(b)])[round(beta*many)] else
right = sort(b[1:many])[round(beta*many)]
c(left,right)
}
edf <- function (x, col='gray', lwd=3, ...) {
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,...)
#lines(s, 1:length(x)/length(x),type = 's',col=col,lwd=lwd,...)
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,...)
}
plotbox <- function(b,new=TRUE,col='blue',lwd=2,xlim=range(b[is.finite(b)]),ylim=c(0,1),xlab='',ylab='Probability',...) {
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)
}
# top event failure probability (tank rupture under pressure) assuming independence
#e = t %|||% (k2 %|||% (s %|&|% (s1 %|||% (k1 %|||% r))))
e = orI(t, orI(k2, andI(s, orI(s1, orI(k1, r)))))
f = function(x1,x2,x3,x4,x5,x6) orI(x1, orI(x2, andI(x3, orI(x4, orI(x5, x6)))))
# c-boxes using hypothetical data in the form of k failures observed out of n trials
Kt = 0; Nt = 2000
Kk2 = 3; Nk2 = 500
Ks = 1; Ns = 150
Ks1 = 0; Ns1 = 460
Kk1 = 7; Nk1 = 10000
Kr = 0; Nr = 380
t = KN(Kt,Nt)
k2 = KN(Kk2,Nk2)
s = KN(Ks,Ns)
s1 = KN(Ks1,Ns1)
k1 = KN(Kk1,Nk1)
r = KN(Kr,Nr)
e = f(t,k2,s,s1,k1,r)
# naive frequentist estimate
freqt = freqKN(Kt,Nt)
freqk2 = freqKN(Kk2,Nk2)
freqs = freqKN(Ks,Ns)
freqs1 = freqKN(Ks1,Ns1)
freqk1 = freqKN(Kk1,Nk1)
freqr = freqKN(Kr,Nr)
freqe = f(freqt,freqk2,freqs,freqs1,freqk1,freqr)
# Bayesian estimate
jefft = jeffKN(0, 2000)
jeffk2 = jeffKN(3, 500)
jeffs = jeffKN(1, 150)
jeffs1 = jeffKN(0, 460)
jeffk1 = jeffKN(7, 10000)
jeffr = jeffKN(0, 380)
jeffe = f(jefft,jeffk2,jeffs,jeffs1,jeffk1,jeffr)
par(mfrow=c(2,4))
showem = function(t,jefft,freqt,tlab) {plotbox(t); edf(jefft,col='gray'); abline(v=freqt); title(tlab)}
showem(t, jefft,freqt+0.5/Nt,'tank T')
showem(k2, jeffk2,freqk2,'relay K2')
showem(s, jeffs,freqs,'pressure switch S')
showem(s1, jeffs1,freqs1+0.5/Ns1,'on-off switch S1')
showem(k1, jeffk1,freqk1,'relay K1')
showem(r, jeffr,freqr+0.5/Nr,'relay K1')
showem(e, jeffe,freqe,'top event E')
# reasonable true probabilities and sample sizes
p = c(freqt+0.5/Nt,freqk2,freqs,freqs1+0.5/Ns1,freqk1,freqr+0.5/Nr) # zeros are bumped up
N = c(Nt,Nk2,Ns,Ns1,Nk1,Nr)
true = f(p[[1]],p[[2]],p[[3]],p[[4]],p[[5]],p[[6]])
cb = function(i) KN(rbinom(1,N[[i]],p[[i]]), N[[i]])
LEFT = 1:many
RIGHT = (many+1):(2*many)
u = w = array(0, some)
for (i in 1:some) {
d = f(cb(1),cb(2),cb(3),cb(4),cb(5),cb(6))
u[i] = sum((d[LEFT] < true)) / many
w[i] = sum((d[RIGHT] < true)) / many
}
plot(c(0,1),c(0,1),col='red',type='l', lwd=2)
edf(u)
edf(w,col='blue')
title('Singh plot')
# also consider various, perhaps unreasonable values for the probabilities and sample sizes
par(mfrow=c(1,1))
K = 6
some = 1000 # increase for more accuracy
many = 1000 # increase for more accuracy
LEFT = 1:many
RIGHT = (many+1):(2*many)
plot(c(0,1),c(0,1),col='red',type='l', lwd=2)
title('Singh plot for various probabilities and sample sizes')
for (m in 1:100) {
p = runif(K)
N = round(200*runif(K))+1
true = f(p[[1]],p[[2]],p[[3]],p[[4]],p[[5]],p[[6]])
cat(N,'\n')
cat(p,true,'\n')
u = w = array(0, some)
for (i in 1:some) {
d = f(cb(1),cb(2),cb(3),cb(4),cb(5),cb(6))
u[i] = sum((d[LEFT] < true)) / many
w[i] = sum((d[RIGHT] < true)) / many
}
edf(u)
edf(w,col='blue')
}
lines(c(0,1),c(0,1),col='red',type='l', lwd=2)
areametric = function(u) sum(abs(sort(u)-((1:length(u))/length(u))))/length(u)
Page updated
Google Sites
Report abuse