Altogether

Post date: Mar 11, 2014 5:4:3 PM

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

# infrastructural functions

many = 10000 # increase for more accuracy

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)

L <- function(x,random=TRUE) if (!random) leftside(x) else leftside(x)[sampleindices()]

R <- function(x,random=TRUE) if (!random) rightside(x) else rightside(x)[sampleindices()]

env <- function(x,y) if ((precise(x) && precise(y))) c(x,y) else stop('env error') # serves as env, but 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]))} # might prefer to shuffle indices rather than sample them

kn <- function(k,n) return(pairsides(env(beta(k, n-k+1), beta(k+1, n-k)))) # c-box for Bernoulli parameter: x[i] ~ Bernoulli(parameter), x[i] is either 0 or 1, n = length(x), k = sum(x)

breadth <- function(b) if (precise(b)) 0 else mean(sort(rightside(b))-sort(leftside(b)))

endpoints <- function(b,side, which=1:many) if (constant(b)) return(b) else if (side) return(rightside(b)[which]) else return(leftside(b)[which])

index1 <- function(mon,side) if (mon==2) return(0) else if (mon==1) return(side) else return(1-side)

index2 <- function(mon,side) if (mon==2) return(1) else if (mon==1) return(side) else return(1-side)

corners <- function(f, pr, a1=0,a2=0,a3=0,a4=0,a5=0,a6=0,a7=0,a8=0,a9=0) {

cat('The function will be evaluated',max(2*cumprod(ifelse(pr==0,1,pr))),'x',many,'times\n')

if (length(pr)<9) pr = c(pr,rep(0,9-length(pr)))

f0 = rep(Inf, many)

f1 = rep(-Inf, many)

for (side in 0:1)

for (i1 in index1(pr[[1]],side):index2(pr[[1]],side)) {A1 = endpoints(a1,i1)

for (i2 in index1(pr[[2]],side):index2(pr[[2]],side)) {A2 = endpoints(a2,i2)

for (i3 in index1(pr[[3]],side):index2(pr[[3]],side)) {A3 = endpoints(a3,i3)

for (i4 in index1(pr[[4]],side):index2(pr[[4]],side)) {A4 = endpoints(a4,i4)

for (i5 in index1(pr[[5]],side):index2(pr[[5]],side)) {A5 = endpoints(a5,i5)

for (i6 in index1(pr[[6]],side):index2(pr[[6]],side)) {A6 = endpoints(a6,i6)

for (i7 in index1(pr[[7]],side):index2(pr[[7]],side)) {A7 = endpoints(a7,i7)

for (i8 in index1(pr[[8]],side):index2(pr[[8]],side)) {A8 = endpoints(a8,i8)

for (i9 in index1(pr[[9]],side):index2(pr[[9]],side)) {A9 = endpoints(a9,i9)

ff = f(A1,A2,A3,A4,A5,A6,A7,A8,A9)

#if (side==0) f0 = pmin(f0,ff)

#if (side==1) f1 = pmax(f1,ff)

f0 = pmin(f0,ff)

f1 = pmax(f1,ff)

}}}}}}}}}

env(f0,f1)

}

gone = 0

decr = 0

incr = 1

both = 2

#cat('first pass, decr:',index1(decr,0), index2(decr,0),' second pass, decr:',index1(decr,1), index2(decr,1), '\n')

#cat('first pass, incr:',index1(incr,0), index2(incr,0),' second pass, incr:',index1(incr,1), index2(incr,1), '\n')

#cat('first pass, both:',index1(both,0), index2(both,0),' second pass, both:',index1(both,1), index2(both,1), '\n')

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)

}

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

# 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 MORE TRIALS approach to estimating possible information gain

# via empirical study of the five components of the bridge system

#

# We measure the improvement of uncertainty in the bridge's reliability

# as its fractional reduction when we "pinch" one component's uncertainty

# by hypothetically increasing the number of trials used to compute its

# c-box. This measure is

#

# (unc(B) - unc(P)) / unc(B) = 1 - unc(P) / unc(B)

#

# where B denotes the baseline calculation that includes the uncertainties

# of all the components, and P denotes the pinched calculation with one or

# some of those uncertainties reduced. The unc function characterizes

# the imprecision uncertainty of a c-box or p-box. We use breadth, i.e.,

# the area between the box's left and right bounds, as the unc function.

# The functions below use a Monte-Carlo approach rather than a rigorous

# one in order to simply certain constructions (such as when k=0 or k=n)

# and also to effect the left-right trick in evaluating uncertainties of

# monotone functions. This use is only a matter of convenience; Monte

# Carlo methods are not required for this analysis, for either reason.

# infrastructural functions

#

#many = 10000 # increase for more accuracy

#

#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)])

#

#env <- function(x,y) if ((precise(x) && precise(y))) c(x,y) else stop('env error') # serves as env, but 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 rbeta(many, v, w)

#

#kn <- function(k,n) return(env(beta(k, n-k+1), beta(k+1, n-k))) # c-box for Bernoulli parameter: x[i] ~ Bernoulli(parameter), x[i] is either 0 or 1, n = length(x), k = sum(x)

