tests/testthat/test-bayesian-bridge.R

In rbiogeme: Interface to 'Biogeme'

test_that("Bayesian estimation uses the native bridge and serialises results", {
  skip_if_not(
    identical(Sys.getenv("RBIOGEME_RUN_INTEGRATION"), "1"),
    "Set RBIOGEME_RUN_INTEGRATION=1 to run Bayesian bridge integration tests"
  )
  if (!nzchar(Sys.getenv("PYTENSOR_FLAGS"))) {
    Sys.setenv(PYTENSOR_FLAGS = paste0("base_compiledir=", tempfile("pytensor-")))
  }

  skip_if_not(
    rbiogeme_test_configure_python(),
    "Set RBIOGEME_PYTHON to a compatible native Biogeme interpreter"
  )

  output <- tempfile("rbiogeme-bayesian-bridge-")
  dir.create(output)
  old_directory <- setwd(output)
  on.exit(setwd(old_directory), add = TRUE)

  database <- biogeme_database(
    "toy_bayesian",
    data.frame(choice = c(1, 2, 1, 2), income = c(10, 20, 15, 25))
  )
  beta <- biogeme_beta("beta_income", start = -0.1)
  model <- logit_model(
    database = database,
    choice = "choice",
    utilities = list(`1` = beta * variable("income"), `2` = 0)
  )

  fit <- bayesian_estimate(
    model,
    model_name = "bayesian_bridge_smoke",
    controls = biogeme_control(
      bayesian_draws = 8,
      warmup = 8,
      chains = 1,
      mcmc_sampling_strategy = "pymc",
      sample_from_prior = FALSE,
      calculate_likelihood = FALSE,
      calculate_loo = FALSE,
      generate_html = FALSE
    )
  )

  expect_s3_class(fit, "biogeme_bayesian_fit")
  expect_equal(fit$posterior_draws, 8L)
  expect_named(fit$beta_values, fit$beta_names)
  expect_true(file.exists(file.path(output, fit$netcdf_file)))
  expect_true(file.exists(file.path(output, fit$yaml_file)))
  expect_s3_class(summary(fit), "summary.biogeme_bayesian_fit")
  expect_named(coef(fit), fit$beta_names)
})

test_that("declarative priors and posterior simulation stay native", {
  skip_if_not(
    identical(Sys.getenv("RBIOGEME_RUN_INTEGRATION"), "1"),
    "Set RBIOGEME_RUN_INTEGRATION=1 to run Bayesian bridge integration tests"
  )
  if (!nzchar(Sys.getenv("PYTENSOR_FLAGS"))) {
    Sys.setenv(PYTENSOR_FLAGS = paste0("base_compiledir=", tempfile("pytensor-")))
  }

  skip_if_not(
    rbiogeme_test_configure_python(),
    "Set RBIOGEME_PYTHON to a compatible native Biogeme interpreter"
  )

  output <- tempfile("rbiogeme-bayesian-prior-")
  dir.create(output)
  old_directory <- setwd(output)
  on.exit(setwd(old_directory), add = TRUE)

  database <- biogeme_database(
    "toy_bayesian_prior",
    data.frame(choice = c(1, 2, 1, 2), x = c(1, 2, 1, 2))
  )
  beta <- biogeme_beta(
    "beta",
    start = -0.1,
    upper = 0,
    prior = biogeme_prior("student_t", sigma = 10, nu = 5)
  )
  model <- logit_model(
    database = database,
    choice = "choice",
    utilities = list(`1` = beta * variable("x"), `2` = 0)
  )
  fit <- bayesian_estimate(
    model,
    model_name = "bayesian_prior_smoke",
    control = biogeme_control(
      bayesian_draws = 8,
      warmup = 8,
      chains = 1,
      mcmc_sampling_strategy = "pymc",
      calculate_likelihood = FALSE,
      calculate_loo = FALSE,
      generate_html = FALSE
    )
  )

  simulation_model <- biogeme_model(
    database,
    simulations = list(
      probability = logit_probability(
        utilities = list(`1` = beta * variable("x"), `2` = 0),
        alternative = 1
      )
    )
  )
  simulated <- simulate_bayesian(
    simulation_model,
    fit,
    percentage_of_draws_to_use = 100,
    control = biogeme_control(
      calculate_likelihood = FALSE,
      calculate_loo = FALSE
    )
  )
  expect_equal(simulated$posterior_draws, 8L)
  expect_equal(nrow(simulated$values), 4L)
  expect_true("probability_mean" %in% names(simulated$values))
})

