inst/examples/bayesian_swissmetro/plot_b01c_logit_simul.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b01c. Traditional and Bayesian simulation of a Swissmetro logit model.
#
# The model is specified here in full.  The first simulation uses posterior
# means; the second uses native Biogeme's posterior-draw simulation and
# returns native mean and quantile summaries for every observation.

library(rbiogeme)

# prepare_swissmetro_example() is defined in ../swissmetro/example_utils.R.
# It makes this script independent of the 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)), "..", "swissmetro", "example_utils.R"))

build_estimation_model <- function(database) {
  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)
  logit_model(
    database = database,
    choice = "CHOICE",
    utilities = list(
      `1` = asc_train + b_time * variable("TRAIN_TT_SCALED") +
        b_cost * variable("TRAIN_COST_SCALED"),
      `2` = asc_sm + b_time * variable("SM_TT_SCALED") +
        b_cost * variable("SM_COST_SCALED"),
      `3` = asc_car + b_time * variable("CAR_TT_SCALED") +
        b_cost * variable("CAR_CO_SCALED")
    ),
    availability = list(
      `1` = variable("TRAIN_AV_SP"),
      `2` = variable("SM_AV"),
      `3` = variable("CAR_AV_SP")
    )
  )
}

build_simulation_model <- function(database) {
  # TRAIN_TT, SM_TT, and CAR_TT appear explicitly because derive() is a
  # symbolic derivative with respect to a named database variable.
  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)
  v <- list(
    `1` = asc_train + b_time * variable("TRAIN_TT") / 100 +
      b_cost * variable("TRAIN_COST_SCALED"),
    `2` = asc_sm + b_time * variable("SM_TT") / 100 +
      b_cost * variable("SM_COST_SCALED"),
    `3` = asc_car + b_time * variable("CAR_TT") / 100 +
      b_cost * variable("CAR_CO_SCALED")
  )
  probability <- logit_probability(
    utilities = v,
    availability = list(
      `1` = variable("TRAIN_AV_SP"),
      `2` = variable("SM_AV"),
      `3` = variable("CAR_AV_SP")
    ),
    alternative = 1
  )
  biogeme_model(
    database = database,
    simulations = list(
      "Prob. train" = probability,
      "train time elasticity" = derive(probability, "TRAIN_TT") *
        variable("TRAIN_TT") / probability,
      "Value of time" = b_time / b_cost
    )
  )
}

prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "b01c_logit_simul"
)
database <- swissmetro_data(prepared$data)

# A clean directory has no b01a NetCDF file, so create the prerequisite native
# result explicitly.  This preserves the native example's model names while
# ensuring the R script never silently consumes an old result.
unlink(file.path(prepared$output, c("b01a_logit.yaml", "b01a_logit.nc", "b01a_logit.html")))
estimation_model <- build_estimation_model(database)
fit <- bayesian_estimate(
  estimation_model,
  model_name = "b01a_logit",
  control = biogeme_control(
    output_directory = prepared$output,generate_html = FALSE, generate_yaml = TRUE)
)

simulation_model <- build_simulation_model(database)
posterior_mean_simulation <- simulate(
  simulation_model,
  beta = fit,
  expressions = simulation_model$simulations
)
cat("Simulation using posterior means\n")
print(head(posterior_mean_simulation$values))

# This operation loads the native NetCDF result and delegates all posterior
# draw evaluation and quantile calculation to Python Biogeme.
bayesian_simulation <- simulate_bayesian(
  simulation_model,
  bayesian_results = fit,
  percentage_of_draws_to_use = 3
)
cat("Bayesian simulation\n")
print(head(bayesian_simulation$values))
invisible(bayesian_simulation)