3 Tank confidence

Post date: Apr 01, 2014 11:56:26 AM

#################################################################################

# Compute confidence intervals for the probability of Vesely's tank rupturing under pressurization

# wipe everything from memory

rm(list=ls())

par(mfrow=c(1,1))

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

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)

# confidence intervals from c-boxes and related structures with the confidence interpretation

ci <- function(b, conf=0.95, alpha=(1-conf)/2, beta=1-(1-conf)/2) {

   left = sort(b[1:many])[min(many,max(1,round(alpha*many)))]

   if (many < length(b)) right = sort(b[(many+1):length(b)])[min(many,max(1,round(beta*many)))] else

   right = sort(b[1:many])[round(beta*many)]

   c(left,right)

   }

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

e1 = orI(t, orI(k2, andI(s, orI(s1, orI(k1, r)))))

showconfidence <- function(b, title='', conf=0.95, alpha=(1-conf)/2, beta=1-(1-conf)/2, wh=-0.005) {

  plotbox(b)

  title(title)

  if (beta < alpha) {save = alpha; alpha = beta; beta = save}

  CI = ci(b,alpha=alpha,beta=beta)

  lines(c(-1,CI[[1]]),rep(alpha,2),col='darkgray',lty=3)

  lines(c(-1,CI[[2]]),rep(beta,2),col='darkgray',lty=3)

  lines(rep(CI[[1]],2),c(alpha,wh),col='darkgray',lty=3)

  lines(rep(CI[[2]],2),c(beta,wh),col='darkgray',lty=3)

  lines(CI,rep(wh,2),col='red',lwd=6)

  text(max(CI)*1.2,-0.1,paste((beta-alpha)*100,'% confidence interval'))

  }

par(mfrow=c(2,1))

showconfidence(e1, 'top event E1')

#showconfidence(e1, alpha=0.25, beta=0.75)

showconfidence(e1, alpha=0.0, beta=0.95)