## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(AlphaPowerHazard)


## ----dist_funcs---------------------------------------------------------------
x_vals <- c(0.5, 1.0, 1.5, 2.0)
alpha <- 1.5
beta <- 0.8

# PDF, CDF, Survival, Hazard, Cumulative Hazard
dalpha_power(x_vals, alpha, beta)
palpha_power(x_vals, alpha, beta)
salpha_power(x_vals, alpha, beta)
halpha_power(x_vals, alpha, beta)
Halpha_power(x_vals, alpha, beta)

# Quantiles
qalpha_power(c(0.25, 0.50, 0.75), alpha, beta)


## ----moments------------------------------------------------------------------
mom <- moments_alpha_power(k = 4, alpha = 1.5, beta = 1.1)
cat("Mean:", mom$mean, "\n")
cat("Variance:", mom$variance, "\n")
cat("Skewness:", mom$skewness, "\n")
cat("Kurtosis:", mom$kurtosis, "\n")

qs <- quantile_stats_alpha_power(alpha = 1.5, beta = 1.1)
cat("Bowley Skewness:", qs$bowley_skewness, "\n")
cat("Moors Kurtosis:", qs$moors_kurtosis, "\n")


## ----hrf_min------------------------------------------------------------------
# Unique minimum location for beta < 1 (bathtub shape)
hrf_m <- hrf_min_alpha_power(alpha = 1.5, beta = 0.8)
cat("HRF Minimum x:", hrf_m$xmin, "h(x):", hrf_m$hmin, "\n")

# Order statistics
order_stats_alpha_power(c(0.5, 1.0), r = 2, n = 5, alpha = 1.5, beta = 1.1, type = "pdf")


## ----estimation_methods-------------------------------------------------------
set.seed(123)
sim_x <- ralpha_power(50, alpha = 1.5, beta = 0.8)

fit_mle  <- fit_alpha_power_base(sim_x, method = "mle")
fit_lse  <- fit_alpha_power_base(sim_x, method = "lse")
fit_wlse <- fit_alpha_power_base(sim_x, method = "wlse")
fit_mpse <- fit_alpha_power_base(sim_x, method = "mpse")
fit_cme  <- fit_alpha_power_base(sim_x, method = "cme")

fit_mle


## ----regression_models--------------------------------------------------------
set.seed(42)
n <- 60
df <- data.frame(
  time   = ralpha_power(n, alpha = 1.5, beta = 0.8),
  status = sample(c(1, 1, 1, 0), n, replace = TRUE),
  age    = rnorm(n),
  bmi    = rnorm(n)
)

fit_m1 <- alpha_power_hazard(survival::Surv(time, status) ~ age + bmi, data = df, model = "M1")
summary(fit_m1)


## ----diagnostics--------------------------------------------------------------
res <- residuals_alpha_power(fit_m1)
head(res$cox_snell)
head(res$martingale)

# Predictions
pred_s <- predict_alpha_power(fit_m1, type = "survival")
head(pred_s$survival[, 1:4])