test_that("Bayesian MCMC estimation works accurately", {
  set.seed(123)
  sim <- rbitl(30, xi1 = 1.2, xi2 = 1.5, delta = 0.3)
  bayes_fit <- bayes_bitl(sim[, "x"], sim[, "y"], prior = "flat",
                          n_iter = 500, burn_in = 100, thin = 1)

  expect_s3_class(bayes_fit, "bayes_bitl")
  expect_equal(nrow(bayes_fit$estimates), 3)
  expect_named(bayes_fit$estimates, c("Parameter", "SELF_Est", "SELF_Risk",
                                     "MQSELF_Est", "MQSELF_Risk",
                                     "PLF_Est", "PLF_Risk"))
  expect_equal(dim(bayes_fit$etci), c(3, 2))
  expect_equal(dim(bayes_fit$hpd), c(3, 2))
  expect_equal(length(bayes_fit$heidelberger_pvalue), 3)
  expect_true(all(bayes_fit$heidelberger_pvalue >= 0 & bayes_fit$heidelberger_pvalue <= 1))
})