nonParBayesSystemInferencePriorSets: Non-parametric Bayesian posterior predictive system survival...
In louisaslett/ReliabilityTheory: Structural Reliability Analysis
View source: R/nonParBayesSystemInferencePriorSets.R
nonParBayesSystemInferencePriorSets 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
computeSystemSurvivalSignature.
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
survival.signature argument.
nLower, nUpper
the reparameterised lower/upper prior parameter n for the Beta distribution, where
n = \alpha+\beta. Each must match in type and can be:
a single scalar for a fixed prior across time and component types;
a vector of parameters of the same length as the at.times
argument, which indicates the time-varying prior parameter at the
corresponding time in at.times. This is therefore time-varying,
but indicates the same time-varying prior for all component types;
or a data frame where each column is named using the same names as for
the survival.signature argument and each row corresponds to the
time-varying prior parameter at the corresponding time in
at.times.
In all cases, nUpper but be elementwise greater than or equal to
nLower.
By default the 'uninformative' (but certain) prior with nLower=2 and
nUpper=1 is used for all components at all times.
yLower, yUpper
the reparameterised lower/upper prior parameter y for the Beta distribution, where
y = \alpha/(\alpha+\beta). Each must match in type and can be:
a single scalar for a fixed prior across time and component types;
a vector of parameters of the same length as the at.times
argument, which indicates the time-varying prior parameter at the
corresponding time in at.times. This is therefore time-varying,
but indicates the same time-varying prior for all component types;
or a data frame where each column is named using the same names as for
the survival.signature argument and each row corresponds to the
time-varying prior parameter at the corresponding time in
at.times.
In all cases, yUpper but be elementwise greater than or equal to
yLower.
By default the 'uninformative' (but certain) prior with yLower=0.5
and yUpper=0.5 is used for all components at all times.
cores
a scalar indicating how many CPU cores on which to execute parallel parts of
the algorithm (uses the parallel library internally).
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)
louisaslett/ReliabilityTheory documentation built on Feb. 22, 2024, 8:02 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
View source: R/nonParBayesSystemInferencePriorSets.R
| nonParBayesSystemInferencePriorSets | 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)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.