inst/examples/indicators/plot_b01expressions.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b01. Mathematical expressions and native derivatives.
#
# This is the R counterpart of plot_b01expressions.py. The R syntax below
# constructs a neutral Biogeme expression. evaluate_biogeme_expression()
# compiles that complete expression in the Python bridge and delegates the
# value, gradient, Hessian, and BHHH calculations to native Biogeme.

library(rbiogeme)

# When an interpreter is supplied, configure it explicitly so the example is
# reproducible from any working directory. The package still supports its
# normal Python-discovery behavior when RBIOGEME_PYTHON is unset.
python <- Sys.getenv("RBIOGEME_PYTHON", unset = "")
if (nzchar(python)) {
  if (!file.exists(python)) {
    stop("The selected Python executable does not exist: ", python, call. = FALSE)
  }
  biogeme_config(python = python)
}

# Beta() creates a free native parameter. Its name and starting value are
# preserved exactly when the expression is compiled by the bridge.
b <- biogeme_beta("b", start = 1)

# Arithmetic on Biogeme expressions builds a symbolic tree; it is not
# evaluated by R. This is the same expression as exp(-b * b + 1) in Python.
expression <- exp(-b * b + 1)

# With beta omitted, the initial evaluation uses the native Beta starting
# value. The returned components contain only ordinary R values and lists.
evaluation <- evaluate_biogeme_expression(expression)
cat("exp(-b * b + 1) = ", evaluation$initial[["function"]], "\n", sep = "")
cat("f = ")
print(evaluation$initial[["function"]])
cat("g = ")
print(evaluation$initial$gradient)
cat("h = ")
print(evaluation$initial$hessian)
cat("BHHH = ")
print(evaluation$initial$bhhh)

# The bridge compiles the expression once and evaluates all requested
# parameter vectors with the same native callable. Each row corresponds to
# one vector of free-parameter values, just as the Python callable accepts a
# vector of Betas.
grid <- c(2, 3, seq(-2, 1.9, by = 0.1))
repeated <- evaluate_biogeme_expression(
  expression,
  points = data.frame(b = grid)
)

as_grid_data_frame <- function(evaluations) {
  do.call(rbind, lapply(evaluations, function(item) {
    data.frame(
      b = as.numeric(item$beta_values[["b"]]),
      f = as.numeric(item[["function"]]),
      g = as.numeric(item$gradient[["b"]]),
      h = as.numeric(item$hessian[["b"]][["b"]])
    )
  }))
}

values <- as_grid_data_frame(repeated$evaluations)
cat("f(2) = ", values$f[[1L]], "\n", sep = "")
cat("g(2) = ", values$g[[1L]], "\n", sep = "")
cat("h(2) = ", values$h[[1L]], "\n", sep = "")
cat("f(3) = ", values$f[[2L]], "\n", sep = "")
cat("g(3) = ", values$g[[2L]], "\n", sep = "")
cat("h(3) = ", values$h[[2L]], "\n", sep = "")

# The native example displays the function and its first two derivatives.
# Rscript is normally non-interactive, so plotting is enabled when the file
# is sourced from an interactive R session.
if (interactive()) {
  plot(
    values$b[-c(1L, 2L)],
    values$f[-c(1L, 2L)],
    type = "l",
    xlab = "b",
    ylab = "value",
    col = "black",
    lty = 1,
    ylim = range(values[-c(1L, 2L), c("f", "g", "h")]),
    main = "Native Biogeme expression and derivatives"
  )
  lines(values$b[-c(1L, 2L)], values$g[-c(1L, 2L)], col = "blue")
  lines(values$b[-c(1L, 2L)], values$h[-c(1L, 2L)], col = "red")
  legend(
    "top",
    legend = c("f(b)", "f'(b)", "f''(b)"),
    col = c("black", "blue", "red"),
    lty = 1,
    bty = "n"
  )
}

invisible(repeated)