Post date: Mar 11, 2014 5:7:24 PM
#################################################################################
#################################################################################
#################################################################################
#################################################################################
# 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)