inst/examples/swissmetro/plot_b27_monte_carlo_diagnostic.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b27. Post-estimation Monte Carlo draw-stability diagnostic
#
# This example mirrors plot_b27_monte_carlo_diagnostic.py. It estimates a
# small mixed-logit model, then evaluates the native objective and gradient at
# fresh Monte Carlo draw designs while keeping the estimated parameters fixed.

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. It does not define the model.
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_b27_monte_carlo_model <- function(database) {
  # Parameter names, starts, and the fixed Swissmetro ASC match the native
  # model. The random time coefficient is integrated by Monte Carlo.
  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_cost <- biogeme_beta("b_cost", start = 0)
  b_time <- biogeme_beta("b_time", start = 0)
  b_time_s <- biogeme_beta("b_time_s", start = 1)
  b_time_draws <- biogeme_draws(
    "b_time_rnd", "NORMAL", number_of_draws = 2000, seed = 1223
  )
  b_time_rnd <- b_time + b_time_s * b_time_draws

  v_train <- asc_train + b_time_rnd * variable("TRAIN_TT_SCALED") +
    b_cost * variable("TRAIN_COST_SCALED")
  v_swissmetro <- asc_sm + b_time_rnd * variable("SM_TT_SCALED") +
    b_cost * variable("SM_COST_SCALED")
  v_car <- asc_car + b_time_rnd * variable("CAR_TT_SCALED") +
    b_cost * variable("CAR_CO_SCALED")
  conditional_probability <- logit_probability(
    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")
    ),
    alternative = variable("CHOICE")
  )
  log_probability <- log(monte_carlo(conditional_probability))

  biogeme_model(
    database = database,
    formula = log_probability,
    draws = b_time_draws,
    control = biogeme_control(
    output_directory = prepared$output,
      model_name = "b27_monte_carlo",
      number_of_draws = 2000,
      seed = 1223,
      calculating_second_derivatives = "never",
      monte_carlo_diagnostic_auto = FALSE,
      monte_carlo_diagnostic_draw_factors = "0.5,1.0,2.0",
      monte_carlo_diagnostic_replications = 1,
      monte_carlo_diagnostic_time_budget = 300,
      monte_carlo_diagnostic_max_draws = 4000,
      user_notes = paste0(
        "Post-estimation Monte Carlo draw-stability diagnostic for a mixed ",
        "logit model using the Swissmetro data."
      ),
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
}

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

# Remove only this example's files. This keeps the diagnostic reproducible and
# prevents an old estimation or diagnostic checkpoint from changing the run.
stale_files <- c(
  "b27_monte_carlo.yaml",
  "b27_monte_carlo_diagnostic.yaml",
  "b27_monte_carlo_diagnostic.md"
)
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_b27_monte_carlo_model(database)

# Estimation is separate from the diagnostic. The latter never re-estimates;
# it evaluates native Biogeme's fixed fitted parameter vector.
fit <- estimate(
  model,
  model_name = "b27_monte_carlo",
  control = model$control
)
diagnostic <- check_monte_carlo_stability(
  model = model,
  fit = fit,
  model_name = "b27_monte_carlo",
  control = model$control,
  output_directory = prepared$output,
  basename = "b27",
  resume = FALSE
)

print(diagnostic)

invisible(diagnostic)