inst/examples/montecarlo/plot_b04normal_mixture_numerical.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b04. Normal mixture of logits with numerical integration.
#
# This is the R counterpart of plot_b04normal_mixture_numerical.py.  The
# parameter values are fixed, as in the native example.  random_variable() and
# integrate_normal() create native expression nodes; the numerical integral is
# evaluated by Biogeme rather than approximated in R.

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)))
source(file.path(example_directory, "example_utils.R"))
source(file.path(example_directory, "swissmetro_one.R"))

prepared <- prepare_swissmetro_one_example(
  commandArgs(trailingOnly = TRUE),
  example_directory = example_directory,
  default_number_of_draws = 10000L
)
inputs <- prepare_swissmetro_one_database(prepared$data)
database <- inputs$database

# Fixed values copied from the native example.
asc_car <- 0.137
asc_train <- -0.402
asc_sm <- 0
b_time <- -2.26
b_time_s <- 1.66
b_cost <- -1.29

# omega is a native random-variable node used by IntegrateNormal().
omega <- random_variable("omega")
b_time_random <- b_time + b_time_s * omega

# Construct the three native utilities and their availability expressions.
v_train <- asc_train + b_time_random * inputs$train_tt_scaled +
  b_cost * inputs$train_cost_scaled
v_swissmetro <- asc_sm + b_time_random * inputs$sm_tt_scaled +
  b_cost * inputs$sm_cost_scaled
v_car <- asc_car + b_time_random * inputs$car_tt_scaled +
  b_cost * inputs$car_co_scaled
utilities <- list(`1` = v_train, `2` = v_swissmetro, `3` = v_car)
availability <- list(
  `1` = inputs$train_av_sp,
  `2` = inputs$sm_av,
  `3` = inputs$car_av_sp
)

# logit_probability() is compiled to Biogeme's native logit probability.
conditional_probability <- logit_probability(
  utilities = utilities,
  availability = availability,
  alternative = inputs$choice
)
numerical_integral <- integrate_normal(conditional_probability, "omega")

model <- biogeme_model(
  database = database,
  simulations = list(Numerical = numerical_integral)
)
simulation <- simulate(
  model,
  beta = empty_beta_values(),
  control = montecarlo_control(
    "b04normal_mixture_numerical",
    prepared$number_of_draws,
    prepared$seed
  )
)
cat("Mixture of logit - numerical integration: ")
cat(as.data.frame(simulation, check.names = FALSE)[["Numerical"]], "\n")

invisible(simulation)