Package {kofn}

Title:	Maximum Likelihood Estimation for k-Out-of-n System Data
Version:	0.4.0
Description:	Maximum likelihood estimation of component lifetime parameters from system-level observations of k-out-of-n systems. Supports exponential and Weibull component distributions under multiple observation schemes: Scheme 0 (system lifetime only), Scheme 1 (periodic inspection), and Scheme 2 (complete monitoring). Provides an EM algorithm for Weibull parallel systems and Fisher information comparison across schemes. The k-out-of-n framework unifies series (k=1) and parallel (k=m) systems as a censoring problem on component lifetimes. Conforms to the 'likelihood.model' generics and returns fitted objects compatible with 'algebraic.mle'. The data-generating process and topology infrastructure (system survival, density, signature, structure function, importance measures) are delegated to the 'dist.structure' package; 'kofn' focuses exclusively on inference for the k-out-of-n family.
License:	MIT + file LICENSE
URL:	https://github.com/queelius/kofn, https://queelius.github.io/kofn/
BugReports:	https://github.com/queelius/kofn/issues
Depends:	R (≥ 4.1.0)
Imports:	stats, numDeriv, likelihood.model, compositional.mle, generics, dist.structure, algebraic.dist, flexhaz (≥ 0.5.2)
Suggests:	algebraic.mle, maskedcauses, testthat (≥ 3.0.0), knitr, rmarkdown
Config/testthat/edition:	3
Encoding:	UTF-8
RoxygenNote:	7.3.3
VignetteBuilder:	knitr
NeedsCompilation:	no
Packaged:	2026-06-11 11:33:51 UTC; spinoza
Author:	Alexander Towell [aut, cre]
Maintainer:	Alexander Towell <lex@metafunctor.com>
Repository:	CRAN
Date/Publication:	2026-06-19 12:00:29 UTC

System CDF for exponential parallel systems

Description

Computes F_{sys}(t) = \prod_j (1 - e^{-\lambda_j t}) using the inclusion-exclusion expansion.

Usage

F_sys_exp(t, par)

Arguments

Value

Scalar CDF value P(T_{sys} \le t).

See Also

S_sys_exp() for the survival function.

Examples

F_sys_exp(1, c(1, 2))  # P(T_sys <= 1) for 2-component parallel

System survival function for exponential parallel systems

Description

Computes S_{sys}(t) = 1 - F_{sys}(t) for a parallel system with exponential components.

Usage

S_sys_exp(t, par)

Arguments

Value

Scalar survival probability P(T_{sys} > t).

See Also

F_sys_exp() for the CDF.

Examples

S_sys_exp(1, c(1, 2))  # P(T_sys > 1) for 2-component parallel

Assumptions for exponential k-out-of-n model

Description

Returns a character vector listing the assumptions made by the exponential k-out-of-n likelihood model.

Usage

## S3 method for class 'exp_kofn'
assumptions(model, ...)

Arguments

Value

Character vector of assumptions.

Examples

model <- kofn(k = 3, m = 3, component = dfr_exponential())
assumptions(model)

Assumptions for Weibull k-out-of-n system model

Description

Returns a list of assumptions made by the model.

Usage

## S3 method for class 'wei_kofn'
assumptions(model, ...)

Arguments

Value

A character vector of assumptions.

Examples

model <- kofn(k = 2, m = 2, component = dfr_weibull())
assumptions(model)

Compare Fisher information across observation schemes

Description

For a given parameter configuration, computes the observed Fisher information under Scheme 0 (system only), Scheme 1 (periodic inspection), and Scheme 2 (complete component data) via simulation.

Usage

compare_fisher_info(
  shapes = NULL,
  scales = NULL,
  rates = NULL,
  n = 200L,
  delta = 1,
  n_rep = 50L,
  component = dfr_exponential()
)

Arguments

Details

The determinant of the observed information matrix is used as a scalar summary. Efficiency ratios indicate relative information content: a ratio less than 1 means the denominator scheme carries more information.

For Scheme 2 (complete data), the Fisher information is computed analytically:

• Exponential: I_{jj} = n / \lambda_j^2 (diagonal).

• Weibull: Block-diagonal with per-component 2x2 Fisher information matrices using the standard Weibull FIM formulas.

