test_that("Moments and correlation functions work correctly", {
  a1 <- 3; b1 <- 2; a2 <- 3; b2 <- 2; d <- 0.4

  # Marginal mean and variance
  m1 <- mean_sps_bxii(a1, b1)
  expect_true(is.numeric(m1) && m1 > 0)

  v1 <- var_sps_bxii(a1, b1)
  expect_true(is.numeric(v1) && v1 > 0)

  # Product moment
  pm11 <- product_moment_sps_bxii(1, 1, a1, b1, a2, b2, d)
  expect_true(is.numeric(pm11) && pm11 > 0)

  # Covariance and correlation
  cov_val <- cov_sps_bxii(a1, b1, a2, b2, d)
  expect_true(is.numeric(cov_val))

  cor_val <- cor_sps_bxii(a1, b1, a2, b2, d)
  expect_true(is.numeric(cor_val))
  expect_true(cor_val >= -1 && cor_val <= 1)

  # Zero delta yields zero covariance and correlation
  expect_equal(cov_sps_bxii(a1, b1, a2, b2, 0), 0)
  expect_equal(cor_sps_bxii(a1, b1, a2, b2, 0), 0)

  # Conditional moment
  cm1 <- cond_moment_sps_bxii(1, 0.5, a1, b1, a2, b2, d)
  expect_true(is.numeric(cm1) && cm1 > 0)
})

test_that("Moment existence conditions are enforced", {
  # Existence requires alpha * beta > order
  expect_error(mean_sps_bxii(0.5, 1))
  expect_error(var_sps_bxii(1, 1))
  expect_error(product_moment_sps_bxii(2, 2, 1, 1, 1, 1, 0.5))
})