inst/examples/montecarlo/plot_b02simple_integral.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b02. Simple integral using several native draw types.
#
# This is the R counterpart of plot_b02simple_integral.py.  The model
# specification is intentionally kept in this script.  The custom HALTON13
# generator is declared as data-only metadata; the Python bridge implements it
# with Biogeme's public get_halton_draws() API, so R never supplies a callback.

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
)

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

# Each draw() node is compiled into a native Biogeme Draws expression.  The
# built-in draw types below are provided directly by native Biogeme.
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)

# HALTON13 is the user-defined generator from the Python example.  This
# descriptor registers the bridge-owned generator for the draw type.
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

# The sample variance and standard error are also native expressions.  This
# preserves the operation order of the native example.
sample_variance <- monte_carlo(integrand * integrand) -
  simulated_integral * simulated_integral
standard_error <- sqrt(sample_variance / prepared$number_of_draws)
error <- simulated_integral - true_integral

sample_variance_halton <- monte_carlo(integrand_halton * integrand_halton) -
  simulated_integral_halton * simulated_integral_halton
standard_error_halton <- sqrt(sample_variance_halton / prepared$number_of_draws)
error_halton <- simulated_integral_halton - true_integral

sample_variance_halton13 <-
  monte_carlo(integrand_halton13 * integrand_halton13) -
  simulated_integral_halton13 * simulated_integral_halton13
standard_error_halton13 <-
  sqrt(sample_variance_halton13 / prepared$number_of_draws)
error_halton13 <- simulated_integral_halton13 - true_integral

sample_variance_mlhs <- monte_carlo(integrand_mlhs * integrand_mlhs) -
  simulated_integral_mlhs * simulated_integral_mlhs
standard_error_mlhs <- sqrt(sample_variance_mlhs / prepared$number_of_draws)
error_mlhs <- simulated_integral_mlhs - true_integral

simulation_expressions <- list(
  `Analytical Integral` = true_integral,
  `Simulated Integral` = simulated_integral,
  `Sample variance   ` = sample_variance,
  `Std Error         ` = standard_error,
  `Error             ` = error,
  `Simulated Integral (Halton)` = simulated_integral_halton,
  `Sample variance (Halton)   ` = sample_variance_halton,
  `Std Error (Halton)         ` = standard_error_halton,
  `Error (Halton)             ` = error_halton,
  `Simulated Integral (Halton13)` = simulated_integral_halton13,
  `Sample variance (Halton13)   ` = sample_variance_halton13,
  `Std Error (Halton13)         ` = standard_error_halton13,
  `Error (Halton13)             ` = error_halton13,
  `Simulated Integral (MLHS)` = simulated_integral_mlhs,
  `Sample variance (MLHS)   ` = sample_variance_mlhs,
  `Std Error (MLHS)         ` = standard_error_mlhs,
  `Error (MLHS)             ` = error_mlhs
)

# The metadata is declarative.  The complete expression tree and this draw
# registration are compiled once before native Biogeme starts simulation.
draw_metadata <- biogeme_draws(
  name = "U_halton13",
  draw_type = "HALTON13",
  number_of_draws = prepared$number_of_draws,
  seed = prepared$seed,
  generator = "HALTON13"
)
model <- biogeme_model(
  database = database,
  simulations = simulation_expressions,
  draws = draw_metadata
)

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

# Keep the same compact comparison table shown by the native example.
row <- values[1L, , drop = FALSE]
cat("Analytical integral: ", format(row[["Analytical Integral"]]), "\n", sep = "")
cat("                 Uniform       Halton     Halton13         MLHS\n")
cat(
  "Simulated       ",
  format(row[["Simulated Integral"]]), " ",
  format(row[["Simulated Integral (Halton)"]]), " ",
  format(row[["Simulated Integral (Halton13)"]]), " ",
  format(row[["Simulated Integral (MLHS)" ]]),
  "\n",
  sep = ""
)

invisible(simulation)