| nonParBayesSystemInferencePriorSets {ReliabilityTheory} | R Documentation |
Non-parametric Bayesian posterior predictive system survival inference using sets of priors
Description
Computes a non-parametric Bayesian posterior predictive survival probability
given the survival signature of a system, test data on each of the
components and a set of priors. This is the methodology described in Walter
et al (2017), which extends the method in nonParBayesSystemInference (Aslett et al, 2015) to allow modelling imperfect prior knowledge.
Usage
nonParBayesSystemInferencePriorSets(at.times, survival.signature, test.data,
nLower=2, nUpper=2, yLower=0.5, yUpper=0.5, cores=NA)
Arguments
at.times |
a vector of times at which the posterior predictive estimate of survival probability should be computed. |
survival.signature |
the survival signature matrix of the system/network for which inference is
performed. This should be in the same format as returned by
|
test.data |
a list of vectors containing the component test data. The elements of the
list should be named identically to the component columns in the
|
nLower, nUpper |
the reparameterised lower/upper prior parameter
In all cases, By default the 'uninformative' (but certain) prior with |
yLower, yUpper |
the reparameterised lower/upper prior parameter
In all cases, By default the 'uninformative' (but certain) prior with |
cores |
a scalar indicating how many CPU cores on which to execute parallel parts of
the algorithm (uses the |
Details
This function implements the technique described in Walter et al (2017), which extends the methodology of Aslett et al (2015) to allow modelling partial or imperfect prior knowledge on component failure distributions.
In brief Aslett et al (2015) consider, at any fixed time t, the functioning of a single component of type k to be Bernoulli(p_t^k) distributed for suitable p_t^k, irrespective of the lifetime distribution of the component. Correspondingly, the distribution of the number of components still functioning at time t in a collection of n_k iid components of type k is Binomial(n_k, p_t^k).
Taking the priors p_t^k \sim Beta(\alpha_t^k, \beta_t^k), Aslett et al (2015) show that this leads to a posterior predictive survival distribution with a nice closed form (see equations 9 and 10 in Section 4).
Walter et al (2017) use the standard reparameterisation (dropping sub/super-scripts for readability) n = \alpha+\beta and y = \alpha/(\alpha+\beta). This allows a more natural interpretation, where n represents the prior strength (it represents a pseudo-count for number of failures informing the prior specification) and y represents the prior expectation for the probability a component functions.
In particular, Walter et al (2017) then enable imprecise prior specification by allowing lower and upper bounds on n and y, which may optionally be time varying. This is then propagated to construct bounds on the posterior predictive distribution for the functioning of a new system containing components exchangeable with those provided in the testing set and used in a system with design specified by the survival signature provided.
Value
A list containing two slots, lower and upper, each of which is a vector of the same length as the at.times argument, where each element is the lower/upper posterior predictive probability of a new system surviving to the corresponding time in at.times.
Note
Please feel free to email louis.aslett@durham.ac.uk with any queries or if you encounter errors when running this function.
Author(s)
Louis J. M. Aslett louis.aslett@durham.ac.uk (https://www.louisaslett.com/)
References
Walter, G., Aslett, L. J. M. and Coolen, F. P. A. (2017), ‘Bayesian nonparametric system reliability using sets of priors’, International Journal of Approximate Reasoning, 80, 67–88. Download paper
Aslett, L. J. M., Coolen, F. P. A. and Wilson, S. P. (2015), ‘Bayesian Inference for Reliability of Systems and Networks using the Survival Signature’, Risk Analysis 39(9), 1640–1651. Download paper
See Also
computeSystemSurvivalSignature, and also nonParBayesSystemInference which is the precise counterpart to this method.
Examples
## Exactly reproduce the toy bridge system example in Section 7.2.1 of Walter et al (2017)
# Produces survival signature matrix for one component of type "name", for use
# in nonParBayesSystemInference()
oneCompSurvSign <- function(name){
res <- data.frame(name=c(0,1), Probability=c(0,1))
names(res)[1] <- name
res
}
# Produces data frame with prior and posterior lower & upper component survival
# function for component of type "name" based on
# nonParBayesSystemInferencePriorSets() inputs for all components except
# survival signature; nLower, nUpper, yLower, yUpper must be data frames where
# each column corresponds to the component type, so there must be a match
oneCompPriorPostSet <- function(name, at.times, test.data, nLower, nUpper, yLower, yUpper){
sig <- oneCompSurvSign(name)
nodata <- list(name=NULL)
names(nodata) <- name
nL <- nLower[, match(name, names(nLower))]
nU <- nUpper[, match(name, names(nUpper))]
yL <- yLower[, match(name, names(yLower))]
yU <- yUpper[, match(name, names(yUpper))]
data <- test.data[match(name, names(test.data))]
# NB limit to 1 core on CRAN due to Windows -- make larger to speed up locally!
prio <- nonParBayesSystemInferencePriorSets(at.times, sig, nodata, nL, nU, yL, yU, cores = 1)
post <- nonParBayesSystemInferencePriorSets(at.times, sig, data, nL, nU, yL, yU, cores = 1)
data.frame(Time=rep(at.times,2),
Lower=c(prio$lower,post$lower),
Upper=c(prio$upper,post$upper),
Item=rep(c("Prior", "Posterior"), each=length(at.times)))
}
# ----------------------------------------------
# System
b3 <- createSystem(s -- 2:3 -- 4 -- 5:6 -- 1 -- t, 2 -- 5, 3 -- 6,
types = list(T1 = c(2,3,5,6), T2 = c(4), T3 = c(1)))
# Data
b3nulldata <- list("T1"=NULL, "T2"=NULL, "T3"=NULL)
b3testdata <- list("T1"=c(2.2, 2.4, 2.6, 2.8),
"T2"=c(3.2, 3.4, 3.6, 3.8),
"T3"=(1:4)/10+4) # T3 late failures
b3testdata <- list("T1"=c(2.2, 2.4, 2.6, 2.8),
"T2"=c(3.2, 3.4, 3.6, 3.8),
"T3"=(1:4)/10+0.5) # T3 early failures
b3testdata <- list("T1"=c(2.2, 2.4, 2.6, 2.8),
"T2"=c(3.2, 3.4, 3.6, 3.8),
"T3"=(1:4)-0.5) # T3 fitting failures
b3dat <- reshape2::melt(b3testdata); names(b3dat) <- c("x", "Part")
b3dat$Part <- ordered(b3dat$Part, levels=c("T1", "T2", "T3", "System"))
# Setup to run
b3sig <- computeSystemSurvivalSignature(b3)
b3t <- seq(0, 5, length.out=301)
b3nL <- data.frame(T1=rep(1,301), T2=rep(1,301), T3=rep(1,301))
b3nU <- data.frame(T1=rep(2,301), T2=rep(2,301), T3=rep(4,301))
b3yL <- data.frame(T1=rep(0.001, 301),
T2=rep(0.001, 301),
T3=c(rep(c(0.625,0.375,0.250,0.125,0.010), each=60), 0.01))
b3yU <- data.frame(T1=rep(0.999, 301),
T2=rep(0.999, 301),
T3=c(rep(c(0.999,0.875,0.500,0.375,0.250), each=60), 0.25))
b3T1 <- oneCompPriorPostSet("T1", b3t, b3testdata, b3nL, b3nU, b3yL, b3yU)
b3T2 <- oneCompPriorPostSet("T2", b3t, b3testdata, b3nL, b3nU, b3yL, b3yU)
b3T3 <- oneCompPriorPostSet("T3", b3t, b3testdata, b3nL, b3nU, b3yL, b3yU)
# Compute prior and posterior sets!!
# NB limit to 1 core on CRAN due to Windows -- make larger to speed up locally!
b3prio <- nonParBayesSystemInferencePriorSets(b3t, b3sig, b3nulldata,
b3nL, b3nU, b3yL, b3yU, cores = 1)
b3post <- nonParBayesSystemInferencePriorSets(b3t, b3sig, b3testdata,
b3nL, b3nU, b3yL, b3yU, cores = 1)
b3df <- rbind(data.frame(b3T1, Part="T1"),
data.frame(b3T2, Part="T2"),
data.frame(b3T3, Part="T3"),
data.frame(Time=rep(b3t,2),
Lower=c(b3prio$lower,b3post$lower),
Upper=c(b3prio$upper,b3post$upper),
Item=rep(c("Prior", "Posterior"), each=length(b3t)), Part="System"))
b3df$Item <- ordered(b3df$Item, levels=c("Prior", "Posterior"))
b3df$Part <- ordered(b3df$Part, levels=c("T1", "T2", "T3", "System"))
library("ggplot2")
library("xtable")
ggplot(b3df, aes(x=Time)) +
scale_fill_manual(values = c("#b2df8a", "#1f78b4")) +
scale_colour_manual(values = c("#b2df8a", "#1f78b4")) +
geom_line(aes(y=Upper, group=Item, colour=Item)) +
geom_line(aes(y=Lower, group=Item, colour=Item)) +
geom_ribbon(aes(ymin=Lower, ymax=Upper, group=Item, colour=Item, fill=Item), alpha=0.5) +
facet_wrap(~Part, nrow=2) +
geom_rug(aes(x=x), data=b3dat) +
xlab("Time") +
ylab("Survival Probability") +
theme_bw() +
theme(legend.title = element_blank())
b3sigtable <- b3sig[b3sig$T3 == 1,]
b3sigtable$T1 <- as.factor(b3sigtable$T1)
b3sigtable$T2 <- as.factor(b3sigtable$T2)
b3sigtable$T3 <- as.factor(b3sigtable$T3)
xtable(b3sigtable)