Post date: Mar 25, 2014 3:14:52 PM
#################################################################################
# wipe everything from memory
rm(list=ls())
many = 500000
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(env(beta(k, n-k+1), beta(k+1, n-k)))
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',cex.axis=1.3,cex.lab=1.5,...) {
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,cex.axis=cex.axis,cex.lab=cex.lab)
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)))))
png("c:/users/visiteur/desktop/sallak/sallak pictures/tank4.png")
par(mfrow=c(2,3))
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')
dev.off()
png("c:/users/visiteur/desktop/sallak/sallak pictures/tank5.png")
par(mfrow=c(1,1))
plotbox(e1); title('top event E1')
dev.off()
# remember e1 for later confidence plots
E1 = e1
png("c:/users/visiteur/desktop/sallak/sallak pictures/tank6.png")
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')
dev.off()
# file c:/users/visiteur/desktop/sallak/sallak pictures/tank7.png
# comes from slides "2 interval.ppt"
#################################################################################
# Compute reliability of Vesely's pressurized tank system from sparse sample data
# without making any assumptions about dependence among the components
# the Frechet logical operators, where x = P(X), y = P(Y)
# P(X & Y) = [max(0, x + y – 1), min(x, y)]
# P(X or Y) = [max(x, y), min(1, x + y)]
andF <- function(x,y) env(pmax(0,leftside(x)+leftside(y)-1), pmin(rightside(x),rightside(y)))
orF <- function(x,y) env(pmax(leftside(x), leftside(y)), pmin(1, rightside(x)+rightside(y)))
################
# compute failure probability of top event (tank rupturing under pressurization)
# without making any assumptions about the dependence among components
#e2 = t %|% (k2 %|% (s %&% (s1 %|% (k1 %|% r))))
e1 = orF(t, orF(k2, andF(s, orF(s1, orF(k1, r)))))
png("c:/users/visiteur/desktop/sallak/sallak pictures/tank8.png")
par(mfrow=c(1,1))
plotbox(e1); title('top event E1')
dev.off()
png("c:/users/visiteur/desktop/sallak/sallak pictures/tank9.png")
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')
dev.off()
# 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)
}
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'))
}
png("c:/users/visiteur/desktop/sallak/sallak pictures/tank10.png")
par(mfrow=c(1,1))
showconfidence(E1, 'top event E1')
dev.off()
png("c:/users/visiteur/desktop/sallak/sallak pictures/tank11.png")
showconfidence(E1, alpha=0.0, beta=0.95)
dev.off()
png("c:/users/visiteur/desktop/sallak/sallak pictures/tank12.png")
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)
dev.off()
# independent
mean(leftside(E1))
mean(rightside(E1))
# frechet
mean(leftside(e1))
mean(rightside(e1))
# independent
ci(E1)
ci(E1,alpha=0,beta=0.95)
# frechet
ci(e1)
ci(e1,alpha=0,beta=0.95)