## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
oldopts <- options(digits = 4)

## ----setup--------------------------------------------------------------------
library(algebraic.mle)

## ----mle-direct---------------------------------------------------------------
# Exponential rate from a published study: lambda_hat = 0.5, se = 0.07, n = 200
fit_pub <- mle(
  theta.hat = c(lambda = 0.5),
  sigma = matrix(0.07^2),
  nobs = 200L
)
params(fit_pub)
se(fit_pub)

## ----mle-numerical------------------------------------------------------------
set.seed(42)
x <- rexp(80, rate = 2)

loglik <- function(rate) {
  if (rate <= 0) return(-Inf)
  sum(dexp(x, rate = rate, log = TRUE))
}

result <- optim(
  par = c(lambda = 1),
  fn = loglik,
  method = "Brent", lower = 0.01, upper = 20,
  hessian = TRUE,
  control = list(fnscale = -1)
)

fit_num <- mle_numerical(result, options = list(nobs = length(x)))
params(fit_num)
se(fit_num)

## ----mle-boot-----------------------------------------------------------------
set.seed(42)
y <- rexp(30, rate = 2)

boot_result <- boot::boot(
  data = y,
  statistic = function(d, i) c(lambda = 1 / mean(d[i])),
  R = 999
)

fit_boot <- mle_boot(boot_result)
params(fit_boot)
se(fit_boot)
confint(fit_boot, type = "perc")

## ----joint-example------------------------------------------------------------
# Lab A: component 1 failure rate
fit_A <- mle(
  theta.hat = c(lambda1 = 0.02),
  sigma = matrix(0.001^2),
  nobs = 150L
)

# Lab B: component 2 failure rate
fit_B <- mle(
  theta.hat = c(lambda2 = 0.05),
  sigma = matrix(0.003^2),
  nobs = 80L
)

# Joint MLE: block-diagonal covariance
fit_joint <- joint(fit_A, fit_B)
params(fit_joint)
vcov(fit_joint)

## ----joint-methods------------------------------------------------------------
confint(fit_joint)
se(fit_joint)

## ----joint-marginal-----------------------------------------------------------
# Recover component 2 parameters
fit_B_recovered <- marginal(fit_joint, 2)
params(fit_B_recovered)
se(fit_B_recovered)

## ----joint-three--------------------------------------------------------------
fit_C <- mle(
  theta.hat = c(lambda3 = 0.01),
  sigma = matrix(0.0005^2),
  nobs = 200L
)

fit_system <- joint(fit_A, fit_B, fit_C)
params(fit_system)

## ----combine-example----------------------------------------------------------
lab1 <- mle(theta.hat = c(lambda = 0.050), sigma = matrix(0.004^2), nobs = 50L)
lab2 <- mle(theta.hat = c(lambda = 0.048), sigma = matrix(0.002^2), nobs = 200L)
lab3 <- mle(theta.hat = c(lambda = 0.053), sigma = matrix(0.003^2), nobs = 100L)

combined <- combine(lab1, lab2, lab3)
params(combined)
se(combined)

## ----combine-precision--------------------------------------------------------
cat("Lab SEs:     ", se(lab1), se(lab2), se(lab3), "\n")
cat("Combined SE: ", se(combined), "\n")

## ----combine-equal------------------------------------------------------------
eq1 <- mle(theta.hat = c(mu = 10.0), sigma = matrix(1))
eq2 <- mle(theta.hat = c(mu = 12.0), sigma = matrix(1))
eq3 <- mle(theta.hat = c(mu = 11.0), sigma = matrix(1))

eq_combined <- combine(eq1, eq2, eq3)
params(eq_combined)  # (10 + 12 + 11) / 3 = 11

## ----rmap-univariate----------------------------------------------------------
fit_rate <- mle(
  theta.hat = c(lambda = 0.02),
  sigma = matrix(0.001^2),
  nobs = 150L
)

