inst/examples/montecarlo/plot_b01simple_integral.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b01. Simple integral using native Monte Carlo integration.
#
# This is the R counterpart of plot_b01simple_integral.py.  The expressions
# below are neutral R nodes; simulate() compiles the complete tree once and
# native Biogeme evaluates every draw and arithmetic operation.

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 = 2000L,
  default_multiplier = 100000L
)

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

# Draws and MonteCarlo() are symbolic nodes.  No random numbers are generated
# in R and no R function is called from a native likelihood evaluation.
integrand <- exp(draw("U", "UNIFORM"))
simulated_integral <- monte_carlo(integrand)

# These are the same derived expressions as in the native example.
true_integral <- exp(1.0) - 1.0
sample_variance <- monte_carlo(integrand * integrand) -
  simulated_integral * simulated_integral
standard_error <- sqrt(sample_variance / prepared$number_of_draws)
error <- simulated_integral - true_integral

simulation_expressions <- list(
  `Analytical Integral` = true_integral,
  `Simulated Integral` = simulated_integral,
  `Sample variance   ` = sample_variance,
  `Std Error         ` = standard_error,
  `Error             ` = error
)

# The model is simulation-only: its named expressions are evaluated by native
# Biogeme at an empty, but explicitly named, Beta vector.
model <- biogeme_model(
  database = database,
  simulations = simulation_expressions
)

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

# Repeat the same native expression tree with the larger draw count, as in the
# second BIOGEME object in the Python example.
second_number_of_draws <- prepared$number_of_draws * prepared$multiplier
second_simulation <- simulate(
  model,
  beta = empty_beta_values(),
  control = montecarlo_control(
    paste0("01simpleIntegral_", second_number_of_draws),
    second_number_of_draws,
    prepared$seed
  )
)
cat("Number of draws: ", second_simulation$number_of_draws, "\n", sep = "")
print(as.data.frame(second_simulation, check.names = FALSE))

invisible(list(first = first_simulation, second = second_simulation))