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

## ----hazard-sum---------------------------------------------------------------
library(serieshaz)

sys <- dfr_dist_series(list(
    dfr_weibull(shape = 2, scale = 100),
    dfr_exponential(0.01),
    dfr_gompertz(a = 0.001, b = 0.05)
))

h_sys <- hazard(sys)
h1 <- component_hazard(sys, 1)
h2 <- component_hazard(sys, 2)
h3 <- component_hazard(sys, 3)

t_vals <- c(10, 50, 100, 200)
for (t0 in t_vals) {
    lhs <- h_sys(t0)
    rhs <- h1(t0) + h2(t0) + h3(t0)
    cat(sprintf("t=%3d: h_sys=%.6f, sum h_j=%.6f, diff=%.2e\n",
                t0, lhs, rhs, abs(lhs - rhs)))
}

## ----surv-product-------------------------------------------------------------
S_sys <- surv(sys)
comp1 <- component(sys, 1)
comp2 <- component(sys, 2)
comp3 <- component(sys, 3)
S1 <- surv(comp1)
S2 <- surv(comp2)
S3 <- surv(comp3)

for (t0 in t_vals) {
    lhs <- S_sys(t0)
    rhs <- S1(t0) * S2(t0) * S3(t0)
    cat(sprintf("t=%3d: S_sys=%.6f, prod S_j=%.6f, diff=%.2e\n",
                t0, lhs, rhs, abs(lhs - rhs)))
}

## ----cumhaz-sum---------------------------------------------------------------
# The system's cumulative hazard is the sum of component cumulative hazards
H_sys <- cum_haz(sys)
H1 <- cum_haz(comp1)
H2 <- cum_haz(comp2)
H3 <- cum_haz(comp3)

for (t0 in t_vals) {
    lhs <- H_sys(t0)
    rhs <- H1(t0) + H2(t0) + H3(t0)
    cat(sprintf("t=%3d: H_sys=%.6f, sum H_j=%.6f, diff=%.2e\n",
                t0, lhs, rhs, abs(lhs - rhs)))
}

## ----density-formula----------------------------------------------------------
f_sys <- density(sys)
for (t0 in t_vals) {
    from_density <- f_sys(t0)
    from_hS <- h_sys(t0) * S_sys(t0)
    cat(sprintf("t=%3d: f(t)=%.8f, h(t)*S(t)=%.8f, diff=%.2e\n",
                t0, from_density, from_hS, abs(from_density - from_hS)))
}

## ----layout-demo--------------------------------------------------------------
sys <- dfr_dist_series(list(
    dfr_weibull(shape = 2, scale = 100),    # 2 params
    dfr_exponential(0.05),            # 1 param
    dfr_gompertz(a = 0.01, b = 0.1)         # 2 params
))

# Full parameter vector (5 elements)
params(sys)

# Layout maps global indices to components
param_layout(sys)

# Component 1 uses par[1:2], component 2 uses par[3], component 3 uses par[4:5]

## ----analytical-check---------------------------------------------------------
# All standard distributions provide analytical cumulative hazard
sys_analytical <- dfr_dist_series(list(
    dfr_weibull(shape = 2, scale = 100),
    dfr_exponential(0.05)
))
cat("Has analytical cum_haz:", !is.null(sys_analytical$cum_haz_rate), "\n")

# A custom dfr_dist without cum_haz_rate forces numerical fallback
custom <- dfr_dist(
    rate = function(t, par, ...) par[1] * t^(par[1] - 1),
    par = c(1.5)
)
sys_numerical <- dfr_dist_series(list(
    custom,
    dfr_exponential(0.05)
))
cat("Has analytical cum_haz:", !is.null(sys_numerical$cum_haz_rate), "\n")

## ----score-decomposed---------------------------------------------------------
sys <- dfr_dist_series(list(
    dfr_weibull(shape = 2, scale = 100),
    dfr_exponential(0.01),
    dfr_gompertz(a = 0.001, b = 0.05)
))

