## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 4
)

# Heavy simulations are disabled by default for CRAN.
# Set run_long <- TRUE to reproduce the full comparison tables.
run_long <- FALSE

## ----load-package-------------------------------------------------------------
library(kofn)
library(flexhaz)

## ----weibull-hazard-plot, fig.cap = "Weibull hazard functions for different shape parameters."----
t_grid <- seq(0.01, 5, length.out = 200)
alphas <- c(0.5, 1.0, 1.5, 2.5)
beta <- 2.0

plot(NULL, xlim = c(0, 5), ylim = c(0, 3),
     xlab = "Time", ylab = "Hazard rate h(t)",
     main = "Weibull hazard: shape controls failure mechanism")

cols <- c("steelblue", "grey40", "darkorange", "firebrick")
for (i in seq_along(alphas)) {
  a <- alphas[i]
  h <- (a / beta) * (t_grid / beta)^(a - 1)
  lines(t_grid, h, col = cols[i], lwd = 2)
}
legend("topright",
       legend = paste0("alpha = ", alphas),
       col = cols, lwd = 2, bty = "n")

## ----direct-mle-demo----------------------------------------------------------
# Generate data from a Weibull parallel system
set.seed(42)
model_mle <- kofn(k = 2, m = 2, component = dfr_weibull(), method = "mle")
gen <- rdata(model_mle)

true_par <- c(1.5, 2.0,   # component 1: shape=1.5, scale=2.0
              2.0, 3.0)    # component 2: shape=2.0, scale=3.0
df <- gen(true_par, n = 50)

# Fit with direct MLE
fit_mle_fn <- fit(model_mle)
result_mle <- fit_mle_fn(df)

cat("Direct MLE estimates:\n")
cat(sprintf("  Component 1: shape = %.3f, scale = %.3f\n",
            result_mle$shapes[1], result_mle$scales[1]))
cat(sprintf("  Component 2: shape = %.3f, scale = %.3f\n",
            result_mle$shapes[2], result_mle$scales[2]))
cat(sprintf("  Converged: %s\n", result_mle$converged))

## ----truncated-moments-demo---------------------------------------------------
# E[T^alpha | T < t] for T ~ Weibull(alpha=1.5, beta=2.0)
alpha <- 1.5
beta <- 2.0
t_trunc <- 3.0

# Scalar interface
pow_moment <- trunc_pow_moment(k = alpha, t = t_trunc,
                               alpha = alpha, beta = beta)
cat(sprintf("E[T^%.1f | T < %.1f] = %.4f\n", alpha, t_trunc, pow_moment))
cat(sprintf("  (compare: t^alpha = %.4f)\n", t_trunc^alpha))

# Vectorized interface (used in EM inner loop)
v_max <- (t_trunc / beta)^alpha
pow_vec <- trunc_pow_moment_vec(k = alpha, v_max_vec = v_max,
                                alpha = alpha, beta = beta)
cat(sprintf("  Vectorized result: %.4f (same)\n", pow_vec))

# Truncated log-moment
log_moment <- trunc_log_moment_vec(v_max_vec = v_max,
                                   alpha = alpha, beta = beta)
cat(sprintf("E[log T | T < %.1f] = %.4f\n", t_trunc, log_moment))
cat(sprintf("  (compare: log t = %.4f)\n", log(t_trunc)))

## ----em-vs-mle----------------------------------------------------------------
set.seed(123)

# Create models
model_em  <- kofn(k = 2, m = 2, component = dfr_weibull(), method = "em")
model_mle <- kofn(k = 2, m = 2, component = dfr_weibull(), method = "mle")

# True parameters: shapes (1.5, 2.0), scales (2.0, 3.0)
true_par <- c(1.5, 2.0, 2.0, 3.0)

# Generate data
gen <- rdata(model_em)
df <- gen(true_par, n = 50)

# Fit with EM
fit_em_fn <- fit(model_em)
result_em <- fit_em_fn(df)

# Fit with direct MLE
fit_mle_fn <- fit(model_mle)
result_mle <- fit_mle_fn(df)

# Compare
cat("True parameters:\n")
cat(sprintf("  Component 1: shape = %.2f, scale = %.2f\n", 1.5, 2.0))
cat(sprintf("  Component 2: shape = %.2f, scale = %.2f\n", 2.0, 3.0))
cat("\nEM estimates:\n")
cat(sprintf("  Component 1: shape = %.3f, scale = %.3f\n",
            result_em$shapes[1], result_em$scales[1]))
cat(sprintf("  Component 2: shape = %.3f, scale = %.3f\n",
            result_em$shapes[2], result_em$scales[2]))
cat(sprintf("  Converged: %s, Iterations: %d\n",
            result_em$converged, result_em$iterations))
cat("\nDirect MLE estimates:\n")
cat(sprintf("  Component 1: shape = %.3f, scale = %.3f\n",
            result_mle$shapes[1], result_mle$scales[1]))
cat(sprintf("  Component 2: shape = %.3f, scale = %.3f\n",
            result_mle$shapes[2], result_mle$scales[2]))
cat(sprintf("  Converged: %s\n", result_mle$converged))