#

#data = function(p, trials) sum(runif(trials) <= p)

#

#L <- function(x,random=FALSE) if (!random) leftside(x) else leftside(x)[round(runif(many)*(many-1))+1]

#

#R <- function(x,random=FALSE) if (!random) rightside(x) else rightside(x)[round(runif(many)*(many-1))+1]

#

#breadth <- function(b) if (precise(b)) 0 else mean(sort(rightside(b))-sort(leftside(b)))

#

#plotbox <- function(b,new=TRUE,col='blue',lwd=2,xlim=range(b[is.finite(b)]),ylim=c(0,1),xlab='',ylab='Prob',...) {

# edf <- function (x, col, lwd, ...) {

# n <- length(x)

# s <- sort(x)

# 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)

# }

#

# functions specifically for the bridge calculations

#

#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)

#

#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))))

# }

control <- function(p) { # these functions are not suitable for use with intervals or boxes because they have repeated variables

a = p[[1]]

b = p[[2]]

c = p[[3]]

d = p[[4]]

e = p[[5]]

da = c + d*e - b*c*e - c*d*e - b*d*e - b*c*d + 2*b*c*d*e

db = d + c*e - a*c*e - a*d*e - c*d*e - a*c*d + 2*a*c*d*e

dc = a + b*e - a*b*e - a*d*e - b*d*e - a*b*d + 2*a*b*d*e

dd = b + a*e - a*c*e - a*b*e - b*c*e - a*b*c + 2*a*b*c*e

de = a*d + b*c - a*b*c - a*c*d - a*b*d - b*c*d + 2*a*b*c*d

c(da,db,dc,dd,de)

}

improve <- function(p,n,m) {

pa = p[[1]]; na = n[[1]]; ma = m[[1]]

pb = p[[2]]; nb = n[[2]]; mb = m[[2]]

pc = p[[3]]; nc = n[[3]]; mc = m[[3]]

pd = p[[4]]; nd = n[[4]]; md = m[[4]]

pe = p[[5]]; ne = n[[5]]; me = m[[5]]

ka = data(pa,na)

kb = data(pb,nb)

kc = data(pc,nc)

kd = data(pd,nd)

ke = data(pe,ne)

a = kn(ka,na) # to avoid making the c-boxes perfectly dependent, do not use a construction like "a = b = c = d = e = kn(2,60)"

b = kn(kb,nb)

c = kn(kc,nc)

d = kn(kd,nd)

e = kn(kd,ne)

BB = breadth(bridgeI(a,b,c,d,e))

better <- function(k,n,p,m) kn(k+data(p,m), n+m)

BBi = breadth(bridgeI(better(ka,na,pa,ma),better(kb,nb,pb,mb),better(kc,nc,pc,mc),better(kd,nd,pd,md),better(ke,ne,pe,me)))

i = 100*(1-BBi/BB)

cat('investing', sum(m),'samples yields a',i,'% improvement\n')

i

}

study <- function(p,n,effort,new=TRUE) {

N = length(p)

# compute the true Birnbaum importance measures, based on true probabilities

c = control(p)

s = rep(0,N+1)

# allocate the samples to the each component, in turn, to be studied, and to a sixth strategy to share the samples among all components

d = diag(rep(effort,N))

for (i in 1:N) s[[i]] = improve(p,n,d[i,])

sharedeffort = rep(round(effort/N),N)

# in case there are too many samples because effort is not divisible by N, reduce the samples randomly until there's the right number

while (effort < sum(sharedeffort)) { i=round(1+runif(1)*N); sharedeffort[[i]] = sharedeffort[[i]]-1; }

# compute the improvement for each sampling strategy

s[[N+1]] = improve(p,n,sharedeffort)

# output the results as plots and printouts

names = c(letters[1:5], 'share')

myplot <- function(x, y, txt, ...) { if (new) plot(NULL,xlim=range(x),ylim=range(y),...); text(x, y, txt, adj = c(0.5, 0.5), cex=2, ...) }

#myplot(c(p,mean(p)),s,names)

myplot(c(c,mean(c)),s,names,col=c(1:5,40))

cat(rev(names[order(s)]),'\n')

}

studies <- function(p,n,m,xlim=c(0,1),ylim=c(-100,100),count=20) {

plot(NULL,xlim=xlim,ylim=ylim,xlab='True Birnbaum importance',ylab='Uncertainty reduction (%)')

t = 'p ='; for (i in p) t=paste(t,i);title(t)

for (i in 1:count) study(p, n, m, FALSE)

}

samples = 60

moresamples = 60

par(mfrow=c(2,3))

studies(c(0.564, 0.234, 0.685, 0.199, 0.292),rep(samples,5), moresamples, xlim=c(0.05,0.7),ylim=c(-20,30))

studies(c(0.9, 0.9, 0.8, 0.8, 0.99),rep(samples,5), moresamples,xlim=c(0,.25),ylim=c(0,50))

studies(c(0.9, 0.2, 0.2, 0.9, 0.5),rep(samples,5), moresamples,xlim=c(.4,.55),ylim=c(-10,30))

