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).
- author:
Michel Bierlaire, EPFL
- date:
Thu Jul 13 21:31:54 2023
import biogeme.biogeme as bio
import biogeme.biogeme_logging as blog
from biogeme import models
from biogeme.expressions import Expression, Beta
from biogeme.models import boxcox
from biogeme.catalog import Catalog
from biogeme.results import compile_estimation_results, pareto_optimal
from biogeme.data.swissmetro import (
read_data,
CHOICE,
SM_AV,
CAR_AV_SP,
TRAIN_AV_SP,
TRAIN_TT_SCALED,
TRAIN_COST_SCALED,
SM_TT_SCALED,
SM_COST_SCALED,
CAR_TT_SCALED,
CAR_CO_SCALED,
)
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.
ell_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, ell_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, ell_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, ell_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.
V1 = ASC_TRAIN + B_TIME * train_tt_catalog + B_COST * TRAIN_COST_SCALED
V2 = B_TIME * sm_tt_catalog + B_COST * SM_COST_SCALED
V3 = ASC_CAR + B_TIME * car_tt_catalog + B_COST * CAR_CO_SCALED
Associate utility functions with the numbering of alternatives.
V = {1: V1, 2: V2, 3: V3}
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.
logprob = models.loglogit(V, av, CHOICE)
Read the data
database = read_data()
Create the Biogeme object.
the_biogeme = bio.BIOGEME(database, logprob)
the_biogeme.modelName = 'b02nonlinear'
the_biogeme.generate_html = False
the_biogeme.generate_pickle = False
Biogeme parameters read from biogeme.toml.
Estimate the parameters.
dict_of_results = the_biogeme.estimate_catalog()
Estimating 3 models.
Biogeme parameters read from biogeme.toml.
As the model is not too complex, we activate the calculation of second derivatives. If you want to change it, change the name of the algorithm in the TOML file from "automatic" to "simple_bounds"
*** Initial values of the parameters are obtained from the file __b02nonlinear_000000.iter
Cannot read file __b02nonlinear_000000.iter. Statement is ignored.
As the model is not too complex, we activate the calculation of second derivatives. If you want to change it, change the name of the algorithm in the TOML file from "automatic" to "simple_bounds"
Optimization algorithm: hybrid Newton/BFGS with simple bounds [simple_bounds]
** Optimization: Newton with trust region for simple bounds
Iter. ASC_CAR ASC_TRAIN B_COST B_TIME cube_tt_coef square_tt_coef Function Relgrad Radius Rho
0 0 0 0 0 0 0 1.1e+04 0.26 0.5 0.031 -
1 0.012 -0.5 -0.3 -0.5 0 0 9.2e+03 4 5 1 ++
2 0.012 -0.5 -0.3 -0.5 0 0 9.2e+03 4 2.5 -11 -
3 0.012 -0.5 -0.3 -0.5 0 0 9.2e+03 4 1.2 -8.7 -
4 0.012 -0.5 -0.3 -0.5 0 0 9.2e+03 4 0.62 -6.7 -
5 0.012 -0.5 -0.3 -0.5 0 0 9.2e+03 4 0.31 -2.4 -
6 0.012 -0.5 -0.3 -0.5 0 0 9.2e+03 4 0.16 -1.3 -
7 0.012 -0.5 -0.3 -0.5 0 0 9.2e+03 4 0.078 -1.2 -
8 0.012 -0.5 -0.3 -0.5 0 0 9.2e+03 4 0.039 -1.5 -
9 0.012 -0.5 -0.3 -0.5 0 0 9.2e+03 4 0.02 -2.3 -
10 0.012 -0.5 -0.3 -0.5 0 0 9.2e+03 4 0.0098 -3.2 -
11 0.012 -0.5 -0.3 -0.5 0 0 9.2e+03 4 0.0049 -4 -
12 0.012 -0.5 -0.3 -0.5 0 0 9.2e+03 4 0.0024 -4.6 -
13 0.012 -0.5 -0.3 -0.5 0 0 9.2e+03 4 0.0012 -2.5 -
14 0.012 -0.5 -0.3 -0.5 0 0 9.2e+03 4 0.00061 -1.3 -
15 0.012 -0.5 -0.3 -0.5 0 0 9.2e+03 4 0.00031 -0.28 -
16 0.012 -0.5 -0.3 -0.5 -0.00031 0.00031 9.2e+03 3 0.00031 0.64 +
17 0.012 -0.5 -0.3 -0.5 -0.00023 0.00061 9.2e+03 1.2 0.00031 0.81 +
18 0.012 -0.5 -0.3 -0.5 -0.00026 0.00092 9.2e+03 0.15 0.0031 0.99 ++
19 0.011 -0.5 -0.3 -0.5 -0.00026 0.004 9.2e+03 0.27 0.031 1 ++
20 0.0074 -0.52 -0.3 -0.53 -0.0004 0.034 9.1e+03 0.32 0.31 1 ++
21 -0.0046 -0.65 -0.43 -0.83 -0.0014 0.27 8.8e+03 2.4 0.31 0.6 +
22 0.022 -0.74 -0.73 -0.95 -0.0003 0.021 8.7e+03 8.2 0.31 0.68 +
23 -0.051 -0.75 -0.78 -1.3 -0.00023 -0.032 8.6e+03 8.2 3.1 0.96 ++
24 -0.051 -0.75 -0.78 -1.3 -0.00023 -0.032 8.6e+03 8.2 1.5 -87 -
25 -0.051 -0.75 -0.78 -1.3 -0.00023 -0.032 8.6e+03 8.2 0.76 -27 -
26 -0.051 -0.75 -0.78 -1.3 -0.00023 -0.032 8.6e+03 8.2 0.38 -4 -
27 -0.051 -0.75 -0.78 -1.3 -0.00023 -0.032 8.6e+03 8.2 0.19 -1.4 -
28 -0.051 -0.75 -0.78 -1.3 -0.00023 -0.032 8.6e+03 8.2 0.095 -1 -
29 -0.051 -0.75 -0.78 -1.3 -0.00023 -0.032 8.6e+03 8.2 0.048 -1.1 -
30 -0.051 -0.75 -0.78 -1.3 -0.00023 -0.032 8.6e+03 8.2 0.024 -1.3 -
31 -0.051 -0.75 -0.78 -1.3 -0.00023 -0.032 8.6e+03 8.2 0.012 -1.5 -
32 -0.051 -0.75 -0.78 -1.3 -0.00023 -0.032 8.6e+03 8.2 0.006 -1.4 -
33 -0.051 -0.75 -0.78 -1.3 -0.00023 -0.032 8.6e+03 8.2 0.003 -1.2 -
34 -0.051 -0.75 -0.78 -1.3 -0.00023 -0.032 8.6e+03 8.2 0.0015 -1.1 -
35 -0.051 -0.75 -0.78 -1.3 -0.00023 -0.032 8.6e+03 8.2 0.00075 -1 -
36 -0.051 -0.75 -0.78 -1.3 -0.00023 -0.032 8.6e+03 8.2 0.00037 -1 -
37 -0.051 -0.75 -0.78 -1.3 -0.00023 -0.032 8.6e+03 8.2 0.00019 -1 -
38 -0.051 -0.75 -0.78 -1.3 -0.00023 -0.032 8.6e+03 8.2 9.3e-05 -0.058 -
39 -0.051 -0.75 -0.78 -1.3 -0.00014 -0.032 8.6e+03 4.2 0.00093 1 ++
40 -0.051 -0.75 -0.78 -1.3 -0.00012 -0.031 8.6e+03 2.6 0.00093 0.67 +
41 -0.051 -0.75 -0.78 -1.3 -0.00013 -0.031 8.6e+03 0.094 0.0093 1 ++
42 -0.049 -0.75 -0.78 -1.3 -0.00015 -0.027 8.6e+03 0.0093 0.093 1 ++
43 -0.03 -0.73 -0.78 -1.4 -3.5e-05 -0.053 8.6e+03 0.29 0.93 0.93 ++
44 0.12 -0.49 -0.8 -2 0.00026 -0.12 8.6e+03 23 0.93 0.1 +
45 0.19 -0.38 -0.81 -2.2 0.00015 -0.094 8.6e+03 12 0.93 0.74 +
46 0.15 -0.46 -0.81 -2.1 0.00016 -0.097 8.6e+03 3.2 9.3 0.98 ++
47 0.15 -0.46 -0.8 -2.1 0.00016 -0.097 8.6e+03 0.031 93 1 ++
48 0.15 -0.46 -0.8 -2.1 0.00016 -0.097 8.6e+03 2.4e-05 93 1 ++
Biogeme parameters read from biogeme.toml.
As the model is not too complex, we activate the calculation of second derivatives. If you want to change it, change the name of the algorithm in the TOML file from "automatic" to "simple_bounds"
*** Initial values of the parameters are obtained from the file __b02nonlinear_000001.iter
Cannot read file __b02nonlinear_000001.iter. Statement is ignored.
As the model is not too complex, we activate the calculation of second derivatives. If you want to change it, change the name of the algorithm in the TOML file from "automatic" to "simple_bounds"
Optimization algorithm: hybrid Newton/BFGS with simple bounds [simple_bounds]
** Optimization: Newton with trust region for simple bounds
Iter. ASC_CAR ASC_TRAIN B_COST B_TIME Function Relgrad Radius Rho
0 0.0067 -0.64 -0.71 -1.1 8.7e+03 0.028 1 0.77 +
1 0.013 -0.65 -0.78 -1.3 8.7e+03 0.0017 10 1 ++
2 0.013 -0.65 -0.78 -1.3 8.7e+03 6.5e-06 10 1 ++
Biogeme parameters read from biogeme.toml.
As the model is not too complex, we activate the calculation of second derivatives. If you want to change it, change the name of the algorithm in the TOML file from "automatic" to "simple_bounds"
*** Initial values of the parameters are obtained from the file __b02nonlinear_000002.iter
Cannot read file __b02nonlinear_000002.iter. Statement is ignored.
As the model is not too complex, we activate the calculation of second derivatives. If you want to change it, change the name of the algorithm in the TOML file from "automatic" to "simple_bounds"
Optimization algorithm: hybrid Newton/BFGS with simple bounds [simple_bounds]
** Optimization: Newton with trust region for simple bounds
Iter. ASC_CAR ASC_TRAIN B_COST B_TIME lambda_travel_t Function Relgrad Radius Rho
0 0.2 -0.37 -0.8 -1.8 0.39 8.6e+03 0.0014 1 0.65 +
1 0.14 -0.47 -0.79 -1.6 0.52 8.6e+03 0.0011 10 1.1 ++
2 0.14 -0.47 -0.79 -1.6 0.52 8.6e+03 2.9e-05 10 1 ++
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
)
compiled_results
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
compiled_pareto_results
Glossary.
for short_name, spec in pareto_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
Total running time of the script: (0 minutes 1.902 seconds)