test_that("Bivariate Teissier functions work correctly", {
  xi1 <- 2.5
  xi2 <- 1.8
  delta <- 0.6
  z1 <- 0.4
  z2 <- 0.7
  
  # Joint PDF
  pdf_val <- dbivteissier(z1, z2, xi1, xi2, delta)
  expect_true(is.numeric(pdf_val) && pdf_val > 0)
  expect_equal(dbivteissier(z1, z2, xi1, xi2, delta, log = TRUE), log(pdf_val))
  expect_equal(dbivteissier(0, z2, xi1, xi2, delta), 0)
  expect_equal(dbivteissier(z1, 0, xi1, xi2, delta), 0)
  
  # Independence case delta = 0
  pdf_indep <- dbivteissier(z1, z2, xi1, xi2, 0)
  expect_equal(pdf_indep, dteissier(z1, xi1) * dteissier(z2, xi2))
  
  # Joint CDF and Survival
  cdf_val <- pbivteissier(z1, z2, xi1, xi2, delta)
  surv_val <- sbivteissier(z1, z2, xi1, xi2, delta)
  expect_true(cdf_val > 0 && cdf_val < 1)
  expect_true(surv_val > 0 && surv_val < 1)
  expect_equal(pbivteissier(0, z2, xi1, xi2, delta), 0)
  expect_equal(pbivteissier(z1, 0, xi1, xi2, delta), 0)
  expect_equal(pbivteissier(z1, z2, xi1, xi2, delta, lower.tail = FALSE), surv_val)
  
  # Random number generation
  set.seed(123)
  n <- 500
  sim <- rbivteissier(n, xi1, xi2, delta)
  expect_s3_class(sim, "data.frame")
  expect_equal(nrow(sim), n)
  expect_equal(ncol(sim), 2)
  expect_true(all(sim$z1 > 0))
  expect_true(all(sim$z2 > 0))
  
  # Missing values and parameter validation
  expect_true(is.na(dbivteissier(NA, z2, xi1, xi2, delta)))
  expect_true(is.na(pbivteissier(z1, NA, xi1, xi2, delta)))
  expect_warning(res_bad <- dbivteissier(z1, z2, -1, xi2, delta))
  expect_true(is.nan(res_bad))
  expect_warning(res_bad2 <- dbivteissier(z1, z2, xi1, xi2, 1.5))
  expect_true(is.nan(res_bad2))
})