test_that("Bivariate BITL functions work accurately", {
  xi1 <- 1.5
  xi2 <- 1.2
  delta <- 0.5

  x <- c(0.2, 0.5, 1.0)
  y <- c(0.3, 0.6, 1.2)

  # Joint PDF
  pdf_val <- dbitl(x, y, xi1 = xi1, xi2 = xi2, delta = delta)
  expect_true(all(pdf_val > 0))
  expect_equal(length(pdf_val), length(x))

  # Joint CDF and Survival
  cdf_val <- pbitl(x, y, xi1 = xi1, xi2 = xi2, delta = delta)
  surv_val <- sbitl(x, y, xi1 = xi1, xi2 = xi2, delta = delta)
  expect_true(all(cdf_val >= 0 & cdf_val <= 1))
  expect_true(all(surv_val >= 0 & surv_val <= 1))

  # Joint Hazard & Reversed Hazard
  h_val <- hbitl(x, y, xi1 = xi1, xi2 = xi2, delta = delta)
  rh_val <- rhbitl(x, y, xi1 = xi1, xi2 = xi2, delta = delta)
  expect_true(all(h_val > 0))
  expect_true(all(rh_val > 0))

  # Hazard gradient
  hgrad <- hgradbitl(x, y, xi1 = xi1, xi2 = xi2, delta = delta)
  expect_s3_class(hgrad, "data.frame")
  expect_named(hgrad, c("eta1", "eta2"))
  expect_true(all(hgrad > 0))
  expect_true(all(hgrad > 0))

  # Conditional distributions
  cond_pdf <- cond_bitl(x, y, xi1 = xi1, xi2 = xi2, delta = delta, given = "X", type = "pdf")
  cond_cdf <- cond_bitl(x, y, xi1 = xi1, xi2 = xi2, delta = delta, given = "X", type = "cdf")
  expect_true(all(cond_pdf > 0))
  expect_true(all(cond_cdf >= 0 & cond_cdf <= 1))

  # Random bivariate generation
  set.seed(123)
  sim <- rbitl(n = 50, xi1 = xi1, xi2 = xi2, delta = delta)
  expect_equal(dim(sim), c(50, 2))
  expect_true(all(sim[, "x"] > 0))
  expect_true(all(sim[, "y"] > 0))

  # Input validations
  expect_error(dbitl(x, y, xi1 = -1, xi2 = 1, delta = 0))
  expect_error(dbitl(x, y, xi1 = 1, xi2 = 1, delta = 2))
})