studies(rep(0.5,5), rep(samples,5), moresamples,xlim=c(0.1,0.4),ylim=c(-10,30))

studies(rep(0.99,5),rep(samples,5), moresamples,xlim=c(0.0,0.015),ylim=c(-20,50))

studies(signif(runif(5),3),rep(samples,5), moresamples,xlim=c(0.0,1),ylim=c(-20,50))

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

# Reliability of a noncoherent system

# Consider a pump, whose reliability is p, pumping liquid into

# a vessel monitored by a level sensor whose reliability is s,

# which informs a level controller whose reliability is c.

# The reliability of this system is nonunate. Below is an

# elaboration of the interval calculation of the system's

# reliability

F <- function(a,b,c) return(a * c + b - b * c)

f <- function(a,b,c, u,v,w,x,y,z) return(F(a,b,c))

sc = corners(f, c(both,both,both), a,b,c)

plotbox(sc,new=FALSE,col='red')

# translate input from experts into parameters

q <- function(lambda,mu) return(pexp(500, lambda+mu) / (1 + mu/lambda)) # translate lambda and mu to characterizations of the component unavailability

# do the analysis

pump = q1 = q(1e-3, 0.1)

sensor = q2 = q(2e-3, 5e-2)

controller = q3 = q(3e-3, 1/60)

pump; sensor; controller # 0.00990099, 0.03846154, 0.1525342

f(pump,sensor,controller) # 0.03410508

//////////////////////////////////////////////////////////////////////////////////////

// characterize the uncertainty arbitrarily...plus or minus 30%

pump = measure(q1, 30%)

sensor = measure(q2, 30%)

controller = measure(q3, 30%)

range pump; range sensor; range controller

[ 0.003960396, 0.01584159]

[ 0.01538461, 0.06153847]

[ 0.06101367, 0.2440548]

n = ncs(pump,sensor,controller)

c = corners(pump,sensor,controller)

clear; s = sir(pump,sensor,controller, 20)

m = mca(pump,sensor,controller, 1000)

show s in blue; show n in black; show c in green; show m in red

range(n); range(s); range(c); range(m)

[ 0.0006075025, 0.06446601]

[ 0.01199703, 0.05903612]

[ 0.01259648, 0.05875033]

[ 0.01516103, 0.05562977]

//////////////////////////////////////////////////////////////////////////////////////

// natural (minimal) uncertainty, from significant digits

pump = q1 = q([1e-3 per hour], [0.1 per hour])

sensor = q2 = q([2e-3 per hour], [5e-2 per hour])

controller = q3 = q([3e-3 per hour], 1/[60 hours]) + 0

pump; sensor; controller

[ 0.003322259, 0.02912622]

[ 0.02654867, 0.05263158]

[ 0.1208633, 0.1853325]

[1e-3 per hour]; [0.1 per hour]

[ 0.0005, 0.0015] per hour

[ 0.05, 0.15] per hour

[2e-3 per hour]; [5e-2 per hour]

[ 0.0015, 0.0025] per hour

[ 0.045, 0.055] per hour

[3e-3 per hour]; [60 hours]

[ 0.0025, 0.0035] per hour

[ 55, 65] hours

n = ncs(pump,sensor,controller)

c = corners(pump,sensor,controller)

clear; s = sir(pump,sensor,controller, 20)

m = mca(pump,sensor,controller, 1000)

clear; show s in blue; show n in black; show c in green; show m in red

range(n); range(s); range(c); range(m)

[ 0.01719587, 0.05482085]

[ 0.02199165, 0.05004216]

[ 0.02224406, 0.04979065]

[ 0.0235281, 0.04926684]

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

# Vesely's risk of the tank rupturing under pressurization

# test c-boxes with tank calculation

fixed = TRUE

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

# original 'data' from Vesely et al.

t = 5e-6

k2 = 3e-5

s = 1e-4

s1 = 3e-5

k1 = 3e-5

r = 1e-4

m = rep(5000,6) # number of trials for each of the components above

pfixed = fixed

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

# random failure probability for each component

t = runif(1)

k2 = runif(1)

s = runif(1)

s1 = runif(1)

k1 = runif(1)

r = runif(1)

m = rep(20,6) # number of trials for each of the components above

pfixed = fixed

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

# fixed, small failure probability for each component

t = k2 = s = s1 = k1 = r = 0.01

m = rep(20,6) # number of trials for each of the components above

pfixed = fixed

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

# original 'data' from Vesely et al.

t = 5e-6

k2 = 3e-5

s = 1e-4

s1 = 3e-5

k1 = 3e-5

r = 1e-4

m = rep(5000,6) # number of trials for each of the components above

pfixed = fixed

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

# details

plotting = FALSE

many = 10000 # increase for more accuracy

lots = 5000

al = c(0.005, 0.025, 0.050, 0.125, 0.00, 0.00, 0.00, 0.00)

be = c(0.995, 0.975, 0.950, 0.875, 0.99, 0.95, 0.90, 0.75)

checks = length(al)

cov = rep(0,checks)

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

# overhead

count = function(m) round(rexp(1,1/m))

