inst/examples/swissmetro/plot_b01d_logit_simul.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b01d. Simulation of probabilities and direct elasticities
#
# This example estimates the b01a model in a clean run, then simulates choice
# probabilities and travel-time elasticities. It compares the general
# symbolic Derive formulation with the closed-form logit formula.

library(rbiogeme)

# The shared helper contains command-line parsing and data preparation. The
# estimation and simulation expressions are specified in this script.
script_path <- commandArgs(trailingOnly = FALSE)
script_path <- sub("^--file=", "", script_path[startsWith(script_path, "--file=")][[1L]])
source(file.path(dirname(normalizePath(script_path)), "example_utils.R"))

build_b01d_model <- function(database) {
  # Define symbolic parameters. The Swissmetro ASC is fixed at zero for
  # identification, just as in b01a.
  asc_car <- biogeme_beta("asc_car", start = 0)
  asc_train <- biogeme_beta("asc_train", start = 0)
  asc_sm <- biogeme_beta("asc_sm", start = 0, fixed = TRUE)
  b_time <- biogeme_beta("b_time", start = 0)
  b_cost <- biogeme_beta("b_cost", start = 0)

  # Keep TRAIN_TT, SM_TT, and CAR_TT explicitly in the utilities. Derive()
  # differentiates with respect to these named data variables; differentiating
  # only through a precomputed scaled column would hide that dependency.
  v_train <- asc_train + b_time * variable("TRAIN_TT") / 100 +
    b_cost * variable("TRAIN_COST_SCALED")
  v_swissmetro <- asc_sm + b_time * variable("SM_TT") / 100 +
    b_cost * variable("SM_COST_SCALED")
  v_car <- asc_car + b_time * variable("CAR_TT") / 100 +
    b_cost * variable("CAR_CO_SCALED")
  utilities <- list(`1` = v_train, `2` = v_swissmetro, `3` = v_car)
  availability <- list(
    `1` = variable("TRAIN_AV_SP"),
    `2` = variable("SM_AV"),
    `3` = variable("CAR_AV_SP")
  )

  # logit_probability() is a symbolic wrapper for native biogeme.models.logit.
  # It returns the probability of one selected alternative.
  prob_train <- logit_probability(utilities, availability, alternative = 1)
  prob_swissmetro <- logit_probability(utilities, availability, alternative = 2)
  prob_car <- logit_probability(utilities, availability, alternative = 3)

  # General elasticity: Derive() creates a native symbolic derivative node.
  general_time_elasticity_train <-
    Derive(prob_train, "TRAIN_TT") * variable("TRAIN_TT") / prob_train
  general_time_elasticity_swissmetro <-
    Derive(prob_swissmetro, "SM_TT") * variable("SM_TT") / prob_swissmetro
  general_time_elasticity_car <-
    Derive(prob_car, "CAR_TT") * variable("CAR_TT") / prob_car

  # Closed-form direct elasticity for the multinomial logit model.
  logit_time_elasticity_train <-
    variable("TRAIN_AV_SP") * (1 - prob_train) * variable("TRAIN_TT") * b_time / 100
  logit_time_elasticity_swissmetro <-
    variable("SM_AV") * (1 - prob_swissmetro) * variable("SM_TT") * b_time / 100
  logit_time_elasticity_car <-
    variable("CAR_AV_SP") * (1 - prob_car) * variable("CAR_TT") * b_time / 100

  # The named list becomes the columns of the native simulation data frame.
  database_model <- logit_model(
    database = database,
    choice = "CHOICE",
    utilities = utilities,
    availability = availability
  )
  database_model$simulations <- list(
    `Prob. train` = prob_train,
    `Prob. Swissmetro` = prob_swissmetro,
    `Prob. car` = prob_car,
    `logit elas. 1` = logit_time_elasticity_train,
    `generic elas. 1` = general_time_elasticity_train,
    `logit elas. 2` = logit_time_elasticity_swissmetro,
    `generic elas. 2` = general_time_elasticity_swissmetro,
    `logit elas. 3` = logit_time_elasticity_car,
    `generic elas. 3` = general_time_elasticity_car
  )
  database_model
}

describe_simulation <- function(values) {
  # This is presentation-only post-processing, equivalent to pandas
  # DataFrame.describe(); the simulation values themselves come from Biogeme.
  columns <- lapply(values, function(column) {
    c(
      count = sum(!is.na(column)),
      mean = mean(column, na.rm = TRUE),
      sd = stats::sd(column, na.rm = TRUE),
      min = min(column, na.rm = TRUE),
      `25%` = unname(stats::quantile(column, 0.25, na.rm = TRUE)),
      `50%` = unname(stats::quantile(column, 0.50, na.rm = TRUE)),
      `75%` = unname(stats::quantile(column, 0.75, na.rm = TRUE)),
      max = max(column, na.rm = TRUE)
    )
  })
  result <- do.call(cbind, columns)
  colnames(result) <- names(values)
  as.data.frame(result, check.names = FALSE)
}

# prepare_swissmetro_example() is defined in example_utils.R. It parses the
# command line, validates the data/Python paths, configures the bridge, reads
# the data, and creates a fresh output directory.
prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "b01d_logit_simul"
)

# b01d's native Python version loads b01a estimates from a previous YAML file.
# For a clean, self-contained R example, estimate those same b01a parameters
# freshly here rather than silently loading an old result or iteration file.
stale_files <- c("b01a_logit.yaml", "__b01a_logit.iter")
stale_files <- file.path(prepared$output, stale_files)
stale_files <- stale_files[file.exists(stale_files)]
if (length(stale_files) > 0L) unlink(stale_files, force = TRUE)

database <- swissmetro_data(prepared$data)
model <- build_b01d_model(database)
fit <- estimate(
  model,
  model_name = "b01a_logit",
  control = biogeme_control(
    output_directory = prepared$output,
    model_name = "b01a_logit",
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
)

# simulate() compiles the complete simulation tree once, then calls native
# Biogeme with the fixed estimates. No R callback is run during evaluation.
simulated <- simulate(
  model,
  beta = fit,
  control = biogeme_control(
    output_directory = prepared$output,model_name = "b01d_logit_simul")
)
print(describe_simulation(as.data.frame(simulated)))

invisible(simulated)