For Schemes 0 and 1, the observed information is computed numerically via numDeriv::hessian of the negative log-likelihood evaluated at the true parameter values.

Value

A list with components:

scheme0_det

Numeric vector of Scheme 0 information determinants (length n_rep).

scheme1_det

Numeric vector of Scheme 1 information determinants.

scheme2_det

Numeric vector of Scheme 2 (complete data) information determinants.

median_det

Named numeric vector of median determinants across schemes.

efficiency_01

Ratio of median Scheme 0 to Scheme 1 determinant (< 1 means Scheme 1 is more informative).

efficiency_02

Ratio of median Scheme 0 to Scheme 2 determinant.

efficiency_12

Ratio of median Scheme 1 to Scheme 2 determinant.

n_valid

Named integer vector of valid (non-NA) replicate counts per scheme.

Examples

set.seed(1)
result <- compare_fisher_info(rates = c(1, 2), n = 50, delta = 1.0, n_rep = 5)
result$median_det
result$efficiency_01  # Scheme 0 vs Scheme 1

Compute default initial values for exponential rates

Description

Uses the method-of-moments estimator based on the harmonic number relationship for i.i.d. parallel exponentials: ET = H_m / lambda. Returns a spread of initial values to break permutation symmetry.

Usage

default_init_exp(df, m, lifetime, omega, lifetime_upper)

Arguments

Value

Numeric vector of length m with initial rate estimates.

MLE fitting for exponential k-out-of-n model

Description

Returns a closure ⁠function(df, par0, n_starts)⁠ that computes the maximum likelihood estimate of the component rate parameters.

Usage

## S3 method for class 'exp_kofn'
fit(object, ...)

Arguments

Details

Uses multi-start optimization with L-BFGS-B as primary solver and Nelder-Mead on the log-scale as fallback. Standard errors are computed from the observed Fisher information (negative Hessian) at the MLE.

Value

A function ⁠function(df, par0 = NULL, n_starts = 5L)⁠ returning a fisher_mle object (from likelihood.model).

Examples

model <- kofn(k = 2, m = 2, component = dfr_exponential())
set.seed(42)
df <- rdata(model)(c(1, 2), n = 50)
result <- fit(model)(df)
coef(result)

Fit Weibull k-out-of-n system model

Description

Returns a closure that fits the model to data. Two methods are available:

Usage

## S3 method for class 'wei_kofn'
fit(object, ...)

Arguments

Details

"em"

EM algorithm for parallel systems (k = m only). Treats the identity of the last-failing component as latent data. Uses truncated Weibull moments for the E-step and profile optimization over shape with closed-form scale for the M-step.

"mle"

Direct MLE via L-BFGS-B with Nelder-Mead fallback. Works for any k-out-of-n structure.

The returned solver accepts these additional arguments:

par0

Initial parameter vector (length 2*m, interleaved shape/scale). If NULL, method-of-moments initialization is used.

n_starts

Number of multi-start attempts (default: 5).

tol

Convergence tolerance for EM (default: 1e-8).

maxiter

Maximum EM iterations (default: 2000).

shape_bounds

Bounds for shape optimization in EM (default: c(0.05, 15)).

verbose

Print iteration progress (default: FALSE).

Value

A function function(df, par0, ...) that returns a fisher_mle object (from the likelihood.model package).

Examples

model <- kofn(k = 2, m = 2, component = dfr_weibull())
set.seed(42)
df <- rdata(model)(c(1.5, 2.0, 2.0, 3.0), n = 50)
result <- fit(model)(df)
coef(result)

Fit k-out-of-n system MLE under Scheme 1 (periodic inspection)

Description

Returns a closure that fits the model to Scheme 1 data using multi-start optimization with L-BFGS-B and Nelder-Mead fallback.

Usage

fit_scheme1(model, ...)

Arguments

Details

The solver uses L-BFGS-B as the primary optimization method with positivity constraints, falling back to Nelder-Mead on the log-parameter scale if L-BFGS-B fails to converge.

Standard errors are computed from the numerical Hessian at the MLE.

Value

A function function(df, par0 = NULL, n_starts = 5L) that returns a fisher_mle object (from the likelihood.model package).

