Post date: Mar 11, 2014 5:13:32 PM
#################################################################################
#################################################################################
#################################################################################
#################################################################################
# functions specifically for the bridge calculations
#
# These calculation concern the reliability of a bridge with structure
#
# +---a----+----c---+
# | | |
# ------| e |-------
# | | |
# +---b----+----d---+
#
# where the reliabilities of the five (unrepairable) components are
# given as independent uncertain numbers. This calculation uses
# a combination of monotonicity and Monte Carlo methods to
# take account of independence assumptions and sidestep
# the problem of the repeated variables.
#
# The Boolean function is (a&c)|(b&d)|(a&e&d)|(b&e&c) which
# implies an arithmetic function 1+(ac-1)(1-bd)(1-ade)(1-bec),
# after simplification. Replacing squares (because a*a = a in
# Boolean systmes), yields the function
#
# ac+bd+ade+bce-abce-acde-abde-bcde-abcd+2abcde
#
# which gives the reliability of the system in terms of probabilities
# a, b, c, d, e of its component parts.
#
#################################################################################
BridgeI = function(a,b,c,d,e) (a*c+b*d+a*d*e+b*c*e-a*b*c*e-a*c*d*e-a*b*d*e-b*c*d*e-a*b*c*d+2*a*b*c*d*e)
#################################################################################
# use Monte Carlo to handle independence among variables and
# assume function is monotone increasing in all 5 variables
bridgeI = function(a,b,c,d,e) return(env(BridgeI(L(a),L(b),L(c),L(d),L(e)), BridgeI(R(a),R(b),R(c),R(d),R(e))))
#################################################################################
# compute Birnbaum measures generically, without requiring knowledge about monotone variables
birnbaum <- function(a,b,c,d,e, F) {
rr = sampleindices()
A0 = endpoints(a, 0, which=rr)
A1 = endpoints(a, 1, which=rr)
rr = sampleindices()
B0 = endpoints(b, 0, which=rr)
B1 = endpoints(b, 1, which=rr)
rr = sampleindices()
C0 = endpoints(c, 0, which=rr)
C1 = endpoints(c, 1, which=rr)
rr = sampleindices()
D0 = endpoints(d, 0, which=rr)
D1 = endpoints(d, 1, which=rr)
rr = sampleindices()
E0 = endpoints(e, 0, which=rr)
E1 = endpoints(e, 1, which=rr)
# do all 32 corners so we don't have to figure out which variables the functions are monotone in; in general, you'd have to compute 2^(#dims) corners
F0 = F(A0,B0,C0,D0,E0)
F1 = F(A0,B0,C0,D0,E1)
F2 = F(A0,B0,C0,D1,E0)
F3 = F(A0,B0,C0,D1,E1)
F4 = F(A0,B0,C1,D0,E0)
F5 = F(A0,B0,C1,D0,E1)
F6 = F(A0,B0,C1,D1,E0)
F7 = F(A0,B0,C1,D1,E1)
F8 = F(A0,B1,C0,D0,E0)
F9 = F(A0,B1,C0,D0,E1)
F10 = F(A0,B1,C0,D1,E0)
F11 = F(A0,B1,C0,D1,E1)
F12 = F(A0,B1,C1,D0,E0)
F13 = F(A0,B1,C1,D0,E1)
F14 = F(A0,B1,C1,D1,E0)
F15 = F(A0,B1,C1,D1,E1)
F16 = F(A1,B0,C0,D0,E0)
F17 = F(A1,B0,C0,D0,E1)
F18 = F(A1,B0,C0,D1,E0)
F19 = F(A1,B0,C0,D1,E1)
F20 = F(A1,B0,C1,D0,E0)
F21 = F(A1,B0,C1,D0,E1)
F22 = F(A1,B0,C1,D1,E0)
F23 = F(A1,B0,C1,D1,E1)
F24 = F(A1,B1,C0,D0,E0)
F25 = F(A1,B1,C0,D0,E1)
F26 = F(A1,B1,C0,D1,E0)
F27 = F(A1,B1,C0,D1,E1)
F28 = F(A1,B1,C1,D0,E0)
F29 = F(A1,B1,C1,D0,E1)
F30 = F(A1,B1,C1,D1,E0)
F31 = F(A1,B1,C1,D1,E1)
f0 = pmin(F0,F1,F2,F3,F4,F5,F6,F7,F8,F9,F10,F11,F12,F13,F14,F15,F16,F17,F18,F19,F20,F21,F22,F23,F24,F25,F26,F27,F28,F29,F30,F31)
f1 = pmax(F0,F1,F2,F3,F4,F5,F6,F7,F8,F9,F10,F11,F12,F13,F14,F15,F16,F17,F18,F19,F20,F21,F22,F23,F24,F25,F26,F27,F28,F29,F30,F31)
env(f0,f1)
}
#################################################################################
# inputs for the bridge example
a = kn(58,60) # Mohamed, what should these guys be?
b = kn(59,60)
c = kn(57,60)
d = kn(60,60)
e = kn(54,60)
a = kn(158,160) # Mohamed, what should these guys be?
b = kn(59,60)
c = kn(5,10)
d = kn(6,6)
e = kn(5,6)
#################################################################################
# compute bridge reliability, making use of the unateness of the variables
par(mfrow=c(2,3))
B = bridgeI(a,b,c,d,e)
plotbox(B)
title('Bridge reliability')
# redundant calculation using more general 'corners' function
f <- function(a,b,c,d,e, w,x,y,z) return(BridgeI(a,b,c,d,e))
sc = corners(f, rep(incr,5), a,b,c,d,e) # assumes function is monotone increasing in all 5 variables
plotbox(sc,new=FALSE,col='red')
#################################################################################
# compute Birnbaum measures generically, without requiring knowledge about monotone variables
# Birnbaum function for variable a
F <- function(a,b,c,d,e) return(c-b*c*e-b*c*d+d*e-c*d*e-b*d*e+2*b*c*d*e)
fa = birnbaum(a,b,c,d,e,F)
plotbox(fa)
title('Birnbaum for variable a')
# redundant calculation using more general 'corners' function
f <- function(a,b,c,d,e, w,x,y,z) return(F(a,b,c,d,e))
sc = corners(f, rep(both,5), a,b,c,d,e)
plotbox(sc,new=FALSE,col='red')
# Birnbaum function for variable b
F = function(a,b,c,d,e) (d+c*e-a*c*e-a*d*e-c*d*e-a*c*d+2*a*c*d*e)
fb = birnbaum(a,b,c,d,e,F)
plotbox(fb)
title('Birnbaum for variable b')
# redundant calculation using more general 'corners' function
f <- function(a,b,c,d,e, w,x,y,z) return(F(a,b,c,d,e))
sc = corners(f, rep(both,5), a,b,c,d,e)
plotbox(sc,new=FALSE,col='red')
# Birnbaum function for variable c
F = function(a,b,c,d,e) (a+b*e-a*b*e-a*d*e-b*d*e-a*b*d+2*a*b*d*e)
fc = birnbaum(a,b,c,d,e,F)
plotbox(fc)
title('Birnbaum for variable c')
# redundant calculation using more general 'corners' function
f <- function(a,b,c,d,e, w,x,y,z) return(F(a,b,c,d,e))
sc = corners(f, rep(both,5), a,b,c,d,e)
plotbox(sc,new=FALSE,col='red')
# Birnbaum function for variable d
F = function(a,b,c,d,e) (b+a*e-a*c*e-a*b*e-b*c*e-a*b*c+2*a*b*c*e)
fd = birnbaum(a,b,c,d,e,F)
plotbox(fd)
title('Birnbaum for variable d')
# redundant calculation using more general 'corners' function
f <- function(a,b,c,d,e, w,x,y,z) return(F(a,b,c,d,e))
sc = corners(f, rep(both,5), a,b,c,d,e)
plotbox(sc,new=FALSE,col='red')
# Birnbaum function for variable e
F = function(a,b,c,d,e) (a*d+b*c-a*b*c-a*c*d-a*b*d-b*c*d+2*a*b*c*d)
fe = birnbaum(a,b,c,d,e,F)
plotbox(fe)
title('Birnbaum for variable e')
# redundant calculation using more general 'corners' function
f <- function(a,b,c,d,e, w,x,y,z) return(F(a,b,c,d,e))
sc = corners(f, rep(both,5), a,b,c,d,e)
plotbox(sc,new=FALSE,col='red')
#################################################################################
# compute Birnbaum measure for variable a, but taking into account monotonicity of some variables
F <- function(A,B,C,D,E) return(C - B * C * E - B * C * D + D * E - C * D * E - B * D * E + 2 * B * C * D * E)
for (s in 0:1) {
A = 0
B = endpoints(b, !s, sampleindices()) # monotone decreasing
C = endpoints(c, s, sampleindices()) # monotone increasing
rr = sampleindices()
D0 = endpoints(d, 0, rr)
D1 = endpoints(d, 1, rr)
rr = sampleindices()
E0 = endpoints(e, 0, rr)
E1 = endpoints(e, 1, rr)
# do all four corners; in general, you'd have to compute 2^(#dims) corners
F0 = F(A,B,C,D0,E0)
F1 = F(A,B,C,D0,E1)
F2 = F(A,B,C,D1,E0)
F3 = F(A,B,C,D1,E1)
if (s==0) f0 = pmin(F0,F1,F2,F3)
if (s==1) f1 = pmax(F0,F1,F2,F3)
}
fa2 = env(f0,f1)
plotbox(fa2, col='red')
# redundant calculation using more general 'corners' function
f <- function(a,b,c,d,e, w,x,y,z) return(F(a,b,c,d,e))
sc = corners(f, c(gone,decr,incr,both,both), a,b,c,d,e)
plotbox(sc,new=FALSE,col='red')
#################################################################################
# compute Birnbaum measure for variable b
F = function(a,b,c,d,e) (d+c*e-a*c*e-a*d*e-c*d*e-a*c*d+2*a*c*d*e)
for (s in 0:1) {
A = endpoints(a, !s, sampleindices()) # monotone decreasing
B = 0
rr = sampleindices()
C0 = endpoints(c, 0, rr)
C1 = endpoints(c, 1, rr)
D = endpoints(d, s, sampleindices()) # monotone increasing
rr = sampleindices()
E0 = endpoints(e, 0, rr)
E1 = endpoints(e, 1, rr)
# do all four corners; in general, you'd have to compute 2^(#dims) corners
F0 = F(A,B,C0,D,E0)
F1 = F(A,B,C0,D,E1)
F2 = F(A,B,C1,D,E0)
F3 = F(A,B,C1,D,E1)
if (s==0) f0 = pmin(F0,F1,F2,F3)
if (s==1) f1 = pmax(F0,F1,F2,F3)
}
fb2 = env(f0,f1)
plotbox(fb2,new=FALSE,col='green')
title('a (red) and b (green)')
# redundant calculation using more general 'corners' function
f <- function(a,b,c,d,e, w,x,y,z) return(F(a,b,c,d,e))
sc = corners(f, c(decr,gone,both,incr,both), a,b,c,d,e)
plotbox(sc,new=FALSE,col='green')
###############################################################################
# The Boolean function (a&c)|(b&d)|(a&e&d)|(b&e&c) has repetitions, but it
# cannot be evaluated using the same paired-sides approach that we use for
# the arithmetic function ac+bd+ade+bce-abce-acde-abde-bcde-abcd+2abcde.
# The plots below illustrate why. The addends of the arithmetic expression
# correspond to non-overlapping regions of the Venn diagram, so they are
# mutually exclusive. The disjunctands of the Boolean expression, however,
# are not independent of one another, but in a way that the Monte Carlo
# calculation does not fix.
bb = function(a,b,c,d,e) a*c+b*d+a*d*e+b*c*e-a*b*c*e-a*c*d*e-a*b*d*e-b*c*d*e-a*b*c*d+2*a*b*c*d*e
orI <- function(x,y) return(1-(1-x)*(1-y))
andI <- function(x,y) return(x*y)
BB = function(a,b,c,d,e) orI(orI(andI(a,c), andI(b,d)), orI(andI(a, andI(e,d)), andI(b, andI(e,c))))
BB = function(a,b,c,d,e) orI(orI(a*c, b*d), orI(a*e*d, b*e*c))
BB2 = function(a,b,c,d,e) {
Many = 1000
n = length(a)
w = rep(0, n)
for (i in 1:Many) {
A = runif(n) < a
B = runif(n) < b
C = runif(n) < c
D = runif(n) < d
E = runif(n) < e
w = w + orI(orI(A*C, B*D), orI(A*E*D, B*E*C))
}
w / Many
}
bridgeI = function(a,b,c,d,e,bb) {
L = leftside
R = rightside
return(env(bb(L(a),L(b),L(c),L(d),L(e)), bb(R(a),R(b),R(c),R(d),R(e))))
}
z = bridgeI(a,b,c,d,e,bb)
y = bridgeI(a,b,c,d,e,BB)
x = bridgeI(a,b,c,d,e,BB2)
par(mfrow=c(1,1))
plotbox(z)
plotbox(y, new = FALSE, col='gray')
plotbox(x, new = FALSE, col='red')
title('bridge (assuming independence)')