inst/examples/mdcev_no_outside_good/plot_non_monotonic_forecasting.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# Non-monotonic MDCEV forecasting.
#
# This script is self-contained and includes the complete non-monotonic
# specification, including the mu utility expressions. The native Python
# example loads a saved YAML fit; this R version estimates a fresh fit first
# so no old YAML or iteration file can silently change the forecast.

library(rbiogeme)

# prepare_mdcev_example() is defined in example_utils.R next to this script.
# It reads data.csv by default and accepts --data=, --python=, and --output=
# when the script is run from any current working directory.
script_arguments <- commandArgs(trailingOnly = FALSE)
script_file <- script_arguments[startsWith(script_arguments, "--file=")]
if (length(script_file) != 1L) {
  stop("This example must be run as an R script.", call. = FALSE)
}
script_file <- sub("^--file=", "", script_file)
source(file.path(dirname(normalizePath(script_file)), "example_utils.R"))

prepared <- prepare_mdcev_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "non_monotonic_forecasting"
)

stale_files <- c("non_monotonic.yaml", "non_monotonic.html", "__non_monotonic.iter")
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 <- prepare_mdcev_database(prepared$data)

# The WESML multiplier and all parameter names match the native specification.
weight <- variable("weight") * 1.7718243289995812
cte_shopping <- biogeme_beta("cte_shopping", start = 0)
cte_socializing <- biogeme_beta("cte_socializing", start = 0)
cte_recreation <- biogeme_beta("cte_recreation", start = 0)
number_members_socializing <- biogeme_beta("number_members_socializing", start = 0)
number_members_recreation <- biogeme_beta("number_members_recreation", start = 0)
metropolitan_shopping <- biogeme_beta("metropolitan_shopping", start = 0)
male_shopping <- biogeme_beta("male_shopping", start = 0)
male_socializing <- biogeme_beta("male_socializing", start = 0)
male_recreation <- biogeme_beta("male_recreation", start = 0)
age_15_40_shopping <- biogeme_beta("age_15_40_shopping", start = 0)
age_15_40_recreation <- biogeme_beta("age_15_40_recreation", start = 0)
age_41_60_socializing <- biogeme_beta("age_41_60_socializing", start = 0)
age_41_60_personal <- biogeme_beta("age_41_60_personal", start = 0)
bachelor_socializing <- biogeme_beta("bachelor_socializing", start = 0)
bachelor_personal <- biogeme_beta("bachelor_personal", start = 0)
white_personal <- biogeme_beta("white_personal", start = 0)
spouse_shopping <- biogeme_beta("spouse_shopping", start = 0)
spouse_recreation <- biogeme_beta("spouse_recreation", start = 0)
employed_shopping <- biogeme_beta("employed_shopping", start = 0)
sunday_socializing <- biogeme_beta("sunday_socializing", start = 0)
sunday_personal <- biogeme_beta("sunday_personal", start = 0)

shopping <- cte_shopping +
  metropolitan_shopping * variable("metro") +
  male_shopping * variable("male") +
  age_15_40_shopping * variable("age15_40") +
  spouse_shopping * variable("spousepr") +
  employed_shopping * variable("employed")
socializing <- cte_socializing +
  number_members_socializing * variable("hhsize") +
  male_socializing * variable("male") +
  age_41_60_socializing * variable("age41_60") +
  bachelor_socializing * variable("bachigher") +
  sunday_socializing * variable("Sunday")
recreation <- cte_recreation +
  number_members_recreation * variable("hhsize") +
  male_recreation * variable("male") +
  age_15_40_recreation * variable("age15_40") +
  spouse_recreation * variable("spousepr")
personal <- age_41_60_personal * variable("age41_60") +
  bachelor_personal * variable("bachigher") +
  white_personal * variable("white") +
  sunday_personal * variable("Sunday")

baseline_utilities <- list(
  `1` = shopping,
  `2` = socializing,
  `3` = recreation,
  `4` = personal
)
consumed_quantities <- list(
  `1` = variable("t1") / 60,
  `2` = variable("t2") / 60,
  `3` = variable("t3") / 60,
  `4` = variable("t4") / 60
)
lowest_positive_value <- 0.0001
gamma_parameters <- list(
  `1` = biogeme_beta("gamma_shopping", start = 1, lower = lowest_positive_value),
  `2` = biogeme_beta("gamma_socializing", start = 1, lower = lowest_positive_value),
  `3` = biogeme_beta("gamma_recreation", start = 1, lower = lowest_positive_value),
  `4` = biogeme_beta("gamma_personal", start = 1, lower = lowest_positive_value)
)
alpha_parameters <- list(
  `1` = biogeme_beta("alpha_shopping", start = 0.5, lower = lowest_positive_value, upper = 1 - lowest_positive_value),
  `2` = biogeme_beta("alpha_socializing", start = 0.5, lower = lowest_positive_value, upper = 1 - lowest_positive_value),
  `3` = biogeme_beta("alpha_recreation", start = 0.5, lower = lowest_positive_value, upper = 1 - lowest_positive_value),
  `4` = biogeme_beta("alpha_personal", start = 0.5, lower = lowest_positive_value, upper = 1 - lowest_positive_value)
)
scale_parameter <- biogeme_beta("scale", start = 1, lower = lowest_positive_value)

