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

## ----device-------------------------------------------------------------------
library(maskedhaz)

device <- dfr_series_md(components = list(
  battery    = dfr_weibull(shape = 3, scale = 10),
  electronics = dfr_exponential(0.02),
  lead_wire  = dfr_gompertz(a = 0.02, b = 0.3),
  sensor     = dfr_loglogistic(alpha = 8, beta = 2)
))
device

## ----layout-------------------------------------------------------------------
param_layout(device$series)

## ----simulate-----------------------------------------------------------------
true_par <- c(3, 10,        # battery: shape, scale
              0.02,         # electronics: rate
              0.02, 0.3,    # lead wire: a, b
              8, 2)         # sensor: alpha, beta

set.seed(2026)
rdata_fn <- rdata(device)
df <- rdata_fn(theta = true_par, n = 500, tau = 8, p = 0.3)

table(df$omega)

## ----fit----------------------------------------------------------------------
solver <- fit(device)

# Initialize reasonably far from the truth to exercise the optimizer
init <- c(2, 8, 0.03, 0.03, 0.2, 6, 1.5)
result <- solver(df, par = init)

## ----compare------------------------------------------------------------------
comparison <- data.frame(
  parameter = c("battery shape", "battery scale",
                "electronics rate",
                "lead wire a", "lead wire b",
                "sensor alpha", "sensor beta"),
  truth    = true_par,
  estimate = round(coef(result), 4),
  se       = round(sqrt(diag(vcov(result))), 4)
)
comparison

## ----ccp, fig.width = 7, fig.height = 4---------------------------------------
ccp_fn <- conditional_cause_probability(device)
t_grid <- seq(0.5, 15, length.out = 60)
probs <- ccp_fn(t_grid, par = true_par)

matplot(t_grid, probs, type = "l", lwd = 2,
        xlab = "system failure time (years)",
        ylab = "P(component j caused failure | T = t)",
        main = "Which component fails when?")
legend("topright",
       legend = c("battery", "electronics", "lead wire", "sensor"),
       col = 1:4, lty = 1:4, lwd = 2, bty = "n")

## ----custom-dist, eval = FALSE------------------------------------------------
# library(flexhaz)
# 
# bathtub <- dfr_dist(
#   rate      = function(t, par, ...) par[1]/sqrt(t) + par[2] + par[3]*t,
#   cum_haz_rate = function(t, par, ...) 2*par[1]*sqrt(t) + par[2]*t + par[3]*t^2/2,
#   n_par = 3L,
#   par = c(0.1, 0.05, 0.01),
#   par_lower = c(1e-6, 1e-6, 1e-6)
# )
# 
# # Plug it into a series system like any other dfr_dist
# mixed <- dfr_series_md(components = list(
#   bathtub,
#   dfr_exponential(0.02)
# ))