inst/examples/swissmetro/plot_b05c_normal_mixture_simul.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b05c. Simulation of a normal mixture
#
# This example first estimates the b05a normal-mixture model and then uses
# native Biogeme simulation formulas to calculate the integration error and
# the individual time coefficients.  Estimation and simulation both use the
# same symbolic model specification; no R callback is evaluated by Biogeme.

library(rbiogeme)

# The shared helper contains command-line parsing and data preparation. The
# complete parameter, utility, probability, and simulation specification is
# written below 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_b05c_normal_mixture_components <- function(
    database,
    number_of_draws = 10000L,
    seed = 1223L
) {
  # These parameter names and starting values match native b05a. 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)

  # draw() creates a named native Draws node. The random time coefficient is
  # b_time + b_time_s * NORMAL draw, exactly as in the Python example.
  b_time_rnd <- b_time + b_time_s * draw("b_time_rnd", "NORMAL")

  # Utilities and availability are symbolic expressions. They are compiled
  # once into native Biogeme expressions before estimation or simulation.
  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 on the random coefficient, this is the native logit kernel.
  # The observed-choice selector is compiled to models.logit(..., i=CHOICE).
  conditional_probability <- logit_probability(
    utilities = utilities,
    availability = availability,
    alternative = variable("CHOICE")
  )

  # These are the four named expressions simulated by native b05c. Monte
  # Carlo integration, multiplication, subtraction, and square root remain
  # native Biogeme operations; only the final presentation arithmetic below
  # is performed on the returned R data frame.
  integral <- monte_carlo(conditional_probability)
  integral_square <- monte_carlo(conditional_probability * conditional_probability)
  variance <- integral_square - integral * integral
  integration_error <- sqrt(variance / 2.0)
  simulations <- list(
    Numerator = monte_carlo(b_time_rnd * conditional_probability),
    Denominator = integral,
    Integral = integral,
    `Integration error` = integration_error
  )

  draws <- biogeme_draws(
    name = "b_time_rnd",
    draw_type = "NORMAL",
    number_of_draws = number_of_draws,
    seed = seed
  )

  list(
    # Estimation uses only the log likelihood. The simulation expressions are
    # passed separately after fresh estimation, mirroring native b05c.
    model = biogeme_model(
      database = database,
      formula = log(integral),
      draws = draws
    ),
    simulations = simulations,
    draws = draws
  )
}

# 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, --draws, --seed, and --plot options work from any working
# directory.
prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "b05c_normal_mixture_simul"
)

number_of_draws <- if (!is.null(prepared$options$draws) && nzchar(prepared$options$draws)) {
  example_integer(prepared$options$draws, "draws")
} else {
  10000L
}
seed <- if (!is.null(prepared$options$seed) && nzchar(prepared$options$seed)) {
  example_integer(prepared$options$seed, "seed")
} else {
  1223L
}
make_plot <- example_flag(prepared$options$plot, default = TRUE)

# Both estimation and simulation are fresh. Remove only artifacts belonging
# to these exact native model names, so an old YAML or iteration file cannot
# silently affect this self-contained example.
stale_files <- c(
  "b05a_normal_mixture.yaml",
  "__b05a_normal_mixture.iter",
  "b05a_normal_mixture.html",
  "b05normal_mixture_simul.yaml",
  "__b05normal_mixture_simul.iter",
  "b05normal_mixture_simul.html",
  "b05c_normal_mixture_simul.png"
)
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)
components <- build_b05c_normal_mixture_components(
  database,
  number_of_draws = number_of_draws,
  seed = seed
)

estimation_control <- biogeme_control(
    output_directory = prepared$output,
  model_name = "b05a_normal_mixture",
  user_notes = paste0(
    "Example of a mixture of logit models with three alternatives, ",
    "approximated using Monte-Carlo integration."
  ),
  number_of_draws = number_of_draws,
  seed = seed,
  analytical_hessian_mode = "automatic",
  generate_html = TRUE,
  generate_yaml = FALSE,
  save_iterations = FALSE
)

cat(sprintf("Number of draws: %s\n", format(number_of_draws, big.mark = "_")))

# This is a deliberate fresh estimation of b05a. Native b05c reads the b05a
# YAML result, but this R example remains runnable from a clean directory.
fit <- estimate(
  components$model,
  model_name = "b05a_normal_mixture",
  control = estimation_control
)
print(summary(fit))
print(coef(fit))

simulation_control <- biogeme_control(
    output_directory = prepared$output,
  model_name = "b05normal_mixture_simul",
  number_of_draws = number_of_draws,
  seed = seed,
  generate_html = FALSE,
  generate_yaml = FALSE,
  save_iterations = FALSE
)

# simulate() compiles the complete named simulation expression dictionary
# once, then native Biogeme evaluates all rows and draws at the fixed beta
# values. The returned object is an R data frame for presentation only.
simulation <- simulate(
  components$model,
  expressions = components$simulations,
  beta = fit,
  control = simulation_control
)
simulation_values <- as.data.frame(simulation)

# Native b05c computes these post-estimation quantities from the simulation
# columns. The probability integration and its error were already computed
# by native Biogeme; this is only scalar/vector presentation arithmetic.
simulation_values$left <- log(
  simulation_values$Integral - 1.96 * simulation_values[["Integration error"]]
)
simulation_values$right <- log(
  simulation_values$Integral + 1.96 * simulation_values[["Integration error"]]
)
simulation_values$Beta <- simulation_values$Numerator / simulation_values$Denominator

log_likelihood <- sum(log(simulation_values$Integral))
total_integration_error <- sum(simulation_values[["Integration error"]])
average_integration_error <- mean(simulation_values[["Integration error"]])
confidence_interval <- c(
  sum(simulation_values$left),
  sum(simulation_values$right)
)
cat(sprintf("Log likelihood: %.12f\n", log_likelihood))
cat(sprintf(
  "Integration error for %s draws: %.12f\n",
  format(number_of_draws, big.mark = "_"),
  total_integration_error
))
cat(sprintf("In average %.12f per observation.\n", average_integration_error))
cat(sprintf(
  "95%% confidence interval: [%.12f - %.12f]\n",
  confidence_interval[[1L]],
  confidence_interval[[2L]]
))

# The native example displays a histogram of individual coefficients and
# overlays their estimated normal distribution. Save the equivalent plot so
# the script is also useful in a non-interactive clean working directory.
if (make_plot) {
  beta_values <- coef(fit)[c("b_time", "b_time_s")]
  normalpdf <- function(value, mean = 0, standard_deviation = 1) {
    exp(-((value - mean)^2) / (2 * standard_deviation^2)) /
      (standard_deviation * sqrt(2 * pi))
  }
  png(file.path(prepared$output, "b05c_normal_mixture_simul.png"), width = 1000, height = 700)
  hist(
    simulation_values$Beta,
    probability = TRUE,
    breaks = 20,
    main = "Individual random time coefficients",
    xlab = "Beta"
  )
  x <- seq(
    min(simulation_values$Beta),
    max(simulation_values$Beta),
    by = 0.01
  )
  if (length(x) > 1L) {
    lines(x, normalpdf(x, beta_values[["b_time"]], beta_values[["b_time_s"]]))
  }
  dev.off()
}

invisible(list(fit = fit, simulation = simulation_values))