inst/examples/bayesian_swissmetro/plot_b04_validation.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b04. Out-of-sample validation after Bayesian estimation.
#
# Bayesian estimation is performed by native Biogeme. Native cross-validation
# then re-estimates each training fold and evaluates its validation fold. The
# R wrapper returns the native fold records without implementing a second
# validation or likelihood engine.

library(rbiogeme)

# prepare_swissmetro_example() is defined in ../swissmetro/example_utils.R;
# it reads the data and Python paths and creates a clean example output area.
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_model <- function(database) {
  # This is the same complete specification as native b04validation.
  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")
    )
  )
}

prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "b04validation"
)
database <- swissmetro_data(prepared$data)
unlink(file.path(prepared$output, c("b04validation.yaml", "b04validation.nc", "b04validation.html")))

model <- build_model(database)
fit <- bayesian_estimate(
  model,
  model_name = "b04validation",
  control = biogeme_control(
    output_directory = prepared$output,generate_html = FALSE, generate_yaml = TRUE)
)

# validate() accepts the Bayesian result wrapper and uses its posterior means
# as the native starting values, matching the native b04 workflow.
folds <- validate(model, fit, folds = 5L)
for (fold in folds) {
  values <- fold$simulated_values
  cat(
    "Log likelihood for ", nrow(values), " validation data: ",
    sum(values[[1L]]), "\n", sep = ""
  )
}
invisible(folds)