samp = function(p, trials) sum(runif(trials) <= p)

alphabeta = function(conf=0.95,alpha=(1-conf)/2, beta=1-(1-conf)/2) sort(c(alpha, beta))

env <- function(x,y) c(x,y)

vary <- function(w) if (!w) return(' but these values varied within the simulation') else return('')

# 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)))

# confidence intervals

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

if (alpha==0) left = min(b) else 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)

}

# logical operators

orI <- function(x,y) return(1-(1-x)*(1-y))

andI <- function(x,y) return(x*y)

# distribution constructors

beta <- function(v,w) if (v==0) rep(0,many) else if (w==0) rep(1,many) else rbeta(many, v, w)

bernoulli <- function(p) runif(many) < p

betabinomial <- function(size,v,w) rbinom(many, size, beta(v,w))

# plotting

plotbox <- function(b,new=TRUE,col='blue',lwd=2,xlim=range(b[is.finite(b)]),ylim=c(0,1),xlab='',ylab='Prob',...) {

edf <- function (x, col, lwd, ...) {

n <- length(x)

s <- sort(x)

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)

}

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

# compute failure probability of top event (tank rupturing under pressurization)

#e2 = t %|||% (k2 %|||% (s %|&|% (s1 %|||% (k1 %|||% r))))

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

# how could we simulation the Frechet calculation??

#plot(e1,xlim=c(0,1))

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

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

# generate hypothetical trial sizes

nt = count(m[[1]])

nk2 = count(m[[2]])

ns = count(m[[3]])

ns1 = count(m[[4]])

nk1 = count(m[[5]])

nr = count(m[[6]])

nfixed = fixed

nt = nk2 = ns = ns1 = nk1 = nr = 2000

nfixed = fixed

if (nfixed) cat('n-values: ',nt,nk2,ns,ns1,nk1,nr,'\n')

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

# iterate to compute coverages

for (iteration in 1:lots) {

fixed = FALSE

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

# generate hypothetical sample data

kt = samp(t,nt)

kk2 = samp(k2,nk2)

ks = samp(s,ns)

ks1 = samp(s1,ns1)

kk1 = samp(k1,nk1)

kr = samp(r,nr)

kfixed = fixed

if (kfixed) cat('k-values: ',kt,kk2,ks,ks1,kk1,kr,'\n')

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

# make c-boxes from the hypothetical sample data

T = kn(kt,nt)

K2 = kn(kk2,nk2)

S = kn(ks,ns)

S1 = kn(ks1,ns1)

K1 = kn(kk1,nk1)

R = kn(kr,nr)

if (plotting) {

par(mfrow=c(2,4))

plotbox(T,xlab='T'); title(paste(kt,'out of',nt,'~',signif(t,3)))

plotbox(K2,xlab='K2'); title(paste(kk2,'out of',nk2,'~',signif(k2,3)))

plotbox(S,xlab='S'); title(paste(ks,'out of',ns,'~',signif(s,3)))

plotbox(S1,xlab='S1'); title(paste(ks1,'out of',ns1,'~',signif(s1,3)))

plotbox(K1,xlab='K1'); title(paste(kk1,'out of',nk1,'~',signif(k1,3)))

plotbox(R,xlab='R'); title(paste(kr,'out of',nr,'~',signif(r,3)))

}

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

# project c-boxes through to estimate risk of top event

#E2 = T %|||% (K2 %|||% (S %|&|% (S1 %|||% (K1 %|||% R))))

E1 = orI(T, orI(K2, andI(S, orI(S1, orI(K1, R)))))

E1 = andI(T, andI(K2, andI(S, andI(S1, andI(K1, R)))))

if (plotting) {

plotbox(E1)

lines(c(e1,e1),c(0,1),lwd=4)

}

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

# assess performance

# in theory, we could check infinitely many confidence statements, but we'll limit ourselves to checking just eight:

# confidence intervals at levels of 99%, 95%, 90%, and 75%, in both the symmetric two-sided and the left one-sided forms (((the "left" one is the one with all the outside mass on the right, right?)))

for (ch in 1:checks) {

answer = ci(E1,alpha=al[[ch]],beta=be[[ch]])

# lines(c(Sci95@lo,Sci95@hi), rep(iteration/lots,2), col='red')

if ((left(answer)<=e1) & (e1<=right(answer))) cov[[ch]]=cov[[ch]]+1

}

if (plotting) cat('iterate\n')

}

for (ch in 1:checks) cat(ch,100*cov[[ch]]/lots,'%\n')

cat('number of Monte Carlo replicates per distribution: ',many,'\n')

cat('number of Monte Carlo iterations in the study: ',lots,'\n')

cat('probabilities: ', t, k2, s, s1, k1, r,vary(pfixed),'\n')

cat('n-values: ',nt,nk2,ns,ns1,nk1,nr,vary(nfixed),'\n')

cat('k-values: ',kt,kk2,ks,ks1,kk1,kr,vary(kfixed),'\n')

cat('planned confidence levels: ', be-al,'\n')

cat('actual observed coverages: ', cov/lots ,'\n')

number of Monte Carlo replicates per distribution: 10000

