test_that("MLE estimation and S3 methods work and match paper results", {
  data(burr_data)
  expect_equal(nrow(burr_data), 50)
  expect_equal(ncol(burr_data), 2)
  
  # Fit MLE
  fit <- mlebivteissier(burr_data)
  expect_s3_class(fit, "mlebivteissier")
  
  # Check MLE parameter values match Table 2 of paper (5.9549, 6.3414, 0.4950)
  expect_equal(as.numeric(fit$par["xi1"]), 5.9549, tolerance = 1e-3)
  expect_equal(as.numeric(fit$par["xi2"]), 6.3414, tolerance = 1e-3)
  expect_equal(as.numeric(fit$par["delta"]), 0.4950, tolerance = 1e-3)
  
  # Check Standard Errors match Table 2 (0.3206, 0.3414, 0.3404)
  expect_equal(as.numeric(fit$se["xi1"]), 0.3206, tolerance = 1e-2)
  expect_equal(as.numeric(fit$se["xi2"]), 0.3414, tolerance = 1e-2)
  expect_equal(as.numeric(fit$se["delta"]), 0.3404, tolerance = 1e-2)
  
  # Check Information Criteria match Table 4 (LogL = 114.7026, AIC = -223.4052, BIC = -217.6691, AICc = -222.8835)
  expect_equal(fit$loglik, 114.7026, tolerance = 1e-2)
  expect_equal(fit$aic, -223.4052, tolerance = 1e-2)
  expect_equal(fit$bic, -217.6691, tolerance = 1e-2)
  expect_equal(fit$aicc, -222.8835, tolerance = 1e-2)
  
  # S3 methods
  expect_equal(length(coef(fit)), 3)
  expect_s3_class(logLik(fit), "logLik")
  expect_equal(dim(vcov(fit)), c(3, 3))
  ci <- confint(fit)
  expect_equal(dim(ci), c(3, 2))
  
  # Output printing methods
  expect_output(print(fit))
  expect_output(summary(fit))
})

test_that("Bayesian estimation and MCMC diagnostics work", {
  data(burr_data)
  set.seed(42)
  bayes_fit <- bayesbivteissier(burr_data, prior = "vague", n.iter = 1000, burn.in = 200, thin = 1)
  expect_s3_class(bayes_fit, "bayesbivteissier")
  
  # Loss function estimates
  expect_true(is.list(bayes_fit$estimates))
  expect_true(all(c("SELF", "MQSELF", "PLF") %in% names(bayes_fit$estimates)))
  expect_equal(dim(bayes_fit$hpd_ci), c(3, 2))
  expect_equal(dim(bayes_fit$equal_tail_ci), c(3, 2))
  
  # Convergence diagnostic
  expect_true(is.data.frame(bayes_fit$heidel_diag))
  expect_equal(nrow(bayes_fit$heidel_diag), 3)
  
  # S3 methods
  expect_equal(length(coef(bayes_fit, loss = "SELF")), 3)
  expect_equal(dim(confint(bayes_fit, type = "HPD")), c(3, 2))
  expect_output(print(bayes_fit))
  expect_output(summary(bayes_fit))
})