Note
Go to the end to download the full example code.
Mixture of logit models
Example of a mixture of logit models, using numerical integration. The mixing distribution is uniform.
- author:
Michel Bierlaire, EPFL
- date:
Sun Apr 9 17:52:52 2023
import biogeme.biogeme_logging as blog
import biogeme.biogeme as bio
from biogeme import models
from biogeme.expressions import (
Beta,
Integrate,
RandomVariable,
exp,
log,
)
See the data processing script: Data preparation for Swissmetro.
from swissmetro_data import (
database,
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)
logger.info('Example b06unif_mixture_integral.py')
Example b06unif_mixture_integral.py
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_COST = Beta('B_COST', 0, None, None, 0)
Define a random parameter, normally distributed, designed to be used for numerical integration
B_TIME = Beta('B_TIME', 0, None, None, 0)
B_TIME_S = Beta('B_TIME_S', 1, None, None, 0)
omega = RandomVariable('omega')
As the numerical integration ranges from -∞ to +∞, we need to perform a change of variable in order to integrate between -1 and 1.
LOWER_BND = -1
UPPER_BND = 1
x = LOWER_BND + (UPPER_BND - LOWER_BND) / (1 + exp(-omega))
dx = (UPPER_BND - LOWER_BND) * exp(-omega) * (1 + exp(-omega)) ** (-2)
B_TIME_RND = B_TIME + B_TIME_S * x
Definition of the utility functions.
V1 = ASC_TRAIN + B_TIME_RND * TRAIN_TT_SCALED + B_COST * TRAIN_COST_SCALED
V2 = ASC_SM + B_TIME_RND * SM_TT_SCALED + B_COST * SM_COST_SCALED
V3 = ASC_CAR + B_TIME_RND * CAR_TT_SCALED + 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}
Conditional on omega, we have a logit model (called the kernel).
condprob = models.logit(V, av, CHOICE)
We integrate over omega using numerical integration.
logprob = log(Integrate(condprob * dx / (UPPER_BND - LOWER_BND), 'omega'))
Create the Biogeme object.
the_biogeme = bio.BIOGEME(database, logprob)
the_biogeme.modelName = '06unif_mixture_integral'
Biogeme parameters read from biogeme.toml.
Estimate the parameters
results = the_biogeme.estimate()
As the model is rather complex, we cancel the calculation of second derivatives. If you want to control the parameters, 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 __06unif_mixture_integral.iter
Cannot read file __06unif_mixture_integral.iter. Statement is ignored.
As the model is rather complex, we cancel the calculation of second derivatives. If you want to control the parameters, 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: BFGS with trust region for simple bounds
Iter. ASC_CAR ASC_TRAIN B_COST B_TIME B_TIME_S Function Relgrad Radius Rho
0 -1 -1 -1 -1 2 5.6e+03 0.09 1 0.39 +
1 -1 -1 -1 -1 2 5.6e+03 0.09 0.5 -0.11 -
2 -0.5 -1.5 -1.5 -1.3 1.5 5.4e+03 0.074 0.5 0.32 +
3 -0.5 -1.5 -1.5 -1.3 1.5 5.4e+03 0.074 0.25 -0.18 -
4 -0.25 -1.2 -1.2 -1 1.8 5.3e+03 0.028 0.25 0.47 +
5 -0.5 -1 -1 -1.3 1.5 5.3e+03 0.04 0.25 0.16 +
6 -0.5 -1 -1 -1.3 1.5 5.3e+03 0.04 0.12 -0.058 -
7 -0.38 -0.88 -1.1 -1.2 1.5 5.3e+03 0.024 0.12 0.63 +
8 -0.25 -1 -1.2 -1.3 1.4 5.3e+03 0.017 0.12 0.33 +
9 -0.24 -0.88 -1.1 -1.4 1.5 5.3e+03 0.0093 0.12 0.88 +
10 -0.11 -0.75 -1.3 -1.5 1.7 5.2e+03 0.012 0.12 0.71 +
11 -0.079 -0.68 -1.1 -1.6 1.7 5.2e+03 0.0073 0.12 0.71 +
12 -0.043 -0.56 -1.2 -1.7 1.9 5.2e+03 0.0094 1.2 0.91 ++
13 -0.043 -0.56 -1.2 -1.7 1.9 5.2e+03 0.0094 0.62 -2.1 -
14 -0.043 -0.56 -1.2 -1.7 1.9 5.2e+03 0.0094 0.31 -1.5 -
15 -0.043 -0.56 -1.2 -1.7 1.9 5.2e+03 0.0094 0.16 -0.15 -
16 0.05 -0.59 -1.2 -1.9 2 5.2e+03 0.01 0.16 0.38 +
17 -0.011 -0.52 -1.2 -1.9 2.2 5.2e+03 0.0047 0.16 0.81 +
18 -0.011 -0.52 -1.2 -1.9 2.2 5.2e+03 0.0047 0.078 -0.45 -
19 0.067 -0.51 -1.2 -1.9 2.2 5.2e+03 0.008 0.078 0.3 +
20 0.039 -0.48 -1.2 -2 2.3 5.2e+03 0.0038 0.78 0.94 ++
21 0.039 -0.48 -1.2 -2 2.3 5.2e+03 0.0038 0.39 -1.4 -
22 0.039 -0.48 -1.2 -2 2.3 5.2e+03 0.0038 0.2 -0.89 -
23 0.039 -0.48 -1.2 -2 2.3 5.2e+03 0.0038 0.098 0.017 -
24 0.062 -0.5 -1.3 -2.1 2.4 5.2e+03 0.0066 0.098 0.47 +
25 0.095 -0.45 -1.2 -2.1 2.5 5.2e+03 0.0047 0.098 0.83 +
26 0.11 -0.44 -1.2 -2.2 2.6 5.2e+03 0.0032 0.98 0.94 ++
27 0.11 -0.44 -1.2 -2.2 2.6 5.2e+03 0.0032 0.49 -22 -
28 0.11 -0.44 -1.2 -2.2 2.6 5.2e+03 0.0032 0.24 -3.8 -
29 0.11 -0.44 -1.2 -2.2 2.6 5.2e+03 0.0032 0.12 -0.74 -
30 0.098 -0.41 -1.3 -2.2 2.7 5.2e+03 0.0028 0.12 0.23 +
31 0.098 -0.41 -1.3 -2.2 2.7 5.2e+03 0.0028 0.061 -0.22 -
32 0.12 -0.41 -1.3 -2.2 2.7 5.2e+03 0.00093 0.061 0.77 +
33 0.12 -0.41 -1.3 -2.2 2.7 5.2e+03 0.00093 0.031 -0.029 -
34 0.13 -0.39 -1.3 -2.3 2.8 5.2e+03 0.0015 0.031 0.61 +
35 0.13 -0.4 -1.3 -2.3 2.8 5.2e+03 0.00073 0.31 1 ++
36 0.13 -0.4 -1.3 -2.3 2.8 5.2e+03 0.00073 0.15 -4.9 -
37 0.13 -0.4 -1.3 -2.3 2.8 5.2e+03 0.00073 0.076 -0.56 -
38 0.13 -0.4 -1.3 -2.3 2.8 5.2e+03 0.00073 0.038 -0.47 -
39 0.13 -0.4 -1.3 -2.3 2.8 5.2e+03 0.00073 0.019 0.00081 -
40 0.14 -0.39 -1.3 -2.3 2.8 5.2e+03 0.0016 0.019 0.43 +
41 0.14 -0.39 -1.3 -2.3 2.8 5.2e+03 0.00058 0.19 0.92 ++
42 0.14 -0.39 -1.3 -2.3 2.8 5.2e+03 0.00058 0.095 -3.4 -
43 0.14 -0.39 -1.3 -2.3 2.8 5.2e+03 0.00058 0.048 -0.51 -
44 0.14 -0.4 -1.3 -2.3 2.9 5.2e+03 0.00047 0.048 0.17 +
45 0.14 -0.4 -1.3 -2.3 2.9 5.2e+03 0.00047 0.015 -3.6 -
46 0.15 -0.38 -1.3 -2.3 2.9 5.2e+03 0.00032 0.015 0.73 +
47 0.15 -0.38 -1.3 -2.3 2.9 5.2e+03 0.00032 0.0076 -0.71 -
48 0.15 -0.38 -1.3 -2.3 2.9 5.2e+03 0.00032 0.0038 0.07 -
49 0.14 -0.38 -1.3 -2.3 2.9 5.2e+03 0.00032 0.0038 0.39 +
50 0.14 -0.39 -1.3 -2.3 2.9 5.2e+03 0.00019 0.0038 0.47 +
51 0.14 -0.39 -1.3 -2.3 2.9 5.2e+03 8.4e-05 0.0038 0.7 +
Results saved in file 06unif_mixture_integral.html
Results saved in file 06unif_mixture_integral.pickle
print(results.short_summary())
Results for model 06unif_mixture_integral
Nbr of parameters: 5
Sample size: 6768
Excluded data: 3960
Final log likelihood: -5215.073
Akaike Information Criterion: 10440.15
Bayesian Information Criterion: 10474.24
pandas_results = results.get_estimated_parameters()
pandas_results
Total running time of the script: (0 minutes 27.385 seconds)