test_that("posterior means by observation stay native", {
  skip_if_not(
    identical(Sys.getenv("RBIOGEME_RUN_INTEGRATION"), "1"),
    "Set RBIOGEME_RUN_INTEGRATION=1 to run Bayesian bridge integration tests"
  )
  if (!nzchar(Sys.getenv("PYTENSOR_FLAGS"))) {
    Sys.setenv(PYTENSOR_FLAGS = paste0("base_compiledir=", tempfile("pytensor-")))
  }
  skip_if_not(
    rbiogeme_test_configure_python(),
    "Set RBIOGEME_PYTHON to a compatible native Biogeme interpreter"
  )

  output <- tempfile("rbiogeme-bayesian-observation-")
  dir.create(output)
  old_directory <- setwd(output)
  on.exit(setwd(old_directory), add = TRUE)

  database <- biogeme_database(
    "toy_bayesian_observation",
    data.frame(choice = c(1, 2, 1, 2), x = c(1, 2, 1, 2))
  )
  beta <- biogeme_beta("beta_observation", start = -0.1)
  random_beta <- distributed_parameter(
    "beta_observation_rnd",
    beta + draw("beta_observation_eps", "NORMAL")
  )
  model <- biogeme_model(
    database,
    formula = logit_log_probability(
      utilities = list(`1` = random_beta * variable("x"), `2` = 0),
      alternative = variable("choice")
    )
  )

  fit <- bayesian_estimate(
    model,
    model_name = "bayesian_observation_smoke",
    controls = biogeme_control(
      bayesian_draws = 4,
      warmup = 4,
      chains = 1,
      mcmc_sampling_strategy = "pymc",
      calculate_likelihood = FALSE,
      calculate_loo = FALSE,
      generate_html = FALSE,
      generate_yaml = FALSE,
      generate_netcdf = TRUE
    )
  )
  table <- bayesian_posterior_mean_by_observation(fit, "beta_observation_rnd")
  expect_s3_class(table, "data.frame")
  expect_equal(nrow(table), 4L)
  expect_identical(names(table), "beta_observation_rnd")
  expect_true(all(is.finite(table[[1L]])))
})

test_that("triangular Bayesian draws use the native PyMC factory", {
  skip_if_not(
    identical(Sys.getenv("RBIOGEME_RUN_INTEGRATION"), "1"),
    "Set RBIOGEME_RUN_INTEGRATION=1 to run Bayesian bridge integration tests"
  )
  if (!nzchar(Sys.getenv("PYTENSOR_FLAGS"))) {
    Sys.setenv(PYTENSOR_FLAGS = paste0("base_compiledir=", tempfile("pytensor-")))
  }
  skip_if_not(
    rbiogeme_test_configure_python(),
    "Set RBIOGEME_PYTHON to a compatible native Biogeme interpreter"
  )

  output <- tempfile("rbiogeme-bayesian-triangular-")
  dir.create(output)
  old_directory <- setwd(output)
  on.exit(setwd(old_directory), add = TRUE)

  database <- biogeme_database(
    "toy_bayesian_triangular",
    data.frame(choice = c(1, 2, 1, 2), x = c(1, 2, 1, 2))
  )
  beta <- biogeme_beta("beta_triangular", start = -0.1)
  random_beta <- distributed_parameter(
    "beta_triangular_rnd",
    beta + draw("beta_triangular_eps", "TRIANGULAR")
  )
  model <- biogeme_model(
    database,
    formula = logit_log_probability(
      utilities = list(`1` = random_beta * variable("x"), `2` = 0),
      alternative = variable("choice")
    )
  )
  fit <- bayesian_estimate(
    model,
    model_name = "bayesian_triangular_smoke",
    controls = biogeme_control(
      bayesian_draws = 4,
      warmup = 4,
      chains = 1,
      mcmc_sampling_strategy = "pymc",
      calculate_likelihood = FALSE,
      calculate_loo = FALSE,
      generate_html = FALSE,
      generate_yaml = FALSE,
      generate_netcdf = TRUE
    )
  )
  expect_s3_class(fit, "biogeme_bayesian_fit")
  expect_true(file.exists(fit$netcdf_file))
})