inst/examples/swissmetro/plot_b06c_unif_mixture_integral.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b06c. Uniform mixture with numerical integration
#
# Biogeme's IntegrateNormal integrates over (-Inf, Inf). This example uses a
# logistic change of variable to map that domain to the uniform coefficient's
# support [-1, 1], then includes the Jacobian and density correction in the
# native integrand.

library(rbiogeme)

# The shared helper contains command-line parsing and data preparation. The
# complete change-of-variable and model specification remains in this script.
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_b06c_unif_mixture_integral_model <- function(
    database,
    number_of_quadrature_points = 60L
) {
  # These parameter names and starting values match native b06c. The
  # Swissmetro ASC is fixed at zero to identify the utility scale.
  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)
  omega <- random_variable("omega")

  # Map the standard-normal integration variable omega to x in [-1, 1].
  # The derivatives and density terms are symbolic native expressions.
  lower_bound <- -1.0
  upper_bound <- 1.0
  x <- lower_bound + (upper_bound - lower_bound) / (1 + exp(-omega))
  dx <- (upper_bound - lower_bound) * exp(-omega) /
    ((1 + exp(-omega)) ^ 2)
  b_time_rnd <- b_time + b_time_s * x

  utilities <- list(
    `1` = asc_train + b_time_rnd * variable("TRAIN_TT_SCALED") +
      b_cost * variable("TRAIN_COST_SCALED"),
    `2` = asc_sm + b_time_rnd * variable("SM_TT_SCALED") +
      b_cost * variable("SM_COST_SCALED"),
    `3` = asc_car + b_time_rnd * 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")
  )
  conditional_probability <- logit_probability(
    utilities = utilities,
    availability = availability,
    alternative = variable("CHOICE")
  )

  # The uniform density is 1/(upper-lower). IntegrateNormal expects a normal
  # density, so divide by native normal_pdf(omega) after applying dx and the
  # uniform density. No R calculation occurs during integration or estimation.
  pdf_uniform <- 1 / (upper_bound - lower_bound)
  new_integrand <- conditional_probability * dx * pdf_uniform /
    normal_pdf(omega)
  log_probability <- log(integrate_normal(
    new_integrand,
    name = "omega",
    number_of_quadrature_points = number_of_quadrature_points
  ))

  biogeme_model(
    database = database,
    formula = log_probability,
    control = biogeme_control(
    output_directory = prepared$output,
      model_name = "b06c_unif_mixture_integral",
      generate_html = TRUE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
}

# 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,
# --output, and --quadrature-points options work from any current working
# directory.
prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "b06c_unif_mixture_integral"
)

number_of_quadrature_points <- if (
    !is.null(prepared$options$quadrature_points) &&
      nzchar(prepared$options$quadrature_points)
) {
  example_integer(prepared$options$quadrature_points, "quadrature-points")
} else {
  60L
}

# estimate() always performs fresh native estimation. Remove only exact b06c
# artifacts so an old YAML or iteration file cannot silently be reused.
stale_files <- c(
  "b06c_unif_mixture_integral.yaml",
  "__b06c_unif_mixture_integral.iter",
  "b06c_unif_mixture_integral.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_b06c_unif_mixture_integral_model(
  database,
  number_of_quadrature_points = number_of_quadrature_points
)

cat(sprintf("Number of quadrature points: %d\n", number_of_quadrature_points))

# The complete change-of-variable integration graph is compiled once. Native
# Biogeme performs quadrature, differentiation, optimization, and reporting.
fit <- estimate(
  model,
  model_name = "b06c_unif_mixture_integral",
  control = model$control
)

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