number of Monte Carlo iterations in the study: 5000

probabilities: 0.8061411 0.9631644 0.1903984 0.2276681 0.5869775 0.4176342

n-values: 13 65 6 4 3 3

planned confidence levels: 0.95 0.9 0.75 0.95 0.9 0.75

actual observed coverages: 0.9932 0.9842 0.9300 0.9942 0.9826 0.9384

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

# Coverage simulations to check performance

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

# utility to set alpha and beta, possibly symmetrically from a confidence level

alphabeta = function(confidence=0.95,alpha=(1-confidence)/2,beta=1-(1-confidence)/2) sort(c(alpha, beta))

# Monte Carlo simulation of the coverage performance of the normal mean c-box

coveragenormal.mu = function(n,MU,SIGMA,many=1e4,lots=1e3, ... ) {

ab = alphabeta(...)

mu = seq((MU-5*SIGMA),(MU+5*SIGMA),length.out=many)

coverage = 0

for (i in 1:lots) {

x = rnorm(n, MU, SIGMA)

h = pt(sqrt(n)*(mu-mean(x))/sd(x),n-1) # h = pcnorm.mu(mu, x)

ci = range(mu[(ab[1]<=h) & (h<=ab[2])]) # ci = cinorm.mu(x, alpha=ab[1], beta=ab[2])

if ((ci[1]<=MU)&(MU<=ci[2])) coverage=coverage+1

}

cat('With n=',n,' normal samples at (µ,s)=(',MU,',',SIGMA,'), ',

diff(ab)*100,'% intervals performed at ', 100*coverage/lots,'%\n',sep='')

coverage/lots

}

# Confidence intervals for a normal distribution’s mean

cinorm.mu = function(x, ...) # same as ci(parameter.normal.mean(x))

mean(x)+qt(alphabeta(...),df=length(x)-1)*sd(x)/sqrt(length(x))

# Confidence box for a normal distribution's mean (distribution as a cumulative function of mu-values)

pcnorm.mu = function(mu, x) pt(sqrt(length(x))*(mu-mean(x))/sd(x),length(x)-1)

# examples

coveragenormal.mu(4, 20, 1)

for (s in c(0.1,1,2,10)) coveragenormal.mu(4, 20, s)

for (s in c(0.1,1,2,10)) coveragenormal.mu(n=4, MU=20, SIGMA=s, alpha=0.2, beta=0.8)

# Confidence box for a normal distribution's mean (distribution as a collection of deviates)

rcnorm.mu = function(many, z) # same as parameter.normal.mu(z)

mean(z)+sd(z)*rt(many, length(z)-1)/sqrt(length(z))

# Confidence intervals for the difference of normal means

cinormdiff.mu <- function(x, y, many=2000, ...) {

d = sort(rcnorm.mu(many, x) - rcnorm.mu(many, y))

range(d[round(alphabeta(...)*many)])

}

# Monte Carlo simulation of the coverage performance of the c-box for the difference of normal means

coveragenormal.mudiff=function(n,MU,SIGMA,many=1e4,lots=1e3,...){

ab = alphabeta(...)

truediff = MU[1] - MU[2]

coverage = 0

for (i in 1:lots) {

x = rnorm(n[1], MU[1], SIGMA[1])

y = rnorm(n[2], MU[2], SIGMA[2])

ci=range(sort(rcnorm.mu(many,x)-

rcnorm.mu(many,y))[round(many*ab)])

if ((ci[1] <= truediff) & (truediff <= ci[2])) coverage = coverage + 1

}

cat(' Intended',diff(ab)*100,'%\n','Observed',100*coverage/lots,'%\n')

coverage/lots

}

# examples

plotbox(rcnorm.mu(1000,rnorm(4,5,1)))

plotbox(rcnorm.mu(1000,rnorm(1000,5,1)))

coveragenormal.mudiff(n=c(5,6), MU=c(10,14), SIGMA=c(1,2), confidence=0.95)

# C-box for a constant Y masked by error e of unknown distribution and bias, given independent random samples of Y and e

nonparametric.debiasing <- function(x, error) {# i.e., the c-box for Y such that x[i] = Y + e[i]

z = NULL

for(jj in 1:length(error)) z = c(z, y - error[jj])

z = sort(z)

Q = Get_Q(length(y), length(error))

w = Q / sum( Q )

env(mixture(z,w), mixture(c(z[-1],Inf),w))

}

# Monte Carlo simulation of the coverage performance of the debiased constant c-box

coveragenonparametric.debias <- function(truth=0, err=error, few=4, some=9, many=1000, ...) {

ab = alphabeta(...)

cc = NULL

for (n in 1:many) {

e = err(some+few) # error structure

y = truth + e[(few+1):(some+few)]

e = e[1:few]

z = NULL

for(jj in 1:few) z = c(z,y - e[jj])

z = sort(z)

Q = Get_Q(some, few)

c = cumsum(Q / sum(Q))

j = length(c[c<ab[1]])+1

k = length(c[c<=ab[2]])

if ((c(-Inf , z)[j] <= truth) & (truth <= c(z , Inf)[k])) col=3 else col=4

cc = c(cc,col)

}

length(cc[cc==3])/length(cc) # coverage performance

}

