Note
Go to the end to download the full example code.
Catalog of nonlinear specificationsΒΆ
Investigate of nonlinear specifications for the travel time variables:
linear specification,
Box-Cox transform,
power series,
for a total of 3 specifications. See Bierlaire and Ortelli (2023).
Michel Bierlaire, EPFL Sun Apr 27 2025, 15:47:18
from IPython.core.display_functions import display
import biogeme.biogeme_logging as blog
from biogeme.biogeme import BIOGEME
from biogeme.catalog import Catalog
from biogeme.data.swissmetro 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,
read_data,
)
from biogeme.expressions import Beta, Expression
from biogeme.models import boxcox, loglogit
from biogeme.results_processing import compile_estimation_results, pareto_optimal
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)
b_time = Beta('b_time', 0, None, 0, 0)
b_cost = Beta('b_cost', 0, None, 0, 0)
Non-linear specifications for the travel time.
Parameter of the Box-Cox transform.
lambda_travel_time = Beta('lambda_travel_time', 1, -10, 10, 0)
Coefficients of the power series.
square_tt_coef = Beta('square_tt_coef', 0, None, None, 0)
cube_tt_coef = Beta('cube_tt_coef', 0, None, None, 0)
Function calculation the power series.
def power_series(the_variable: Expression) -> Expression:
"""Generate the expression of a polynomial of degree 3
:param the_variable: variable of the polynomial
"""
return (
the_variable
+ square_tt_coef * the_variable**2
+ cube_tt_coef * the_variable * the_variable**3
)
Train travel time
Linear specification.
linear_train_tt = TRAIN_TT_SCALED
Box-Cox transform.
boxcox_train_tt = boxcox(TRAIN_TT_SCALED, lambda_travel_time)
Power series.
power_train_tt = power_series(TRAIN_TT_SCALED)
Definition of the catalog.
train_tt_catalog = Catalog.from_dict(
catalog_name='train_tt_catalog',
dict_of_expressions={
'linear': linear_train_tt,
'boxcox': boxcox_train_tt,
'power': power_train_tt,
},
)
Swissmetro travel time
Linear specification.
linear_sm_tt = SM_TT_SCALED
Box-Cox transform.
boxcox_sm_tt = boxcox(SM_TT_SCALED, lambda_travel_time)
Power series.
power_sm_tt = power_series(SM_TT_SCALED)
Definition of the catalog. Note that the controller is the same as for train.
sm_tt_catalog = Catalog.from_dict(
catalog_name='sm_tt_catalog',
dict_of_expressions={
'linear': linear_sm_tt,
'boxcox': boxcox_sm_tt,
'power': power_sm_tt,
},
controlled_by=train_tt_catalog.controlled_by,
)
Car travel time
Linear specification.
linear_car_tt = CAR_TT_SCALED
Box-Cox transform.
boxcox_car_tt = boxcox(CAR_TT_SCALED, lambda_travel_time)
Power series.
power_car_tt = power_series(CAR_TT_SCALED)
Definition of the catalog. Note that the controller is the same as for train.
car_tt_catalog = Catalog.from_dict(
catalog_name='car_tt_catalog',
dict_of_expressions={
'linear': linear_car_tt,
'boxcox': boxcox_car_tt,
'power': power_car_tt,
},
controlled_by=train_tt_catalog.controlled_by,
)
Definition of the utility functions.
v_train = asc_train + b_time * train_tt_catalog + b_cost * TRAIN_COST_SCALED
v_swissmetro = b_time * sm_tt_catalog + b_cost * SM_COST_SCALED
v_car = asc_car + b_time * car_tt_catalog + 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)
Read the data
database = read_data()
Create the Biogeme object.
the_biogeme = BIOGEME(
database, log_probability, generate_html=False, generate_yaml=False
)
the_biogeme.model_name = 'b02nonlinear'
Biogeme parameters read from biogeme.toml.
Estimate the parameters.
dict_of_results = the_biogeme.estimate_catalog()
Estimating 3 models.
Biogeme parameters provided by the user.
*** Initial values of the parameters are obtained from the file __b02nonlinear_000000.iter
Cannot read file __b02nonlinear_000000.iter. Statement is ignored.
Starting values for the algorithm: {}
As the model is not too complex, we activate the calculation of second derivatives. To change this behavior, modify the algorithm to "simple_bounds" in the TOML file.
Optimization algorithm: hybrid Newton/BFGS with simple bounds [simple_bounds]
** Optimization: Newton with trust region for simple bounds
Iter. asc_train b_time square_tt_coef cube_tt_coef b_cost asc_car Function Relgrad Radius Rho
0 0 0 0 0 0 0 1.1e+04 0.26 0.5 0.031 -
1 -0.5 -0.5 0 0 -0.3 0.012 9.2e+03 4 5 1 ++
2 -0.5 -0.5 0 0 -0.3 0.012 9.2e+03 4 2.5 -11 -
3 -0.5 -0.5 0 0 -0.3 0.012 9.2e+03 4 1.2 -8.7 -
4 -0.5 -0.5 0 0 -0.3 0.012 9.2e+03 4 0.62 -6.7 -
5 -0.5 -0.5 0 0 -0.3 0.012 9.2e+03 4 0.31 -2.4 -
6 -0.5 -0.5 0 0 -0.3 0.012 9.2e+03 4 0.16 -1.3 -
7 -0.5 -0.5 0 0 -0.3 0.012 9.2e+03 4 0.078 -1.2 -
8 -0.5 -0.5 0 0 -0.3 0.012 9.2e+03 4 0.039 -1.5 -
9 -0.5 -0.5 0 0 -0.3 0.012 9.2e+03 4 0.02 -2.3 -
10 -0.5 -0.5 0 0 -0.3 0.012 9.2e+03 4 0.0098 -3.2 -
11 -0.5 -0.5 0 0 -0.3 0.012 9.2e+03 4 0.0049 -4 -
12 -0.5 -0.5 0 0 -0.3 0.012 9.2e+03 4 0.0024 -4.6 -
13 -0.5 -0.5 0 0 -0.3 0.012 9.2e+03 4 0.0012 -2.5 -
14 -0.5 -0.5 0 0 -0.3 0.012 9.2e+03 4 0.00061 -1.3 -
15 -0.5 -0.5 0 0 -0.3 0.012 9.2e+03 4 0.00031 -0.28 -
16 -0.5 -0.5 0.00031 -0.00031 -0.3 0.012 9.2e+03 3 0.00031 0.64 +
17 -0.5 -0.5 0.00061 -0.00023 -0.3 0.012 9.2e+03 1.2 0.00031 0.81 +
18 -0.5 -0.5 0.00092 -0.00026 -0.3 0.012 9.2e+03 0.15 0.0031 0.99 ++
19 -0.5 -0.5 0.004 -0.00026 -0.3 0.011 9.2e+03 0.27 0.031 1 ++
20 -0.52 -0.53 0.034 -0.0004 -0.3 0.0074 9.1e+03 0.32 0.31 1 ++
21 -0.65 -0.83 0.27 -0.0014 -0.43 -0.0046 8.8e+03 2.4 0.31 0.6 +
22 -0.74 -0.95 0.021 -0.0003 -0.73 0.022 8.7e+03 8.2 0.31 0.68 +
23 -0.75 -1.3 -0.032 -0.00023 -0.78 -0.051 8.6e+03 8.2 3.1 0.96 ++
24 -0.75 -1.3 -0.032 -0.00023 -0.78 -0.051 8.6e+03 8.2 1.5 -89 -
25 -0.75 -1.3 -0.032 -0.00023 -0.78 -0.051 8.6e+03 8.2 0.76 -27 -
26 -0.75 -1.3 -0.032 -0.00023 -0.78 -0.051 8.6e+03 8.2 0.38 -4 -
27 -0.75 -1.3 -0.032 -0.00023 -0.78 -0.051 8.6e+03 8.2 0.19 -1.4 -
28 -0.75 -1.3 -0.032 -0.00023 -0.78 -0.051 8.6e+03 8.2 0.095 -1 -
29 -0.75 -1.3 -0.032 -0.00023 -0.78 -0.051 8.6e+03 8.2 0.048 -1.1 -
30 -0.75 -1.3 -0.032 -0.00023 -0.78 -0.051 8.6e+03 8.2 0.024 -1.3 -
31 -0.75 -1.3 -0.032 -0.00023 -0.78 -0.051 8.6e+03 8.2 0.012 -1.4 -
32 -0.75 -1.3 -0.032 -0.00023 -0.78 -0.051 8.6e+03 8.2 0.006 -1.3 -
33 -0.75 -1.3 -0.032 -0.00023 -0.78 -0.051 8.6e+03 8.2 0.003 -1.1 -
34 -0.75 -1.3 -0.032 -0.00023 -0.78 -0.051 8.6e+03 8.2 0.0015 -1 -
35 -0.75 -1.3 -0.032 -0.00023 -0.78 -0.051 8.6e+03 8.2 0.00075 -0.96 -
36 -0.75 -1.3 -0.032 -0.00023 -0.78 -0.051 8.6e+03 8.2 0.00037 -0.94 -
37 -0.75 -1.3 -0.032 -0.00023 -0.78 -0.051 8.6e+03 8.2 0.00019 -0.93 -
38 -0.75 -1.3 -0.032 -0.00023 -0.78 -0.051 8.6e+03 8.2 9.3e-05 -0.025 -
39 -0.75 -1.3 -0.032 -0.00014 -0.78 -0.051 8.6e+03 4.4 0.00093 1 ++
40 -0.75 -1.3 -0.031 -0.00012 -0.78 -0.051 8.6e+03 3.3 0.00093 0.58 +
41 -0.75 -1.3 -0.031 -0.00013 -0.78 -0.051 8.6e+03 0.1 0.0093 1 ++
42 -0.75 -1.3 -0.027 -0.00015 -0.78 -0.049 8.6e+03 0.0093 0.093 1 ++
43 -0.73 -1.4 -0.053 -3.5e-05 -0.78 -0.03 8.6e+03 0.29 0.93 0.93 ++
44 -0.73 -1.4 -0.053 -3.5e-05 -0.78 -0.03 8.6e+03 0.29 0.34 0.098 -
45 -0.64 -1.7 -0.088 0.00011 -0.78 0.018 8.6e+03 3.2 3.4 1 ++
46 -0.47 -2 -0.096 0.00015 -0.8 0.14 8.6e+03 1.3 34 1.1 ++
47 -0.46 -2.1 -0.097 0.00016 -0.8 0.14 8.6e+03 0.09 3.4e+02 1 ++
48 -0.46 -2.1 -0.097 0.00016 -0.8 0.15 8.6e+03 0.00014 3.4e+03 1 ++
49 -0.46 -2.1 -0.097 0.00016 -0.8 0.15 8.6e+03 0.00014 3.4e+04 1 ++
50 -0.46 -2.1 -0.097 0.00016 -0.8 0.15 8.6e+03 4.4e-06 3.4e+04 1 ++
Optimization algorithm has converged.
Relative gradient: 4.389040213358555e-06
Cause of termination: Relative gradient = 4.4e-06 <= 6.1e-06
Number of function evaluations: 92
Number of gradient evaluations: 41
Number of hessian evaluations: 20
Algorithm: Newton with trust region for simple bound constraints
Number of iterations: 51
Proportion of Hessian calculation: 20/20 = 100.0%
Optimization time: 0:00:00.740818
Calculate second derivatives and BHHH
Biogeme parameters provided by the user.
*** Initial values of the parameters are obtained from the file __b02nonlinear_000001.iter
Cannot read file __b02nonlinear_000001.iter. Statement is ignored.
Starting values for the algorithm: {}
Optimization algorithm: hybrid Newton/BFGS with simple bounds [simple_bounds]
** Optimization: Newton with trust region for simple bounds
Iter. asc_train b_time b_cost asc_car Function Relgrad Radius Rho
0 -0.76 -0.77 -0.7 -0.29 8.8e+03 0.04 10 1.1 ++
1 -0.66 -1.2 -0.77 -0.0015 8.7e+03 0.0064 1e+02 1.1 ++
2 -0.65 -1.3 -0.79 0.016 8.7e+03 0.00012 1e+03 1 ++
3 -0.65 -1.3 -0.79 0.016 8.7e+03 4e-08 1e+03 1 ++
Optimization algorithm has converged.
Relative gradient: 3.954408093567457e-08
Cause of termination: Relative gradient = 4e-08 <= 6.1e-06
Number of function evaluations: 13
Number of gradient evaluations: 9
Number of hessian evaluations: 4
Algorithm: Newton with trust region for simple bound constraints
Number of iterations: 4
Proportion of Hessian calculation: 4/4 = 100.0%
Optimization time: 0:00:00.303418
Calculate second derivatives and BHHH
Biogeme parameters provided by the user.
*** Initial values of the parameters are obtained from the file __b02nonlinear_000002.iter
Cannot read file __b02nonlinear_000002.iter. Statement is ignored.
Starting values for the algorithm: {}
Optimization algorithm: hybrid Newton/BFGS with simple bounds [simple_bounds]
** Optimization: Newton with trust region for simple bounds
Iter. asc_train b_time lambda_travel_t b_cost asc_car Function Relgrad Radius Rho
0 -0.66 -1 1.6 -0.57 -0.29 8.9e+03 0.045 1 0.83 +
1 -0.57 -1.7 0.55 -0.94 -0.0006 8.7e+03 0.016 1 0.89 +
2 -0.48 -1.6 0.55 -0.78 0.14 8.6e+03 0.00087 10 0.97 ++
3 -0.48 -1.6 0.55 -0.78 0.14 8.6e+03 3.5e-06 10 1 ++
Optimization algorithm has converged.
Relative gradient: 3.525451557160213e-06
Cause of termination: Relative gradient = 3.5e-06 <= 6.1e-06
Number of function evaluations: 13
Number of gradient evaluations: 9
Number of hessian evaluations: 4
Algorithm: Newton with trust region for simple bound constraints
Number of iterations: 4
Proportion of Hessian calculation: 4/4 = 100.0%
Optimization time: 0:00:01.037716
Calculate second derivatives and BHHH
Number of estimated models.
print(f'A total of {len(dict_of_results)} models have been estimated')
A total of 3 models have been estimated
All estimation results
compiled_results, specs = compile_estimation_results(
dict_of_results, use_short_names=True
)
display('All estimated models')
display(compiled_results)
All estimated models
Model_000000 ... Model_000002
Number of estimated parameters 6 ... 5
Sample size 10719 ... 10719
Final log likelihood -8581.314 ... -8626.22
Akaike Information Criterion 17174.63 ... 17262.44
Bayesian Information Criterion 17218.31 ... 17298.84
asc_train (t-test) -0.459 (-10.1) ... -0.474 (-10.1)
b_time (t-test) -2.09 (-26.2) ... -1.62 (-28.9)
square_tt_coef (t-test) -0.0971 (-18.9) ...
cube_tt_coef (t-test) 0.000159 (6.09) ...
b_cost (t-test) -0.8 (-15.6) ... -0.787 (-15.7)
asc_car (t-test) 0.148 (4.56) ... 0.137 (4.04)
lambda_travel_time (t-test) ... 0.548 (8.77)
[12 rows x 3 columns]
Glossary
for short_name, spec in specs.items():
print(f'{short_name}\t{spec}')
Model_000000 train_tt_catalog:power
Model_000001 train_tt_catalog:linear
Model_000002 train_tt_catalog:boxcox
Estimation results of the Pareto optimal models.
pareto_results = pareto_optimal(dict_of_results)
compiled_pareto_results, pareto_specs = compile_estimation_results(
pareto_results, use_short_names=True
)
No Pareto file has been provided
display('Non dominated models')
display(compiled_pareto_results)
Non dominated models
Model_000000 ... Model_000002
Number of estimated parameters 5 ... 6
Sample size 10719 ... 10719
Final log likelihood -8626.22 ... -8581.314
Akaike Information Criterion 17262.44 ... 17174.63
Bayesian Information Criterion 17298.84 ... 17218.31
asc_train (t-test) -0.474 (-10.1) ... -0.459 (-10.1)
b_time (t-test) -1.62 (-28.9) ... -2.09 (-26.2)
lambda_travel_time (t-test) 0.548 (8.77) ...
b_cost (t-test) -0.787 (-15.7) ... -0.8 (-15.6)
asc_car (t-test) 0.137 (4.04) ... 0.148 (4.56)
square_tt_coef (t-test) ... -0.0971 (-18.9)
cube_tt_coef (t-test) ... 0.000159 (6.09)
[12 rows x 3 columns]
Glossary.
for short_name, spec in pareto_specs.items():
print(f'{short_name}\t{spec}')
Model_000000 train_tt_catalog:boxcox
Model_000001 train_tt_catalog:linear
Model_000002 train_tt_catalog:power
Total running time of the script: (0 minutes 4.301 seconds)