Examples

model <- kofn(k = 2, m = 2, component = dfr_exponential())
set.seed(42)
df <- rdata_scheme1(model)(c(1, 2), n = 50, delta = 1.0)
result <- fit_scheme1(model)(df)
coef(result)

Hessian of the log-likelihood for exponential k-out-of-n model

Description

Returns a closure ⁠function(df, par)⁠ that computes the Hessian matrix of the log-likelihood with respect to the rate parameters.

Usage

## S3 method for class 'exp_kofn'
hess_loglik(model, ...)

Arguments

Details

Uses numerical differentiation via numDeriv::hessian() applied to the log-likelihood closure.

Value

A function ⁠function(df, par)⁠ returning an ⁠m x m⁠ Hessian matrix.

Examples

model <- kofn(k = 2, m = 2, component = dfr_exponential())
H <- hess_loglik(model)
set.seed(1)
df <- rdata(model)(c(1, 2), n = 30)
H(df, c(1, 2))  # 2x2 Hessian matrix

Hessian of log-likelihood for Weibull k-out-of-n system

Description

Returns a closure that computes the Hessian matrix of the log-likelihood via numerical differentiation using numDeriv::hessian.

Usage

## S3 method for class 'wei_kofn'
hess_loglik(model, ...)

Arguments

Value

A function function(df, par) returning a square matrix (the Hessian of the log-likelihood).

Examples

model <- kofn(k = 2, m = 2, component = dfr_weibull())
H <- hess_loglik(model)
set.seed(1)
df <- rdata(model)(c(1.5, 2.0, 2.0, 3.0), n = 30)
H(df, c(1.5, 2.0, 2.0, 3.0))  # 4x4 Hessian matrix

Inclusion-exclusion expansion of a product of CDFs

Description

For a set of exponential rates lam, computes the inclusion-exclusion expansion of \prod_{i} (1 - e^{-\lambda_i t}).

Usage

ie_expand(lam)

Arguments

Details

Returns sign and rate-sum vectors such that:

\prod_i (1 - e^{-\lambda_i t}) = \sum_k \text{sign}_k \cdot e^{-\text{rate\_sum}_k \cdot t}

The number of terms is 2^m where m is length(lam).

Value

A list with elements:

sign

Integer vector of signs (+1 or -1).

rate_sum

Numeric vector of summed rates for each term.

Examples

ie <- ie_expand(c(1, 2))
# Product: (1 - exp(-t)) * (1 - exp(-2t))
# = 1 - exp(-t) - exp(-2t) + exp(-3t)
ie$sign      # c(1, -1, -1, 1)
ie$rate_sum  # c(0, 1, 2, 3)

Test if an object is a kofn model

Description

Test if an object is a kofn model

Usage

is_kofn(x)

Arguments

Value

Logical indicating whether x inherits from "kofn".

Examples

is_kofn(kofn(k = 3, m = 3, component = dfr_exponential()))
is_kofn(42)

Create a k-out-of-n system estimation model

Description

Constructs a likelihood model for component lifetime estimation from k-out-of-n system data. The system fails when k components have failed: k=1 is a series system (one failure kills it), k=m is a parallel system (all must fail).

Usage

kofn(
  k = 1L,
  m = 2L,
  component = dfr_exponential(),
  method = "mle",
  lifetime = "t",
  omega = "omega",
  lifetime_upper = "t_upper"
)

Arguments

Details

This model satisfies the likelihood_model concept from the likelihood.model package by providing methods for loglik, score, and hess_loglik.

The class hierarchy is:

• "exp_kofn" or "wei_kofn" (distribution-specific dispatch)

• "kofn" (shared methods)

• "likelihood_model" (generic inference infrastructure)

The concrete subclass is determined by the type of component. The current closed-form machinery is homogeneous: all m components share the same distribution family, and component acts as a prototype for that family.

For non-k-of-n topologies (bridges, arbitrary coherent systems), use coherent_dist or one of the topology shortcuts in dist.structure directly. kofn is exclusively for the k-out-of-n family.

Value

An S3 object of class c("exp_kofn"/"wei_kofn", "kofn", "likelihood_model").

Examples

