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

## ----example-setup------------------------------------------------------------
library(kofn)
library(flexhaz)
set.seed(42)

model <- kofn(k = 1, m = 3, component = dfr_exponential())
gen   <- rdata(model)
df    <- gen(theta = c(0.3, 0.8, 1.5), n = 100)
head(df, 3)

## ----model-class--------------------------------------------------------------
class(model)

## ----likelihood-generics------------------------------------------------------
ll <- loglik(model)                  # ll: function(df, par) -> numeric
sc <- score(model)                   # sc: function(df, par) -> gradient
H  <- hess_loglik(model)             # H:  function(df, par) -> Hessian

par0 <- c(0.3, 0.8, 1.5)
ll(df, par0)
sc(df, par0)

## ----rdata-assumptions--------------------------------------------------------
assumptions(model)[1:3]              # first 3 assumptions

## ----fit-default--------------------------------------------------------------
fit_fn <- fit(model)
res <- fit_fn(df, n_starts = 1L)
round(coef(res), 3)

## ----compose-solver-----------------------------------------------------------
library(compositional.mle)

# Build the problem yourself to use compositional.mle directly
prob <- mle_problem(
  loglike    = function(par) ll(df, par),
  constraint = mle_constraint(support = function(par) all(par > 0))
)

# Sequential: grid search as a warm start, then L-BFGS-B for polish
strategy <- grid_search(lower = rep(0.05, 3),
                         upper = rep(2.0, 3), n = 3) %>>%
            lbfgsb(lower = rep(1e-10, 3))
res_chain <- strategy(prob, theta0 = rep(0.5, 3))
round(coef(res_chain), 3)

## ----parallel-race------------------------------------------------------------
race <- lbfgsb(lower = rep(1e-10, 3)) %|% nelder_mead()
res_race <- race(prob, theta0 = rep(0.5, 3))
round(coef(res_race), 3)

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

round(coef(res), 3)                  # point estimates
round(se(res), 3)                    # standard errors
round(confint(res), 3)[1:3, ]        # 95% CIs, first three params
logLik(res)                          # log-likelihood
AIC(res); BIC(res)                   # information criteria
nobs(res); nparams(res)              # bookkeeping

## ----mle-advanced-------------------------------------------------------------
round(bias(res), 4)                  # O(1/n) bias (zero for MLE at truth)
round(observed_fim(res), 1)[1:3, 1:3]  # observed Fisher information

## ----combine-estimates--------------------------------------------------------
df2  <- gen(theta = c(0.3, 0.8, 1.5), n = 100)
res1 <- fit_fn(df,  n_starts = 1L)
res2 <- fit_fn(df2, n_starts = 1L)

pooled <- combine(res1, res2)        # inverse-variance weighted
round(coef(pooled), 3)
round(se(pooled), 3)
nobs(pooled)                         # 200 = 100 + 100

## ----as-dist------------------------------------------------------------------
d <- as_dist(res)
class(d)

## ----rmap-delta---------------------------------------------------------------
# g(lambda) = 1/sum(lambda): mean time to first failure (series system)
g <- function(lam) 1 / sum(lam)
res_mttff <- rmap(res, g)
round(coef(res_mttff), 3)            # point estimate
round(se(res_mttff), 4)              # delta-method SE