| Title: | Likelihood Models for Systems with Masked Component Cause of Failure |
|---|---|
| Description: | Maximum likelihood estimation for series systems where the component cause of failure is masked. Implements analytical log-likelihood, score, and Hessian functions for exponential, homogeneous Weibull, and heterogeneous Weibull component lifetimes under masked cause conditions (C1, C2, C3). Supports exact, right-censored, left-censored, and interval-censored observations via composable observation functors. Provides random data generation, model fitting, and Fisher information for asymptotic inference. See Lin, Loh, and Bai (1993) <doi:10.1109/24.257799> and Craiu and Reiser (2006) <doi:10.1111/j.1541-0420.2005.00498.x>. |
| Authors: | Alexander Towell [aut, cre] (ORCID: <https://orcid.org/0000-0001-6443-9897>) |
| Maintainer: | Alexander Towell <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 0.10.0 |
| Built: | 2026-05-24 09:29:37 UTC |
| Source: | https://github.com/queelius/maskedcauses |
exp_series_md_c1_c2_c3 model.Returns the assumptions made by this likelihood model.
## S3 method for class 'exp_series_md_c1_c2_c3' assumptions(model, ...)## S3 method for class 'exp_series_md_c1_c2_c3' assumptions(model, ...)
model |
the likelihood model object |
... |
additional arguments (ignored) |
character vector of assumptions
assumptions(exp_series_md_c1_c2_c3())assumptions(exp_series_md_c1_c2_c3())
wei_series_homogeneous_md_c1_c2_c3 model.Returns the assumptions made by this likelihood model.
## S3 method for class 'wei_series_homogeneous_md_c1_c2_c3' assumptions(model, ...)## S3 method for class 'wei_series_homogeneous_md_c1_c2_c3' assumptions(model, ...)
model |
the likelihood model object |
... |
additional arguments (ignored) |
character vector of assumptions
assumptions(wei_series_homogeneous_md_c1_c2_c3())assumptions(wei_series_homogeneous_md_c1_c2_c3())
wei_series_md_c1_c2_c3 model.Returns the assumptions made by this likelihood model.
## S3 method for class 'wei_series_md_c1_c2_c3' assumptions(model, ...)## S3 method for class 'wei_series_md_c1_c2_c3' assumptions(model, ...)
model |
the likelihood model object |
... |
additional arguments (ignored) |
character vector of assumptions
assumptions(wei_series_md_c1_c2_c3())assumptions(wei_series_md_c1_c2_c3())
Returns a closure computing P(K=j | theta) for all components j, marginalized over the system failure time T. By Theorem 5 of the foundational paper, this equals E_T[P(K=j | T, theta)].
cause_probability(model, ...) ## S3 method for class 'series_md' cause_probability(model, ...)cause_probability(model, ...) ## S3 method for class 'series_md' cause_probability(model, ...)
model |
a likelihood model object |
... |
additional arguments passed to the returned closure |
The default method uses Monte Carlo integration via rdata().
a function with signature function(par, ...) returning an m-vector
where element j gives P(K=j | theta)
cause_probability(series_md): Default for series_md models via Monte Carlo
integration (Theorem 5)
model <- exp_series_md_c1_c2_c3() cp <- cause_probability(model) cp(par = c(0.5, 0.3, 0.2))model <- exp_series_md_c1_c2_c3() cp <- cause_probability(model) cp(par = c(0.5, 0.3, 0.2))
Returns a closure computing the hazard function h_j(t; theta) for the j-th
component. The returned function takes the full parameter vector par and
extracts the relevant component parameters internally.
component_hazard(model, j, ...)component_hazard(model, j, ...)
model |
a likelihood model object |
j |
component index (integer from 1 to m) |
... |
additional arguments passed to the returned closure (e.g., covariates for proportional hazards extensions) |
a function with signature function(t, par, ...) computing h_j(t)
model <- exp_series_md_c1_c2_c3() h1 <- component_hazard(model, j = 1) h1(t = 1.0, par = c(0.5, 0.3, 0.2))model <- exp_series_md_c1_c2_c3() h1 <- component_hazard(model, j = 1) h1(t = 1.0, par = c(0.5, 0.3, 0.2))
Returns a closure computing P(K=j | T=t, theta) for all components j, conditional on a specific failure time t. By Theorem 6 of the foundational paper, this equals h_j(t; theta) / sum_l h_l(t; theta).
conditional_cause_probability(model, ...) ## S3 method for class 'series_md' conditional_cause_probability(model, ...)conditional_cause_probability(model, ...) ## S3 method for class 'series_md' conditional_cause_probability(model, ...)
model |
a likelihood model object |
... |
additional arguments passed to the returned closure |
a function with signature function(t, par, ...) returning an
n x m matrix where n = length(t) and column j gives P(K=j | T=t, theta)
conditional_cause_probability(series_md): Default for series_md models via
component hazard ratios (Theorem 6)
model <- exp_series_md_c1_c2_c3() ccp <- conditional_cause_probability(model) ccp(t = c(1, 2), par = c(0.5, 0.3, 0.2))model <- exp_series_md_c1_c2_c3() ccp <- conditional_cause_probability(model) ccp(t = c(1, 2), par = c(0.5, 0.3, 0.2))
Thin wrapper. Equivalent to
algebraic.dist::density(dist.structure::exp_series(rates))(t, log = log).
dexp_series(t, rates, log = FALSE)dexp_series(t, rates, log = FALSE)
t |
series system lifetime |
rates |
rate parameters for exponential component lifetimes |
log |
return the log of the pdf |
Density values for the exponential series distribution.
Preserved for backwards compatibility. New code: use
dist.structure::exp_series() then algebraic.dist::density().
rates <- c(0.5, 0.3, 0.2) dexp_series(1.0, rates) dexp_series(c(0.5, 1.0, 2.0), rates, log = TRUE)rates <- c(0.5, 0.3, 0.2) dexp_series(1.0, rates) dexp_series(c(0.5, 1.0, 2.0), rates, log = TRUE)
exp_series_md_c1_c2_c3.Likelihood model for exponential series systems with masked component cause of failure with candidate sets that satisfy conditions C1, C2, and C3, described below.
exp_series_md_c1_c2_c3( rates = NULL, lifetime = "t", lifetime_upper = "t_upper", omega = "omega", candset = "x" )exp_series_md_c1_c2_c3( rates = NULL, lifetime = "t", lifetime_upper = "t_upper", omega = "omega", candset = "x" )
rates |
rate parameters for exponential component lifetimes (optional, used as initial values for MLE if provided) |
lifetime |
column name for system lifetime, defaults to |
lifetime_upper |
column name for interval upper bound, defaults to
|
omega |
column name for observation type, defaults to |
candset |
column prefix for candidate set indicators, defaults to |
This model satisfies the concept of a likelihood_model in the
likelihood.model package by providing the following methods:
(1) loglik.exp_series_md_c1_c2_c3
(2) score.exp_series_md_c1_c2_c3
(3) hess_loglik.exp_series_md_c1_c2_c3
These are useful for doing maximum likelihood estimation, hypothesis testing (e.g., likelihood ratio test), estimation of asymptotic sampling distribution given data from the DGP according to the specified model, etc.
It is designed to work well with the likelihood_model R package. In
particular, it is intended to be used with the likelihood_contr_model
object, which is a likelihood_model object that allows likelihood
contributions to be added for whatever data model you have in mind.
In this likelihood model, masked component data approximately satisfies the following conditions:
C1: Pr{K[i] in C[i]) = 1
C2: Pr{C[i]=c[i] | K[i]=j, T[i]=t[i]) = Pr(C[i]=c[i] | K[i]=j', T[i]=t[i])
for any j, j' in c[i].
C3: masking probabilities are independent of theta
As a special case, this model also includes exact component cause of failure data where the candidate set is a singleton.
likelihood model object with class
c("exp_series_md_c1_c2_c3", "series_md", "likelihood_model")
model <- exp_series_md_c1_c2_c3() # Generate data and evaluate log-likelihood set.seed(123) gen <- rdata(model) df <- gen(theta = c(0.5, 0.3, 0.2), n = 50, tau = 10, p = 0.3) ll <- loglik(model) ll(df, par = c(0.5, 0.3, 0.2))model <- exp_series_md_c1_c2_c3() # Generate data and evaluate log-likelihood set.seed(123) gen <- rdata(model) df <- gen(theta = c(0.5, 0.3, 0.2), n = 50, tau = 10, p = 0.3) ll <- loglik(model) ll(df, par = c(0.5, 0.3, 0.2))
Thin wrapper. The series of m independent exponentials has constant
system hazard sum(rates). Equivalent to
algebraic.dist::hazard(dist.structure::exp_series(rates))(t).
hazard_exp_series(t, rates, log.p = FALSE)hazard_exp_series(t, rates, log.p = FALSE)
t |
Vector of series system lifetimes. |
rates |
Vector of rate parameters for exponential component lifetimes. |
log.p |
Logical; if TRUE, return the log of the hazard function. |
The hazard function evaluated at the specified lifetimes.
Preserved for backwards compatibility. New code: use
dist.structure::exp_series() then algebraic.dist::hazard().
rates <- c(0.5, 0.3, 0.2) hazard_exp_series(1.0, rates) hazard_exp_series(c(0.5, 1.0), rates, log.p = TRUE)rates <- c(0.5, 0.3, 0.2) hazard_exp_series(1.0, rates) hazard_exp_series(c(0.5, 1.0), rates, log.p = TRUE)
exp_series_md_c1_c2_c3 model.Returns the Hessian (second derivative matrix) of the log-likelihood for an
exponential series system with respect to parameter theta for masked data
with candidate sets that satisfy conditions C1, C2, and C3.
## S3 method for class 'exp_series_md_c1_c2_c3' hess_loglik(model, ...)## S3 method for class 'exp_series_md_c1_c2_c3' hess_loglik(model, ...)
model |
the likelihood model object |
... |
additional arguments (passed to returned function) |
All four observation types have closed-form Hessian contributions.
hessian function that takes the following arguments:
df: masked data frame
par: rate parameters of exponential component lifetime distributions
lifetime: system lifetime column name (default from model)
lifetime_upper: interval upper bound column name (default from model)
omega: observation type column name (default from model)
candset: prefix of Boolean matrix encoding candidate sets
model <- exp_series_md_c1_c2_c3() set.seed(1) df <- rdata(model)(theta = c(0.5, 0.3, 0.2), n = 30, tau = 10, p = 0.3) H <- hess_loglik(model) H(df, par = c(0.5, 0.3, 0.2))model <- exp_series_md_c1_c2_c3() set.seed(1) df <- rdata(model)(theta = c(0.5, 0.3, 0.2), n = 30, tau = 10, p = 0.3) H <- hess_loglik(model) H(df, par = c(0.5, 0.3, 0.2))
wei_series_homogeneous_md_c1_c2_c3.Returns the Hessian (second derivative matrix) of the log-likelihood for a Weibull series system with homogeneous shape. Computed numerically via the Jacobian of the score.
## S3 method for class 'wei_series_homogeneous_md_c1_c2_c3' hess_loglik(model, ...)## S3 method for class 'wei_series_homogeneous_md_c1_c2_c3' hess_loglik(model, ...)
model |
the likelihood model object |
... |
additional arguments (passed to returned function) |
hessian function that takes the following arguments:
df: masked data frame
par: parameter vector (shape, scale1, scale2, ...)
lifetime: system lifetime column name (default from model)
lifetime_upper: interval upper bound column name (default from model)
omega: observation type column name (default from model)
candset: prefix of Boolean matrix encoding candidate sets
model <- wei_series_homogeneous_md_c1_c2_c3() set.seed(1) theta <- c(1.2, 1000, 900, 850) df <- rdata(model)(theta = theta, n = 30, tau = 500, p = 0.3) H <- hess_loglik(model) H(df, par = theta)model <- wei_series_homogeneous_md_c1_c2_c3() set.seed(1) theta <- c(1.2, 1000, 900, 850) df <- rdata(model)(theta = theta, n = 30, tau = 500, p = 0.3) H <- hess_loglik(model) H(df, par = theta)
wei_series_md_c1_c2_c3 model.Returns the Hessian (second derivative matrix) of the log-likelihood for a Weibull series system. Computed numerically via the Jacobian of the score.
## S3 method for class 'wei_series_md_c1_c2_c3' hess_loglik(model, ...)## S3 method for class 'wei_series_md_c1_c2_c3' hess_loglik(model, ...)
model |
the likelihood model object |
... |
additional arguments (passed to returned function) |
hessian function that takes the following arguments:
df: masked data frame
par: parameter vector (shape1, scale1, shape2, scale2, ...)
lifetime: system lifetime column name (default from model)
lifetime_upper: interval upper bound column name (default from model)
omega: observation type column name (default from model)
candset: prefix of Boolean matrix encoding candidate sets
model <- wei_series_md_c1_c2_c3() set.seed(1) theta <- c(1.2, 1000, 0.8, 900) df <- rdata(model)(theta = theta, n = 30, tau = 500, p = 0.3) H <- hess_loglik(model) H(df, par = theta)model <- wei_series_md_c1_c2_c3() set.seed(1) theta <- c(1.2, 1000, 0.8, 900) df <- rdata(model)(theta = theta, n = 30, tau = 500, p = 0.3) H <- hess_loglik(model) H(df, par = theta)
Given a hazard rate function , returns a function computing
via numerical integration.
integrate_hazard(haz)integrate_hazard(haz)
haz |
hazard rate function |
A function that computes the cumulative hazard at time t
# Exponential hazard h(t) = lambda haz <- function(t, ...) rep(0.5, length(t)) H <- integrate_hazard(haz) H(2) # Should be 1.0 (0.5 * 2)# Exponential hazard h(t) = lambda haz <- function(t, ...) rep(0.5, length(t)) H <- integrate_hazard(haz) H(2) # Should be 1.0 (0.5 * 2)
exp_series_md_c1_c2_c3 model.Returns a log-likelihood function for an exponential series system with respect to rate parameters for masked data with candidate sets that satisfy conditions C1, C2, and C3.
## S3 method for class 'exp_series_md_c1_c2_c3' loglik(model, ...)## S3 method for class 'exp_series_md_c1_c2_c3' loglik(model, ...)
model |
the likelihood model object |
... |
additional arguments (passed to returned function) |
Supports four observation types: exact failures, right-censored, left-censored, and interval-censored. All have closed-form likelihood contributions for the exponential model.
log-likelihood function that takes the following arguments:
df: masked data frame
par: rate parameters of exponential component lifetime distributions
lifetime: system lifetime column name (default from model)
lifetime_upper: interval upper bound column name (default from model)
omega: observation type column name (default from model)
candset: prefix of Boolean matrix encoding candidate sets
model <- exp_series_md_c1_c2_c3() set.seed(1) df <- rdata(model)(theta = c(0.5, 0.3, 0.2), n = 30, tau = 10, p = 0.3) ll <- loglik(model) ll(df, par = c(0.5, 0.3, 0.2))model <- exp_series_md_c1_c2_c3() set.seed(1) df <- rdata(model)(theta = c(0.5, 0.3, 0.2), n = 30, tau = 10, p = 0.3) ll <- loglik(model) ll(df, par = c(0.5, 0.3, 0.2))
wei_series_homogeneous_md_c1_c2_c3 model.Returns a log-likelihood function for a Weibull series system with homogeneous shape parameter. The parameter vector is (k, scale_1, ..., scale_m) for masked data with candidate sets that satisfy conditions C1, C2, and C3.
## S3 method for class 'wei_series_homogeneous_md_c1_c2_c3' loglik(model, ...)## S3 method for class 'wei_series_homogeneous_md_c1_c2_c3' loglik(model, ...)
model |
the likelihood model object |
... |
additional arguments (passed to returned function) |
Supports four observation types. Left-censored and interval-censored have closed-form likelihood contributions because all shapes are equal.
log-likelihood function that takes the following arguments:
df: masked data frame
par: parameter vector (shape, scale1, scale2, ...)
lifetime: system lifetime column name (default from model)
lifetime_upper: interval upper bound column name (default from model)
omega: observation type column name (default from model)
candset: prefix of Boolean matrix encoding candidate sets
model <- wei_series_homogeneous_md_c1_c2_c3() set.seed(1) # theta: (shape, scale1, scale2, scale3) theta <- c(1.2, 1000, 900, 850) df <- rdata(model)(theta = theta, n = 30, tau = 500, p = 0.3) ll <- loglik(model) ll(df, par = theta)model <- wei_series_homogeneous_md_c1_c2_c3() set.seed(1) # theta: (shape, scale1, scale2, scale3) theta <- c(1.2, 1000, 900, 850) df <- rdata(model)(theta = theta, n = 30, tau = 500, p = 0.3) ll <- loglik(model) ll(df, par = theta)
wei_series_md_c1_c2_c3 model.Returns a log-likelihood function for a Weibull series system with respect to parameter vector (shape_1, scale_1, ..., shape_m, scale_m) for masked data with candidate sets that satisfy conditions C1, C2, and C3.
## S3 method for class 'wei_series_md_c1_c2_c3' loglik(model, ...)## S3 method for class 'wei_series_md_c1_c2_c3' loglik(model, ...)
model |
the likelihood model object |
... |
additional arguments (passed to returned function) |
Supports four observation types. Left-censored and interval-censored observations require numerical integration (via stats::integrate) because heterogeneous shapes prevent a closed-form solution.
log-likelihood function that takes the following arguments:
df: masked data frame
par: parameter vector (shape1, scale1, shape2, scale2, ...)
lifetime: system lifetime column name (default from model)
lifetime_upper: interval upper bound column name (default from model)
omega: observation type column name (default from model)
candset: prefix of Boolean matrix encoding candidate sets
model <- wei_series_md_c1_c2_c3() set.seed(1) theta <- c(1.2, 1000, 0.8, 900) df <- rdata(model)(theta = theta, n = 30, tau = 500, p = 0.3) ll <- loglik(model) ll(df, par = theta)model <- wei_series_md_c1_c2_c3() set.seed(1) theta <- c(1.2, 1000, 0.8, 900) df <- rdata(model)(theta = theta, n = 30, tau = 500, p = 0.3) ll <- loglik(model) ll(df, par = theta)
Bernoulli candidate set model is a particular type of uninformed model.
Note that we do not generate candidate sets with this function. See
md_cand_sampler for that.
md_bernoulli_cand_c1_c2_c3( df, p, prob = "q", comp = "t", right_censoring_indicator = "delta" )md_bernoulli_cand_c1_c2_c3( df, p, prob = "q", comp = "t", right_censoring_indicator = "delta" )
df |
masked data. |
p |
a vector of probabilities ( |
prob |
column prefix for component probabilities, defaults to
|
comp |
column prefix for component lifetimes, defaults to |
right_censoring_indicator |
right-censoring indicator column name.
if TRUE, then the system lifetime is right-censored, otherwise it is
observed. If NULL, then no right-censoring is assumed. Defaults to
|
This model satisfies conditions C1, C2, and C3.
The failed component will be in the corresponding candidate set with
probability 1, and the remaining components will be in the candidate set
with probability p (the same probability for each component). p
may be different for each system, but it is assumed to be the same for
each component within a system, so p can be a vector such that the
length of p is the number of systems in the data set (with recycling
if necessary).
Data frame with added candidate set probability columns
(e.g., q1, q2, ..., qm).
# Generate component lifetimes and system data mat <- matrix(rexp(9, rate = rep(c(0.5, 0.3, 0.2), each = 3)), nrow = 3, ncol = 3) df <- md_encode_matrix(mat, "t") df$t <- apply(mat, 1, min) df$delta <- TRUE md_bernoulli_cand_c1_c2_c3(df, p = 0.5)# Generate component lifetimes and system data mat <- matrix(rexp(9, rate = rep(c(0.5, 0.3, 0.2), each = 3)), nrow = 3, ncol = 3) df <- md_encode_matrix(mat, "t") df$t <- apply(mat, 1, min) df$delta <- TRUE md_bernoulli_cand_c1_c2_c3(df, p = 0.5)
Replaces Boolean matrix columns (e.g., x1, x2, x3) with a single
character column showing set notation like {1, 3}.
md_boolean_matrix_to_charsets(df, setvar = "x", cname = NULL, drop_set = FALSE)md_boolean_matrix_to_charsets(df, setvar = "x", cname = NULL, drop_set = FALSE)
df |
data frame containing Boolean matrix columns |
setvar |
column prefix for the Boolean matrix (default |
cname |
name for the new character column (default: same as |
drop_set |
if TRUE, remove the original Boolean columns (default FALSE) |
data frame with character set column added
df <- data.frame(x1 = c(TRUE, FALSE, TRUE), x2 = c(TRUE, TRUE, FALSE), x3 = c(FALSE, TRUE, TRUE)) md_boolean_matrix_to_charsets(df)df <- data.frame(x1 = c(TRUE, FALSE, TRUE), x2 = c(TRUE, TRUE, FALSE), x3 = c(FALSE, TRUE, TRUE)) md_boolean_matrix_to_charsets(df)
Candidate set generator. Requires columns for component probabilities
e.g., q1,...,qm where qj is the probability that the jth component
will be in the corresponding candidate set generated for that observation
in the md table.
md_cand_sampler(df, prob = "q", candset = "x")md_cand_sampler(df, prob = "q", candset = "x")
df |
(masked) data frame |
prob |
column prefix for component probabilities, defaults to
|
candset |
column prefix for candidate sets (as Boolean matrix),
defaults to |
Data frame with added Boolean candidate set columns
(e.g., x1, x2, ..., xm).
# Generate component lifetimes set.seed(42) mat <- matrix(rexp(9, rate = rep(c(0.5, 0.3, 0.2), each = 3)), nrow = 3, ncol = 3) df <- md_encode_matrix(mat, "t") df$t <- apply(mat, 1, min) df$delta <- TRUE df <- md_bernoulli_cand_c1_c2_c3(df, p = 0.5) md_cand_sampler(df)# Generate component lifetimes set.seed(42) mat <- matrix(rexp(9, rate = rep(c(0.5, 0.3, 0.2), each = 3)), nrow = 3, ncol = 3) df <- md_encode_matrix(mat, "t") df$t <- apply(mat, 1, min) df$delta <- TRUE df <- md_bernoulli_cand_c1_c2_c3(df, p = 0.5) md_cand_sampler(df)
Converts a matrix to a data frame with columns named var1, var2, ....
md_encode_matrix(mat, var)md_encode_matrix(mat, var)
mat |
matrix to encode |
var |
character prefix for the column names |
a data frame with named columns
mat <- matrix(1:6, nrow = 2, ncol = 3) md_encode_matrix(mat, "x")mat <- matrix(1:6, nrow = 2, ncol = 3) md_encode_matrix(mat, "x")
Generates right-censored system failure times and right-censoring indicators for a series system with the given data frame of component lifetimes.
md_series_lifetime_right_censoring( df, tau = Inf, comp = "t", lifetime = "t", right_censoring_indicator = "delta" )md_series_lifetime_right_censoring( df, tau = Inf, comp = "t", lifetime = "t", right_censoring_indicator = "delta" )
df |
a data frame with the indicated component lifetimes |
tau |
vector of right-censoring times, defaults to
|
comp |
component lifetime prefix variable name, defaults to |
lifetime |
system lifetime variable name, defaults to |
right_censoring_indicator |
right-censoring indicator variable, defaults
to |
(masked) data frame with masked data as described in the paper
mat <- matrix(rexp(9, rate = 0.5), nrow = 3, ncol = 3) df <- md_encode_matrix(mat, "t") md_series_lifetime_right_censoring(df, tau = 5)mat <- matrix(rexp(9, rate = 0.5), nrow = 3, ncol = 3) df <- md_encode_matrix(mat, "t") md_series_lifetime_right_censoring(df, tau = 5)
Computes the expected value of a series system with exponentially
distributed component lifetimes. For a series system with component
rates , the system lifetime is
exponential with rate , so .
## S3 method for class 'exp_series' mean(x, ...)## S3 method for class 'exp_series' mean(x, ...)
x |
An |
... |
Additional arguments (ignored, for S3 generic compatibility). |
This method is polymorphic: it accepts both the legacy
maskedcauses shape (a numeric vector of rates with class attribute
"exp_series") and the current dist.structure::exp_series()
shape (a list with a $total_rate field). This prevents a breaking
S3 dispatch collision when both packages are attached.
The mean of the exponential series distribution (1/sum of rates).
Preserved for backwards compatibility. New code: construct
with dist.structure::exp_series(rates) and call mean() on it.
# Legacy shape rates <- structure(c(0.5, 0.3, 0.2), class = "exp_series") mean(rates) # Modern shape (dist.structure list) works via the same method mean(dist.structure::exp_series(c(0.5, 0.3, 0.2)))# Legacy shape rates <- structure(c(0.5, 0.3, 0.2), class = "exp_series") mean(rates) # Modern shape (dist.structure list) works via the same method mean(dist.structure::exp_series(c(0.5, 0.3, 0.2)))
Returns the number of components m in the series system. If the model was constructed without parameter values, returns NULL.
ncomponents(model, ...)ncomponents(model, ...)
model |
a likelihood model object |
... |
additional arguments (currently ignored) |
integer number of components, or NULL if not determinable
model <- exp_series_md_c1_c2_c3(rates = c(0.5, 0.3, 0.2)) ncomponents(model)model <- exp_series_md_c1_c2_c3(rates = c(0.5, 0.3, 0.2)) ncomponents(model)
Creates an observation functor for a single-inspection design. If the system
has already failed by inspection time tau, we know it failed before
tau but not exactly when (left-censored). If it is still running, we
know it survived past tau (right-censored).
observe_left_censor(tau)observe_left_censor(tau)
tau |
inspection time (positive numeric) |
A function with signature function(t_true) returning a list
with components:
inspection time tau
"left" if failed before tau, "right" otherwise
NA (not used for this scheme)
obs <- observe_left_censor(tau = 100) obs(50) # left-censored: list(t = 100, omega = "left", t_upper = NA) obs(150) # right-censored: list(t = 100, omega = "right", t_upper = NA)obs <- observe_left_censor(tau = 100) obs(50) # left-censored: list(t = 100, omega = "left", t_upper = NA) obs(150) # right-censored: list(t = 100, omega = "right", t_upper = NA)
Creates an observation functor that randomly selects from a set of observation schemes for each observation. This models heterogeneous monitoring environments where different units are observed differently.
observe_mixture(..., weights = NULL)observe_mixture(..., weights = NULL)
... |
observation functors (created by |
weights |
mixing probabilities (numeric vector). If NULL, uniform weights are used. Weights are normalized to sum to 1. |
A function with signature function(t_true) returning a list
from one of the constituent schemes, selected randomly according to
weights.
obs <- observe_mixture( observe_right_censor(tau = 100), observe_left_censor(tau = 50), weights = c(0.7, 0.3) ) set.seed(42) obs(30) # randomly selects one of the two schemesobs <- observe_mixture( observe_right_censor(tau = 100), observe_left_censor(tau = 50), weights = c(0.7, 0.3) ) set.seed(42) obs(30) # randomly selects one of the two schemes
Creates an observation functor for periodic inspections at intervals of
delta. Failures are bracketed between the last inspection before
failure and the first inspection after failure (interval-censored). Systems
surviving past tau are right-censored.
observe_periodic(delta, tau = Inf)observe_periodic(delta, tau = Inf)
delta |
inspection interval (positive numeric) |
tau |
study end time (positive numeric or Inf for no right-censoring) |
A function with signature function(t_true) returning a list
with components:
lower bound of interval (or tau if right-censored)
"interval" or "right"
upper bound of interval (NA if right-censored)
obs <- observe_periodic(delta = 10, tau = 100) obs(25) # interval: list(t = 20, omega = "interval", t_upper = 30) obs(150) # right-censored: list(t = 100, omega = "right", t_upper = NA)obs <- observe_periodic(delta = 10, tau = 100) obs(25) # interval: list(t = 20, omega = "interval", t_upper = 30) obs(150) # right-censored: list(t = 100, omega = "right", t_upper = NA)
Creates an observation functor that applies right-censoring at time tau.
Systems that fail before tau are observed exactly; systems surviving
past tau are right-censored.
observe_right_censor(tau)observe_right_censor(tau)
tau |
censoring time (positive numeric) |
A function with signature function(t_true) returning a list
with components:
observed time
"exact" or "right"
NA (not used for this scheme)
obs <- observe_right_censor(tau = 100) obs(50) # exact: list(t = 50, omega = "exact", t_upper = NA) obs(150) # right-censored: list(t = 100, omega = "right", t_upper = NA)obs <- observe_right_censor(tau = 100) obs(50) # exact: list(t = 50, omega = "exact", t_upper = NA) obs(150) # right-censored: list(t = 100, omega = "right", t_upper = NA)
Thin wrapper. Equivalent to
algebraic.dist::cdf(dist.structure::exp_series(rates))(t).
pexp_series(t, rates, lower.tail = TRUE, log.p = FALSE)pexp_series(t, rates, lower.tail = TRUE, log.p = FALSE)
t |
Vector of series system lifetimes. |
rates |
Vector of rate parameters for exponential component lifetimes. |
lower.tail |
Logical; if TRUE (default), probabilities are P(X <= x), otherwise, P(X > x). |
log.p |
Logical; if TRUE, return the log of the cdf. |
The cumulative probabilities evaluated at the specified lifetimes.
Preserved for backwards compatibility. New code: use
dist.structure::exp_series() then algebraic.dist::cdf().
rates <- c(0.5, 0.3, 0.2) pexp_series(1.0, rates) pexp_series(c(0.5, 1.0, 2.0), rates)rates <- c(0.5, 0.3, 0.2) pexp_series(1.0, rates) pexp_series(c(0.5, 1.0, 2.0), rates)
Finds the time t such that S(t) = p using root finding. The survival function S(t) is assumed to be monotonically decreasing from S(0) = 1 to S(inf) = 0.
qcomp( p, surv, theta, t_lower = .Machine$double.eps, t_upper = .Machine$double.xmax^0.5, ... )qcomp( p, surv, theta, t_lower = .Machine$double.eps, t_upper = .Machine$double.xmax^0.5, ... )
p |
probability (quantile level), must be in (0, 1) |
surv |
survival function S(t, theta, ...) |
theta |
parameter vector passed to surv |
t_lower |
lower bound for search interval |
t_upper |
upper bound for search interval (sqrt to avoid overflow) |
... |
additional arguments passed to surv |
time t such that S(t) = p
# Exponential survival function surv_exp <- function(t, theta) exp(-theta * t) # Median lifetime (50th percentile) for rate = 2 qcomp(0.5, surv = surv_exp, theta = 2.0)# Exponential survival function surv_exp <- function(t, theta) exp(-theta * t) # Median lifetime (50th percentile) for rate = 2 qcomp(0.5, surv = surv_exp, theta = 2.0)
Thin wrapper around stats::qexp() with rate = sum(rates). Equivalent
to algebraic.dist::inv_cdf(dist.structure::exp_series(rates))(p).
qexp_series(p, rates, lower.tail = TRUE, log.p = FALSE)qexp_series(p, rates, lower.tail = TRUE, log.p = FALSE)
p |
Vector of quantiles. |
rates |
Vector of rate parameters for exponential component lifetimes. |
lower.tail |
Logical, if TRUE (default), probabilities are P(X <= x), otherwise, P(X > x). |
log.p |
Logical, if TRUE, vector of probabilities |
Quantiles corresponding to the given probabilities p.
Preserved for backwards compatibility. New code: use
dist.structure::exp_series() then algebraic.dist::inv_cdf().
rates <- c(0.5, 0.3, 0.2) qexp_series(0.5, rates) qexp_series(c(0.25, 0.5, 0.75), rates)rates <- c(0.5, 0.3, 0.2) qexp_series(0.5, rates) qexp_series(c(0.25, 0.5, 0.75), rates)
Generates random samples using inverse transform sampling.
rcomp(n, surv, theta)rcomp(n, surv, theta)
n |
number of samples to generate |
surv |
survival function S(t, theta) |
theta |
parameter vector passed to surv |
vector of n random samples
# Exponential survival function surv_exp <- function(t, theta) exp(-theta * t) # Generate 10 random samples with rate = 2 set.seed(123) rcomp(10, surv = surv_exp, theta = 2.0)# Exponential survival function surv_exp <- function(t, theta) exp(-theta * t) # Generate 10 random samples with rate = 2 set.seed(123) rcomp(10, surv = surv_exp, theta = 2.0)
exp_series_md_c1_c2_c3 model.Returns a function that generates random masked data from the exponential series system DGP at a given parameter value.
## S3 method for class 'exp_series_md_c1_c2_c3' rdata(model, ...)## S3 method for class 'exp_series_md_c1_c2_c3' rdata(model, ...)
model |
the likelihood model object |
... |
additional arguments (passed to returned function) |
function that takes (theta, n, tau, p, observe, ...) and returns a data frame with columns: t, omega, x1, x2, ..., xm
model <- exp_series_md_c1_c2_c3() gen <- rdata(model) set.seed(42) df <- gen(theta = c(0.5, 0.3, 0.2), n = 10, tau = 5, p = 0.5) head(df)model <- exp_series_md_c1_c2_c3() gen <- rdata(model) set.seed(42) df <- gen(theta = c(0.5, 0.3, 0.2), n = 10, tau = 5, p = 0.5) head(df)
wei_series_homogeneous_md_c1_c2_c3 model.Returns a function that generates random masked data from the homogeneous Weibull series system DGP at a given parameter value.
## S3 method for class 'wei_series_homogeneous_md_c1_c2_c3' rdata(model, ...)## S3 method for class 'wei_series_homogeneous_md_c1_c2_c3' rdata(model, ...)
model |
the likelihood model object |
... |
additional arguments (passed to returned function) |
function that takes (theta, n, tau, p, observe, ...) and returns a data frame with columns: t, omega, x1, x2, ..., xm
model <- wei_series_homogeneous_md_c1_c2_c3() gen <- rdata(model) set.seed(42) df <- gen(theta = c(1.2, 1000, 900, 850), n = 10, tau = 500, p = 0.5) head(df)model <- wei_series_homogeneous_md_c1_c2_c3() gen <- rdata(model) set.seed(42) df <- gen(theta = c(1.2, 1000, 900, 850), n = 10, tau = 500, p = 0.5) head(df)
wei_series_md_c1_c2_c3 model.Returns a function that generates random masked data from the Weibull series system DGP at a given parameter value.
## S3 method for class 'wei_series_md_c1_c2_c3' rdata(model, ...)## S3 method for class 'wei_series_md_c1_c2_c3' rdata(model, ...)
model |
the likelihood model object |
... |
additional arguments (passed to returned function) |
function that takes (theta, n, tau, p, observe, ...) and returns a data frame with columns: t, omega, x1, x2, ..., xm
model <- wei_series_md_c1_c2_c3() gen <- rdata(model) set.seed(42) # theta: (shape1, scale1, shape2, scale2) df <- gen(theta = c(1.2, 1000, 0.8, 900), n = 10, tau = 500, p = 0.5) head(df)model <- wei_series_md_c1_c2_c3() gen <- rdata(model) set.seed(42) # theta: (shape1, scale1, shape2, scale2) df <- gen(theta = c(1.2, 1000, 0.8, 900), n = 10, tau = 500, p = 0.5) head(df)
Generates random variates from an exponential series distribution.
The default path (keep_latent = FALSE) is equivalent to
algebraic.dist::sampler(dist.structure::exp_series(rates))(n). The
keep_latent = TRUE path additionally returns the latent component
lifetimes (sampled independently from each component).
rexp_series(n, rates, keep_latent = FALSE)rexp_series(n, rates, keep_latent = FALSE)
n |
Integer, number of observations to generate. |
rates |
Vector of rate parameters for exponential component lifetimes. |
keep_latent |
Logical; if TRUE, returns a matrix with system lifetimes in the first column and individual component lifetimes in subsequent columns. If FALSE (default), returns only system lifetimes. |
If keep_latent = FALSE, a vector of random variates from the
exponential series distribution. If keep_latent = TRUE, a matrix
with system lifetime in the first column and component lifetimes
in columns 2 through m+1.
Preserved for backwards compatibility. New code: use
dist.structure::exp_series() with algebraic.dist::sampler().
set.seed(123) rexp_series(5, rates = c(0.5, 0.3, 0.2)) rexp_series(3, rates = c(0.5, 0.3, 0.2), keep_latent = TRUE)set.seed(123) rexp_series(5, rates = c(0.5, 0.3, 0.2)) rexp_series(3, rates = c(0.5, 0.3, 0.2), keep_latent = TRUE)
exp_series_md_c1_c2_c3 model.Returns a score (gradient) function for an exponential series system with
respect to parameter theta for masked component failure with candidate
sets that satisfy conditions C1, C2, and C3.
## S3 method for class 'exp_series_md_c1_c2_c3' score(model, ...)## S3 method for class 'exp_series_md_c1_c2_c3' score(model, ...)
model |
the likelihood model object |
... |
additional arguments (passed to returned function) |
All four observation types (exact, right, left, interval) have closed-form score contributions.
score function that takes the following arguments:
df: masked data frame
par: rate parameters of exponential component lifetime distributions
lifetime: system lifetime column name (default from model)
lifetime_upper: interval upper bound column name (default from model)
omega: observation type column name (default from model)
candset: prefix of Boolean matrix encoding candidate sets
model <- exp_series_md_c1_c2_c3() set.seed(1) df <- rdata(model)(theta = c(0.5, 0.3, 0.2), n = 30, tau = 10, p = 0.3) s <- score(model) s(df, par = c(0.5, 0.3, 0.2))model <- exp_series_md_c1_c2_c3() set.seed(1) df <- rdata(model)(theta = c(0.5, 0.3, 0.2), n = 30, tau = 10, p = 0.3) s <- score(model) s(df, par = c(0.5, 0.3, 0.2))
wei_series_homogeneous_md_c1_c2_c3 model.Returns a score (gradient) function for a Weibull series system with homogeneous shape parameter. The parameter vector is (k, scale_1, ..., scale_m) for masked data with candidate sets that satisfy conditions C1, C2, and C3.
## S3 method for class 'wei_series_homogeneous_md_c1_c2_c3' score(model, ...)## S3 method for class 'wei_series_homogeneous_md_c1_c2_c3' score(model, ...)
model |
the likelihood model object |
... |
additional arguments (passed to returned function) |
Uses a hybrid approach: analytical score for exact and right-censored observations, numerical gradient (via numDeriv) for left-censored and interval-censored observations.
score function that takes the following arguments:
df: masked data frame
par: parameter vector (shape, scale1, scale2, ...)
lifetime: system lifetime column name (default from model)
lifetime_upper: interval upper bound column name (default from model)
omega: observation type column name (default from model)
candset: prefix of Boolean matrix encoding candidate sets
model <- wei_series_homogeneous_md_c1_c2_c3() set.seed(1) theta <- c(1.2, 1000, 900, 850) df <- rdata(model)(theta = theta, n = 30, tau = 500, p = 0.3) s <- score(model) s(df, par = theta)model <- wei_series_homogeneous_md_c1_c2_c3() set.seed(1) theta <- c(1.2, 1000, 900, 850) df <- rdata(model)(theta = theta, n = 30, tau = 500, p = 0.3) s <- score(model) s(df, par = theta)
wei_series_md_c1_c2_c3 model.Returns a score (gradient) function for a Weibull series system with respect to parameter vector (shape_1, scale_1, ..., shape_m, scale_m) for masked data with candidate sets that satisfy conditions C1, C2, and C3.
## S3 method for class 'wei_series_md_c1_c2_c3' score(model, ...)## S3 method for class 'wei_series_md_c1_c2_c3' score(model, ...)
model |
the likelihood model object |
... |
additional arguments (passed to returned function) |
Uses a hybrid approach: analytical score for exact and right-censored observations, numerical gradient (via numDeriv) for left-censored and interval-censored observations.
score function that takes the following arguments:
df: masked data frame
par: parameter vector (shape1, scale1, shape2, scale2, ...)
lifetime: system lifetime column name (default from model)
lifetime_upper: interval upper bound column name (default from model)
omega: observation type column name (default from model)
candset: prefix of Boolean matrix encoding candidate sets
model <- wei_series_md_c1_c2_c3() set.seed(1) theta <- c(1.2, 1000, 0.8, 900) df <- rdata(model)(theta = theta, n = 30, tau = 500, p = 0.3) s <- score(model) s(df, par = theta)model <- wei_series_md_c1_c2_c3() set.seed(1) theta <- c(1.2, 1000, 0.8, 900) df <- rdata(model)(theta = theta, n = 30, tau = 500, p = 0.3) s <- score(model) s(df, par = theta)
Thin wrapper. Equivalent to
algebraic.dist::surv(dist.structure::exp_series(rates))(t).
surv.exp_series(t, rates, log.p = FALSE)surv.exp_series(t, rates, log.p = FALSE)
t |
Vector of series system lifetimes. |
rates |
Vector of rate parameters for exponential component lifetimes. |
log.p |
Logical; if TRUE, return the log of the survival function. |
The survival function evaluated at the specified lifetimes.
Preserved for backwards compatibility. New code: use
dist.structure::exp_series() then algebraic.dist::surv().
rates <- c(0.5, 0.3, 0.2) surv.exp_series(1.0, rates) surv.exp_series(c(0.5, 1.0, 2.0), rates)rates <- c(0.5, 0.3, 0.2) surv.exp_series(1.0, rates) surv.exp_series(c(0.5, 1.0, 2.0), rates)
wei_series_homogeneous_md_c1_c2_c3.Likelihood model for Weibull series systems with homogeneous shape parameter and masked component cause of failure with candidate sets that satisfy conditions C1, C2, and C3.
wei_series_homogeneous_md_c1_c2_c3( shape = NULL, scales = NULL, lifetime = "t", lifetime_upper = "t_upper", omega = "omega", candset = "x" )wei_series_homogeneous_md_c1_c2_c3( shape = NULL, scales = NULL, lifetime = "t", lifetime_upper = "t_upper", omega = "omega", candset = "x" )
shape |
common shape parameter for all Weibull component lifetimes |
scales |
scale parameters for Weibull component lifetimes (optional) |
lifetime |
column name for system lifetime, defaults to |
lifetime_upper |
column name for interval upper bound, defaults to
|
omega |
column name for observation type, defaults to |
candset |
column prefix for candidate set indicators, defaults to |
This is a reduced model where all components share a common shape parameter k, while retaining individual scale parameters. The parameter vector is (k, scale_1, ..., scale_m), giving m+1 parameters instead of 2m.
A key property of this model is that the series system lifetime is itself
Weibull distributed with shape k and scale .
This model satisfies the concept of a likelihood_model in the
likelihood.model package by providing the following methods:
(1) loglik.wei_series_homogeneous_md_c1_c2_c3
(2) score.wei_series_homogeneous_md_c1_c2_c3
(3) hess_loglik.wei_series_homogeneous_md_c1_c2_c3
In this likelihood model, masked component data approximately satisfies:
C1: Pr{K[i] in C[i]} = 1
C2: Pr{C[i]=c[i] | K[i]=j, T[i]=t[i]} = Pr{C[i]=c[i] | K[i]=j', T[i]=t[i]}
for any j, j' in c[i].
C3: masking probabilities are independent of theta
likelihood model object with class
c("wei_series_homogeneous_md_c1_c2_c3", "series_md", "likelihood_model")
wei_series_md_c1_c2_c3() for the full model with heterogeneous shapes
# Create model and fit to data using generic dispatch model <- wei_series_homogeneous_md_c1_c2_c3() # solver <- fit(model) # theta: (shape, scale1, scale2, ...) # mle <- solver(data, par = c(1.2, 1000, 900, 850))# Create model and fit to data using generic dispatch model <- wei_series_homogeneous_md_c1_c2_c3() # solver <- fit(model) # theta: (shape, scale1, scale2, ...) # mle <- solver(data, par = c(1.2, 1000, 900, 850))
wei_series_md_c1_c2_c3.Likelihood model for Weibull series systems with masked component cause of failure with candidate sets that satisfy conditions C1, C2, and C3.
wei_series_md_c1_c2_c3( shapes = NULL, scales = NULL, lifetime = "t", lifetime_upper = "t_upper", omega = "omega", candset = "x" )wei_series_md_c1_c2_c3( shapes = NULL, scales = NULL, lifetime = "t", lifetime_upper = "t_upper", omega = "omega", candset = "x" )
shapes |
shape parameters for Weibull component lifetimes (optional) |
scales |
scale parameters for Weibull component lifetimes (optional) |
lifetime |
column name for system lifetime, defaults to |
lifetime_upper |
column name for interval upper bound, defaults to
|
omega |
column name for observation type, defaults to |
candset |
column prefix for candidate set indicators, defaults to |
This model satisfies the concept of a likelihood_model in the
likelihood.model package by providing the following methods:
(1) loglik.wei_series_md_c1_c2_c3
(2) score.wei_series_md_c1_c2_c3
(3) hess_loglik.wei_series_md_c1_c2_c3
The Weibull series system has 2m parameters: (shape_1, scale_1, ..., shape_m, scale_m).
In this likelihood model, masked component data approximately satisfies:
C1: Pr{K[i] in C[i]} = 1
C2: Pr{C[i]=c[i] | K[i]=j, T[i]=t[i]} = Pr{C[i]=c[i] | K[i]=j', T[i]=t[i]}
for any j, j' in c[i].
C3: masking probabilities are independent of theta
likelihood model object with class
c("wei_series_md_c1_c2_c3", "series_md", "likelihood_model")
# Create model and fit to data using generic dispatch model <- wei_series_md_c1_c2_c3() # solver <- fit(model) # theta: (shape1, scale1, shape2, scale2, ...) # mle <- solver(data, par = c(1, 1000, 1, 1000, 1, 1000))# Create model and fit to data using generic dispatch model <- wei_series_md_c1_c2_c3() # solver <- fit(model) # theta: (shape1, scale1, shape2, scale2, ...) # mle <- solver(data, par = c(1, 1000, 1, 1000, 1, 1000))
For a series system with Weibull components sharing shape k but with individual scales, the system lifetime is itself Weibull with shape k and this computed scale.
wei_series_system_scale(k, scales)wei_series_system_scale(k, scales)
k |
common shape parameter |
scales |
vector of component scale parameters |
system scale parameter
# 3-component system with common shape 1.2 wei_series_system_scale(k = 1.2, scales = c(1000, 900, 850))# 3-component system with common shape 1.2 wei_series_system_scale(k = 1.2, scales = c(1000, 900, 850))