# Parallel system with 3 exponential components (k = m)
model <- kofn(k = 3, m = 3, component = dfr_exponential())
print(model)

# Series system (k = 1)
model_series <- kofn(k = 1, m = 4, component = dfr_exponential())

# Weibull parallel system with EM estimation
model_wei <- kofn(k = 2, m = 2, component = dfr_weibull(), method = "em")

Log-likelihood for exponential k-out-of-n model

Description

Returns a closure ⁠function(df, par)⁠ that computes the log-likelihood of exponential component rates given system-level data.

Usage

## S3 method for class 'exp_kofn'
loglik(model, ...)

Arguments

Details

For parallel systems (k = m), the closed-form IE expansion is used, supporting all four observation types: exact, right, left, and interval censored. For general k-out-of-n systems, the per-observation contributions are computed via dist.structure::exp_kofn(k, par) and the algebraic.dist generics density, surv, and cdf.

Value

A function ⁠function(df, par)⁠ where df is a data frame with columns for lifetime, observation type, and (for interval-censored) upper bounds, and par is a numeric vector of m component rates.

Examples

model <- kofn(k = 3, m = 3, component = dfr_exponential())
ll <- loglik(model)
set.seed(1)
df <- rdata(model)(c(1, 2, 3), n = 30)
ll(df, c(1, 2, 3))

Log-likelihood for Weibull k-out-of-n system

Description

Returns a closure function(df, par) that computes the log-likelihood for a Weibull k-out-of-n system given data and parameters.

Usage

## S3 method for class 'wei_kofn'
loglik(model, ...)

Arguments

Details

For parallel systems (k = m), uses the direct Weibull system density f_{sys}(t) = \sum_j f_j(t) \prod_{i \ne j} F_i(t). Supports "exact" and "right" observation types. Left and interval censoring are not supported for Weibull (no IE expansion); use the exponential model or Scheme 1 for those cases.

For general k, delegates per-observation to dist.structure::wei_kofn(k, shapes, scales) via the algebraic.dist generics density, surv, and cdf.

Value

A function function(df, par) returning a scalar log-likelihood.

Examples

model <- kofn(k = 2, m = 2, component = dfr_weibull())
ll <- loglik(model)
set.seed(1)
df <- rdata(model)(c(1.5, 2.0, 2.0, 3.0), n = 30)
ll(df, c(1.5, 2.0, 2.0, 3.0))

Masked k-out-of-n log-likelihood

Description

Computes the log-likelihood for a k-out-of-n system with candidate failed sets (generalized masking). The data frame must contain columns for system lifetime, observation type, and Boolean candidate set indicators (c1, c2, ..., cm).

Usage

loglik_masked(model, candset = "c", ...)

Arguments

Details

At system failure, k components have failed. The candidate set C_i is a superset of the true failed set. The likelihood marginalizes over all \binom{|C_i|}{k} possible k-subsets.

For series systems (k = 1), this reduces to the standard masked cause-of-failure likelihood \sum_{j \in C_i} w_j(t_i). For parallel systems (k = m), the candidate set must be {1, ..., m} and the likelihood equals the system density.

The function uses the exponential or Weibull distribution functions directly (not the IE expansion), so it works for any family supported by parse_params.

Value

A function function(df, par) returning a scalar log-likelihood.

Examples

# 2-out-of-4 system with candidate failed sets
model <- kofn(k = 2, m = 4, component = dfr_exponential())
ll <- loglik_masked(model)

# Manually construct masked data: k=2 failed, C = {1,2,3}
df <- data.frame(
  t = 1.5, omega = "exact",
  c1 = TRUE, c2 = TRUE, c3 = TRUE, c4 = FALSE
)
ll(df, c(1, 0.8, 0.6, 0.4))

Log-likelihood for Scheme 1 (periodic inspection) parallel system

Description

Returns a closure that computes the log-likelihood under Scheme 1 observation. The likelihood combines the system density at the exact system failure time with interval-censored component contributions.

Usage

loglik_scheme1(model, ...)

Arguments

Details

The log-likelihood for observation i is:

\log L_i(\theta) = \log f_{sys}(t_i) + \sum_j \log[F_j(a_{ij}^+) - F_j(a_{ij}^-)]

