inst/examples/montecarlo/plot_b03antithetic_explicit.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b03 explicit. Antithetic pairs generated inside each integrand.
#
# This is the R counterpart of plot_b03antithetic_explicit.py.  Unlike the
# companion example, the antithetic pair is written explicitly as exp(U) plus
# exp(1 - U), and the Monte Carlo result is divided by two.  This demonstrates
# that the antithetic construction is part of the symbolic model expression.

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)
)

# The complement 1 - U is symbolic.  R does not evaluate any draw locally.
uniform_draw <- draw("U", "UNIFORM")
integrand <- exp(uniform_draw) + exp(1 - uniform_draw)
simulated_integral <- monte_carlo(integrand) / 2.0

halton13_draw <- draw("U_halton13", "HALTON13")
integrand_halton13 <- exp(halton13_draw) + exp(1 - halton13_draw)
simulated_integral_halton13 <- monte_carlo(integrand_halton13) / 2.0

mlhs_draw <- draw("U_mlhs", "UNIFORM_MLHS")
integrand_mlhs <- exp(mlhs_draw) + exp(1 - mlhs_draw)
simulated_integral_mlhs <- monte_carlo(integrand_mlhs) / 2.0

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
)

# HALTON13 is a bridge-owned custom generator matching the native example.
model <- biogeme_model(
  database = database,
  simulations = simulation_expressions,
  draws = biogeme_draws(
    name = "U_halton13",
    draw_type = "HALTON13",
    number_of_draws = number_of_draws,
    seed = prepared$seed,
    generator = "HALTON13"
  )
)

simulation <- simulate(
  model,
  beta = empty_beta_values(),
  control = montecarlo_control(
    "b03antithetic_explicit",
    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)