test_that("Distribution functions work properly and return valid values", {
  a1 <- 2; b1 <- 3; a2 <- 2; b2 <- 3; d <- 0.5

  # Joint PDF
  dens <- dSPSBXII(1, 1, a1, b1, a2, b2, d)
  expect_true(is.numeric(dens))
  expect_true(dens > 0)

  log_dens <- dSPSBXII(1, 1, a1, b1, a2, b2, d, log = TRUE)
  expect_equal(log_dens, log(dens))

  # Non-positive values return 0 / -Inf
  expect_equal(dSPSBXII(-1, 1, a1, b1, a2, b2, d), 0)
  expect_equal(dSPSBXII(-1, 1, a1, b1, a2, b2, d, log = TRUE), -Inf)

  # Joint CDF
  prob <- pSPSBXII(1, 1, a1, b1, a2, b2, d)
  expect_true(prob >= 0 && prob <= 1)

  # Joint Survival Function
  surv <- sSPSBXII(1, 1, a1, b1, a2, b2, d)
  expect_true(surv >= 0 && surv <= 1)

  # Joint HRF
  hrf <- hSPSBXII(1, 1, a1, b1, a2, b2, d)
  expect_true(hrf > 0)
  expect_equal(hrf, dens / surv)

  # Conditional densities
  cyx <- cond_pdf_y_given_x(1, 1, a1, b1, a2, b2, d)
  cxy <- cond_pdf_x_given_y(1, 1, a1, b1, a2, b2, d)
  expect_true(cyx > 0)
  expect_true(cxy > 0)

  # Random variate generation
  set.seed(42)
  samp <- rSPSBXII(100, a1, b1, a2, b2, d)
  expect_equal(nrow(samp), 100)
  expect_equal(ncol(samp), 2)
  expect_true(all(samp$x > 0))
  expect_true(all(samp$y > 0))
  expect_false(any(is.na(samp)))
  expect_false(any(is.nan(as.matrix(samp))))
})

test_that("Parameter boundary validation stops on invalid inputs", {
  expect_error(dSPSBXII(1, 1, -1, 2, 2, 2, 0.5))
  expect_error(dSPSBXII(1, 1, 1, -2, 2, 2, 0.5))
  expect_error(dSPSBXII(1, 1, 1, 2, 2, 2, 1.5))
  expect_error(pSPSBXII(1, 1, 1, 2, 2, 2, -1.5))
  expect_error(rSPSBXII(-5, 1, 2, 2, 2, 0.5))
})