where f_{sys} is the parallel system density and [a_{ij}^-, a_{ij}^+) is the inspection interval containing component j's failure time.

Value

A function function(df, par) returning a scalar log-likelihood.

Examples

model <- kofn(k = 2, m = 2, component = dfr_exponential())
ll <- loglik_scheme1(model)
set.seed(1)
df <- rdata_scheme1(model)(c(1, 2), n = 30, delta = 1.0)
ll(df, c(1, 2))

Expected parameter count for a homogeneous kofn system.

Description

Expected parameter count for a homogeneous kofn system.

Usage

n_par_kofn(m, component)

Arguments

Value

Integer, the total parameter count for the flat par vector.

Number of components in a kofn model

Description

Returns the number of components m in the k-out-of-n system.

Usage

## S3 method for class 'kofn'
ncomponents(x, ...)

Arguments

Value

Integer number of components.

Examples

ncomponents(kofn(k = 5, m = 5, component = dfr_exponential()))

Exact observation scheme (no censoring)

Description

Creates an observation functor that records the exact failure time with no censoring. This is the trivial case — included for completeness and for explicit use in Fisher information comparisons.

Usage

observe_exact()

Value

Observation functor: function(t_true) returning a list with t (exact time), omega ("exact"), and t_upper (NA).

Examples

obs <- observe_exact()
obs(42.5)  # list(t = 42.5, omega = "exact", t_upper = NA)

Interval-censoring observation scheme

Description

Creates an observation functor that places all failure times into a fixed window [a, b). Systems failing within the window are interval-censored; systems failing outside are observed exactly.

Usage

observe_interval_censor(a, b)

Arguments

Details

For regular grids, use observe_periodic instead.

Value

Observation functor: function(t_true) returning a list with t, omega ("interval" or "exact"), and t_upper.

Examples

obs <- observe_interval_censor(a = 5, b = 10)
obs(7)   # interval: list(t = 5, omega = "interval", t_upper = 10)
obs(3)   # exact:    list(t = 3, omega = "exact", t_upper = NA)
obs(12)  # exact:    list(t = 12, omega = "exact", t_upper = NA)

Left-censoring observation scheme

Description

Creates an observation functor that applies left-censoring at time tau. Systems that fail after tau are observed exactly; systems that fail before tau are left-censored (we only know T < tau).

Usage

observe_left_censor(tau)

Arguments

Value

Observation functor: function(t_true) returning a list with t, omega ("left" or "exact"), and t_upper (NA).

Examples

obs <- observe_left_censor(tau = 10)
obs(5)   # left:  list(t = 10, omega = "left", t_upper = NA)
obs(15)  # exact: list(t = 15, omega = "exact", t_upper = NA)

Mixture of observation schemes

Description

Creates an observation functor that randomly selects from a set of observation schemes for each observation. This models heterogeneous monitoring environments where different units may be observed under different protocols.

Usage

observe_mixture(..., weights = NULL)

Arguments

Value

Observation functor: function(t_true) that randomly selects a constituent scheme and returns its output.

Examples

obs <- observe_mixture(
  observe_right_censor(tau = 100),
  observe_periodic(delta = 10, tau = 100),
  weights = c(0.7, 0.3)
)
set.seed(42)
obs(30)

Periodic inspection observation scheme

Description

Creates an observation functor for periodic inspections at intervals of width delta. The failure time is known only to lie within the inspection interval containing it (interval-censored). Systems surviving past tau are right-censored.

Usage

observe_periodic(delta, tau = Inf)

Arguments

Value

Observation functor: function(t_true) returning a list with t (interval lower bound or tau), omega ("interval" or "right"), and t_upper (interval upper bound or NA).

Examples

obs <- observe_periodic(delta = 10, tau = 100)
obs(25)   # interval: list(t = 20, omega = "interval", t_upper = 30)
obs(150)  # right:    list(t = 100, omega = "right", t_upper = NA)

Right-censoring observation scheme

Description

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.

Usage

observe_right_censor(tau = Inf)

Arguments

Value

Observation functor: function(t_true) returning a list with t, omega ("exact" or "right"), and t_upper (NA).

Examples

obs <- observe_right_censor(tau = 100)
obs(50)   # exact:  list(t = 50, omega = "exact", t_upper = NA)
obs(150)  # right:  list(t = 100, omega = "right", t_upper = NA)