#examples

coveragenonparametric.debias(truth=5, function(many) return(rnorm(many)) ) # unbiased

coveragenonparametric.debias(truth=5, function(many) return(rnorm(many,10)) ) # biased

error <- function(many) return(rbeta(many, 2.0, 8.0) * 4 + 3 ) # this error structure could be anything

coveragenonparametric.debias()

coveragenonparametric.debias(alpha=0.9, beta=0.1)

coveragenonparametric.debias(alpha=0.9, beta=0.0)

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

///////////////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////

// make up some interval values to play with

clear

function rand() {r__=random; return env(r__,random); }

a=rand(); b=rand(); c=rand(); d=rand(); e=rand()

//a;b;c;d;e

a = [ 0.03375244, 0.4723816]

b = [ 0.6900725, 0.9585024]

c = [ 0.7711181, 0.9777833]

d = [ 0.3543395, 0.6956177]

e = [ 0.09057617, 0.6003418]

////////////////////////////////////////////////////////////////////////////////

// "naive" calculation of the Birnbaum measure for the variables a

func F() return C - B * C * E - B * C * D + D * E - C * D * E - B * D * E + 2 * B * C * D * E

B = b

C = c

D = d

E = e

naive = F()

////////////////////////////////////////////////////////////////////////////////

// subinterval reconstitution is slow but rigorous

run brute

//Brute force interval library has been loaded

func map2() return C * (1 - B * ($1 + $2)) + $1 * $2 * (1 - C - B + 2 * B * C)

B = left(b)

C = right(c)

R1 = sir2(d,e,50)

B = right(b)

C = left(c)

R2 = sir2(d,e,50)

R = env(R1,R2)

R // [ 0.11755425, 0.7023882]

[ 0.1175542, 0.7023882]

////////////////////////////////////////////////////////////////////////////////

// Monte Carlo searching yields an inner estimate

function try() {

A = specific(a)

B = specific(b)

C = specific(c)

D = specific(d)

E = specific(e)

return C - B*C*E - B*C*D + D * E * (1 - C - B + 2 * B * C)

}

mc = try()

for i = 1 to 1000 do mc = env(mc, try())

////////////////////////////////////////////////////////////////////////////////

// algebraic rearrangements are rigorous, but limited

a0 = tri(c + d*e - b*c*e - c*d*e - b*d*e - b*c*d + 2*b*c*d*e)

a1 = tri(c - b*c*e - b*c*d + d * e * (1 -c - b + 2 * b * c))

a2 = tri(c * (1 - b * e - b*d) + d * e * (1 - c - b + 2 * b * c))

a3 = tri(c * (1 - b * (d + e)) + d * e * (1 - c - b + 2 * b * c))

a4 = tri(c - b*c*e + d * (e * (1 -c - b + 2 * b * c) - b*c))

a5 = tri(c - b*c*d + e * (d * (1 -c - b + 2 * b * c) - b*c))

a;b;c;d;e

[ 0.03375244, 0.4723816]

[ 0.6900725, 0.9585024]

[ 0.7711181, 0.9777833]

[ 0.3543395, 0.6956177]

[ 0.09057617, 0.6003418]

////////////////////////////////////////////////////////////////////////////////

// checking endpoints seems to work, but I don't quite understand why,

// nor can I believe it will work generally

func endpoint() if ($2==0) return left($1) else return right($1)

procedure assignall() begin

A = endpoint(a,$1)

B = endpoint(b,$2)

C = endpoint(c,$3)

D = endpoint(d,$4)

E = endpoint(e,$5)

end

// get an anchor for the env function

assignall(0,0,0,0,0)

f = F()

// search corners, but don't repeat for the variable a

i = 0 // for i=0 to 0 do

for j=0 to 1 do for k=0 to 1 do for l=0 to 1 do for m=0 to 1 do begin

assignall(i,j,k,l,m)

ff = F()

f = env(f,ff)

end

ends = f

////////////////////////////////////////////////////////////////////////////////

// combine endpoint checking with monotonicity

assignall(0,0,0,0,0)

f = F()

for n = 0 to 1 do begin

B = endpoint(b, n) // monotone increasing

C = endpoint(c, not n) // monotone decreasing

for l = 0 to 1 do for m = 0 to 1 do begin

D = endpoint(d,l)

E = endpoint(e,m)

ff = F()

f = env(f,ff)

end

monends = f

////////////////////////////////////////////////////////////////////////////////

// summarize and compare output intervals

show mc in red

show R

show a0; show a1 in olive; show a2 in blue; show a3 in green; show a4 in teal; show a5 in gray

ff = tri(monends)

show ff in green

range(monends)

[ 0.1258789, 0.699457]

range(ends)

[ 0.1258789, 0.699457]

R

[ 0.1175542, 0.7023882]

mc

[ 0.1605074, 0.6508642]

range(a5);range(a4);range(a3);range(a2);range(a1);range(a0)

[ -0.4162429, 1.059946]

[ -0.4354022, 1.149603]

[ -0.2326928, 1.267755]

