inst/examples/indicators/plot_b03simulation.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b03. Simulation of a nested-logit choice model.
#
# This is the R counterpart of plot_b03simulation.py. The model equations are
# kept in this file so the example remains readable and self-contained. The
# data and output helpers are intentionally separate from the specification.

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

# prepare_indicator_estimation_example() is defined in indicator_utils.R. It
# selects optima.dat by default, configures Python, and creates a fresh output
# directory. read_optima_database() is defined in optima.R and reproduces the
# native Optima filters and derived variables.
source(file.path(example_directory, "indicator_utils.R"))
source(file.path(example_directory, "optima.R"))

prepared <- prepare_indicator_estimation_example(
  commandArgs(trailingOnly = TRUE),
  example_directory = example_directory,
  default_model = "b03simulation",
  default_bootstrap_samples = 1L
)
database <- read_optima_database(
  prepared$data_path,
  name = "b03simulation_optima"
)

# Parameter definitions match scenarios.py and plot_b02estimation.R. The
# public-transport ASC is fixed, while all other parameters are estimated.
asc_car <- biogeme_beta("asc_car", start = 0)
asc_pt <- biogeme_beta("asc_pt", start = 0, fixed = TRUE)
asc_sm <- biogeme_beta("asc_sm", start = 0)
beta_time_fulltime <- biogeme_beta("beta_time_fulltime", start = 0)
beta_time_other <- biogeme_beta("beta_time_other", start = 0)
beta_dist_male <- biogeme_beta("beta_dist_male", start = 0)
beta_dist_female <- biogeme_beta("beta_dist_female", start = 0)
beta_dist_unreported <- biogeme_beta("beta_dist_unreported", start = 0)
beta_cost <- biogeme_beta("beta_cost", start = 0)
mu_no_car <- biogeme_beta("mu_no_car", start = 1, lower = 1, upper = 2)

# Scale the data and construct native indicator expressions. Comparisons are
# compiled to Biogeme and are not evaluated by R.
time_pt_scaled <- variable("TimePT") / 200
time_car_scaled <- variable("TimeCar") / 200
cost_car_scaled <- variable("CostCarCHF") / 10
distance_scaled <- variable("distance_km") / 5
male <- variable("Gender") == 1
female <- variable("Gender") == 2
unreported_gender <- variable("Gender") == -1
fulltime <- variable("OccupStat") == 1
not_fulltime <- variable("OccupStat") != 1
marginal_cost_pt_scaled <- variable("MarginalCostPT") / 10

# Utility equations use the native alternative coding: 0 = PT, 1 = car,
# and 2 = slow modes.
v_pt <- asc_pt + beta_time_fulltime * time_pt_scaled * fulltime +
  beta_time_other * time_pt_scaled * not_fulltime +
  beta_cost * marginal_cost_pt_scaled
v_car <- asc_car + beta_time_fulltime * time_car_scaled * fulltime +
  beta_time_other * time_car_scaled * not_fulltime +
  beta_cost * cost_car_scaled
v_sm <- asc_sm + beta_dist_male * distance_scaled * male +
  beta_dist_female * distance_scaled * female +
  beta_dist_unreported * distance_scaled * unreported_gender
utilities <- list(`0` = v_pt, `1` = v_car, `2` = v_sm)
nests <- nested_nests(
  choice_set = c(0, 1, 2),
  nests = list(
    nested_nest(mu_no_car, alternatives = c(0, 2), name = "no_car"),
    nested_nest(1, alternatives = 1, name = "car")
  )
)
prob_pt <- nested_probability(utilities, NULL, nests, 0)
prob_car <- nested_probability(utilities, NULL, nests, 1)
prob_sm <- nested_probability(utilities, NULL, nests, 2)
log_probability <- nested_log_probability(
  utilities,
  availability = NULL,
  nests = nests,
  alternative = variable("Choice")
)

# Fit the same likelihood as b02 in a separate model object. Keeping the
# simulation formulas out of this object prevents the weight expression below
# from becoming an estimation weight.
estimation_model <- biogeme_model(
  database = database,
  formula = log_probability
)
fit <- estimate(
  estimation_model,
  model_name = "b02estimation",
  control = biogeme_control(
    output_directory = prepared$output,
    model_name = "b02estimation",
    bootstrap_samples = prepared$bootstrap_samples,
    generate_html = TRUE,
    generate_yaml = TRUE,
    save_iterations = FALSE,
    user_notes = "Nested-logit estimation and direct simulation for the Optima indicators example."
  ),
  run_bootstrap = FALSE
)

# The first route mirrors the explicit native get_value_c calls. Every
# expression is compiled before native evaluation, and the returned vectors
# are ordinary R values only after Biogeme has evaluated them.
simulation_formulas <- list(
  weight = variable("normalized_weight"),
  `Utility PT` = v_pt,
  `Utility car` = v_car,
  `Utility SM` = v_sm,
  `Prob. PT` = prob_pt,
  `Prob. car` = prob_car,
  `Prob. SM` = prob_sm
)
simulation_model <- biogeme_model(
  database = database,
  simulations = simulation_formulas
)
direct_values <- as.data.frame(do.call(
  cbind,
  lapply(simulation_formulas, function(expression) {
    evaluate_biogeme_expression_c(
      model = simulation_model,
      expression = expression,
      beta = fit,
      aggregation = FALSE,
      number_of_draws = 1000L,
      numerically_safe = FALSE,
      use_jit = TRUE
    )
  })
), check.names = FALSE)
names(direct_values) <- names(simulation_formulas)

# The second route sends all named expressions to native BIOGEME.simulate at
# once. Native Biogeme can recycle common subexpressions across these formulas.
biogeme_values <- as.data.frame(simulate(
  simulation_model,
  beta = fit,
  control = biogeme_control(
    output_directory = prepared$output,model_name = "b03simulation")
), check.names = FALSE)

cat("Explicit native get_value_c evaluation:\n")
print(utils::head(direct_values))
cat("Native BIOGEME.simulate evaluation:\n")
print(utils::head(biogeme_values))

invisible(list(
  fit = fit,
  direct_values = direct_values,
  biogeme_values = biogeme_values
))