---
title: "A Likelihood Framework for Masked Series Systems"
output:
    rmarkdown::html_vignette:
        toc: true
vignette: >
  %\VignetteIndexEntry{A Likelihood Framework for Masked Series Systems}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
header-includes:
- \renewcommand{\v}[1]{\boldsymbol{#1}}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  echo = TRUE,
  comment = "#>")

# Set to TRUE to regenerate long-running simulation results
run_long <- FALSE

library(maskedcauses)
old_opts <- options(digits = 4)
```


Introduction {#introduction}
============================

In reliability engineering, many systems are designed as **series systems**:
the system functions only when all of its $m$ components function. When
the system fails, the failure time is observed, but the *component* that
caused the failure is often unknown or only partially identified. This
situation --- where the failed component is hidden behind incomplete
diagnostic information --- is called **masked failure cause data**.

This vignette develops a self-contained mathematical framework for
estimating component lifetime parameters from masked series system data.
Everything is derived from first principles; no external paper is required.
The framework is then instantiated for three distribution families
implemented in `maskedcauses`:

1. **Exponential** (`exp_series_md_c1_c2_c3`) --- constant failure rates
2. **Homogeneous Weibull** (`wei_series_homogeneous_md_c1_c2_c3`) ---
   shared shape, individual scales
3. **Heterogeneous Weibull** (`wei_series_md_c1_c2_c3`) --- individual
   shapes and scales

Each section pairs the mathematical derivation with live R code that
demonstrates the corresponding package functions.


Series System Model {#series-system}
=====================================

## Definition

Consider a system of $m$ independent components. The lifetime of
component $j$ in the $i$-th system is $T_{ij}$, with reliability
(survival) function $R_j(t) = P(T_{ij} > t)$, pdf $f_j(t)$, and hazard
function
$$
  h_j(t) = \frac{f_j(t)}{R_j(t)}.
$$

In a series configuration ("weakest link"), the system fails as soon as
any component fails:
$$
  T_i = \min(T_{i1}, \ldots, T_{im}).
$$

**Assumption (independence).**
Component lifetimes $T_{i1}, \ldots, T_{im}$ are mutually independent
given the parameter vector $\theta$.

## System reliability

The system reliability is the product of component reliabilities:
\begin{equation}
\label{eq:sys_reliability}
R(t \mid \theta) = P(T_i > t \mid \theta) = \prod_{j=1}^m R_j(t \mid \theta_j).
\end{equation}
\begin{proof}
$P(\min_j T_{ij} > t) = P(T_{i1} > t, \ldots, T_{im} > t) = \prod_j P(T_{ij} > t)$,
where the last equality uses independence.
\end{proof}

## Additive hazards

The system hazard is the sum of component hazards:
\begin{equation}
\label{eq:sys_hazard}
h(t \mid \theta) = \sum_{j=1}^m h_j(t \mid \theta_j).
\end{equation}
\begin{proof}
By definition, $h(t) = f(t) / R(t)$. The system pdf satisfies
$f(t) = -R'(t) = R(t) \sum_j h_j(t)$, from differentiating the log of
Eq.\ \eqref{eq:sys_reliability}. Dividing by $R(t)$ yields the result.
\end{proof}

The additive hazard decomposition is the central structural result for
series systems. It means that each component contributes independently
to the instantaneous system failure risk at every time $t$.

The system pdf follows immediately:
$$
  f(t \mid \theta) = h(t \mid \theta) \cdot R(t \mid \theta)
    = \left[\sum_{j=1}^m h_j(t)\right]
      \prod_{l=1}^m R_l(t).
$$

### R code: Exponential example

We demonstrate with a 3-component exponential series system. Each
component has a constant hazard (failure rate) $\lambda_j$:

```{r sys-demo}
model_exp <- exp_series_md_c1_c2_c3()

# Component failure rates
theta_exp <- c(1.0, 1.1, 0.95)
m <- length(theta_exp)

# Extract hazard closures
h1 <- component_hazard(model_exp, 1)
h2 <- component_hazard(model_exp, 2)
h3 <- component_hazard(model_exp, 3)

t_grid <- seq(0.01, 3, length.out = 200)
h_sys <- h1(t_grid, theta_exp) + h2(t_grid, theta_exp) + h3(t_grid, theta_exp)

cat("System hazard (constant):", h_sys[1], "\n")
cat("Sum of component rates:  ", sum(theta_exp), "\n")
```

For exponential components, the system hazard is constant and equals the
sum of the individual rates --- consistent with the additive hazard
theorem.


Component Cause of Failure {#cause-of-failure}
===============================================

## Definition

Let $K_i$ denote the index of the component that caused the $i$-th
system failure:
$$
  K_i = \arg\min_{j \in \{1,\ldots,m\}} T_{ij}.
$$

$K_i$ is a discrete random variable on $\{1, \ldots, m\}$ that, in
general, depends on the failure time $T_i$.

## Conditional cause probability

The probability that component $j$ caused the system failure, given that
the failure occurred at time $t$, is:
\begin{equation}
\label{eq:cond_cause}
P(K_i = j \mid T_i = t, \theta)
  = \frac{h_j(t \mid \theta_j)}{\sum_{l=1}^m h_l(t \mid \theta_l)}
  = \frac{h_j(t)}{h(t)}.
\end{equation}
\begin{proof}
The joint density of $(K_i, T_i)$ is $f_{K,T}(j, t) = f_j(t) \prod_{l \neq j} R_l(t)$.
The marginal density of $T_i$ is $f_T(t) = \sum_j f_{K,T}(j, t) = R(t) \sum_j h_j(t)$.
Then $P(K = j \mid T = t) = f_{K,T}(j,t) / f_T(t) = h_j(t) / \sum_l h_l(t)$.
\end{proof}

**Intuition.** At any instant $t$, the component with the highest hazard
is the most likely cause. The cause probability is the component's
*share* of the total instantaneous risk --- a proportional allocation.

### R code: Exponential vs Weibull cause probabilities

For exponential components, the cause probability is *constant* in time
(because all hazards are constant). For heterogeneous Weibull components,
the cause probability *varies* with time:

```{r cause-prob-demo, fig.width=7, fig.height=4.5}
# Exponential: constant cause probabilities
ccp_exp <- conditional_cause_probability(model_exp)
probs_exp <- ccp_exp(t_grid, theta_exp)

# Heterogeneous Weibull: time-varying
theta_wei <- c(0.7, 200,   # DFR electronics
               1.0, 150,   # CFR seals
               2.0, 100)   # IFR bearing
model_wei <- wei_series_md_c1_c2_c3(
  shapes = theta_wei[seq(1, 6, 2)],
  scales = theta_wei[seq(2, 6, 2)]
)
ccp_wei <- conditional_cause_probability(model_wei)

t_wei <- seq(0.5, 120, length.out = 200)
probs_wei <- ccp_wei(t_wei, theta_wei)

oldpar <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
on.exit(par(oldpar))

# Exponential
plot(t_grid, probs_exp[, 1], type = "l", col = "steelblue", lwd = 2,
     ylim = c(0, 0.6), xlab = "Time", ylab = "P(K = j | T = t)",
     main = "Exponential (constant)")
lines(t_grid, probs_exp[, 2], col = "forestgreen", lwd = 2, lty = 2)
lines(t_grid, probs_exp[, 3], col = "firebrick", lwd = 2, lty = 3)
legend("right", paste("Comp", 1:3),
       col = c("steelblue", "forestgreen", "firebrick"),
       lty = 1:3, lwd = 2, cex = 0.8)

# Heterogeneous Weibull
plot(t_wei, probs_wei[, 1], type = "l", col = "steelblue", lwd = 2,
     ylim = c(0, 1), xlab = "Time", ylab = "P(K = j | T = t)",
     main = "Heterogeneous Weibull (time-varying)")
lines(t_wei, probs_wei[, 2], col = "forestgreen", lwd = 2, lty = 2)
lines(t_wei, probs_wei[, 3], col = "firebrick", lwd = 2, lty = 3)
legend("right", c("DFR (k=0.7)", "CFR (k=1.0)", "IFR (k=2.0)"),
       col = c("steelblue", "forestgreen", "firebrick"),
       lty = 1:3, lwd = 2, cex = 0.8)
```

In the right panel, the DFR component (high infant-mortality hazard)
dominates at early times, but the IFR component (wear-out) takes over as
the system ages. This crossover cannot be captured by models that force
a common shape parameter.


The Observational Model {#observational-model}
===============================================

In practice, we do not observe the latent component lifetimes
$(T_{i1}, \ldots, T_{im})$ directly. Instead, for each of $n$ systems
we observe a triple $(s_i, \omega_i, c_i)$:

- $s_i$: a time value (possibly censored)
- $\omega_i$: the observation type
- $c_i \subseteq \{1, \ldots, m\}$: the candidate set of possibly-failed
  components

## Four observation types

The package supports four types, each arising from a different monitoring
scheme:

| Type | $\omega_i$ | What we know | Likelihood contribution |
|------|-----------|--------------|------------------------|
| Exact | `"exact"` | System failed at $s_i$ | $R(s_i) \cdot h_c(s_i)$ |
| Right-censored | `"right"` | System survived past $s_i$ | $R(s_i)$ |
| Left-censored | `"left"` | System failed before $s_i$ | $\int_0^{s_i} h_c(u) R(u)\, du$ |
| Interval-censored | `"interval"` | Failure in $(s_i, s_i^{\mathrm{upper}})$ | $\int_{s_i}^{s_i^{\mathrm{upper}}} h_c(u) R(u)\, du$ |

Here $h_c(t) = \sum_{j \in c_i} h_j(t)$ is the candidate-set hazard.

## Masking

When $|c_i| = 1$, the failed component is known exactly. When
$|c_i| > 1$, the cause of failure is **masked**: the diagnostic
identifies a set of possible culprits but not the specific one.

## Observe functors

The package provides composable observation functors for generating data
under different monitoring schemes:

```{r observe-demo}
# Continuous monitoring with right-censoring at tau
obs_right <- observe_right_censor(tau = 5)
obs_right(3.2)   # failure before tau -> exact
obs_right(7.1)   # survival past tau  -> right-censored

# Single inspection at tau (left-censoring)
obs_left <- observe_left_censor(tau = 5)
obs_left(3.2)    # failed before inspection -> left-censored
obs_left(7.1)    # surviving at inspection  -> right-censored

# Periodic inspections every delta time units
obs_periodic <- observe_periodic(delta = 1, tau = 5)
obs_periodic(3.2)  # failure in (3, 4) -> interval-censored

# Mixture: random assignment to schemes
obs_mixed <- observe_mixture(
  observe_right_censor(tau = 5),
  observe_left_censor(tau = 3),
  observe_periodic(delta = 0.5, tau = 5),
  weights = c(0.5, 0.2, 0.3)
)
```

All `rdata()` methods accept an `observe` argument to generate data under
a specified monitoring scheme:

```{r observe-rdata}
gen_exp <- rdata(model_exp)

set.seed(42)
df_demo <- gen_exp(theta_exp, n = 500, p = 0.3,
                   observe = observe_mixture(
                     observe_right_censor(tau = 5),
                     observe_left_censor(tau = 3),
                     weights = c(0.7, 0.3)
                   ))
cat("Observation type counts:\n")
print(table(df_demo$omega))
```


The C1--C2--C3 Likelihood {#c1-c2-c3}
=======================================

This is the core contribution: three conditions on the masking process
that reduce the full joint likelihood to a tractable form that depends
only on the lifetime parameters $\theta$.

## The three conditions

**C1 (containment).** The candidate set always contains the true failed
component:
$$
  P(K_i \in c_i) = 1.
$$
This is a minimal accuracy requirement on the diagnostic. Without C1,
the candidate set could be misleading, and the likelihood would need to
model diagnostic errors explicitly.

*Example:* An on-board diagnostic (OBD) system reports a set of
possibly-failed modules. C1 requires that the truly failed module is
always in this set --- the OBD may include false positives but never a
false negative.

**C2 (symmetric masking).** The masking probability is the same for all
components in the candidate set:
$$
  P(c_i \mid K_i = j, T_i = t_i, \theta)
    = P(c_i \mid K_i = j', T_i = t_i, \theta)
    \quad \text{for all } j, j' \in c_i.
$$
Intuitively, the diagnostic cannot distinguish among the candidates ---
if it could, it would narrow the candidate set.

*Example:* A line-replaceable unit (LRU) grouping replaces all
components in a subsystem together. The grouping is determined by
physical layout, not by which component failed, so C2 holds by
construction.

**C3 (parameter independence).** The masking probabilities do not depend
on the lifetime parameters $\theta$:
$$
  P(c_i \mid K_i = j, T_i = t_i) \text{ is not a function of } \theta.
$$

*Example:* Diagnostic accuracy depends on sensor calibration and
technician skill, not on component failure rates. If the diagnostic
equipment was calibrated before the study, C3 holds.

## Deriving the likelihood

We derive the likelihood step by step. Consider observation $i$ with
failure time $t_i$ and candidate set $c_i$.

**Step 1: Full joint density.** The observed-data density includes a sum
over the latent cause $K_i$:
$$
  L_i(\theta) = \sum_{j=1}^m f_{K,T}(j, t_i) \cdot P(c_i \mid K_i = j, T_i = t_i, \theta).
$$

**Step 2: Apply C1.** Since $P(K_i \in c_i) = 1$, contributions from
$j \notin c_i$ are zero:
$$
  L_i(\theta) = \sum_{j \in c_i} f_{K,T}(j, t_i) \cdot P(c_i \mid K_i = j, T_i = t_i, \theta).
$$

**Step 3: Apply C2.** Symmetric masking means
$P(c_i \mid K_i = j, T_i = t_i, \theta) = \beta_i$ for all $j \in c_i$:
$$
  L_i(\theta) = \beta_i \sum_{j \in c_i} f_{K,T}(j, t_i).
$$

**Step 4: Apply C3.** Since $\beta_i$ does not depend on $\theta$, it is
a constant with respect to the parameters and drops from the
maximization:
$$
  L_i(\theta) \propto \sum_{j \in c_i} f_{K,T}(j, t_i)
    = \sum_{j \in c_i} h_j(t_i) \cdot R(t_i)
    = R(t_i) \cdot h_c(t_i).
$$

The **exact-failure likelihood contribution** is therefore:
\begin{equation}
\label{eq:exact_contrib}
  \ell_i^{\text{exact}}(\theta)
    = \log R(t_i \mid \theta) + \log h_c(t_i \mid \theta),
\end{equation}
where $h_c(t) = \sum_{j \in c_i} h_j(t)$ is the candidate-set hazard.

## Censored observations

**Right-censored.** The system survived past $t_i$. No failure was
observed, so there is no candidate set and the contribution is simply:
$$
  \ell_i^{\text{right}}(\theta) = \log R(t_i \mid \theta).
$$

**Left-censored.** The system was found failed at inspection time
$\tau_i$, but the exact failure time is unknown. Applying C1--C3 to the
integrated likelihood:
$$
  \ell_i^{\text{left}}(\theta)
    = \log \int_0^{\tau_i} h_c(u) \, R(u) \, du.
$$

**Interval-censored.** The failure occurred in $(a_i, b_i)$:
$$
  \ell_i^{\text{interval}}(\theta)
    = \log \int_{a_i}^{b_i} h_c(u) \, R(u) \, du.
$$

## Combined log-likelihood

The full log-likelihood decomposes as:
\begin{equation}
\label{eq:full_loglik}
  \ell(\theta) = \ell_E + \ell_R + \ell_L + \ell_I,
\end{equation}
where
\begin{align}
  \ell_E &= \sum_{i : \omega_i = \text{exact}}
    \bigl[\log R(t_i) + \log h_{c_i}(t_i)\bigr], \\
  \ell_R &= \sum_{i : \omega_i = \text{right}} \log R(t_i), \\
  \ell_L &= \sum_{i : \omega_i = \text{left}}
    \log \int_0^{\tau_i} h_{c_i}(u) R(u) \, du, \\
  \ell_I &= \sum_{i : \omega_i = \text{interval}}
    \log \int_{a_i}^{b_i} h_{c_i}(u) R(u) \, du.
\end{align}

### R code: Log-likelihood evaluation

The `loglik()` generic returns a closure that evaluates this
log-likelihood for any model:

```{r loglik-demo}
ll_exp <- loglik(model_exp)

# Evaluate at true parameters on the mixed-censoring demo data
ll_val <- ll_exp(df_demo, theta_exp)
cat("Log-likelihood at true theta:", round(ll_val, 4), "\n")

# Perturb parameters and compare
theta_bad <- c(2, 2, 2)
cat("Log-likelihood at wrong theta:", round(ll_exp(df_demo, theta_bad), 4), "\n")
```


Distribution Families {#families}
==================================

The C1--C2--C3 framework applies to *any* component lifetime distribution.
The following table lists common families with their hazard, reliability,
and pdf:

| Family | Parameters | $h_j(t)$ | $R_j(t)$ |
|--------|-----------|-----------|-----------|
| Exponential | $\lambda_j$ | $\lambda_j$ | $e^{-\lambda_j t}$ |
| Weibull | $k_j, \beta_j$ | $\frac{k_j}{\beta_j}\left(\frac{t}{\beta_j}\right)^{k_j-1}$ | $\exp\!\left(-\left(\frac{t}{\beta_j}\right)^{k_j}\right)$ |
| Pareto | $\alpha_j, x_{m,j}$ | $\frac{\alpha_j}{t}$ | $\left(\frac{x_{m,j}}{t}\right)^{\alpha_j}$ |
| Log-normal | $\mu_j, \sigma_j$ | $\frac{\phi(z_j)}{t\sigma_j \Phi(-z_j)}$ | $\Phi(-z_j)$ |
| Gamma | $\alpha_j, \beta_j$ | $f_j(t)/R_j(t)$ | $1 - \frac{\gamma(\alpha_j, t/\beta_j)}{\Gamma(\alpha_j)}$ |

where $z_j = (\log t - \mu_j)/\sigma_j$.
The package currently implements the first two families (Exponential and
Weibull). The framework extends to any family by substituting the
appropriate $h_j$ and $R_j$ into the general log-likelihood.


Worked Example: Exponential Components {#exponential}
=====================================================

The exponential model is the simplest instantiation: all hazards are
constant, yielding closed-form expressions for every quantity.

## Setup

```{r exp-setup}
theta <- c(1.0, 1.1, 0.95, 1.15, 1.1)
m <- length(theta)
model <- exp_series_md_c1_c2_c3()

cat("True rates:", theta, "\n")
cat("System rate:", sum(theta), "\n")
```

## Specializing the likelihood

For exponential components with $h_j(t) = \lambda_j$ and
$R_j(t) = e^{-\lambda_j t}$:

\begin{align}
\text{Exact:}     \quad & \log\!\left(\sum_{j \in c_i} \lambda_j\right)
  - \left(\sum_{j=1}^m \lambda_j\right) t_i, \\
\text{Right:}     \quad & -\lambda_{\text{sys}} t_i, \\
\text{Left:}      \quad & \log \lambda_c + \log(1 - e^{-\lambda_{\text{sys}} \tau_i})
  - \log \lambda_{\text{sys}}, \\
\text{Interval:}  \quad & \log \lambda_c - \lambda_{\text{sys}} a_i
  + \log(1 - e^{-\lambda_{\text{sys}}(b_i - a_i)}) - \log \lambda_{\text{sys}},
\end{align}
where $\lambda_c = \sum_{j \in c_i} \lambda_j$ and
$\lambda_{\text{sys}} = \sum_j \lambda_j$.

All four observation types have **closed-form** log-likelihood, score,
*and* Hessian --- the exponential model is unique in this respect.

## Data generation and fitting

```{r exp-generate}
gen <- rdata(model)
set.seed(7231)

df <- gen(theta, n = 2000, p = 0.3,
          observe = observe_right_censor(tau = 3))

cat("Observation types:\n")
print(table(df$omega))
```

```{r exp-fit, warning=FALSE}
solver <- fit(model)
theta0 <- rep(1, m)
estimate <- solver(df, par = theta0, method = "Nelder-Mead")
```

```{r exp-results}
recovery <- data.frame(
  Component = 1:m,
  True = theta,
  MLE = estimate$par,
  SE = sqrt(diag(estimate$vcov)),
  Rel_Error_Pct = 100 * (estimate$par - theta) / theta
)

knitr::kable(recovery, digits = 4,
             caption = "Exponential MLE: parameter recovery",
             col.names = c("Component", "True", "MLE", "SE", "Rel. Error %"))
```

## Score and Hessian verification

The analytical score should vanish at the MLE. We also verify it
against numerical differentiation:

```{r exp-score-verify}
ll_fn <- loglik(model)
scr_fn <- score(model)
hess_fn <- hess_loglik(model)

scr_at_mle <- scr_fn(df, estimate$par)
cat("Score at MLE:", round(scr_at_mle, 4), "\n")

scr_num <- numDeriv::grad(function(th) ll_fn(df, th), estimate$par)
cat("Max |analytical - numerical| score:",
    formatC(max(abs(scr_at_mle - scr_num)), format = "e", digits = 2), "\n")

# Hessian eigenvalues (should all be negative for concavity)
H <- hess_fn(df, estimate$par)
cat("Hessian eigenvalues:", round(eigen(H)$values, 2), "\n")
```

## Fisher information and confidence intervals

The observed Fisher information is $I(\hat\theta) = -H(\hat\theta)$.
Asymptotic $(1-\alpha)$% Wald confidence intervals are:
$$
  \hat\theta_j \pm z_{1-\alpha/2} \sqrt{I(\hat\theta)^{-1}_{jj}}.
$$

```{r exp-ci}
alpha <- 0.05
z <- qnorm(1 - alpha / 2)
se <- sqrt(diag(estimate$vcov))

ci_df <- data.frame(
  Component = 1:m,
  Lower = estimate$par - z * se,
  MLE = estimate$par,
  Upper = estimate$par + z * se,
  True = theta,
  Covered = (estimate$par - z * se <= theta) & (theta <= estimate$par + z * se)
)

knitr::kable(ci_df, digits = 4,
             caption = "95% Wald confidence intervals",
             col.names = c("Comp.", "Lower", "MLE", "Upper", "True", "Covered?"))
```


Homogeneous Weibull {#homogeneous-weibull}
==========================================

The **homogeneous Weibull** model constrains all $m$ components to share
a common shape parameter $k$, while allowing individual scale parameters
$\beta_1, \ldots, \beta_m$. The parameter vector is
$\theta = (k, \beta_1, \ldots, \beta_m) \in \mathbb{R}^{m+1}_{> 0}$.

## Key properties

**Property 1: The system lifetime is Weibull.**
Because all shapes are equal, the system reliability simplifies:
$$
  R(t) = \prod_{j=1}^m \exp\!\left(-\left(\frac{t}{\beta_j}\right)^k\right)
       = \exp\!\left(-\left(\frac{t}{\beta_{\text{sys}}}\right)^k\right),
$$
where
$\beta_{\text{sys}} = \bigl(\sum_{j=1}^m \beta_j^{-k}\bigr)^{-1/k}$.
This closure under the minimum is the structural advantage of the
homogeneous model.

**Property 2: Constant cause probabilities.**
The time dependence $t^{k-1}$ cancels in the hazard ratio:
$$
  P(K = j \mid T = t) = \frac{\beta_j^{-k}}{\sum_l \beta_l^{-k}}
    \eqqcolon w_j.
$$
Just as in the exponential case, the conditional cause probability does
not depend on $t$.

**Property 3: Closed-form censored contributions.**
Left- and interval-censored terms factor as $\log w_{c_i} + \log [\text{Weibull CDF difference}]$ --- no numerical integration needed.

## R code: Setup and hazard visualization

```{r hom-setup, fig.width=7, fig.height=4.5}
theta_hom <- c(k = 1.5, beta1 = 100, beta2 = 150, beta3 = 200)
k <- theta_hom[1]
scales <- theta_hom[-1]
m_hom <- length(scales)

model_hom <- wei_series_homogeneous_md_c1_c2_c3()

# System scale
beta_sys <- wei_series_system_scale(k, scales)
cat("System scale:", round(beta_sys, 2), "\n")

# Plot component hazards
t_grid <- seq(0.1, 200, length.out = 300)
cols <- c("steelblue", "forestgreen", "firebrick")

plot(NULL, xlim = c(0, 200), ylim = c(0, 0.06),
     xlab = "Time", ylab = "Hazard h_j(t)",
     main = "Homogeneous Weibull: component hazards (k = 1.5)")
for (j in seq_len(m_hom)) {
  h_j <- component_hazard(model_hom, j)
  lines(t_grid, h_j(t_grid, theta_hom), col = cols[j], lwd = 2)
}
legend("topleft", paste0("Comp ", 1:m_hom, " (beta=", scales, ")"),
       col = cols, lwd = 2, cex = 0.9)
grid()
```

All hazards are increasing ($k > 1$), sharing the same shape but
differing in magnitude. Component 1 (smallest $\beta$) dominates at all
times.

## Cause probabilities

```{r hom-cause-prob}
# Analytical cause weights
w <- scales^(-k) / sum(scales^(-k))
names(w) <- paste0("Component ", 1:m_hom)
cat("Time-invariant cause weights:\n")
print(round(w, 4))

# Verify with package function (pass scales so model knows m)
ccp_fn <- conditional_cause_probability(
  wei_series_homogeneous_md_c1_c2_c3(scales = scales)
)
probs <- ccp_fn(c(10, 50, 100, 150), theta_hom)
knitr::kable(probs, digits = 4,
             caption = "P(K=j | T=t) at four time points (should be constant)",
             col.names = paste0("Comp ", 1:m_hom))
```

## Data generation and MLE

```{r hom-generate}
gen_hom <- rdata(model_hom)
set.seed(2024)

df_hom <- gen_hom(theta_hom, n = 500, p = 0.3,
                  observe = observe_periodic(delta = 20, tau = 250))

cat("Observation types:\n")
print(table(df_hom$omega))
```

```{r hom-fit, warning=FALSE}
solver_hom <- fit(model_hom)
theta0_hom <- c(1.2, 110, 140, 180)

est_hom <- solver_hom(df_hom, par = theta0_hom, method = "Nelder-Mead")

recovery_hom <- data.frame(
  Parameter = c("k", paste0("beta_", 1:m_hom)),
  True = theta_hom,
  MLE = est_hom$par,
  SE = sqrt(diag(est_hom$vcov)),
  Rel_Error_Pct = 100 * (est_hom$par - theta_hom) / theta_hom
)

knitr::kable(recovery_hom, digits = 3,
             caption = paste0("Homogeneous Weibull MLE (n = 500)"),
             col.names = c("Parameter", "True", "MLE", "SE", "Rel. Error %"))
```


Heterogeneous Weibull {#heterogeneous-weibull}
===============================================

The **heterogeneous Weibull** model allows each component to have its own
shape $k_j$ and scale $\beta_j$, giving a $2m$-dimensional parameter
vector $\theta = (k_1, \beta_1, \ldots, k_m, \beta_m)$.

## Differences from the homogeneous model

| Property | Homogeneous | Heterogeneous |
|----------|-------------|---------------|
| Parameters | $m + 1$ | $2m$ |
| System lifetime | Weibull($k$, $\beta_{\text{sys}}$) | **Not** Weibull |
| Cause probabilities | Constant in $t$ | Time-varying |
| Left/interval likelihood | Closed form | Numerical integration |

The additional flexibility comes at a cost: the system lifetime
$T = \min(T_1, \ldots, T_m)$ is no longer Weibull when shapes differ,
and left-censored or interval-censored likelihood contributions require
numerical integration (`stats::integrate`).

## R code: Mixed failure regimes

We construct a 3-component system with DFR (infant mortality), CFR
(random), and IFR (wear-out) components:

```{r het-setup, fig.width=7, fig.height=4.5}
theta_het <- c(0.7, 200,   # DFR electronics
               1.0, 150,   # CFR seals
               2.0, 100)   # IFR bearing
m_het <- length(theta_het) / 2

model_het <- wei_series_md_c1_c2_c3(
  shapes = theta_het[seq(1, 6, 2)],
  scales = theta_het[seq(2, 6, 2)]
)

# Component hazard functions
t_het <- seq(0.5, 120, length.out = 300)

h_fns <- lapply(1:m_het, function(j) component_hazard(model_het, j))
h_vals <- sapply(h_fns, function(hf) hf(t_het, theta_het))

plot(t_het, h_vals[, 1], type = "l", col = "steelblue", lwd = 2,
     ylim = c(0, 0.05),
     xlab = "Time", ylab = "Hazard h_j(t)",
     main = "Heterogeneous Weibull: mixed failure regimes")
lines(t_het, h_vals[, 2], col = "forestgreen", lwd = 2, lty = 2)
lines(t_het, h_vals[, 3], col = "firebrick", lwd = 2, lty = 3)
legend("topright",
       c("Comp 1: DFR (k=0.7)", "Comp 2: CFR (k=1.0)", "Comp 3: IFR (k=2.0)"),
       col = c("steelblue", "forestgreen", "firebrick"),
       lty = 1:3, lwd = 2, cex = 0.85)
grid()
```

## Time-varying cause probabilities

```{r het-cause-prob, fig.width=7, fig.height=4.5}
ccp_het <- conditional_cause_probability(model_het)
probs_het <- ccp_het(t_het, theta_het)

plot(t_het, probs_het[, 1], type = "l", col = "steelblue", lwd = 2,
     ylim = c(0, 1), xlab = "Time", ylab = "P(K = j | T = t)",
     main = "Heterogeneous Weibull: time-varying cause probabilities")
lines(t_het, probs_het[, 2], col = "forestgreen", lwd = 2, lty = 2)
lines(t_het, probs_het[, 3], col = "firebrick", lwd = 2, lty = 3)
legend("right",
       c("DFR (k=0.7)", "CFR (k=1.0)", "IFR (k=2.0)"),
       col = c("steelblue", "forestgreen", "firebrick"),
       lty = 1:3, lwd = 2, cex = 0.85)
grid()
```

The crossover pattern is characteristic of bathtub-curve behavior at the
system level: early failures are dominated by infant-mortality mechanisms,
while late failures are dominated by wear-out.

## Data generation and MLE

```{r het-generate}
gen_het <- rdata(model_het)
set.seed(7231)

df_het <- gen_het(theta_het, n = 300, tau = 120, p = 0.3,
                  observe = observe_right_censor(tau = 120))

cat("Observation types:\n")
print(table(df_het$omega))
```

We use right-censoring for speed --- exact and right-censored observations
have closed-form likelihood contributions. Left- and interval-censored
observations require numerical integration per observation per
optimization iteration.

```{r het-fit, warning=FALSE}
solver_het <- fit(model_het)
theta0_het <- rep(c(1, 150), m_het)

est_het <- solver_het(df_het, par = theta0_het, method = "Nelder-Mead")

par_names <- paste0(
  rep(c("k", "beta"), m_het), "_",
  rep(1:m_het, each = 2)
)

recovery_het <- data.frame(
  Parameter = par_names,
  True = theta_het,
  MLE = est_het$par,
  SE = sqrt(diag(est_het$vcov)),
  Rel_Error_Pct = 100 * (est_het$par - theta_het) / theta_het
)

knitr::kable(recovery_het, digits = 3,
             caption = paste0("Heterogeneous Weibull MLE (n = 300)"),
             col.names = c("Parameter", "True", "MLE", "SE", "Rel. Error %"))
```

## Model comparison: heterogeneous vs homogeneous

When the true DGP has heterogeneous shapes, the homogeneous model is
misspecified. We compare the two on the same data:

```{r het-comparison, warning=FALSE}
# Refit on larger dataset for sharper comparison
set.seed(999)
df_comp <- gen_het(theta_het, n = 800, tau = 150, p = 0.2,
                   observe = observe_right_censor(tau = 150))

# Heterogeneous (6 params)
fit_het <- solver_het(df_comp, par = theta0_het, method = "Nelder-Mead")

# Homogeneous (4 params)
model_hom2 <- wei_series_homogeneous_md_c1_c2_c3()
solver_hom2 <- fit(model_hom2)
fit_hom2 <- solver_hom2(df_comp, par = c(1, 150, 150, 150),
                        method = "Nelder-Mead")

# Likelihood ratio test
LRT <- 2 * (fit_het$loglik - fit_hom2$loglik)
df_lrt <- 6 - 4
p_val <- pchisq(LRT, df = df_lrt, lower.tail = FALSE)

comp_df <- data.frame(
  Model = c("Heterogeneous (2m = 6)", "Homogeneous (m+1 = 4)"),
  Params = c(6, 4),
  LogLik = c(fit_het$loglik, fit_hom2$loglik),
  AIC = c(-2 * fit_het$loglik + 2 * 6, -2 * fit_hom2$loglik + 2 * 4)
)

knitr::kable(comp_df, digits = 2,
             caption = "Model comparison",
             col.names = c("Model", "# Params", "Log-Lik", "AIC"))

cat(sprintf("\nLRT statistic: %.2f (df = %d, p = %.4f)\n",
            LRT, df_lrt, p_val))
```

Both AIC and the likelihood ratio test favor the heterogeneous model when
the true DGP has unequal shapes. The common-shape estimate from the
homogeneous fit represents a compromise value between 0.7, 1.0, and 2.0,
leading to biased estimates of the scale parameters.


Monte Carlo Assessment {#monte-carlo}
======================================

We briefly assess finite-sample properties of the MLE for each model
family. Full Monte Carlo studies with sensitivity analyses are provided in
the dedicated vignettes (`exponential_series`, `weibull_series`,
`weibull_homogeneous_series`).

```{r load-precomputed-mc, include=FALSE, eval=!run_long}
list2env(readRDS("precomputed_framework.rds"), envir = environment())
```

```{r mc-setup, eval=run_long}
set.seed(42)
B <- 100
alpha <- 0.05

# Exponential: 5 components
theta_mc_exp <- c(1.0, 1.1, 0.95, 1.15, 1.1)
m_mc_exp <- length(theta_mc_exp)
n_mc <- 1000

model_mc_exp <- exp_series_md_c1_c2_c3()
gen_mc_exp <- rdata(model_mc_exp)
solver_mc_exp <- fit(model_mc_exp)

est_exp <- matrix(NA, B, m_mc_exp)
conv_exp <- logical(B)

for (b in seq_len(B)) {
  df_b <- gen_mc_exp(theta_mc_exp, n = n_mc, p = 0.3,
                     observe = observe_right_censor(tau = 3))
  tryCatch({
    fit_b <- solver_mc_exp(df_b, par = rep(1, m_mc_exp),
                           method = "Nelder-Mead")
    est_exp[b, ] <- fit_b$par
    conv_exp[b] <- fit_b$converged
  }, error = function(e) conv_exp[b] <<- FALSE)
}

# Homogeneous Weibull: 3 components
theta_mc_hom <- c(1.5, 100, 150, 200)
m_mc_hom <- 3
model_mc_hom <- wei_series_homogeneous_md_c1_c2_c3()
gen_mc_hom <- rdata(model_mc_hom)
solver_mc_hom <- fit(model_mc_hom)

beta_sys_mc <- wei_series_system_scale(theta_mc_hom[1], theta_mc_hom[-1])
tau_mc_hom <- qweibull(0.75, shape = theta_mc_hom[1], scale = beta_sys_mc)

est_hom <- matrix(NA, B, m_mc_hom + 1)
conv_hom <- logical(B)

for (b in seq_len(B)) {
  df_b <- gen_mc_hom(theta_mc_hom, n = n_mc, p = 0.3,
                     observe = observe_right_censor(tau = tau_mc_hom))
  tryCatch({
    fit_b <- solver_mc_hom(df_b, par = c(1, 120, 120, 120),
                           method = "L-BFGS-B",
                           lower = rep(1e-6, m_mc_hom + 1))
    est_hom[b, ] <- fit_b$par
    conv_hom[b] <- fit_b$converged
  }, error = function(e) conv_hom[b] <<- FALSE)
}

# Heterogeneous Weibull: 3 components
theta_mc_het <- c(0.8, 150, 1.5, 120, 2.0, 100)
m_mc_het <- 3
model_mc_het <- wei_series_md_c1_c2_c3()
gen_mc_het <- rdata(model_mc_het)
solver_mc_het <- fit(model_mc_het)

est_het <- matrix(NA, B, 2 * m_mc_het)
conv_het <- logical(B)

for (b in seq_len(B)) {
  df_b <- gen_mc_het(theta_mc_het, n = n_mc, tau = 200, p = 0.2,
                     observe = observe_right_censor(tau = 200))
  tryCatch({
    fit_b <- solver_mc_het(df_b, par = rep(c(1, 130), m_mc_het),
                           method = "L-BFGS-B",
                           lower = rep(1e-6, 2 * m_mc_het))
    est_het[b, ] <- fit_b$par
    conv_het[b] <- fit_b$converged
  }, error = function(e) conv_het[b] <<- FALSE)
}
```

```{r mc-save, eval=run_long, include=FALSE}
saveRDS(list(
  est_exp = est_exp, conv_exp = conv_exp, theta_mc_exp = theta_mc_exp,
  est_hom = est_hom, conv_hom = conv_hom, theta_mc_hom = theta_mc_hom,
  est_het = est_het, conv_het = conv_het, theta_mc_het = theta_mc_het,
  B = B, n_mc = n_mc, alpha = alpha,
  m_mc_exp = m_mc_exp, m_mc_hom = m_mc_hom, m_mc_het = m_mc_het
), "precomputed_framework.rds")
```

```{r mc-summary}
mc_summary <- function(estimates, converged, theta, par_names) {
  valid <- converged & !is.na(estimates[, 1])
  ev <- estimates[valid, , drop = FALSE]
  bias <- colMeans(ev) - theta
  rmse <- sqrt(bias^2 + apply(ev, 2, var))
  data.frame(
    Parameter = par_names,
    True = theta,
    Mean_Est = colMeans(ev),
    Bias = bias,
    RMSE = rmse,
    Rel_Bias_Pct = 100 * bias / theta,
    stringsAsFactors = FALSE
  )
}

# Exponential
cat("=== Exponential ===\n")
cat("Convergence:", mean(conv_exp), "\n\n")
exp_table <- mc_summary(est_exp, conv_exp, theta_mc_exp,
                        paste0("lambda_", 1:m_mc_exp))
knitr::kable(exp_table, digits = 4,
             caption = "Exponential MC (B=100, n=1000)",
             col.names = c("Parameter", "True", "Mean Est.", "Bias",
                          "RMSE", "Rel. Bias %"))
```

```{r mc-hom-table}
cat("\n=== Homogeneous Weibull ===\n")
cat("Convergence:", mean(conv_hom), "\n\n")
hom_table <- mc_summary(est_hom, conv_hom, theta_mc_hom,
                        c("k", paste0("beta_", 1:m_mc_hom)))
knitr::kable(hom_table, digits = 4,
             caption = "Homogeneous Weibull MC (B=100, n=1000)",
             col.names = c("Parameter", "True", "Mean Est.", "Bias",
                          "RMSE", "Rel. Bias %"))
```

```{r mc-het-table}
cat("\n=== Heterogeneous Weibull ===\n")
cat("Convergence:", mean(conv_het), "\n\n")
het_names <- paste0(rep(c("k", "beta"), m_mc_het), "_",
                    rep(1:m_mc_het, each = 2))
het_table <- mc_summary(est_het, conv_het, theta_mc_het, het_names)
knitr::kable(het_table, digits = 4,
             caption = "Heterogeneous Weibull MC (B=100, n=1000)",
             col.names = c("Parameter", "True", "Mean Est.", "Bias",
                          "RMSE", "Rel. Bias %"))
```


Practical Considerations {#practical}
======================================

**Starting values.** For the exponential model, starting at $\lambda_j = 1$
for all $j$ generally works. For Weibull models, start with shapes near 1
(the exponential case) and scales estimated from sample quantiles. An
exponential fit can provide good initial scale values via
$\hat\beta_j = 1/\hat\lambda_j$.

**Identifiability.** If two components are *always* masked together (i.e.,
they appear in every candidate set as a pair), their individual parameters
are not identifiable --- only their sum (exponential) or combined
contribution (Weibull) can be estimated. This is a structural limitation
of the observational design, not the estimation method.

**Computational cost.** The three models differ substantially in
computation time:

- **Exponential**: $O(n)$ per log-likelihood evaluation. All four
  observation types have closed-form loglik, score, and Hessian.
- **Homogeneous Weibull**: $O(n)$ for loglik (closed form for all types).
  Score uses `numDeriv::grad` for left/interval observations.
- **Heterogeneous Weibull**: $O(n_{\text{exact}} + n_{\text{right}})$ for
  closed-form terms, plus $O(n_{\text{left}} + n_{\text{interval}})$ calls
  to `stats::integrate` for numerical integration terms.

For datasets with many left- or interval-censored observations under the
heterogeneous model, each optimization iteration involves hundreds of
numerical integrations. Strategies to manage this include using the
homogeneous model as a warm start, or extending the study period to
convert more observations to exact/right-censored.

**Bootstrap confidence intervals.** As an alternative to asymptotic Wald
intervals (based on the observed Fisher information), parametric bootstrap
confidence intervals can be constructed by repeatedly simulating from the
fitted model and refitting. This is more computationally expensive but
avoids reliance on large-sample normality approximations, which can be
inaccurate when the sample size is small or the likelihood surface is
asymmetric.

```{r cleanup, include=FALSE}
options(old_opts)
```
