inst/examples/montecarlo/plot_b06estimation_integral.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b06. Estimation of a normal mixture of logits with numerical integration.
#
# This is the R counterpart of plot_b06estimation_integral.py. The complete
# likelihood specification is written below. R creates only a neutral
# expression tree; native Biogeme performs quadrature, differentiation,
# optimization, and report generation after the tree has been compiled once.

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)))

# prepare_swissmetro_estimation_example() is defined in example_utils.R. It
# parses --data, --python, --output, and --quadrature-points, reads the data,
# and creates an output directory. swissmetro_data() then applies the exact
# native filter: PURPOSE is 1 or 3 and CHOICE is not 0, followed by the
# native derived cost, availability, and /100-scaled variables.
source(file.path(example_directory, "example_utils.R"))

prepared <- prepare_swissmetro_estimation_example(
  commandArgs(trailingOnly = TRUE),
  example_directory = example_directory,
  default_model = "06estimation_integral"
)

# estimate() always performs a fresh native estimation. Remove only artifacts
# for this model so an old YAML or iteration file cannot affect this run.
stale_files <- c(
  "06estimation_integral.yaml",
  "06estimation_integral.html",
  "__06estimation_integral.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 starting values and the fixed ASC match the native example. The
# parameter names are preserved exactly because they are part of equivalence.
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_time_s <- biogeme_beta("b_time_s", start = 1)
b_cost <- biogeme_beta("b_cost", start = 0)

# random_variable() is a named native RandomVariable node. It is not an R
# random draw and no R function is called during native likelihood evaluation.
b_time_random <- b_time + b_time_s * random_variable("omega")

# Build the three utilities and availability conditions symbolically. The
# alternative keys 1, 2, and 3 are the native Train, Swissmetro, and Car codes.
utilities <- list(
  `1` = asc_train + b_time_random * variable("TRAIN_TT_SCALED") +
    b_cost * variable("TRAIN_COST_SCALED"),
  `2` = asc_sm + 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() and integrate_normal() compile to Biogeme's native
# logit and deterministic normal-integration nodes. log() makes this the
# observation-level log likelihood estimated by native BIOGEME.
conditional_probability <- logit_probability(
  utilities = utilities,
  availability = availability,
  alternative = variable("CHOICE")
)
log_probability <- log(integrate_normal(
  conditional_probability,
  name = "omega",
  number_of_quadrature_points = prepared$number_of_quadrature_points
))

model <- biogeme_model(
  database = database,
  formula = log_probability,
  control = biogeme_control(
    output_directory = prepared$output,
    model_name = "06estimation_integral",
    generate_html = TRUE,
    generate_yaml = TRUE,
    save_iterations = FALSE
  )
)

cat("Number of quadrature points: ", prepared$number_of_quadrature_points, "\n", sep = "")
fit <- estimate(
  model,
  model_name = "06estimation_integral",
  control = model$control
)

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