test_that("Bayesian MCMC estimation produces valid loss function estimates and HPD intervals", {
  data("burr_data", package = "BivTrigBurr")

  set.seed(42)
  bayes_out <- bayes_sps_bxii(x = burr_data$x, y = burr_data$y,
                              n_iter = 300, burn_in = 50, thin = 1)

  expect_s3_class(bayes_out, "sps_bxii_bayes")
  expect_equal(length(bayes_out$estimates_self), 5)
  expect_equal(length(bayes_out$estimates_elf), 5)
  expect_equal(length(bayes_out$estimates_plf), 5)
  expect_equal(length(bayes_out$posterior_se), 5)
  expect_equal(nrow(bayes_out$hpd_interval), 5)
  expect_equal(ncol(bayes_out$hpd_interval), 2)
  expect_equal(nrow(bayes_out$heidel_test), 5)

  # Check that estimates are reasonable and finite
  expect_true(all(is.finite(bayes_out$estimates_self)))
  expect_true(all(is.finite(bayes_out$estimates_elf)))
  expect_true(all(is.finite(bayes_out$estimates_plf)))
  expect_true(all(bayes_out$hpd_aw >= 0))
  expect_true(all(bayes_out$hpd_interval[, 2] >= bayes_out$hpd_interval[, 1]))

  expect_output(print(bayes_out))
})