Post date: Apr 07, 2014 2:33:33 PM
###############################################################################
# Compute the reliability of a bridge with structure whose Boolean expression necessarily has repeated variables
# wipe everything from memory
rm(list=ls())
# The bridge system is
#
# +---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 has irreducible repetitions
# of every variable. Simplification to an arithmetic function yields
# ac+bd+ade+bce-abce-acde-abde-bcde-abcd+2abcde.
# 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
See Scott's confusion and Adolphus' derivation of this state expression
Scott
To simplify the logical expression
(a&c)|(b&d)|(a&e&d)|(b&e&c)
to its arithmetic form, let
W = ac
X = bd
Y = aed
Z = bec
then, because the conjunctions under independence are just the products W, X, Y, and Z, the logical expression is
W | X | Y | Z
but a four-fold application of the probabilistic sum (disjunction is the sum minus the product) means that this quadruple disjunction is
W + X + Y + Z - WX - WY - XY - WZ - XZ - YZ + WXZ + WYZ + XYZ - WXYZ
which means the logical expression is equal to this arithmetic expression
ac + bd + aed + bec - acbd - acaed - bdaed - acbec - bdbec - aedbec + acbdbec + acaedbec + bdaedbec - acbdaedbec
Remembering that, for Boolean values, xx = x and reordering the factors within terms yields
ac + bd + ade + bce - abcd - acde - abde - abce - bcde - abcde + abcde + abcde + abcde - abcde
of which four of the last five terms cancel, yielding
ac + bd + ade + bce - abcd - acde - abde - abce - bcde + abcde
Adolphus says this is WRONG, and it's easy to show it's wrong because the final expression doesn't evaluate to
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
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)
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))))
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
bridgeI = function(a,b,c,d,e,many=100000) {
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))))
}
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)
}
a = kn(23,24)
b = kn(23,24) # a and b have the same shape, but don't be tempted to assign them as "a = b = kn(23,24)" which would create or imply perfect dependence
c = kn(14,17)
d = kn(14,17)
e = kn(12,12)
z = bridgeI(a,b,c,d,e)
# mean of the z structure?
c(mean(z[1:many]), mean(z[(many+1):(2*many)]))
par(mfrow=c(2,3))
plotbox(a); title('a')
plotbox(b); title('b')
plotbox(c); title('c')
plotbox(d); title('d')
plotbox(e); title('e')
plotbox(z); title('bridge (assuming independence)')
# reproduce the interval calculation from <<citation>>
a = c(rep(0.7,many), rep(0.8,many))
b = c(rep(0.75,many), rep(0.9,many))
c = c(rep(0.9,many), rep(0.95,many))
d = 0.95
e = 0.82
zz = bridgeI(a,b,c,d,e) # Interval: [0.91556, 0.975327]
plotbox(zz)
bb(.7,.75,.9,.95,.82) # 91556
bb(.8,.9,.95,.95,.82) # 0.975327
# reproduce output picture here <<is the result z correct?>>
###############################################################################
# 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 dependent 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)')