[ -0.439358, 1.331206]

[ -1.185823, 1.894515]

naive

[ -1.185823, 1.894515]

///////////////////////////////////////////////////////////////////////////////////////////

// Reliability of a noncoherent system

// Consider a pump, whose reliability is p, pumping liquid into

// a vessel monitored by a level sensor whose reliability is s,

// which informs a level controller whose reliability is c.

// The reliability of this system is nonunate. Below is an

// elaboration of the interval calculation of the system's

// reliability

//////////////////////////////////////////////////////////////////////////////////////

// naive calculation

func ncs() return $1 * $3 + $2 - $2 * $3

//////////////////////////////////////////////////////////////////////////////////////

// corners calculation

func lr() if ($2) return right($1) else return left($1)

func corners() begin

d = lr($3,0)

e = lr($1,0)*d+lr($2,0)*(1-d)

for i = 0 to 1 do for j = 0 to 1 do for k = 0 to 1 do begin

d = lr($3,k)

e = env(e, lr($1,i)*d+lr($2,j)*(1-d))

end

return e

end

//////////////////////////////////////////////////////////////////////////////////////

// subinterval reconstitution

run brute

//Brute force interval library has been loaded

func map3() return ncs($1,$2,$3)

clear

func sir() return sir3($1, $2, $3, $4)

//////////////////////////////////////////////////////////////////////////////////////

// Monte Carlo inner bounds

func mca() begin

mc_a = ncs(specific($1),specific($2),specific($3)) // this line should be replaceable by "mca = empty"

for mc_i = 1 to ($4-1) do begin

mc_a = env(mc_a,ncs(specific($1),specific($2),specific($3)))

// P = specific($1)

// S = specific($2)

// C = specific($3)

// mc_a = env(mc_a,mc__a)

// if (!contains(c,mc__a)) say P,S,C

end

return mc_a

end

//////////////////////////////////////////////////////////////////////////////////////

// translate input from experts into parameters

t = 500 hours

func pexp() return 1 - exp(-$2 * $1) // pexp(x, lambda) yield cumulative probability at x for exponential distribution with mean 1/lambda

func q() return pexp(t, $1+$2) / (1 + $2/$1) // q(lambda, mu) translates lambda and mu values ($1 and $2) to characterizations of the component unavailability

//////////////////////////////////////////////////////////////////////////////////////

// do the analysis

pump = q1 = q(1e-3 per hour, 0.1 per hour)

sensor = q2 = q(2e-3 per hour, 5e-2 per hour)

controller = q3 = q(3e-3 per hour, 1/60 hours) + 0

pump; sensor; controller

0.00990099

0.03846154

0.1525342

ncs(pump,sensor,controller)

0.03410508

//////////////////////////////////////////////////////////////////////////////////////

// characterize the uncertainty arbitrarily...plus or minus 30%

pump = measure(q1, 30%)

sensor = measure(q2, 30%)

controller = measure(q3, 30%)

range pump; range sensor; range controller

[ 0.003960396, 0.01584159]

[ 0.01538461, 0.06153847]

[ 0.06101367, 0.2440548]

n = ncs(pump,sensor,controller)

c = corners(pump,sensor,controller)

clear; s = sir(pump,sensor,controller, 20)

m = mca(pump,sensor,controller, 1000)

show s in blue; show n in black; show c in green; show m in red

range(n); range(s); range(c); range(m)

[ 0.0006075025, 0.06446601]

[ 0.01199703, 0.05903612]

[ 0.01259648, 0.05875033]

[ 0.01516103, 0.05562977]

//////////////////////////////////////////////////////////////////////////////////////

// natural (minimal) uncertainty, from significant digits

pump = q1 = q([1e-3 per hour], [0.1 per hour])

sensor = q2 = q([2e-3 per hour], [5e-2 per hour])

controller = q3 = q([3e-3 per hour], 1/[60 hours]) + 0

pump; sensor; controller

[ 0.003322259, 0.02912622]

[ 0.02654867, 0.05263158]

[ 0.1208633, 0.1853325]

[1e-3 per hour]; [0.1 per hour]

[ 0.0005, 0.0015] per hour

[ 0.05, 0.15] per hour

[2e-3 per hour]; [5e-2 per hour]

[ 0.0015, 0.0025] per hour

[ 0.045, 0.055] per hour

[3e-3 per hour]; [60 hours]

[ 0.0025, 0.0035] per hour

[ 55, 65] hours

n = ncs(pump,sensor,controller)

c = corners(pump,sensor,controller)

clear; s = sir(pump,sensor,controller, 20)

m = mca(pump,sensor,controller, 1000)

clear; show s in blue; show n in black; show c in green; show m in red

range(n); range(s); range(c); range(m)

[ 0.01719587, 0.05482085]

[ 0.02199165, 0.05004216]

[ 0.02224406, 0.04979065]

[ 0.0235281, 0.04926684]

//////////////////////////////////////////////////////////////////////////////////////

// the q() function is has repetitions; are we evaluating it correctly?

// pump

lambda = [1e-3 per hour]

mu = [0.1 per hour]

