inst/examples/bayesian_swissmetro/plot_b01b_logit.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b01b. Bayesian logit estimation with custom declarative priors.
#
# The prior descriptions are data-only rbiogeme objects.  They are translated
# by the bridge into native PyMC distributions before sampling; no R callback
# is called from a likelihood, gradient, or MCMC evaluation.

library(rbiogeme)

# prepare_swissmetro_example() is defined in ../swissmetro/example_utils.R.
# It handles command-line paths and creates the output directory.
script_path <- commandArgs(trailingOnly = FALSE)
script_path <- sub("^--file=", "", script_path[startsWith(script_path, "--file=")][[1L]])
source(file.path(dirname(normalizePath(script_path)), "..", "swissmetro", "example_utils.R"))

build_model <- function(database) {
  # A larger normal prior standard deviation for the two ASCs is equivalent to
  # native Beta(..., sigma_prior = 30).
  asc_car <- biogeme_beta("asc_car", start = 0, sigma_prior = 30)
  asc_train <- biogeme_beta("asc_train", start = 0, sigma_prior = 30)

  # biogeme_prior() describes the custom Student-t factory used by native
  # Biogeme.  The Beta upper bound truncates the prior at zero.
  student_prior <- biogeme_prior("student_t", sigma = 10, nu = 5)
  b_time <- biogeme_beta(
    "b_time", start = -1, upper = 0, prior = student_prior
  )
  b_cost <- biogeme_beta(
    "b_cost", start = -1, upper = 0, prior = student_prior
  )

  # The Swissmetro ASC is normalized to zero.  It need not occur in the
  # utilities, just as in the native b01b example.
  utilities <- list(
    `1` = asc_train + b_time * variable("TRAIN_TT_SCALED") +
      b_cost * variable("TRAIN_COST_SCALED"),
    `2` = b_time * variable("SM_TT_SCALED") +
      b_cost * variable("SM_COST_SCALED"),
    `3` = asc_car + b_time * variable("CAR_TT_SCALED") +
      b_cost * variable("CAR_CO_SCALED")
  )
  logit_model(
    database = database,
    choice = "CHOICE",
    utilities = utilities,
    availability = list(
      `1` = variable("TRAIN_AV_SP"),
      `2` = variable("SM_AV"),
      `3` = variable("CAR_AV_SP")
    )
  )
}

prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "b01b_logit"
)
database <- swissmetro_data(prepared$data)
unlink(file.path(prepared$output, c("b01b_logit.yaml", "b01b_logit.nc", "b01b_logit.html")))

model <- build_model(database)
fit <- bayesian_estimate(
  model,
  model_name = "b01b_logit",
  control = biogeme_control(
    output_directory = prepared$output,
    mcmc_sampling_strategy = "pymc",
    generate_yaml = TRUE,
    generate_html = TRUE,
    generate_netcdf = TRUE
  )
)

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