Note
Go to the end to download the full example code.
Triangular mixture of logit
Example of a mixture of logit models, using Monte-Carlo integration. The mixing distribution is specified by the user. Here, a triangular distribution.
- author:
Michel Bierlaire, EPFL
- date:
Wed Apr 12 18:24:18 2023
import numpy as np
import biogeme.biogeme_logging as blog
import biogeme.biogeme as bio
from biogeme import models
from biogeme.expressions import Beta, bioDraws, log, MonteCarlo
from biogeme.native_draws import RandomNumberGeneratorTuple
from biogeme.parameters import Parameters
See the data processing script: Data preparation for Swissmetro.
from swissmetro_data import (
database,
CHOICE,
CAR_AV_SP,
TRAIN_AV_SP,
TRAIN_TT_SCALED,
TRAIN_COST_SCALED,
SM_TT_SCALED,
SM_COST_SCALED,
CAR_TT_SCALED,
CAR_CO_SCALED,
SM_AV,
)
logger = blog.get_screen_logger(level=blog.INFO)
logger.info('Example b25triangular_mixture.py')
Example b25triangular_mixture.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 with a triangular distribution, designed to be used for Monte-Carlo simulation. The triangular distribution is not directly available from Biogeme. The draws have to be generated by a function provided by the user.
Mean of the distribution.
B_TIME = Beta('B_TIME', 0, None, None, 0)
Scale of the distribution. It is advised not to use 0 as starting value for the following parameter.
B_TIME_S = Beta('B_TIME_S', 1, None, None, 0)
Function generating the draws.
def the_triangular_generator(sample_size: int, number_of_draws: int) -> np.ndarray:
"""
User-defined random number generator to the database.
See the numpy.random documentation to obtain a list of other distributions.
"""
return np.random.triangular(-1, 0, 1, (sample_size, number_of_draws))
Associate the function with a name.
myRandomNumberGenerators = {
'TRIANGULAR': RandomNumberGeneratorTuple(
generator=the_triangular_generator,
description='Draws from a triangular distribution',
)
}
Submit the generator to the database.
database.set_random_number_generators(myRandomNumberGenerators)
Define a random parameter with a triangular distribution, designed to be used for Monte-Carlo simulation.
B_TIME_RND = B_TIME + B_TIME_S * bioDraws('b_time_rnd', 'TRIANGULAR')
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 to b_time_rnd, we have a logit model (called the kernel)
prob = models.logit(V, av, CHOICE)
We integrate over b_time_rnd using Monte-Carlo
logprob = log(MonteCarlo(prob))
As the objective is to illustrate the syntax, we calculate the Monte-Carlo approximation with a small number of draws.
the_biogeme = bio.BIOGEME(database, logprob, number_of_draws=100, seed=1223)
the_biogeme.modelName = 'b25triangular_mixture'
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 __b25triangular_mixture.iter
Cannot read file __b25triangular_mixture.iter. Statement is ignored.
The number of draws (100) is low. The results may not be meaningful.
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.5e+03 0.092 1 0.39 +
1 -1 -1 -1 -1 2 5.5e+03 0.092 0.5 -0.57 -
2 -0.5 -1.5 -1.5 -0.82 1.9 5.4e+03 0.045 0.5 0.27 +
3 -0.54 -1 -1.2 -1.2 1.9 5.3e+03 0.035 0.5 0.87 +
4 -0.037 -0.92 -1.3 -1.3 2 5.3e+03 0.029 0.5 0.31 +
5 -0.073 -0.52 -1.1 -1.8 2 5.3e+03 0.026 0.5 0.41 +
6 -0.0021 -0.5 -1.4 -1.8 2.5 5.2e+03 0.018 0.5 0.4 +
7 -0.0021 -0.5 -1.4 -1.8 2.5 5.2e+03 0.018 0.25 0.032 -
8 -0.021 -0.5 -1.2 -1.8 2.6 5.2e+03 0.0049 0.25 0.76 +
9 -0.021 -0.5 -1.2 -1.8 2.6 5.2e+03 0.0049 0.12 -3 -
10 -0.021 -0.5 -1.2 -1.8 2.6 5.2e+03 0.0049 0.062 -1.9 -
11 0.042 -0.55 -1.2 -1.8 2.7 5.2e+03 0.006 0.062 0.13 +
12 0.016 -0.49 -1.2 -1.8 2.7 5.2e+03 0.0067 0.062 0.89 +
13 0.044 -0.49 -1.2 -1.9 2.8 5.2e+03 0.0052 0.62 1 ++
14 0.044 -0.49 -1.2 -1.9 2.8 5.2e+03 0.0052 0.31 -0.27 -
15 0.042 -0.52 -1.3 -2 3.1 5.2e+03 0.0077 0.31 0.6 +
16 0.14 -0.42 -1.3 -2 3.4 5.2e+03 0.014 0.31 0.18 +
17 0.14 -0.42 -1.3 -2 3.4 5.2e+03 0.014 0.16 -1.7 -
18 0.14 -0.42 -1.3 -2 3.4 5.2e+03 0.014 0.078 -0.11 -
19 0.11 -0.45 -1.2 -2.1 3.4 5.2e+03 0.0048 0.078 0.58 +
20 0.098 -0.43 -1.3 -2.1 3.4 5.2e+03 0.0028 0.078 0.7 +
21 0.098 -0.43 -1.3 -2.1 3.4 5.2e+03 0.0028 0.039 -0.078 -
22 0.1 -0.42 -1.2 -2.1 3.5 5.2e+03 0.002 0.039 0.63 +
23 0.1 -0.42 -1.2 -2.1 3.5 5.2e+03 0.002 0.02 -0.75 -
24 0.097 -0.44 -1.2 -2.1 3.5 5.2e+03 0.0029 0.02 0.34 +
25 0.1 -0.43 -1.3 -2.1 3.5 5.2e+03 0.0019 0.02 0.89 +
26 0.094 -0.42 -1.3 -2.1 3.5 5.2e+03 0.0018 0.02 0.61 +
27 0.1 -0.43 -1.3 -2.1 3.6 5.2e+03 0.002 0.02 0.9 +
28 0.1 -0.42 -1.3 -2.1 3.6 5.2e+03 0.0015 0.2 0.96 ++
29 0.1 -0.42 -1.3 -2.2 3.8 5.2e+03 0.0055 0.2 0.13 +
30 0.14 -0.4 -1.3 -2.2 4 5.2e+03 0.0025 0.2 0.44 +
31 0.14 -0.4 -1.3 -2.2 4 5.2e+03 0.0025 0.065 -8.9 -
32 0.14 -0.4 -1.3 -2.2 4 5.2e+03 0.0025 0.032 -0.68 -
33 0.14 -0.38 -1.3 -2.3 3.9 5.2e+03 0.00041 0.032 0.61 +
34 0.14 -0.38 -1.3 -2.3 3.9 5.2e+03 0.00041 0.016 -0.62 -
35 0.14 -0.4 -1.3 -2.3 3.9 5.2e+03 0.0016 0.016 0.22 +
36 0.14 -0.4 -1.3 -2.3 3.9 5.2e+03 0.0016 0.0081 0.01 -
37 0.13 -0.39 -1.3 -2.3 3.9 5.2e+03 0.00038 0.0081 0.65 +
38 0.13 -0.4 -1.3 -2.3 3.9 5.2e+03 0.00019 0.081 0.94 ++
39 0.13 -0.4 -1.3 -2.3 3.9 5.2e+03 0.00019 0.014 -0.8 -
40 0.13 -0.4 -1.3 -2.3 3.9 5.2e+03 0.00019 0.0071 -3.5 -
41 0.13 -0.4 -1.3 -2.3 3.9 5.2e+03 0.00019 0.0036 -1.5 -
42 0.13 -0.4 -1.3 -2.3 3.9 5.2e+03 0.00019 0.0018 -0.12 -
43 0.13 -0.4 -1.3 -2.3 3.9 5.2e+03 0.00014 0.0018 0.39 +
44 0.13 -0.4 -1.3 -2.3 3.9 5.2e+03 8.2e-05 0.0018 0.88 +
Results saved in file b25triangular_mixture.html
Results saved in file b25triangular_mixture.pickle
print(results.short_summary())
Results for model b25triangular_mixture
Nbr of parameters: 5
Sample size: 6768
Excluded data: 3960
Final log likelihood: -5215.85
Akaike Information Criterion: 10441.7
Bayesian Information Criterion: 10475.8
pandas_results = results.get_estimated_parameters()
pandas_results
Total running time of the script: (0 minutes 15.622 seconds)