inst/examples/montecarlo/plot_b07estimation_monte_carlo_mlhs.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b07estimation_monte_carlo_mlhs. 10,000 modified Latin hypercube normal draws.
#
# This is the R counterpart of plot_b07estimation_monte_carlo_mlhs.R. The complete
# normal-mixture likelihood is specified in this file. R constructs a neutral
# expression tree once; native Biogeme performs draw generation, Monte Carlo
# integration, differentiation, optimization, and report generation.

library(rbiogeme)

script_arguments <- commandArgs(trailingOnly = FALSE)
script_argument <- script_arguments[startsWith(script_arguments, "--file=")]
if (length(script_argument) != 1L) {
  stop("Run this example as an R script with Rscript.", call. = FALSE)
}
example_directory <- dirname(normalizePath(sub("^--file=", "", script_argument)))

# The launcher helper contains only CLI/data/output handling. The public
# swissmetro_data() function below applies the native PURPOSE/CHOICE filter and
# defines the native derived cost, availability, and /100-scaled variables.
source(file.path(example_directory, "example_utils.R"))

prepared <- prepare_swissmetro_draw_estimation_example(
  commandArgs(trailingOnly = TRUE),
  example_directory = example_directory,
  default_model = "b07estimation_monte_carlo_mlhs",
  default_number_of_draws = 10000L
)

# estimate() always performs fresh estimation. Remove only this model's old
# report/checkpoint files so no prior YAML or iteration file is reused.
stale_files <- c(
  "b07estimation_monte_carlo_mlhs.yaml",
  "b07estimation_monte_carlo_mlhs.html",
  "__b07estimation_monte_carlo_mlhs.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 <- swissmetro_data(prepared$data)

# These are the same starting values and parameter names as the native
# b07estimation_specification.py helper. asc_sm is absent because the
# Swissmetro utility is the reference alternative.
asc_car <- biogeme_beta("asc_car", start = 0)
asc_train <- biogeme_beta("asc_train", start = 0)
b_time <- biogeme_beta("b_time", start = 0)
b_time_s <- biogeme_beta("b_time_s", start = 1)
b_cost <- biogeme_beta("b_cost", start = 0)

# draw() declares a named native Draws node. It does not generate values in R.
# The draw type below is the only difference between the five b07 variants.
draw_type <- "NORMAL_MLHS"
b_time_draws <- draw("b_time_rnd", draw_type)
b_time_random <- b_time + b_time_s * b_time_draws

# Build the utility functions symbolically. Quoted names preserve the native
# alternative codes: 1 = Train, 2 = Swissmetro, and 3 = Car.
utilities <- list(
  "1" = asc_train + b_time_random * variable("TRAIN_TT_SCALED") +
    b_cost * variable("TRAIN_COST_SCALED"),
  "2" = b_time_random * variable("SM_TT_SCALED") +
    b_cost * variable("SM_COST_SCALED"),
  "3" = asc_car + b_time_random * 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")
)

# logit_probability() is native models.logit(..., i=CHOICE). Monte Carlo
# integration wraps the conditional probability before taking its log.
conditional_probability <- logit_probability(
  utilities = utilities,
  availability = availability,
  alternative = variable("CHOICE")
)
log_probability <- log(monte_carlo(conditional_probability))

# Metadata records the draw design declaratively; the actual Draws node above
# remains part of the same expression tree compiled by the Python bridge.
draw_metadata <- biogeme_draws(
  name = "b_time_rnd",
  draw_type = draw_type,
  number_of_draws = prepared$number_of_draws,
  seed = prepared$seed
)
model <- biogeme_model(
  database = database,
  formula = log_probability,
  draws = draw_metadata,
  control = biogeme_control(
    output_directory = prepared$output,
    model_name = "b07estimation_monte_carlo_mlhs",
    number_of_draws = prepared$number_of_draws,
    seed = prepared$seed,
    analytical_hessian_mode = "automatic",
    generate_html = TRUE,
    generate_yaml = TRUE,
    save_iterations = FALSE
  )
)

cat("Draw type: ", draw_type, "\n", sep = "")
cat("Number of draws: ", prepared$number_of_draws, "\n", sep = "")
fit <- estimate(
  model,
  model_name = "b07estimation_monte_carlo_mlhs",
  control = model$control
)
print(summary(fit))
print(coef(fit))
invisible(fit)