inst/examples/montecarlo/plot_b05normal_mixture_monte_carlo.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b05. Normal mixture of logits compared across integration methods.
#
# This is the R counterpart of plot_b05normal_mixture_monte_carlo.py.  The
# helper below is executed only while constructing the symbolic tree.  It is
# not an R callback: after model construction, native Biogeme evaluates every
# logit, draw, and Monte Carlo integral.

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 = 2000000L
)
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

# Build one numerical-integration random variable and five draw-based random
# parameters.  All draw names and types are preserved from the native model.
omega <- random_variable("omega")
b_time_random <- b_time + b_time_s * omega
b_time_random_normal <- b_time + b_time_s * draw("b_normal", "NORMAL")
b_time_random_anti <- b_time + b_time_s * draw("b_anti", "NORMAL_ANTI")
b_time_random_halton <- b_time + b_time_s * draw("b_halton", "NORMAL_HALTON2")
b_time_random_mlhs <- b_time + b_time_s * draw("b_mlhs", "NORMAL_MLHS")
b_time_random_antimlhs <- b_time + b_time_s *
  draw("b_antimlhs", "NORMAL_MLHS_ANTI")

# Construct the conditional logit expression for each random coefficient.
conditional_logit <- function(random_coefficient) {
  v_train <- asc_train + random_coefficient * inputs$train_tt_scaled +
    b_cost * inputs$train_cost_scaled
  v_swissmetro <- asc_sm + random_coefficient * inputs$sm_tt_scaled +
    b_cost * inputs$sm_cost_scaled
  v_car <- asc_car + random_coefficient * inputs$car_tt_scaled +
    b_cost * inputs$car_co_scaled
  logit_probability(
    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
    ),
    alternative = inputs$choice
  )
}

simulation_expressions <- list(
  Numerical = integrate_normal(conditional_logit(b_time_random), "omega"),
  MonteCarlo = monte_carlo(conditional_logit(b_time_random_normal)),
  Antithetic = monte_carlo(conditional_logit(b_time_random_anti)),
  Halton = monte_carlo(conditional_logit(b_time_random_halton)),
  MLHS = monte_carlo(conditional_logit(b_time_random_mlhs)),
  `Antithetic MLHS` = monte_carlo(conditional_logit(b_time_random_antimlhs))
)

model <- biogeme_model(
  database = database,
  simulations = simulation_expressions
)
simulation <- simulate(
  model,
  beta = empty_beta_values(),
  control = montecarlo_control(
    "b05normal_mixture_monte_carlo",
    prepared$number_of_draws,
    prepared$seed
  )
)
cat("Number of draws: ", simulation$number_of_draws, "\n", sep = "")
print(as.data.frame(simulation, check.names = FALSE))

invisible(simulation)