## ----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)

## ----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")

## ----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)

## ----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)
)

## ----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))

## ----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")

## ----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")

## ----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))

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

## ----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 %"))

## ----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")

## ----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?"))

## ----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()

## ----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))

## ----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))

## ----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 %"))

## ----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()

## ----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()

## ----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))

## ----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 %"))

## ----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))

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

## ----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)
# }

## ----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")

## ----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 %"))

## ----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 %"))

## ----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 %"))

## ----cleanup, include=FALSE---------------------------------------------------
options(old_opts)