---
title: "Advanced models with rbiogeme"
author: "rbiogeme contributors"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Advanced models with rbiogeme}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  eval = FALSE,
  collapse = TRUE,
  comment = "#>"
)
```

The advanced interfaces use the same rule as the core interface: R stores a
neutral, data-only specification and the bridge compiles it into native
Biogeme once. Likelihoods, priors, integrations, optimizers, posterior
sampling, and reports remain native operations.

The chunks are intentionally not evaluated during package documentation builds:
they are complete, runnable examples, but advanced estimation and simulation
require a configured Python environment and may create native output files.

## Monte Carlo integration and named draws

`draw()` creates a named draw node inside an expression. `biogeme_draws()` adds
draw metadata such as the draw count, type, seed, or a supplied matrix. A model
with a simulation formula can be used to inspect an integral without writing an
R callback.

```{r monte-carlo}
library(rbiogeme)

database <- biogeme_database(
  "integral_demo",
  data.frame(x = c(-1, 0, 1), id = c(1, 2, 3))
)

z <- random_variable("z")
integrand <- exp(-0.5 * z * z) / sqrt(2 * pi)
quadrature_integral <- integrate_normal(
  expression = integrand,
  name = "z",
  number_of_quadrature_points = 30L
)

u <- draw("u", "UNIFORM_HALTON2")
draw_integral <- monte_carlo(u * variable("x") + 1)

model <- biogeme_model(
  database = database,
  simulations = list(
    quadrature = quadrature_integral,
    monte_carlo = draw_integral
  ),
  draws = biogeme_draws(
    name = "u",
    draw_type = "UNIFORM_HALTON2",
    number_of_draws = 512L,
    seed = 1234L
  )
)

