test_that("Reliability and hazard functions work correctly", {
  xi1 <- 2.0
  xi2 <- 1.5
  delta <- 0.5
  z1 <- 0.6
  z2 <- 0.8
  
  # Joint Hazard
  h_val <- hbivteissier(z1, z2, xi1, xi2, delta)
  f_val <- dbivteissier(z1, z2, xi1, xi2, delta)
  s_val <- sbivteissier(z1, z2, xi1, xi2, delta)
  expect_equal(h_val, f_val / s_val)
  expect_true(h_val > 0)
  
  # Conditional Hazards
  ch1 <- condhbivteissier(z1, z2, xi1, xi2, delta, given = "z2")
  ch2 <- condhbivteissier(z1, z2, xi1, xi2, delta, given = "z1")
  expect_true(ch1 > 0)
  expect_true(ch2 > 0)
  
  # Joint Reversed Hazard
  rhf_val <- rhfbivteissier(z1, z2, xi1, xi2, delta)
  p_val <- pbivteissier(z1, z2, xi1, xi2, delta)
  expect_equal(rhf_val, f_val / p_val)
  expect_true(rhf_val > 0)
  
  # Hazard Gradient
  grad <- hazgradbivteissier(z1, z2, xi1, xi2, delta)
  expect_s3_class(grad, "data.frame")
  expect_equal(nrow(grad), 1)
  expect_true(grad$eta1 > 0)
  expect_true(grad$eta2 > 0)
  
  # When delta = 0, hazard gradient equals marginal hazards
  grad0 <- hazgradbivteissier(z1, z2, xi1, xi2, 0)
  h1_marg <- xi1 * (exp(xi1 * z1) - 1)
  h2_marg <- xi2 * (exp(xi2 * z2) - 1)
  expect_equal(grad0$eta1, h1_marg)
  expect_equal(grad0$eta2, h2_marg)
})