test_that("Distribution functions (d, p, q, r, h, H, s) execute cleanly without NAs, NaNs or warnings", {
  alpha <- 1.5
  beta  <- 0.8
  x_vals <- c(0.1, 0.5, 1.0, 2.0, 5.0)
  
  # Hazard & Cumulative hazard
  h <- halpha_power(x_vals, alpha, beta)
  H <- Halpha_power(x_vals, alpha, beta)
  expect_true(all(is.finite(h) & h > 0))
  expect_true(all(is.finite(H) & H >= 0))
  
  # Density & Survival
  d <- dalpha_power(x_vals, alpha, beta)
  s <- salpha_power(x_vals, alpha, beta)
  expect_true(all(is.finite(d) & d > 0))
  expect_true(all(is.finite(s) & s >= 0 & s <= 1))
  
  # CDF & Quantile
  p <- palpha_power(x_vals, alpha, beta)
  expect_true(all(is.finite(p) & p >= 0 & p <= 1))
  
  q <- qalpha_power(p, alpha, beta)
  expect_equal(q, x_vals, tolerance = 1e-4)
  
  # Random variates
  set.seed(42)
  r <- ralpha_power(50, alpha, beta)
  expect_equal(length(r), 50)
  expect_true(all(is.finite(r) & r > 0))
})

test_that("Moments, quantile stats, HRF min, order stats, and stochastic ordering run cleanly", {
  alpha <- 1.5
  beta  <- 0.8
  
  # Moments
  mom <- moments_alpha_power(k = 4, alpha = alpha, beta = beta)
  expect_true(is.numeric(mom$mean) && is.finite(mom$mean) && mom$mean > 0)
  expect_true(is.numeric(mom$variance) && is.finite(mom$variance) && mom$variance > 0)
  expect_true(is.numeric(mom$skewness) && is.finite(mom$skewness))
  expect_true(is.numeric(mom$kurtosis) && is.finite(mom$kurtosis))
  
  # Quantile stats
  qs <- quantile_stats_alpha_power(alpha = alpha, beta = beta)
  expect_true(all(is.finite(qs$quantiles)))
  expect_true(is.numeric(qs$bowley_skewness) && is.finite(qs$bowley_skewness))
  expect_true(is.numeric(qs$moors_kurtosis) && is.finite(qs$moors_kurtosis))
  
  # HRF minimum location
  hrf_m <- hrf_min_alpha_power(alpha = alpha, beta = beta)
  expect_true(is.numeric(hrf_m$xmin) && is.finite(hrf_m$xmin) && hrf_m$xmin > 0)
  expect_true(is.numeric(hrf_m$hmin) && is.finite(hrf_m$hmin) && hrf_m$hmin > 0)
  
  # Order stats
  os_pdf <- order_stats_alpha_power(c(0.5, 1.0), r = 2, n = 5, alpha = alpha, beta = beta, type = "pdf")
  os_cdf <- order_stats_alpha_power(c(0.5, 1.0), r = 2, n = 5, alpha = alpha, beta = beta, type = "cdf")
  expect_true(all(is.finite(os_pdf) & os_pdf >= 0))
  expect_true(all(is.finite(os_cdf) & os_cdf >= 0 & os_cdf <= 1))
  
  # Stochastic ordering
  so <- stochastic_ordering_alpha_power(alpha1 = 1.2, beta1 = 0.8, alpha2 = 1.5, beta2 = 1.1)
  expect_true(is.logical(so$is_ordered))
})

test_that("5 Classical estimation methods (MLE, LSE, WLSE, MPSE, CME) work without errors", {
  set.seed(123)
  dat <- ralpha_power(40, alpha = 1.5, beta = 0.8)
  
  methods <- c("mle", "lse", "wlse", "mpse", "cme")
  for (m in methods) {
    fit <- fit_alpha_power_base(dat, method = m)
    expect_equal(fit$convergence, 0)
    expect_true(all(is.finite(fit$par) & fit$par > c(1, 0)))
  }
})