test-montecarlo-group1.R
In rbiogeme: Interface to 'Biogeme'

native_montecarlo_b01 <- function(number_of_draws, seed) {
  expressions <- reticulate::import("biogeme.expressions", convert = FALSE)
  database_module <- reticulate::import("biogeme.database", convert = FALSE)
  biogeme_module <- reticulate::import("biogeme.biogeme", convert = FALSE)

  database <- database_module$Database(
    "native_montecarlo_b01",
    reticulate::r_to_py(data.frame(FakeColumn = 1.0))
  )
  integrand <- expressions$exp(expressions$Draws("U", "UNIFORM"))
  simulated_integral <- expressions$MonteCarlo(integrand)
  sample_variance <- expressions$MonteCarlo(integrand * integrand) -
    simulated_integral * simulated_integral
  simulations <- reticulate::dict(
    `Analytical Integral` = expressions$exp(1.0) - 1.0,
    `Simulated Integral` = simulated_integral,
    `Sample variance   ` = sample_variance,
    `Std Error         ` = (sample_variance / as.integer(number_of_draws))^0.5,
    `Error             ` = simulated_integral - (expressions$exp(1.0) - 1.0)
  )
  biogeme <- biogeme_module$BIOGEME(
    database,
    simulations,
    number_of_draws = as.integer(number_of_draws),
    seed = as.integer(seed),
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  reticulate::py_to_r(
    biogeme$simulate(the_beta_values = reticulate::dict())
  )
}

native_montecarlo_b02 <- function(number_of_draws, seed) {
  expressions <- reticulate::import("biogeme.expressions", convert = FALSE)
  database_module <- reticulate::import("biogeme.database", convert = FALSE)
  biogeme_module <- reticulate::import("biogeme.biogeme", convert = FALSE)
  draws_module <- reticulate::import("biogeme.draws", convert = FALSE)
  functools <- reticulate::import("functools", convert = FALSE)

  database <- database_module$Database(
    "native_montecarlo_b02",
    reticulate::r_to_py(data.frame(FakeColumn = 1.0))
  )
  integrand <- expressions$exp(expressions$Draws("U", "UNIFORM"))
  simulated_integral <- expressions$MonteCarlo(integrand)
  integrand_halton <- expressions$exp(
    expressions$Draws("U_halton", "UNIFORM_HALTON2")
  )
  simulated_integral_halton <- expressions$MonteCarlo(integrand_halton)
  integrand_halton13 <- expressions$exp(
    expressions$Draws("U_halton13", "HALTON13")
  )
  simulated_integral_halton13 <- expressions$MonteCarlo(integrand_halton13)
  integrand_mlhs <- expressions$exp(
    expressions$Draws("U_mlhs", "UNIFORM_MLHS")
  )
  simulated_integral_mlhs <- expressions$MonteCarlo(integrand_mlhs)
  true_integral <- expressions$exp(1.0) - 1.0
  sample_variance <- expressions$MonteCarlo(integrand * integrand) -
    simulated_integral * simulated_integral
  sample_variance_halton <- expressions$MonteCarlo(
    integrand_halton * integrand_halton
  ) - simulated_integral_halton * simulated_integral_halton
  sample_variance_halton13 <- expressions$MonteCarlo(
    integrand_halton13 * integrand_halton13
  ) - simulated_integral_halton13 * simulated_integral_halton13
  sample_variance_mlhs <- expressions$MonteCarlo(
    integrand_mlhs * integrand_mlhs
  ) - simulated_integral_mlhs * simulated_integral_mlhs

  simulations <- reticulate::dict(
    `Analytical Integral` = true_integral,
    `Simulated Integral` = simulated_integral,
    `Sample variance   ` = sample_variance,
    `Std Error         ` = (sample_variance / as.integer(number_of_draws))^0.5,
    `Error             ` = simulated_integral - true_integral,
    `Simulated Integral (Halton)` = simulated_integral_halton,
    `Sample variance (Halton)   ` = sample_variance_halton,
    `Std Error (Halton)         ` = (
      sample_variance_halton / as.integer(number_of_draws)
    )^0.5,
    `Error (Halton)             ` = simulated_integral_halton - true_integral,
    `Simulated Integral (Halton13)` = simulated_integral_halton13,
    `Sample variance (Halton13)   ` = sample_variance_halton13,
    `Std Error (Halton13)         ` = (
      sample_variance_halton13 / as.integer(number_of_draws)
    )^0.5,
    `Error (Halton13)             ` = simulated_integral_halton13 - true_integral,
    `Simulated Integral (MLHS)` = simulated_integral_mlhs,
    `Sample variance (MLHS)   ` = sample_variance_mlhs,
    `Std Error (MLHS)         ` = (
      sample_variance_mlhs / as.integer(number_of_draws)
    )^0.5,
    `Error (MLHS)             ` = simulated_integral_mlhs - true_integral
  )
  halton13 <- functools$partial(
    draws_module$get_halton_draws,
    base = 13L,
    skip = 10L
  )
  halton13_generator <- draws_module$RandomNumberGeneratorTuple(
    generator = halton13,
    description = "Halton draws with base 13, skipping 10"
  )
  random_number_generators <- reticulate::dict(
    HALTON13 = halton13_generator
  )
  biogeme <- biogeme_module$BIOGEME(
    database,
    simulations,
    random_number_generators = random_number_generators,
    number_of_draws = as.integer(number_of_draws),
    seed = as.integer(seed),
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  reticulate::py_to_r(
    biogeme$simulate(the_beta_values = reticulate::dict())
  )
}

r_montecarlo_b01 <- function(number_of_draws, seed) {
  database <- biogeme_database(
    "fake_database",
    data.frame(FakeColumn = 1.0)
  )
  integrand <- exp(draw("U", "UNIFORM"))
  simulated_integral <- monte_carlo(integrand)
  true_integral <- exp(1.0) - 1.0
  simulations <- list(
    `Analytical Integral` = true_integral,
    `Simulated Integral` = simulated_integral,
    `Sample variance   ` = monte_carlo(integrand * integrand) -
      simulated_integral * simulated_integral,
    `Std Error         ` = sqrt(
      (
        monte_carlo(integrand * integrand) -
          simulated_integral * simulated_integral
      ) / number_of_draws
    ),
    `Error             ` = simulated_integral - true_integral
  )
  model <- biogeme_model(database = database, simulations = simulations)
  simulation <- simulate(
    model,
    beta = setNames(numeric(0), character(0)),
    control = biogeme_control(
      model_name = "r_montecarlo_b01",
      number_of_draws = as.integer(number_of_draws),
      seed = as.integer(seed),
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
  as.data.frame(simulation, check.names = FALSE)
}

r_montecarlo_b02 <- function(number_of_draws, seed) {
  database <- biogeme_database(
    "fakeDatabase",
    data.frame(FakeColumn = 1.0)
  )
  integrand <- exp(draw("U", "UNIFORM"))
  simulated_integral <- monte_carlo(integrand)
  integrand_halton <- exp(draw("U_halton", "UNIFORM_HALTON2"))
  simulated_integral_halton <- monte_carlo(integrand_halton)
  integrand_halton13 <- exp(draw("U_halton13", "HALTON13"))
  simulated_integral_halton13 <- monte_carlo(integrand_halton13)
  integrand_mlhs <- exp(draw("U_mlhs", "UNIFORM_MLHS"))
  simulated_integral_mlhs <- monte_carlo(integrand_mlhs)
  true_integral <- exp(1.0) - 1.0
  simulations <- list(
    `Analytical Integral` = true_integral,
    `Simulated Integral` = simulated_integral,
    `Sample variance   ` = monte_carlo(integrand * integrand) -
      simulated_integral * simulated_integral,
    `Std Error         ` = sqrt(
      (
        monte_carlo(integrand * integrand) -
          simulated_integral * simulated_integral
      ) / number_of_draws
    ),
    `Error             ` = simulated_integral - true_integral,
    `Simulated Integral (Halton)` = simulated_integral_halton,
    `Sample variance (Halton)   ` = monte_carlo(
      integrand_halton * integrand_halton
    ) - simulated_integral_halton * simulated_integral_halton,
    `Std Error (Halton)         ` = sqrt(
      (
        monte_carlo(integrand_halton * integrand_halton) -
          simulated_integral_halton * simulated_integral_halton
      ) / number_of_draws
    ),
    `Error (Halton)             ` = simulated_integral_halton - true_integral,
    `Simulated Integral (Halton13)` = simulated_integral_halton13,
    `Sample variance (Halton13)   ` = monte_carlo(
      integrand_halton13 * integrand_halton13
    ) - simulated_integral_halton13 * simulated_integral_halton13,
    `Std Error (Halton13)         ` = sqrt(
      (
        monte_carlo(integrand_halton13 * integrand_halton13) -
          simulated_integral_halton13 * simulated_integral_halton13
      ) / number_of_draws
    ),
    `Error (Halton13)             ` = simulated_integral_halton13 - true_integral,
    `Simulated Integral (MLHS)` = simulated_integral_mlhs,
    `Sample variance (MLHS)   ` = monte_carlo(
      integrand_mlhs * integrand_mlhs
    ) - simulated_integral_mlhs * simulated_integral_mlhs,
    `Std Error (MLHS)         ` = sqrt(
      (
        monte_carlo(integrand_mlhs * integrand_mlhs) -
          simulated_integral_mlhs * simulated_integral_mlhs
      ) / number_of_draws
    ),
    `Error (MLHS)             ` = simulated_integral_mlhs - true_integral
  )
  model <- biogeme_model(
    database = database,
    simulations = simulations,
    draws = biogeme_draws(
      name = "U_halton13",
      draw_type = "HALTON13",
      number_of_draws = as.integer(number_of_draws),
      seed = as.integer(seed),
      generator = "HALTON13"
    )
  )
  simulation <- simulate(
    model,
    beta = setNames(numeric(0), character(0)),
    control = biogeme_control(
      model_name = "r_montecarlo_b02",
      number_of_draws = as.integer(number_of_draws),
      seed = as.integer(seed),
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
  as.data.frame(simulation, check.names = FALSE)
}

test_that("Monte Carlo Group 1 example files are syntactically valid", {
  files <- list.files(
    rbiogeme_example_path( "montecarlo"),
    pattern = "\\.R$",
    full.names = TRUE
  )
  expect_true(all(file.exists(files)))
  for (file in files) parse(file)
})

test_that("b01 simple integral matches native Biogeme", {
  skip_if_not(
    identical(Sys.getenv("RBIOGEME_RUN_INTEGRATION"), "1"),
    "Set RBIOGEME_RUN_INTEGRATION=1 to run native Monte Carlo equivalence tests"
  )
  skip_if_not(
    rbiogeme_test_configure_python(),
    "Set RBIOGEME_PYTHON to a compatible native Biogeme interpreter"
  )

  temporary_directory <- tempfile("rbiogeme-montecarlo-group1-b01-")
  dir.create(temporary_directory, recursive = TRUE)
  original_directory <- getwd()
  setwd(temporary_directory)
  on.exit(setwd(original_directory), add = TRUE)

  number_of_draws <- 32L
  seed <- 1223L
  r_values <- r_montecarlo_b01(number_of_draws, seed)
  native_values <- as.data.frame(
    native_montecarlo_b01(number_of_draws, seed),
    check.names = FALSE
  )
  expect_identical(names(r_values), names(native_values))
  expect_equal(
    unname(as.matrix(r_values)),
    unname(as.matrix(native_values)),
    tolerance = 1e-12
  )
})

test_that("b02 draw methods and custom Halton13 match native Biogeme", {
  skip_if_not(
    identical(Sys.getenv("RBIOGEME_RUN_INTEGRATION"), "1"),
    "Set RBIOGEME_RUN_INTEGRATION=1 to run native Monte Carlo equivalence tests"
  )
  skip_if_not(
    rbiogeme_test_configure_python(),
    "Set RBIOGEME_PYTHON to a compatible native Biogeme interpreter"
  )

  temporary_directory <- tempfile("rbiogeme-montecarlo-group1-b02-")
  dir.create(temporary_directory, recursive = TRUE)
  original_directory <- getwd()
  setwd(temporary_directory)
  on.exit(setwd(original_directory), add = TRUE)

  number_of_draws <- 32L
  seed <- 1223L
  r_values <- r_montecarlo_b02(number_of_draws, seed)
  native_values <- as.data.frame(
    native_montecarlo_b02(number_of_draws, seed),
    check.names = FALSE
  )
  expect_identical(names(r_values), names(native_values))
  expect_equal(
    unname(as.matrix(r_values)),
    unname(as.matrix(native_values)),
    tolerance = 1e-12
  )
})