# g: rate -> mean lifetime
fit_lifetime <- rmap(fit_rate, function(theta) c(mean_life = 1 / theta[1]),
  method = "delta")
params(fit_lifetime)
se(fit_lifetime)
confint(fit_lifetime)

## ----rmap-system--------------------------------------------------------------
fit_rates <- joint(
  mle(theta.hat = c(lambda1 = 0.02), sigma = matrix(0.001^2), nobs = 150L),
  mle(theta.hat = c(lambda2 = 0.05), sigma = matrix(0.003^2), nobs = 80L)
)

# System reliability at t = 10 hours
t0 <- 10
R_mle <- rmap(fit_rates,
  function(theta) c(R_t = exp(-(theta[1] + theta[2]) * t0)),
  method = "delta"
)

params(R_mle)
se(R_mle)
confint(R_mle)

## ----as-dist-basic------------------------------------------------------------
fit1 <- mle(theta.hat = c(lambda = 0.02), sigma = matrix(0.001^2))
d1 <- as_dist(fit1)
d1

## ----as-dist-prob-------------------------------------------------------------
# P(lambda > 0.021) under the asymptotic distribution
1 - cdf(d1)(0.021)

## ----as-dist-boot-------------------------------------------------------------
d_boot <- as_dist(fit_boot)
d_boot

## ----pipeline-setup-----------------------------------------------------------
set.seed(123)

# --- Step 1: Independent MLEs from separate experiments ---

# Component 1: 200 observed lifetimes
x1 <- rexp(200, rate = 0.003)
loglik1 <- function(rate) {
  if (rate <= 0) return(-Inf)
  sum(dexp(x1, rate = rate, log = TRUE))
}
fit1 <- mle_numerical(
  optim(par = c(lambda1 = 0.001), fn = loglik1,
        method = "Brent", lower = 1e-6, upper = 1,
        hessian = TRUE, control = list(fnscale = -1)),
  options = list(nobs = length(x1))
)

# Component 2: 120 observed lifetimes
x2 <- rexp(120, rate = 0.008)
loglik2 <- function(rate) {
  if (rate <= 0) return(-Inf)
  sum(dexp(x2, rate = rate, log = TRUE))
}
fit2 <- mle_numerical(
  optim(par = c(lambda2 = 0.001), fn = loglik2,
        method = "Brent", lower = 1e-6, upper = 1,
        hessian = TRUE, control = list(fnscale = -1)),
  options = list(nobs = length(x2))
)

cat("Component 1 rate:", params(fit1), "  SE:", se(fit1), "\n")
cat("Component 2 rate:", params(fit2), "  SE:", se(fit2), "\n")

## ----pipeline-compose---------------------------------------------------------
# --- Step 2: Compose into joint parameter space ---
fit_system <- joint(fit1, fit2)
cat("Joint parameters:", params(fit_system), "\n")
cat("Joint vcov:\n")
vcov(fit_system)

## ----pipeline-transform-------------------------------------------------------
# --- Step 3: Transform to system reliability at t = 50 ---
mission_time <- 50
R_system <- rmap(fit_system,
  function(theta) c(R_sys = exp(-(theta[1] + theta[2]) * mission_time)),
  method = "delta"
)

cat("System reliability R(50):", params(R_system), "\n")
cat("SE of R(50):             ", se(R_system), "\n")

## ----pipeline-inference-------------------------------------------------------
# --- Step 4: Inference ---
cat("95% CI for R(50):\n")
confint(R_system)

## ----pipeline-dist------------------------------------------------------------
# --- Step 5: Bridge to distribution algebra ---
R_dist <- as_dist(R_system)
cat("Asymptotic distribution of R(50):", "\n")
R_dist

# Probability that system reliability exceeds 55%
cat("P(R(50) > 0.55):", 1 - cdf(R_dist)(0.55), "\n")

## ----include = FALSE----------------------------------------------------------
options(oldopts)