inst/examples/swissmetro/plot_b11b_cnl_simul.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b11b. Simulation of a cross-nested logit model
#
# This example estimates the b11a CNL model freshly, then simulates native
# CNL probabilities and travel-time elasticities. The complete specification
# is kept here so the script runs from a clean working directory.

library(rbiogeme)

# 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. The --data, --python, and
# --output options work from any current working directory.
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_b11b_cnl_model <- function(database) {
  # Parameter names, starting values, bounds, and the fixed ASC match b11a.
  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_swissmetro <- biogeme_beta("b_time_swissmetro", start = 0)
  b_time_train <- biogeme_beta("b_time_train", start = 0)
  b_time_car <- biogeme_beta("b_time_car", start = 0)
  b_cost <- biogeme_beta("b_cost", start = 0)
  b_headway_swissmetro <- biogeme_beta("b_headway_swissmetro", start = 0)
  b_headway_train <- biogeme_beta("b_headway_train", start = 0)
  ga_train <- biogeme_beta("ga_train", start = 0)
  ga_swissmetro <- biogeme_beta("ga_swissmetro", start = 0)
  existing_nest_parameter <- biogeme_beta(
    "existing_nest_parameter", start = 1, lower = 1, upper = 5
  )
  public_nest_parameter <- biogeme_beta(
    "public_nest_parameter", start = 1, lower = 1, upper = 5
  )
  alpha_existing <- biogeme_beta(
    "alpha_existing", start = 0.5, lower = 0, upper = 1
  )
  alpha_public <- 1 - alpha_existing

  # Keep raw travel-time variables explicit. Derive() differentiates with
  # respect to these names; using only precomputed scaled columns would make
  # the symbolic derivative zero.
  utilities <- list(
    `1` = asc_train + b_time_train * variable("TRAIN_TT") / 100 +
      b_cost * variable("TRAIN_COST_SCALED") +
      b_headway_train * variable("TRAIN_HE") + ga_train * variable("GA"),
    `2` = asc_sm + b_time_swissmetro * variable("SM_TT") / 100 +
      b_cost * variable("SM_COST_SCALED") +
      b_headway_swissmetro * variable("SM_HE") + ga_swissmetro * variable("GA"),
    `3` = asc_car + b_time_car * variable("CAR_TT") / 100 +
      b_cost * variable("CAR_CO_SCALED")
  )
  availability <- list(
    `1` = variable("TRAIN_AV_SP"),
    `2` = variable("SM_AV"),
    `3` = variable("CAR_AV_SP")
  )

  # Train belongs partly to both nests. Allocation expressions remain native
  # symbolic nodes and are compiled once with the rest of the model.
  nests <- cross_nested_nests(
    choice_set = c(1L, 2L, 3L),
    nests = list(
      cross_nested_nest(
        existing_nest_parameter,
        list(`1` = alpha_existing, `2` = 0, `3` = 1),
        name = "existing"
      ),
      cross_nested_nest(
        public_nest_parameter,
        list(`1` = alpha_public, `2` = 1, `3` = 0),
        name = "public"
      )
    )
  )

  probability_train <- cross_nested_probability(
    utilities, availability, nests, alternative = 1
  )
  probability_swissmetro <- cross_nested_probability(
    utilities, availability, nests, alternative = 2
  )
  probability_car <- cross_nested_probability(
    utilities, availability, nests, alternative = 3
  )

  # Derive() is compiled to native Biogeme differentiation. No R function is
  # called while native Biogeme evaluates these simulation expressions.
  simulations <- list(
    `Prob. train` = probability_train,
    `Prob. Swissmetro` = probability_swissmetro,
    `Prob. car` = probability_car,
    `Elas. 1` = Derive(probability_train, "TRAIN_TT") *
      variable("TRAIN_TT") / probability_train,
    `Elas. 2` = Derive(probability_swissmetro, "SM_TT") *
      variable("SM_TT") / probability_swissmetro,
    `Elas. 3` = Derive(probability_car, "CAR_TT") *
      variable("CAR_TT") / probability_car
  )

  model <- cross_nested_logit_model(
    database = database,
    choice = "CHOICE",
    utilities = utilities,
    availability = availability,
    nests = nests,
    control = biogeme_control(
    output_directory = prepared$output,
      model_name = "b11a_cnl",
      generate_html = TRUE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
  model$simulations <- simulations
  model
}

prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "b11b_cnl_simul"
)

# Always estimate b11a afresh. Remove only exact artifacts belonging to these
# model names so an old YAML or iteration file cannot silently be reused.
stale_files <- c(
  "b11a_cnl.yaml",
  "__b11a_cnl.iter",
  "b11a_cnl.html",
  "b11b_cnl_simul.yaml",
  "__b11b_cnl_simul.iter",
  "b11b_cnl_simul.html"
)
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_b11b_cnl_model(database)

# Estimation, derivatives, and optimization are delegated to native Biogeme.
fit <- estimate(
  model,
  model_name = "b11a_cnl",
  control = model$control
)
print(summary(fit))
print(coef(fit))

# This calls the native CNL nest object's correlation operation.
correlation <- cross_nested_logit_correlation(
  model,
  beta_values = coef(fit),
  alternatives_names = c(`1` = "Train", `2` = "Swissmetro", `3` = "Car")
)
print(correlation)

# simulate() compiles the named probability and elasticity expressions once,
# then evaluates them through native Biogeme at the fixed estimates.
simulation <- simulate(
  model,
  beta = fit,
  control = biogeme_control(
    output_directory = prepared$output,
    model_name = "b11b_cnl_simul",
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
)
simulation_values <- as.data.frame(simulation, check.names = FALSE)
print(utils::head(simulation_values))
cat(sprintf(
  "Aggregate share of train: %.1f%%\n",
  100 * mean(simulation_values[["Prob. train"]])
))

invisible(list(fit = fit, correlation = correlation, simulation = simulation_values))