tests/testthat/test-gamma_pois.R

In bang: Bayesian Analysis, No Gibbs

#context("hef, gamms_pois, in-built prior vs user prior")

my_seed <- 47
my_tol <- 1e-5
my_n <- 10

# ------------------------- Pump failure data ------------------------------- #

# 1. Default (gamma) prior

user_prior_fn <- function(x, hpars) {
  return(stats::dexp(x[1], hpars[1], log = TRUE) +
           stats::dexp(x[2], hpars[2], log = TRUE))
}
user_prior <- set_user_prior(user_prior_fn, model = "gamma_pois",
                             hpars = c(0.01, 0.01))

# (i) sampling on (rotated) (log(mean), log(alpha + beta)) scale

# In-built
set.seed(my_seed)
pump_res_a <- hef(model = "gamma_pois", data = pump, n = my_n)
# User
set.seed(my_seed)
pump_res_b <- hef(model = "gamma_pois", data = pump, n = my_n,
                 prior = user_prior)

test_that("gamma-pois: in-built gamma = user gamma, param = trans", {
  testthat::expect_equal(pump_res_a$sim_vals, pump_res_b$sim_vals,
                         tolerance = my_tol)
})

# (ii) Default prior, sampling on (alpha, beta) scale

# In-built
set.seed(my_seed)
pump_res_a <- hef(model = "gamma_pois", data = pump, n = my_n, prior = "gamma",
                 param = "original")
# User
set.seed(my_seed)
pump_res_b <- hef(model = "gamma_pois", data = pump, n = my_n,
                 param = "original", prior = user_prior)

test_that("gamma-pois: in-built gamma = user gamma, param = original", {
  testthat::expect_equal(pump_res_a$sim_vals, pump_res_b$sim_vals,
                         tolerance = my_tol)
})

# --------------------------- Simulated data -------------------------------- #

# Simulate one draw from the posterior predictive distribution based on
# the rat data

pump_res <- hef(model = "gamma_pois", data = pump)
sim_data <- sim_pred_gamma_pois(pump_res$theta_sim_vals, pump, 1)
sim_data <- cbind(sim_data, pump[, 2])

# 1. Default (gamma) prior

user_prior_fn <- function(x, hpars) {
  return(stats::dexp(x[1], hpars[1], log = TRUE) +
           stats::dexp(x[2], hpars[2], log = TRUE))
}
user_prior <- set_user_prior(user_prior_fn, model = "gamma_pois",
                             hpars = c(0.01, 0.01))

# (i) sampling on (rotated) (log(mean), log(alpha + beta)) scale

# In-built
set.seed(my_seed)
sim_res_a <- hef(model = "gamma_pois", data = sim_data, n = my_n)
# User
set.seed(my_seed)
sim_res_b <- hef(model = "gamma_pois", data = sim_data, n = my_n,
                  prior = user_prior)

test_that("gamma-pois: sim_data, in-built gamma = user gamma, param = trans", {
  testthat::expect_equal(pump_res_a$sim_vals, pump_res_b$sim_vals,
                         tolerance = my_tol)
})

# (ii) Default prior, sampling on (alpha, beta) scale

# In-built
set.seed(my_seed)
sim_res_a <- hef(model = "gamma_pois", data = sim_data, n = my_n, prior = "gamma",
                  param = "original")
# User
set.seed(my_seed)
sim_res_b <- hef(model = "gamma_pois", data = sim_data, n = my_n,
                  param = "original", prior = user_prior)

test_that("gamma-pois: sim_data, in-built gamma = user gamma, param = original", {
  testthat::expect_equal(pump_res_a$sim_vals, pump_res_b$sim_vals,
                         tolerance = my_tol)
})