# These unused definitions are retained because they appear in the native
# specification.py file; only parameters referenced by mu_utilities enter the
# compiled native expression registry.
cte_shopping_mu <- biogeme_beta("cte_shopping_mu", start = 0)
cte_social_mu <- biogeme_beta("cte_social_mu", start = 0)
cte_recreation_mu <- biogeme_beta("cte_recreation_mu", start = 0)
age_15_40_personal_mu <- biogeme_beta("age_15_40_personal_mu", start = 0)

# Non-monotonic utility terms. These equations intentionally match the native
# specification: the four unused mu Betas above are not added to these terms.
holiday_shopping_mu <- biogeme_beta("holiday_shopping_mu", start = 0)
metro_social_mu <- biogeme_beta("metro_social_mu", start = 0)
holiday_recreation_mu <- biogeme_beta("holiday_recreation_mu", start = 0)
male_personal_mu <- biogeme_beta("male_personal_mu", start = 0)
mu_utilities <- list(
  `1` = holiday_shopping_mu * variable("holiday"),
  `2` = metro_social_mu * variable("metro"),
  `3` = holiday_recreation_mu * variable("holiday"),
  `4` = male_personal_mu * variable("male")
)

model <- biogeme_mdcev_model(
  database = database,
  model_type = "non_monotonic",
  baseline_utilities = baseline_utilities,
  gamma_parameters = gamma_parameters,
  alpha_parameters = alpha_parameters,
  mu_utilities = mu_utilities,
  scale_parameter = scale_parameter,
  weights = weight,
  number_of_chosen_alternatives = variable("number_chosen"),
  consumed_quantities = consumed_quantities
)

fit <- mdcev_estimate(
  model = model,
  model_name = "non_monotonic",
  control = biogeme_control(
    output_directory = prepared$output,
    model_name = "non_monotonic",
    tolerance = 0.0004,
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
)

# The native example uses Database.extract_rows([10, 11]). The R helper uses
# one-based row indices, so the equivalent selection is rows 11 and 12.
forecast_database <- biogeme_database_extract_rows(database, c(11L, 12L))
budget_in_hours <- 500

# Validate both native algorithms using the same ten Gumbel draws for each
# observation. The draw generation itself remains native NumPy.
validation_epsilons <- mdcev_generate_epsilons(
  model,
  number_of_observations = biogeme_database_nrow(forecast_database),
  number_of_draws = 10L,
  seed = 12345L
)
stopifnot(mdcev_validate_forecast(
  model = model,
  fit = fit,
  database = forecast_database,
  total_budget = budget_in_hours,
  epsilons = validation_epsilons
))

# Use a larger native draw design for the forecast. Pass identical draw
# matrices to brute-force and analytical native algorithms for a direct
# method comparison. --draws=N supports a quick smoke run; the default follows
# the native example's 2,000 draws.
draws_argument <- prepared$options$draws
number_of_draws <- if (is.null(draws_argument) || !nzchar(draws_argument)) {
  2000L
} else {
  value <- suppressWarnings(as.numeric(draws_argument))
  if (length(value) != 1L || is.na(value) || !is.finite(value) ||
      value < 1 || value != floor(value)) {
    stop("draws must be a positive integer.", call. = FALSE)
  }
  as.integer(value)
}
forecast_epsilons <- mdcev_generate_epsilons(
  model,
  number_of_observations = biogeme_database_nrow(forecast_database),
  number_of_draws = number_of_draws,
  seed = 12345L
)

brute_force_time <- system.time(
  optimal_consumptions_brute_force <- mdcev_forecast(
    model = model,
    fit = fit,
    database = forecast_database,
    total_budget = budget_in_hours,
    epsilons = forecast_epsilons,
    brute_force = TRUE
  )
)
cat(
  "Execution time for ", number_of_draws,
  " draws with brute force algorithm: ", brute_force_time[["elapsed"]],
  " seconds\n",
  sep = ""
)

analytical_time <- system.time(
  optimal_consumptions_analytical <- mdcev_forecast(
    model = model,
    fit = fit,
    database = forecast_database,
    total_budget = budget_in_hours,
    epsilons = forecast_epsilons,
    brute_force = FALSE
  )
)
cat(
  "Execution time for ", number_of_draws,
  " draws with analytical algorithm: ", analytical_time[["elapsed"]],
  " seconds\n",
  sep = ""
)

# mdcev_forecast_describe() returns native pandas DataFrame.describe() output
# for each observation, preserving the post-estimation operation shown in the
# Python example while returning only ordinary R data frames.
print(mdcev_forecast_describe(optimal_consumptions_brute_force)[[1L]])
print(mdcev_forecast_describe(optimal_consumptions_analytical)[[1L]])
print(mdcev_forecast_describe(optimal_consumptions_brute_force)[[2L]])
print(mdcev_forecast_describe(optimal_consumptions_analytical)[[2L]])
invisible(list(
  fit = fit,
  brute_force = optimal_consumptions_brute_force,
  analytical = optimal_consumptions_analytical
))