inst/examples/montecarlo/plot_b03antithetic.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b03. Antithetic draws for a simple integral.
#
# This is the R counterpart of plot_b03antithetic.py.  All three integrals
# are native MonteCarlo() expressions.  The built-in antithetic generators are
# selected by their Biogeme draw names; HALTON13_ANTI is declared below and
# implemented by the bridge with native Biogeme draw utilities.

library(rbiogeme)

script_arguments <- commandArgs(trailingOnly = FALSE)
script_argument <- script_arguments[startsWith(script_arguments, "--file=")]
if (length(script_argument) != 1L) {
  stop("Run this example as an R script with Rscript.", call. = FALSE)
}
example_directory <- dirname(normalizePath(sub("^--file=", "", script_argument)))
source(file.path(example_directory, "example_utils.R"))

prepared <- prepare_montecarlo_example(
  commandArgs(trailingOnly = TRUE),
  default_number_of_draws = 2000000L
)
number_of_draws <- montecarlo_require_even_draws(prepared$number_of_draws)

# A one-row database is required because Biogeme stores draws by observation.
database <- biogeme_database(
  "fake_database",
  data.frame(FakeColumn = 1.0)
)

# Native Biogeme supplies antithetic uniform and antithetic MLHS draws.
integrand <- exp(draw("U", "UNIFORM_ANTI"))
simulated_integral <- monte_carlo(integrand)

# This custom draw type is registered entirely in the Python bridge.  The
# generator returns the first half of the base-13 Halton sequence and then its
# 1 - U complement, exactly as in the native example.
integrand_halton13 <- exp(draw("U_halton13", "HALTON13_ANTI"))
simulated_integral_halton13 <- monte_carlo(integrand_halton13)

integrand_mlhs <- exp(draw("U_mlhs", "UNIFORM_MLHS_ANTI"))
simulated_integral_mlhs <- monte_carlo(integrand_mlhs)

true_integral <- exp(1.0) - 1.0
simulation_expressions <- list(
  `Analytical Integral` = true_integral,
  `Simulated Integral` = simulated_integral,
  `Error             ` = simulated_integral - true_integral,
  `Simulated Integral (Halton13)` = simulated_integral_halton13,
  `Error (Halton13)             ` = simulated_integral_halton13 - true_integral,
  `Simulated Integral (MLHS)` = simulated_integral_mlhs,
  `Error (MLHS)             ` = simulated_integral_mlhs - true_integral
)

draw_metadata <- biogeme_draws(
  name = "U_halton13",
  draw_type = "HALTON13_ANTI",
  number_of_draws = number_of_draws,
  seed = prepared$seed,
  generator = "HALTON13_ANTI"
)
model <- biogeme_model(
  database = database,
  simulations = simulation_expressions,
  draws = draw_metadata
)

simulation <- simulate(
  model,
  beta = empty_beta_values(),
  control = montecarlo_control(
    "b03antithetic",
    number_of_draws,
    prepared$seed
  )
)
values <- as.data.frame(simulation, check.names = FALSE)
cat("Number of draws: ", simulation$number_of_draws, "\n", sep = "")
print(values)

invisible(simulation)