8. Box-Cox transforms

Bayesian estimation of a logit model, with a Box-Cox transform of variables.

Michel Bierlaire, EPFL Mon Nov 03 2025, 13:41:40

from IPython.core.display_functions import display

import biogeme.biogeme_logging as blog
from biogeme.bayesian_estimation import BayesianResults, get_pandas_estimated_parameters
from biogeme.biogeme import BIOGEME
from biogeme.expressions import Beta
from biogeme.models import boxcox, loglogit

See the data processing script: Data preparation for Swissmetro.

from swissmetro_data import (
    CAR_AV_SP,
    CAR_CO_SCALED,
    CAR_TT_SCALED,
    CHOICE,
    SM_AV,
    SM_COST_SCALED,
    SM_TT_SCALED,
    TRAIN_AV_SP,
    TRAIN_COST_SCALED,
    TRAIN_TT_SCALED,
    database,
)

logger = blog.get_screen_logger(level=blog.INFO)

Parameters to be estimated.

asc_car = Beta('asc_car', 0, None, None, 0)
asc_train = Beta('asc_train', 0, None, None, 0)
asc_sm = Beta('asc_sm', 0, None, None, 1)
b_time = Beta('b_time', 0, None, 0, 0)
b_cost = Beta('b_cost', 0, None, 0, 0)
boxcox_parameter = Beta('boxcox_parameter', 1, -10, 10, 0)

Definition of the utility functions.

v_train = (
    asc_train
    + b_time * boxcox(TRAIN_TT_SCALED, boxcox_parameter)
    + b_cost * TRAIN_COST_SCALED
)
v_swissmetro = (
    asc_sm + b_time * boxcox(SM_TT_SCALED, boxcox_parameter) + b_cost * SM_COST_SCALED
)
v_car = (
    asc_car + b_time * boxcox(CAR_TT_SCALED, boxcox_parameter) + b_cost * CAR_CO_SCALED
)

Associate utility functions with the numbering of alternatives.

v = {1: v_train, 2: v_swissmetro, 3: v_car}

Associate the availability conditions with the alternatives.

av = {1: TRAIN_AV_SP, 2: SM_AV, 3: CAR_AV_SP}

Definition of the model. This is the contribution of each observation to the log likelihood function.

log_probability = loglogit(v, av, CHOICE)

Create the Biogeme object.

the_biogeme = BIOGEME(database, log_probability, bayesian_draws=10000, warmup=10000)
the_biogeme.model_name = 'b08_boxcox'

Biogeme parameters read from biogeme.toml.

Estimate the parameters.

try:
    results = BayesianResults.from_netcdf(
        filename=f'saved_results/{the_biogeme.model_name}.nc'
    )
except FileNotFoundError:
    results = the_biogeme.bayesian_estimation()

Loaded NetCDF file size: 511.0 MB
load finished in 5887 ms (5.89 s)

print(results.short_summary())

posterior_predictive_loglike finished in 1561 ms (1.56 s)
expected_log_likelihood finished in 54 ms
best_draw_log_likelihood finished in 54 ms
/Users/bierlair/python_envs/venv313/lib/python3.13/site-packages/arviz/stats/stats.py:1667: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail.
See http://arxiv.org/abs/1507.04544 for details
  warnings.warn(
waic_res finished in 3168 ms (3.17 s)
waic finished in 3168 ms (3.17 s)
/Users/bierlair/python_envs/venv313/lib/python3.13/site-packages/arviz/stats/stats.py:797: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(
loo_res finished in 95121 ms (1.59 min)
loo finished in 95121 ms (1.59 min)
Sample size                                              6768
Sampler                                                  NUTS
Number of chains                                         4
Number of draws per chain                                10000
Total number of draws                                    40000
Acceptance rate target                                   0.9
Run time                                                 0:01:25.194293
Posterior predictive log-likelihood (sum of log mean p)  -5267.39
Expected log-likelihood E[log L(Y|θ)]                    -17340.91
Best-draw log-likelihood (posterior upper bound)         -7689.05
WAIC (Widely Applicable Information Criterion)           -25780955.15
WAIC Standard Error                                      22143194.30
Effective number of parameters (p_WAIC)                  25775687.77
LOO (Leave-One-Out Cross-Validation)                     -37215.61
LOO Standard Error                                       12236.43
Effective number of parameters (p_LOO)                   31948.23

pandas_results = get_pandas_estimated_parameters(estimation_results=results)
display(pandas_results)

Diagnostics computation took 199.0 seconds (cached).
               Name  Value (mean)  ...  ESS (bulk)  ESS (tail)
0         asc_train     -0.142829  ...    4.003203    4.003203
1           asc_car      0.168284  ...    4.003203    4.003203
2            b_time     -1.275373  ...    4.003203    4.003203
3  boxcox_parameter      0.661758  ...    4.003203    4.003203
4            b_cost     -1.377411  ...    4.003203    4.003203

[5 rows x 12 columns]