plot_b08_boxcox.R
In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b08. Box--Cox transformations
#
# This example mirrors plot_b08_boxcox.py. The travel-time variables are
# transformed by a common estimated Box--Cox exponent before entering the
# native multinomial-logit utilities.

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_b08_boxcox_model <- function(database) {
  # These parameter names, starts, bounds, and the fixed Swissmetro ASC match
  # the native Python example exactly.
  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)
  boxcox_parameter <- biogeme_beta(
    "boxcox_parameter",
    start = 1,
    lower = -10,
    upper = 10
  )

  # boxcox() is symbolic: the complete node is compiled to native Biogeme's
  # BoxCox expression, including its limiting behavior at lambda = 0.
  train_time <- boxcox(variable("TRAIN_TT_SCALED"), boxcox_parameter)
  sm_time <- boxcox(variable("SM_TT_SCALED"), boxcox_parameter)
  car_time <- boxcox(variable("CAR_TT_SCALED"), boxcox_parameter)

  utilities <- list(
    `1` = asc_train + b_time * train_time +
      b_cost * variable("TRAIN_COST_SCALED"),
    `2` = asc_sm + b_time * sm_time +
      b_cost * variable("SM_COST_SCALED"),
    `3` = asc_car + b_time * car_time +
      b_cost * variable("CAR_CO_SCALED")
  )
  availability <- list(
    `1` = variable("TRAIN_AV_SP"),
    `2` = variable("SM_AV"),
    `3` = variable("CAR_AV_SP")
  )

  # logit_model() compiles the utility graph to native models.loglogit for
  # estimation. No R function is called during likelihood evaluation.
  logit_model(
    database = database,
    choice = "CHOICE",
    utilities = utilities,
    availability = availability
  )
}

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

# Always estimate from the expression tree. Remove only exact b08 artifacts so
# an old YAML or iteration file cannot silently be recycled.
stale_files <- c(
  "b08_boxcox.yaml",
  "__b08_boxcox.iter",
  "b08_boxcox.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_b08_boxcox_model(database)
control <- biogeme_control(
    output_directory = prepared$output,
  model_name = "b08_boxcox",
  generate_html = TRUE,
  generate_yaml = FALSE,
  save_iterations = FALSE
)

# The Python example explicitly checks derivatives around the starting value,
# where the Box--Cox exponent may approach zero. This delegates the check to
# native BIOGEME.check_derivatives(); no R callback runs during the check.
derivative_check <- check_derivatives(
  model,
  model_name = "b08_boxcox",
  controls = biogeme_control(
    output_directory = prepared$output,
    model_name = "b08_boxcox",
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  ),
  verbose = TRUE
)

cat(sprintf(
  "Maximum derivative error: gradient %.6g, Hessian %.6g\n",
  max(abs(unlist(derivative_check$errors_gradient))),
  max(abs(unlist(derivative_check$errors_hessian)))
))

# estimate() delegates the Box--Cox likelihood, differentiation, optimization,
# and reporting to native Biogeme.
fit <- estimate(
  model,
  model_name = "b08_boxcox",
  control = control
)

print(summary(fit))
print(coef(fit))
invisible(fit)