## ----shape-effect-sim, eval = run_long----------------------------------------
# # Full simulation: vary shape, compare EM vs MLE (takes several minutes).
# # Set run_long <- TRUE in setup chunk to reproduce.
# set.seed(2026)
# alphas <- c(0.5, 1.0, 1.5, 2.0, 3.0)
# R <- 200; n <- 300
# 
# for (a in alphas) {
#   true_par <- c(a, 2.0, a, 3.0)
#   mdl_em  <- kofn(k = 2, m = 2, component = dfr_weibull(), method = "em")
#   mdl_mle <- kofn(k = 2, m = 2, component = dfr_weibull(), method = "mle")
#   gen <- rdata(mdl_em)
#   f_em <- fit(mdl_em); f_mle <- fit(mdl_mle)
# 
#   em_sh <- mle_sh <- matrix(NA, R, 2)
#   for (r in seq_len(R)) {
#     d <- gen(true_par, n = n)
#     re <- tryCatch(f_em(d), error = function(e) NULL)
#     rm <- tryCatch(f_mle(d), error = function(e) NULL)
#     if (!is.null(re) && re$converged) em_sh[r, ] <- sort(re$shapes)
#     if (!is.null(rm) && rm$converged) mle_sh[r, ] <- sort(rm$shapes)
#   }
#   ok_em <- complete.cases(em_sh); ok_mle <- complete.cases(mle_sh)
#   cat(sprintf("alpha=%.1f: EM RMSE=%.3f, MLE RMSE=%.3f\n", a,
#     sqrt(mean((em_sh[ok_em,] - a)^2)), sqrt(mean((mle_sh[ok_mle,] - a)^2))))
# }

## ----mixed-shapes-demo--------------------------------------------------------
set.seed(456)

# Mixed DFR + IFR: shape (0.8, 1.5), scale (2.0, 3.0)
model_em <- kofn(k = 2, m = 2, component = dfr_weibull(), method = "em")
model_mle <- kofn(k = 2, m = 2, component = dfr_weibull(), method = "mle")

true_par <- c(0.8, 2.0, 1.5, 3.0)
gen <- rdata(model_em)
df <- gen(true_par, n = 50)

fit_em_fn  <- fit(model_em)
fit_mle_fn <- fit(model_mle)

result_em  <- fit_em_fn(df)
result_mle <- fit_mle_fn(df)

cat("Mixed DFR + IFR system: alpha = (0.8, 1.5), beta = (2.0, 3.0)\n\n")
cat("EM estimates:\n")
cat(sprintf("  Component 1: shape = %.3f (true 0.8), scale = %.3f (true 2.0)\n",
            result_em$shapes[1], result_em$scales[1]))
cat(sprintf("  Component 2: shape = %.3f (true 1.5), scale = %.3f (true 3.0)\n",
            result_em$shapes[2], result_em$scales[2]))
cat("\nDirect MLE estimates:\n")
cat(sprintf("  Component 1: shape = %.3f (true 0.8), scale = %.3f (true 2.0)\n",
            result_mle$shapes[1], result_mle$scales[1]))
cat(sprintf("  Component 2: shape = %.3f (true 1.5), scale = %.3f (true 3.0)\n",
            result_mle$shapes[2], result_mle$scales[2]))

## ----mixed-shapes-sim, eval = run_long----------------------------------------
# # Full mixed-shape simulation (set run_long <- TRUE to reproduce)
# set.seed(2026)
# scenarios <- list(
#   mild  = list(par = c(0.8, 2.0, 1.5, 3.0), label = "Mild (0.8, 1.5)"),
#   strong = list(par = c(0.5, 2.0, 3.0, 3.0), label = "Strong (0.5, 3.0)"),
#   three = list(par = c(0.7, 1.5, 1.5, 2.0, 2.5, 3.0),
#                label = "Three-comp (0.7, 1.5, 2.5)")
# )
# R <- 200; n <- 300
# for (sc in scenarios) {
#   m <- length(sc$par) / 2
#   mdl <- kofn(k = m, m = m, component = dfr_weibull(), method = "em")
#   gen <- rdata(mdl); f <- fit(mdl)
#   true_sh <- sc$par[seq(1, length(sc$par), by = 2)]
#   errs <- numeric(0)
#   for (r in seq_len(R)) {
#     res <- tryCatch(f(gen(sc$par, n = n)), error = function(e) NULL)
#     if (!is.null(res) && res$converged)
#       errs <- c(errs, (sort(res$shapes) - sort(true_sh))^2)
#   }
#   cat(sprintf("%s: EM shape RMSE = %.2f\n", sc$label, sqrt(mean(errs))))
# }

## ----stats-integration--------------------------------------------------------
set.seed(789)
model_em <- kofn(k = 2, m = 2, component = dfr_weibull(), method = "em")
gen <- rdata(model_em)
true_par <- c(1.5, 2.0, 2.0, 3.0)
df <- gen(true_par, n = 50)

fit_em_fn <- fit(model_em)
result <- fit_em_fn(df)

# Coefficient extraction
cat("coef():\n")
print(coef(result))

# Variance-covariance matrix (from observed Fisher information)
cat("\nvcov():\n")
print(round(vcov(result), 5))

# Log-likelihood with degrees of freedom
cat("\nlogLik():\n")
print(logLik(result))

# Model selection criteria
cat("\nAIC:", AIC(result), "\n")
cat("BIC:", BIC(result), "\n")

# Confidence intervals (Wald-type, based on Hessian)
cat("\nconfint() (95% Wald intervals):\n")
print(confint(result))

# Full summary
cat("\nsummary():\n")
print(summary(result))