set.seed(42)
times <- sampler(sys)(200)
df <- data.frame(t = times, delta = rep(1L, 200))

par0 <- params(sys)
sc_fn <- score(sys)
ll_fn <- loglik(sys)

# Decomposed score (what serieshaz computes internally)
score_val <- sc_fn(df, par = par0)

# Independent finite-difference verification
eps <- 1e-6
fd_score <- numeric(length(par0))
for (k in seq_along(par0)) {
    par_p <- par_m <- par0
    par_p[k] <- par0[k] + eps
    par_m[k] <- par0[k] - eps
    fd_score[k] <- (ll_fn(df, par = par_p) - ll_fn(df, par = par_m)) / (2 * eps)
}

result <- data.frame(
    Parameter = c("shape", "scale", "rate", "a", "b"),
    Decomposed = round(score_val, 6),
    FiniteDiff = round(fd_score, 6),
    AbsDiff = formatC(abs(score_val - fd_score), format = "e", digits = 1)
)
print(result, row.names = FALSE)

## ----score-at-mle-------------------------------------------------------------
sys2 <- dfr_dist_series(list(
    dfr_weibull(shape = 2, scale = 100),
    dfr_exponential(0.01)
))

set.seed(42)
df2 <- data.frame(t = sampler(sys2)(300), delta = rep(1L, 300))

sc2 <- score(sys2)
cat("Score at true params:", round(sc2(df2), 4), "\n")

solver <- fit(sys2)
result <- suppressWarnings(solver(df2, par = c(1.5, 80, 0.02)))

cat("MLE:                ", round(coef(result), 4), "\n")
cat("Score at MLE:       ", round(sc2(df2, par = coef(result)), 6), "\n")

## ----hessian-demo-------------------------------------------------------------
H_fn <- hess_loglik(sys2)
hess_val <- H_fn(df2, par = coef(result))
evals <- round(eigen(hess_val, symmetric = TRUE)$values, 2)
cat("Hessian eigenvalues at MLE:", evals, "\n")
cat("Hessian is symmetric:", all.equal(hess_val, t(hess_val)), "\n")

## ----exp-special--------------------------------------------------------------
rates <- c(0.1, 0.2, 0.3)
sys <- dfr_dist_series(lapply(rates, dfr_exponential))
equiv <- dfr_exponential(sum(rates))

h_sys <- hazard(sys)
h_eq  <- hazard(equiv)
S_sys <- surv(sys)
S_eq  <- surv(equiv)

for (t0 in c(1, 5, 10, 50)) {
    cat(sprintf("t=%2d: h_sys=%.4f h_eq=%.4f | S_sys=%.6f S_eq=%.6f\n",
                t0, h_sys(t0), h_eq(t0), S_sys(t0), S_eq(t0)))
}

## ----weibull-identical--------------------------------------------------------
m <- 3
shape <- 2
scale <- 100

sys <- dfr_dist_series(replicate(
    m, dfr_weibull(shape = shape, scale = scale), simplify = FALSE
))
equiv <- dfr_weibull(shape = shape, scale = scale / m^(1/shape))

S_sys <- surv(sys)
S_eq  <- surv(equiv)

for (t0 in c(10, 30, 50)) {
    cat(sprintf("t=%2d: S_sys=%.6f S_equiv=%.6f diff=%.2e\n",
                t0, S_sys(t0), S_eq(t0), abs(S_sys(t0) - S_eq(t0))))
}

## ----single-component---------------------------------------------------------
single <- dfr_dist_series(list(dfr_weibull(shape = 2, scale = 100)))
direct <- dfr_weibull(shape = 2, scale = 100)

h1 <- hazard(single)
h2 <- hazard(direct)
S1 <- surv(single)
S2 <- surv(direct)

for (t0 in c(10, 50, 100)) {
    cat(sprintf("t=%3d: h=%.6f/%.6f S=%.6f/%.6f\n",
                t0, h1(t0), h2(t0), S1(t0), S2(t0)))
}