simulate(
  model,
  beta = numeric(),
  control = biogeme_control(
    output_directory = tempfile("rbiogeme-mc-"),
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
)
```

For simulation-based likelihoods, place the Monte Carlo average inside the
likelihood formula, for example:

```{r simulated-likelihood}
b <- biogeme_beta("b", start = 0)
conditional_probability <- exp(b * variable("x") * draw("x_draw", "NORMAL"))

simulated_likelihood <- biogeme_model(
  database = database,
  formula = log(monte_carlo(conditional_probability)),
  draws = biogeme_draws(
    name = "x_draw",
    draw_type = "NORMAL_ANTI",
    number_of_draws = 128L,
    seed = 1223L
  )
)
```

`random_variable()` and `integrate_normal()` describe native quadrature.
`monte_carlo()` describes a native draw average. Random-draw results can vary
according to the native draw design and seed; deterministic models should
match native results up to floating-point precision.

## Bayesian estimation

Bayesian priors are data-only descriptors. They do not contain R functions or
callbacks. Attach one to a parameter with `prior = biogeme_prior(...)` and use
`bayesian_estimate()` with native Bayesian controls.

```{r bayesian}
data <- data.frame(
  choice = c(1, 2, 1, 2, 1, 2),
  x = c(1, 2, 1, 3, 2, 1)
)
database <- biogeme_database("bayesian_demo", data)

asc_2 <- biogeme_beta(
  "asc_2",
  start = 0,
  prior = biogeme_prior("normal", sigma = 2)
)
b_x <- biogeme_beta(
  "b_x",
  start = 0,
  prior = biogeme_prior("student_t", sigma = 3, nu = 5)
)

model <- logit_model(
  database,
  choice = "choice",
  utilities = list(
    `1` = 0,
    `2` = asc_2 + b_x * variable("x")
  )
)

bayesian_fit <- bayesian_estimate(
  model,
  model_name = "bayesian_demo",
  controls = list(
    bayesian_draws = 1000L,
    warmup = 500L,
    chains = 2L,
    target_accept = 0.9,
    calculate_likelihood = TRUE,
    calculate_waic = TRUE,
    calculate_loo = TRUE,
    output_directory = tempfile("rbiogeme-bayesian-")
  )
)

summary(bayesian_fit)
coef(bayesian_fit)
bayesian_stored_variables(bayesian_fit)
```

The Bayesian result contains native posterior summaries and paths to native
output files. Posterior draws remain in the NetCDF file managed by native
Biogeme; ordinary R code does not need to create or manage PyMC or ArviZ
objects.

Posterior simulation uses the same model and a named simulation list:

```{r bayesian-simulation}
simulation_model <- biogeme_model(
  database = database,
  formula = logit_log_probability(
    utilities = list(`1` = 0, `2` = asc_2 + b_x * variable("x")),
    alternative = variable("choice")
  ),
  simulations = list(
    probability = logit_probability(
      utilities = list(`1` = 0, `2` = asc_2 + b_x * variable("x")),
      alternative = variable("choice")
    ),
    marginal_index = b_x * variable("x")
  )
)

posterior_simulation <- simulate_bayesian(
  simulation_model,
  bayesian_results = bayesian_fit,
  percentage_of_draws_to_use = 10,
  lower_quantile = 0.025,
  upper_quantile = 0.975
)
as.data.frame(posterior_simulation)
```

## Piecewise, Box--Cox, and symbolic derivatives

Use expression functions when a native model requires a transformation or a
derivative. The result is still a symbolic node and is compiled along with the
rest of the model.

```{r transformations}
x <- variable("x")
lambda <- biogeme_beta("lambda", start = 1)

piecewise_x <- piecewise(
  expression = x,
  thresholds = c(NULL, 0, 10, NULL),
  betas = list(
    biogeme_beta("piece_1", start = 1),
    biogeme_beta("piece_2", start = 1),
    biogeme_beta("piece_3", start = 1)
  ),
  transform = "formula"
)

boxcox_x <- boxcox(x, lambda)
derivative <- derive(boxcox_x, "x")
```

`Elem(mapping, index)` selects one expression from a named integer mapping and
is useful for category-dependent specifications:

```{r indexed-selection}
category_beta <- Elem(
  mapping = list(
    `1` = biogeme_beta("b_category_1", start = 0),
    `2` = biogeme_beta("b_category_2", start = 0)
  ),
  index = variable("category")
)
```

## MDCEV estimation and forecasting

MDCEV model constructors describe the native variants explicitly. Baseline
utilities, gamma parameters, observed quantities, prices, and the number of
chosen alternatives are named by alternative code.

```{r mdcev}
mdcev_data <- data.frame(
  chosen = c(1, 1, 2, 2),
  quantity_1 = c(1, 2, 0, 1),
  quantity_2 = c(0, 1, 2, 1),
  price_1 = c(2, 2, 3, 2),
  price_2 = c(3, 3, 2, 3),
  x = c(1, 2, 1, 3)
)
mdcev_database <- biogeme_database("mdcev_demo", mdcev_data)

gamma_1 <- biogeme_beta("gamma_1", start = 1, lower = 0)
gamma_2 <- biogeme_beta("gamma_2", start = 1, lower = 0)

mdcev_model <- biogeme_mdcev_model(
  database = mdcev_database,
  model_type = "gamma_profile",
  baseline_utilities = list(
    `1` = biogeme_beta("asc_1", start = 0) + variable("x"),
    `2` = biogeme_beta("asc_2", start = 0) + variable("x")
  ),
  gamma_parameters = list(`1` = gamma_1, `2` = gamma_2),
  prices = list(`1` = variable("price_1"), `2` = variable("price_2")),
  number_of_chosen_alternatives = variable("chosen"),
  consumed_quantities = list(
    `1` = variable("quantity_1"),
    `2` = variable("quantity_2")
  )
)

mdcev_fit <- mdcev_estimate(
  mdcev_model,
  model_name = "mdcev_demo",
  control = biogeme_control(
    output_directory = tempfile("rbiogeme-mdcev-"),
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
)
mdcev_short_summary(mdcev_fit)
mdcev_parameter_table(mdcev_fit)
```

Forecasting accepts native Gumbel draws supplied by the user or generates them
through the native bridge. The two native algorithms can be compared with
`brute_force = TRUE` and `brute_force = FALSE`:

```{r mdcev-forecast}
epsilons <- mdcev_generate_epsilons(
  model = mdcev_model,
  number_of_observations = nrow(mdcev_data),
  number_of_draws = 128L,
  seed = 1234L
)

forecast <- mdcev_forecast(
  model = mdcev_model,
  fit = mdcev_fit,
  total_budget = 10,
  epsilons = epsilons,
  brute_force = FALSE
)
mdcev_forecast_describe(forecast)
mdcev_validate_forecast(
  model = mdcev_model,
  fit = mdcev_fit,
  total_budget = 10,
  epsilons = epsilons
)
```

## Catalogs, segmentation, and assisted specification

Catalogs allow native Biogeme to estimate multiple synchronized specifications
without exposing Python controller objects. A catalog is a named mapping of
expressions; shared controllers select the same specification in several
catalogs.

```{r catalogs}
catalog_controller <- biogeme_catalog_controller(
  name = "time_form",
  specification_names = c("linear", "log")
)

time_catalog <- biogeme_catalog(
  name = "time_catalog",
  expressions = list(
    linear = variable("x"),
    log = log(variable("x") + 1)
  ),
  controller = catalog_controller
)

catalog_database <- biogeme_database(
  "catalog_demo",
  data.frame(
    choice = c(1, 2, 1, 2),
    x = c(1, 2, 3, 4),
    category = c(1, 2, 1, 2)
  )
)

catalog_model <- biogeme_model(
  database = catalog_database,
  formula = logit_log_probability(
    utilities = list(`1` = 0, `2` = time_catalog),
    alternative = variable("choice")
  )
)

catalog_fit <- estimate_catalog(
  catalog_model,
  model_name = "catalog_demo",
  force = TRUE,
  control = biogeme_control(output_directory = tempfile("rbiogeme-catalog-"))
)
catalog_fit$summary
catalog_fit$non_dominated
```

For discrete segmentation, define a mapping from observed values to segment
names and use `segment_beta()` to generate the segmented parameter. A reference
segment is omitted from the parameterization when supplied.

```{r segmentation}
segmentation <- biogeme_segmentation(
  variable = variable("category"),
  mapping = c(`1` = "low", `2` = "high"),
  reference = "low"
)

b_income <- biogeme_beta("b_income", start = 0)
segmented_income <- segment_beta(
  beta = b_income,
  segmentations = list(segmentation)
)
```

`assisted_specification()` combines native catalog search, objective values,
validity rules, Pareto persistence, and final estimation. The Pareto file is a
deliberate checkpoint, so give it a fresh path for a reproducible search or set
`force = FALSE` when resuming is the intended behavior.

```{r assisted}
assisted_fit <- assisted_specification(
  model = catalog_model,
  objectives = "loglikelihood_dimension",
  validity = NULL,
  pareto_file_name = file.path(tempdir(), "rbiogeme-assisted.pareto"),
  model_name = "assisted_demo",
  force = TRUE,
  control = biogeme_control(output_directory = tempfile("rbiogeme-assisted-"))
)
assisted_fit$summary
```

## Sampling of alternatives

Sampling-of-alternatives models use two ordinary R data frames: one for the
alternative universe and one for individuals. A partition describes strata
and sample sizes. Native Biogeme regenerates the sampled choice sets and
constructs the sampled likelihood; the R interface does not sample in an R
callback.

```{r sampling}
alternatives <- data.frame(
  alternative_id = 1:6,
  alt_time = c(5, 7, 6, 8, 9, 4),
  alt_cost = c(2, 3, 4, 2, 5, 3),
  nest = c(1, 1, 1, 2, 2, 2)
)
individuals <- data.frame(
  choice = c(1, 5, 3, 6),
  income = c(1, 2, 1, 3)
)

partition <- biogeme_sampling_partition(
  segments = list(c(1, 2, 3), c(4, 5, 6)),
  sample_sizes = c(2L, 2L)
)

b_alt_time <- biogeme_beta("b_alt_time", start = 0)
b_alt_cost <- biogeme_beta("b_alt_cost", start = 0)
sampled_model <- sampled_alternatives_model(
  alternatives = alternatives,
  individuals = individuals,
  choice_column = "choice",
  id_column = "alternative_id",
  utility = b_alt_time * variable("alt_time") + b_alt_cost * variable("alt_cost"),
  partition = partition,
  biogeme_file_name = file.path(tempdir(), "rbiogeme-sampled.dat"),
  model_type = "logit",
  control = biogeme_control(output_directory = tempfile("rbiogeme-sampled-"))
)

sampled_fit <- estimate_sampled_alternatives(
  sampled_model,
  model_name = "sampled_demo",
  control = biogeme_control(output_directory = tempfile("rbiogeme-sampled-fit-"))
)
coef(sampled_fit)
```

`cross_variable()` describes an individual-by-alternative expression when the
utility uses both data frames. `sampling_segment_sizes()` is a helper for
creating nearly equal sample sizes; it does not perform the random sampling.

## Hybrid-choice specifications

Hybrid-choice models combine a choice likelihood with measurement equations
for latent variables. In R, write the latent-variable expressions and the
ordered-logit or ordered-probit measurement likelihood explicitly, then join
the components in a generic `biogeme_model()` formula. The ordered-response
constructors use numeric category codes and cutpoint expressions:

```{r hybrid-choice}
hybrid_database <- biogeme_database(
  "hybrid_demo",
  data.frame(
    choice = c(1, 2, 1, 2),
    indicator = c(1, 2, 3, 2),
    x = c(1, 2, 1, 3)
  )
)

latent <- biogeme_beta("latent", start = 0) +
  biogeme_beta("latent_x", start = 0) * variable("x")
cut_1 <- biogeme_beta("cut_1", start = -1)
cut_2 <- biogeme_beta("cut_2", start = 1)

indicator_log_probability <- ordered_probit_log_probability(
  eta = latent,
  cutpoints = list(cut_1, cut_2),
  alternative = variable("indicator"),
  categories = c(1, 2, 3)
)

choice_log_probability <- logit_log_probability(
  utilities = list(`1` = 0, `2` = latent),
  alternative = variable("choice")
)

hybrid_model <- biogeme_model(
  database = hybrid_database,
  formula = choice_log_probability + indicator_log_probability
)
```

The same pattern supports a simultaneous measurement-and-choice likelihood.
When a native example uses sequential estimation or parameter overrides, keep
the stages explicit with separate models and use `parameter_overrides` in the
generic constructor rather than changing parameter names.

## Equivalence and reproducibility checklist

For a comparison with a native result, keep the following fixed:

```{r equivalence-checklist}
equivalence_directory <- tempfile("rbiogeme-equivalence-")
dir.create(equivalence_directory)

equivalence_control <- biogeme_control(
  output_directory = equivalence_directory,
  seed = 1223L,
  generate_html = FALSE,
  generate_yaml = FALSE,
  save_iterations = FALSE
)
```

Use the same data rows, derived columns, filters, panel ordering, alternative
codes, starting values, parameter names, draw type, draw count, and seed. For
random-draw or sampling models, document the accepted simulation noise. Never
use `estimate_or_load()` with an implicit or shared path in an equivalence
test; use `estimate()` or set `force = TRUE`.

The package's reference pages document each constructor and operation. The
guides here explain the R syntax and workflow without requiring users to look
up a separate implementation manual.