# No censoring
obs_all <- observe_right_censor(tau = Inf)
obs_all(999)  # exact: list(t = 999, omega = "exact", t_upper = NA)

Parse a flat parameter vector into shapes and scales

Description

Unified extraction of Weibull shape/scale parameters from the flat vectors used by optimizers. For exponential components, shapes are all 1 and scales are 1/rate.

Usage

parse_params(par, m, component)

Arguments

Value

A list with components shapes and scales.

Print method for kofn models

Description

Displays a human-readable summary of the k-out-of-n system model, including system type, component distribution, estimation method, and column name conventions.

Usage

## S3 method for class 'kofn'
print(x, ...)

Arguments

Value

The model object, invisibly.

Examples

print(kofn(k = 3, m = 3, component = dfr_exponential()))

Random data generation for exponential k-out-of-n model

Description

Returns a closure that generates random system-level data from the exponential k-out-of-n data-generating process.

Usage

## S3 method for class 'exp_kofn'
rdata(model, ...)

Arguments

Details

Workflow:

1. Generate i.i.d. exponential component lifetimes.

2. Compute system lifetime as the (m - k + 1)-th order statistic.

3. Apply the observation functor (exact observation by default).

Value

A function ⁠function(theta, n, observe = NULL)⁠ returning a data frame with columns t (system lifetime) and omega (observation type). Latent component lifetimes, true critical component, and the true parameters are stored as attributes.

Examples

model <- kofn(k = 3, m = 3, component = dfr_exponential())
gen <- rdata(model)
set.seed(1)
df <- gen(theta = c(1, 2, 3), n = 20)
head(df)

# With right-censoring
df2 <- gen(c(1, 2, 3), n = 20, observe = observe_right_censor(tau = 2))
table(df2$omega)

Generate random data from a Weibull k-out-of-n system

Description

Returns a closure that generates system lifetime data with observation scheme support. Component lifetimes are Weibull-distributed with parameters specified by theta (interleaved shape/scale).

Usage

## S3 method for class 'wei_kofn'
rdata(model, ...)

Arguments

Details

The interface matches rdata.exp_kofn: the returned closure accepts an observe functor for observation schemes (censoring, periodic inspection, etc.). The output includes an omega column.

Value

A function function(theta, n, observe = NULL) returning a data frame with columns t (system lifetimes) and omega (observation type). Attributes:

comp_times

n x m matrix of component lifetimes.

par

The parameter vector used for generation.

Examples

model <- kofn(k = 2, m = 2, component = dfr_weibull())
gen <- rdata(model)
set.seed(42)
df <- gen(theta = c(1.5, 2.0, 2.0, 3.0), n = 50)
head(df)

# With right-censoring
df_cens <- gen(c(1.5, 2.0, 2.0, 3.0), 50,
               observe = observe_right_censor(tau = 3))
table(df_cens$omega)

Generate masked k-out-of-n data

Description

Returns a closure that generates system lifetime data with candidate failed sets. The true failed set (the k components that failed by T_sys) is always included in the candidate set. Additional non-failed components are included independently with probability p_mask.

Usage

rdata_masked(model, candset = "c", ...)

Arguments

Value

A function function(theta, n, p_mask = 0, observe = NULL) returning a data frame with columns for system lifetime, observation type, and Boolean candidate set indicators.

Examples

model <- kofn(k = 2, m = 4, component = dfr_exponential())
gen <- rdata_masked(model)
set.seed(42)
df <- gen(theta = c(1, 0.8, 0.6, 0.4), n = 10, p_mask = 0.3)
head(df)

Generate Scheme 1 (periodic inspection) data

Description

Returns a closure that generates k-out-of-n system data under periodic inspection. Each observation yields the exact system failure time and interval-censored component failure times.

Usage

rdata_scheme1(model, ...)

Arguments

Details

Note: The Scheme 1 likelihood uses a composite approximation that treats the system density and component interval contributions as independent. This works well for k >= 2 but is unreliable for series systems (k = 1), where surviving components' intervals should be conditioned on survival past T_sys. For series systems, use the maskedcauses package instead.

Value