P = q(lambda, mu)

p1 = q(left(lambda),left(mu))

p2 = q(left(lambda),right(mu))

p3 = q(right(lambda),left(mu))

p4 = q(right(lambda),right(mu))

P; env(p1,p2,p3,p4)

[ 0.003322259, 0.02912622]

// sensor

lambda = [2e-3 per hour]

mu = [5e-2 per hour]

P = q(lambda, mu)

p1 = q(left(lambda),left(mu))

p2 = q(left(lambda),right(mu))

p3 = q(right(lambda),left(mu))

p4 = q(right(lambda),right(mu))

P; env(p1,p2,p3,p4)

[ 0.02654867, 0.05263158]

// controller

lambda = [3e-3 per hour]

mu = 1/[60 hour]

P = q(lambda, mu)

p1 = q(left(lambda),left(mu))

p2 = q(left(lambda),right(mu))

p3 = q(right(lambda),left(mu))

p4 = q(right(lambda),right(mu))

P; env(p1,p2,p3,p4)

[ 0.1208633, 0.1853325]

[ 0.1208752, 0.1853214]

//////////////////////////////////////////////////////////////////////////////////////

// the q() function seems monotone in each variable

// (even though it's not) because the exponential cdf

// in the function hardly affects the result at all

t = 500 hours

func pexp() return 1 - exp(-$2 * $1) // pexp(x, lambda) yield cumulative probability at x for exponential distribution with mean 1/lambda

func q() return pexp(t, $1+$2) / (1 + $2/$1) // q(lambda, mu) translates lambda and mu values ($1 and $2) to characterizations of the component unavailability

func qq() return 1 / (1 + $2/$1) // replace cdf value with "1"

pump = q(1e-3 per hour, 0.1 per hour)

sensor = q(2e-3 per hour, 5e-2 per hour)

controller = q(3e-3 per hour, 1/60 hours) + 0

pump;sensor;controller

0.00990099

0.03846154

0.1525342

pump = qq(1e-3 per hour, 0.1 per hour)

sensor = qq(2e-3 per hour, 5e-2 per hour)

controller = qq(3e-3 per hour, 1/60 hours) + 0

pump;sensor;controller

0.00990099

0.03846154

0.1525424

Pump = qq([1e-3 per hour], [0.1 per hour])

Sensor = qq([2e-3 per hour], [5e-2 per hour])

Controller = qq([3e-3 per hour], 1/[60 hours]) + 0

Pump;Sensor;Controller

[ 0.003322259, 0.02912622]

[ 0.02654867, 0.05263158]

[ 0.1208791, 0.1853361]

Pump = q([1e-3 per hour], [0.1 per hour])

Sensor = q([2e-3 per hour], [5e-2 per hour])

Controller = q([3e-3 per hour], 1/[60 hours]) + 0

Pump;Sensor;Controller

[ 0.003322259, 0.02912622]

[ 0.02654867, 0.05263158]

[ 0.1208633, 0.1853325]

+++

// this is a non-coherent system, so the lower bound is not obtained using the lower bounds of the components

// it works if both a and b are working, or if a is working and b isn't, or if both are not working

// and it doesn't work if neither of these three cases is true

a=0.7

b=0.75

a*b + a*(1-b) + (1-a)*(1-b)

0.775

a = [0.7, 0.8]

b = [0.75, 0.90]

a*b + a*(1-b) + (1-a)*(1-b)

[ 0.6149999, 0.9950001]

func logb( ) _return _ln($1)/_ln($2)

_func andc( ) _begin

s_=_tan(_pi*1 radian*(1-$3)/4)

_if (_abs(s_-1)<1e-8) _return ($1+0|0)|*|($2+0|0) _else _if (_abs(s_)<=1e-10) _return _min($1+0|0,$2+0|0) _else _if (_abs(s_+1)<1e-8) _return _max(0,$1+$2-1) _else _return logb(1+(_exp(($1+0|0)*_ln(s_))-1)|*|(_exp(($2+0|0)*_ln(s_))-1)/(s_-1),s_)

_end

_func orc( ) _begin

s_=_tan(_pi*1 radian*(1-$3)/4)

_if (abs(s_-1)<1e-8) _return 1-(1-($1+0|0))|*|(1-($2+0|0)) _else _if (abs(s_)<=1e-10) _return _max($1+0|0,$2+0|0) _else _if (abs(s_+1)<1e-8) _return _min(1,($1+0|0)+($2+0|0)) _else _return 1-logb(1+(_exp((1-($1+0|0))*_ln(s_))-1)|*|(_exp((1-($2+0|0))*_ln(s_))-1)/(s_-1),s_)

_end

_func orx() {s_ = 0; _for i_ = 1 to _args _do s_ = _min(1,s_ + _arg(i_)); _return s_; }

_func andx() _return 0

orx((a&b), (a & !b), (!a & !b))

[ 0.4499999, 1]

a & b

[ 0.4499999, 0.8]

a & !b

[ 0, 0.25]

a & (not b)

[ 0, 0.25]

!a & !b

[ 0, 0.25]

(not a) & (not b)

[ 0, 0.25]