test_that("MLE fits correctly on burr_data and matches paper results", {
  data("burr_data", package = "BivTrigBurr")

  fit <- fit_sps_bxii_mle(x = burr_data$x, y = burr_data$y)

  # Check class and structure
  expect_s3_class(fit, "sps_bxii_mle")
  expect_equal(length(fit$estimates), 5)
  expect_equal(names(fit$estimates), c("alpha1", "beta1", "alpha2", "beta2", "delta"))

  # Check paper values (Table 10: alpha1=2.1479, beta1=39.1249, alpha2=2.0277, beta2=37.4637, delta=0.4083)
  expect_equal(unname(round(fit$estimates["alpha1"], 2)), 2.15, tolerance = 0.05)
  expect_equal(unname(round(fit$estimates["beta1"], 0)), 39, tolerance = 2)
  expect_equal(unname(round(fit$estimates["alpha2"], 2)), 2.03, tolerance = 0.05)
  expect_equal(unname(round(fit$estimates["beta2"], 0)), 37, tolerance = 2)
  expect_equal(unname(round(fit$estimates["delta"], 2)), 0.41, tolerance = 0.05)

  # Check LogLik, AIC, BIC (Table 11: -LogL=-113.931, AIC=-217.862, BIC=-208.302)
  expect_equal(round(fit$loglik, 1), 113.9, tolerance = 0.5)
  expect_equal(round(fit$aic, 1), -217.9, tolerance = 1.0)
  expect_equal(round(fit$bic, 1), -208.3, tolerance = 1.0)

  # Check S3 methods
  expect_output(print(fit))
  sm <- summary(fit)
  expect_s3_class(sm, "summary.sps_bxii_mle")
  expect_output(print(sm))

  ci <- confint(fit)
  expect_equal(nrow(ci), 5)
  expect_equal(ncol(ci), 2)
  expect_true(all(ci[, 2] >= ci[, 1]))

  expect_equal(unname(AIC(fit)), fit$aic)
  expect_equal(unname(BIC(fit)), fit$bic)
  expect_equal(as.numeric(logLik(fit)), fit$loglik)
})

test_that("Score gradient matches numerical derivative", {
  x_sub <- c(0.1, 0.2, 0.3)
  y_sub <- c(0.15, 0.25, 0.35)
  par_test <- c(2.0, 3.0, 2.0, 3.0, 0.4)

  ana_grad <- score_sps_bxii(par_test, x_sub, y_sub)

  eps <- 1e-6
  num_grad <- numeric(5)
  for (i in 1:5) {
    p_up <- par_test
    p_dn <- par_test
    p_up[i] <- p_up[i] + eps
    p_dn[i] <- p_dn[i] - eps
    num_grad[i] <- (loglik_sps_bxii(p_up, x_sub, y_sub) - loglik_sps_bxii(p_dn, x_sub, y_sub)) / (2 * eps)
  }

  expect_equal(unname(ana_grad), num_grad, tolerance = 1e-4)
})