A function function(theta, n, delta) returning a data frame with columns:

t

System failure time (exact).

comp_lower_j

Lower bound of inspection interval for component j.

comp_upper_j

Upper bound of inspection interval for component j.

The data frame has attributes comp_times (true component times), delta (inspection interval), and par (true parameters).

Examples

model <- kofn(k = 2, m = 2, component = dfr_exponential())
gen <- rdata_scheme1(model)
set.seed(1)
df <- gen(theta = c(1, 2), n = 20, delta = 1.0)
head(df)

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

dist.structure

ncomponents

flexhaz

dfr_exponential, dfr_weibull

generics

fit

likelihood.model

assumptions, hess_loglik, loglik, rdata, score

Determinant of a positive-definite Hessian, or NA

Description

Computes the numerical Hessian, checks finite-ness and positive definiteness, and returns det(H) or NA.

Usage

safe_hessian_det(neg_ll, par)

Arguments

Value

Scalar determinant, or NA_real_.

Score function for exponential k-out-of-n model

Description

Returns a closure ⁠function(df, par)⁠ that computes the gradient of the log-likelihood with respect to the rate parameters.

Usage

## S3 method for class 'exp_kofn'
score(model, ...)

Arguments

Details

Uses numerical differentiation via numDeriv::grad() applied to the log-likelihood closure.

Value

A function ⁠function(df, par)⁠ returning a numeric vector of length m (the score vector).

Examples

model <- kofn(k = 2, m = 2, component = dfr_exponential())
sc <- score(model)
set.seed(1)
df <- rdata(model)(c(1, 2), n = 30)
sc(df, c(1, 2))  # gradient at true parameters

Score function for Weibull k-out-of-n system

Description

Returns a closure that computes the score (gradient of log-likelihood) via numerical differentiation using numDeriv::grad.

Usage

## S3 method for class 'wei_kofn'
score(model, ...)

Arguments

Value

A function function(df, par) returning a numeric vector (the gradient of the log-likelihood).

Examples

model <- kofn(k = 2, m = 2, component = dfr_weibull())
sc <- score(model)
set.seed(1)
df <- rdata(model)(c(1.5, 2.0, 2.0, 3.0), n = 30)
sc(df, c(1.5, 2.0, 2.0, 3.0))

Default MLE solver for positive parameters

Description

Multi-start optimization via compositional.mle: L-BFGS-B with positivity constraints as the primary solver, Nelder-Mead on the log-parameter scale as fallback. Runs from n_starts random perturbations of par0 and returns the best result.

Usage

solve_mle(neg_ll, par0, n_par, n_starts = 5L, nobs = NULL)

Arguments

Value

An mle_numerical result object (from algebraic.mle) with coef(), vcov(), logLik(), etc. Returns NULL if all starts fail.

Examples

# Fit a simple exponential rate from positive data
x <- rexp(100, rate = 0.5)
neg_ll <- function(rate) -sum(dexp(x, rate, log = TRUE))
fit <- solve_mle(neg_ll, par0 = 1, n_par = 1, nobs = length(x))
coef(fit)

Vectorized truncated log-moment of the Weibull distribution

Description

Computes E[\log T | T < t_i] for T \sim \text{Weibull}(\alpha, \beta) simultaneously for a vector of truncation points.

Usage

trunc_log_moment_vec(v_max_vec, alpha, beta)

Arguments

Details

Uses the identity:

E[\log T | T < t] = \log \beta + \frac{1}{\alpha} E[\log U | U < v_{\max}]

where U \sim \text{Exp}(1) and v_{\max} = (t/\beta)^\alpha.

The inner expectation E[\log U | U < v_{\max}] is computed via central finite difference of the lower incomplete gamma function at s = 1.

For very small v_{\max} (below 1e-10), the asymptotic approximation E[\log T | T < t] \approx \log t - 1/\alpha is used. For large v_{\max} where numerical issues arise, the untruncated moment E[\log T] = \log \beta + \psi(1)/\alpha is used as fallback, where \psi is the digamma function.

Value

Numeric vector of E[\log T | T < t_i] values, same length as v_max_vec.

Examples

