Heterogeneous Weibull Series Systems: Flexible Hazard Shapes

model <- wei_series_md_c1_c2_c3(
  shapes = theta[seq(1, 6, 2)],
  scales = theta[seq(2, 6, 2)]
)

h1 <- component_hazard(model, 1)
h2 <- component_hazard(model, 2)
h3 <- component_hazard(model, 3)

t_grid <- seq(0.5, 120, length.out = 300)

plot(t_grid, h1(t_grid, theta), type = "l", col = "steelblue", lwd = 2,
     ylim = c(0, 0.05),
     xlab = "Time", ylab = "Hazard rate h(t)",
     main = "Component hazard functions")
lines(t_grid, h2(t_grid, theta), col = "forestgreen", lwd = 2, lty = 2)
lines(t_grid, h3(t_grid, theta), col = "firebrick", lwd = 2, lty = 3)
legend("topright",
       legend = 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()

ccp_fn <- conditional_cause_probability(model)

probs <- ccp_fn(t_grid, theta)

plot(t_grid, probs[, 1], type = "l", col = "steelblue", lwd = 2,
     ylim = c(0, 1),
     xlab = "Time", ylab = "P(K = j | T = t)",
     main = "Conditional cause-of-failure probability")
lines(t_grid, probs[, 2], col = "forestgreen", lwd = 2, lty = 2)
lines(t_grid, probs[, 3], col = "firebrick", lwd = 2, lty = 3)
legend("right",
       legend = c("Comp 1: DFR", "Comp 2: CFR", "Comp 3: IFR"),
       col = c("steelblue", "forestgreen", "firebrick"),
       lty = 1:3, lwd = 2, cex = 0.85)
grid()

shapes <- theta[seq(1, 6, 2)]
scales <- theta[seq(2, 6, 2)]
cind <- c(TRUE, TRUE, FALSE)  # candidate set = {1, 2}

# Evaluate the integrand h_c(t) * S(t) at a few time points
t_eval <- c(10, 50, 100)
vals <- maskedcauses:::wei_series_integrand(
  t_eval, shapes, scales, cind
)
names(vals) <- paste0("t=", t_eval)
vals
#>     t=10     t=50    t=100 
#> 0.012504 0.004574 0.001120

t_fine <- seq(0.5, 150, length.out = 500)

# Full candidate set: all components
vals_full <- maskedcauses:::wei_series_integrand(
  t_fine, shapes, scales, cind = rep(TRUE, 3)
)

plot(t_fine, vals_full, type = "l", col = "black", lwd = 2,
     xlab = "Time", ylab = expression(h[c](t) %.% S(t)),
     main = "Weibull series integrand (full candidate set)")
grid()

set.seed(123)
gen <- rdata(model)
ll_fn <- loglik(model)

# Exact + right only
df_exact <- gen(theta, n = 100, tau = 120, p = 0.3,
                observe = observe_right_censor(tau = 120))

# Mixed: includes left and interval censoring
df_mixed <- gen(theta, n = 100, p = 0.3,
                observe = observe_mixture(
                  observe_right_censor(tau = 120),
                  observe_left_censor(tau = 80),
                  observe_periodic(delta = 20, tau = 120),
                  weights = c(0.5, 0.25, 0.25)
                ))

cat("Observation type counts (exact/right):\n")
#> Observation type counts (exact/right):
print(table(df_exact$omega))
#> 
#> exact right 
#>    97     3

cat("\nObservation type counts (mixed):\n")
#> 
#> Observation type counts (mixed):
print(table(df_mixed$omega))
#> 
#>    exact interval     left    right 
#>       47       31       13        9

t_exact <- system.time(replicate(5, ll_fn(df_exact, theta)))["elapsed"]
t_mixed <- system.time(replicate(5, ll_fn(df_mixed, theta)))["elapsed"]

cat(sprintf("\nLoglik time (exact/right only): %.4f s per eval\n", t_exact / 5))
#> 
#> Loglik time (exact/right only): 0.0008 s per eval
cat(sprintf("Loglik time (mixed censoring):  %.4f s per eval\n", t_mixed / 5))
#> Loglik time (mixed censoring):  0.0286 s per eval
cat(sprintf("Slowdown factor: %.1fx\n", t_mixed / t_exact))
#> Slowdown factor: 35.8x

theta_mc <- results$theta_mc
m_mc <- results$m_mc
alpha <- results$alpha

valid <- results$converged & !is.na(results$estimates[, 1])
est_valid <- results$estimates[valid, , drop = FALSE]

bias <- colMeans(est_valid) - theta_mc
variance <- apply(est_valid, 2, var)
mse <- bias^2 + variance
rmse <- sqrt(mse)

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

results_df <- data.frame(
  Parameter = par_names_mc,
  True = theta_mc,
  Mean_Est = colMeans(est_valid),
  Bias = bias,
  Variance = variance,
  MSE = mse,
  RMSE = rmse,
  Rel_Bias_Pct = 100 * bias / theta_mc
)

knitr::kable(results_df, digits = 4,
             caption = paste0("Monte Carlo results (n = ", results$n_mc,
                              ", B = ", sum(valid), " converged)"),
             col.names = c("Parameter", "True", "Mean Est.",
                          "Bias", "Variance", "MSE", "RMSE",
                          "Rel. Bias %"))

coverage <- numeric(2 * m_mc)
for (j in seq_len(2 * m_mc)) {
  valid_j <- valid & !is.na(results$ci_lower[, j]) &
    !is.na(results$ci_upper[, j])
  covered <- (results$ci_lower[valid_j, j] <= theta_mc[j]) &
    (theta_mc[j] <= results$ci_upper[valid_j, j])
  coverage[j] <- mean(covered)
}

mean_width <- colMeans(
  results$ci_upper[valid, ] - results$ci_lower[valid, ], na.rm = TRUE
)

coverage_df <- data.frame(
  Parameter = par_names_mc,
  True = theta_mc,
  Coverage = coverage,
  Nominal = 1 - alpha,
  Mean_Width = mean_width
)

knitr::kable(coverage_df, digits = 4,
             caption = paste0(100 * (1 - alpha),
                              "% CI coverage (nominal = ",
                              1 - alpha, ")"),
             col.names = c("Parameter", "True", "Coverage",
                          "Nominal", "Mean Width"))

oldpar <- par(mfrow = c(2, 3), mar = c(4, 4, 2, 1))
on.exit(par(oldpar))
for (j in seq_len(2 * m_mc)) {
  hist(est_valid[, j], breaks = 20, probability = TRUE,
       main = par_names_mc[j],
       xlab = expression(hat(theta)),
       col = "lightblue", border = "white")
  abline(v = theta_mc[j], col = "red", lwd = 2, lty = 2)
  abline(v = mean(est_valid[, j]), col = "blue", lwd = 2)
  legend("topright", legend = c("True", "Mean"),
         col = c("red", "blue"), lty = c(2, 1), lwd = 2, cex = 0.7)
}