# E[log T | T < t] for Weibull(shape=2, scale=1) at t = 0.5, 1.0, 2.0
v_max <- (c(0.5, 1.0, 2.0) / 1)^2
trunc_log_moment_vec(v_max_vec = v_max, alpha = 2, beta = 1)

Scalar truncated power moment of the Weibull distribution

Description

Convenience wrapper around trunc_pow_moment_vec() for a single truncation point. Computes E[T^k | T < t] for T \sim \text{Weibull}(\alpha, \beta).

Usage

trunc_pow_moment(k, t, alpha, beta)

Arguments

Value

Numeric scalar E[T^k | T < t].

See Also

trunc_pow_moment_vec() for the vectorized version.

Examples

# E[T^1 | T < 2] for Weibull(shape=1.5, scale=1)
trunc_pow_moment(k = 1, t = 2, alpha = 1.5, beta = 1)

Vectorized truncated power moment of the Weibull distribution

Description

Computes E[T^k | T < t_i] for T \sim \text{Weibull}(\alpha, \beta) simultaneously for a vector of truncation points.

Usage

trunc_pow_moment_vec(k, v_max_vec, alpha, beta)

Arguments

Details

Uses the substitution U = (T/\beta)^\alpha \sim \text{Exp}(1), so that T^k = \beta^k U^{k/\alpha} and the truncation T < t becomes U < v_{\max}. The result is:

E[T^k | T < t] = \beta^k \frac{\gamma(k/\alpha + 1, v_{\max})}{1 - e^{-v_{\max}}}

where \gamma(s, x) is the lower incomplete gamma function.

Computation is done in log-space for numerical stability.

For very small v_{\max} (below 1e-12), the asymptotic approximation E[T^k | T < t] \approx t^k / (k/\alpha + 1) is used to avoid numerical issues.

Value

Numeric vector of E[T^k | T < t_i] values, same length as v_max_vec.

Examples

# E[T^2 | T < t] for Weibull(shape=2, scale=1) at t = 0.5, 1.0, 2.0
v_max <- (c(0.5, 1.0, 2.0) / 1)^2
trunc_pow_moment_vec(k = 2, v_max_vec = v_max, alpha = 2, beta = 1)

Component weight for exponential parallel systems

Description

Computes the contribution of component j to the parallel system density at time t. For component j with rate \lambda_j and other rates \lambda_{-j}:

w_j(t) = \lambda_j e^{-\lambda_j t} \prod_{i \neq j} (1 - e^{-\lambda_i t})

Usage

w_j_exact(t, par, j)

Arguments

Details

Uses the inclusion-exclusion expansion to express this as a finite sum of exponentials.

Value

Scalar value of w_j(t).

See Also

w_j_integral() for the closed-form integral of w_j.

Examples

# Component 1 contribution at t = 0.5, rates = c(1, 2, 3)
w_j_exact(0.5, c(1, 2, 3), j = 1)

Closed-form integral of w_j(t) over an interval

Description

Computes \int_a^b w_j(t)\, dt in closed form using the inclusion-exclusion expansion. Each term integrates as:

\int_a^b e^{-r t}\, dt = \frac{e^{-ra} - e^{-rb}}{r}

Usage

w_j_integral(a, b, par, j)

Arguments

Value

Scalar integral value.

See Also

w_j_exact() for the pointwise evaluation.

Examples

# Integral of w_1(t) from 0 to 1, rates = c(1, 2)
w_j_integral(0, 1, c(1, 2), j = 1)

System survival function for Weibull parallel system (internal)

Description

Computes S_sys(t) = 1 - prod_j F_j(t) for a parallel system.

Usage

weibull_S_sys(t, shapes, scales)

Arguments

Value

Numeric scalar S_sys(t).

System density for Weibull parallel system (internal)

Description

Computes f_{sys}(t) = \sum_j f_j(t) \prod_{i \neq j} F_i(t).

Usage

weibull_f_sys(t, shapes, scales)

Arguments

Value

Numeric scalar f_sys(t).

Component weight for Weibull parallel system (internal)

Description

Computes w_j(t) = f_j(t) \prod_{i \neq j} F_i(t), the contribution of component j to the parallel system density at time t.

Usage

weibull_w_j(t, shapes, scales, j)

Arguments

Value

Numeric scalar w_j(t).