Note
Go to the end to download the full example code.
Measurement equations: discrete indicators
Ordered probit.
- author:
Michel Bierlaire, EPFL
- date:
Thu Apr 13 16:48:55 2023
import biogeme.biogeme_logging as blog
from biogeme.models import piecewise_formula
import biogeme.biogeme as bio
from biogeme.expressions import Beta, log, Elem, bioNormalCdf
from biogeme.data.optima import (
read_data,
age_65_more,
ScaledIncome,
moreThanOneCar,
moreThanOneBike,
individualHouse,
male,
haveChildren,
haveGA,
highEducation,
Envir01,
Envir02,
Envir03,
Mobil11,
Mobil14,
Mobil16,
Mobil17,
)
logger = blog.get_screen_logger(level=blog.INFO)
logger.info('Example b02one_latent_ordered.py')
Example b02one_latent_ordered.py
Parameters to be estimated
coef_intercept = Beta('coef_intercept', 0.0, None, None, 0)
coef_age_65_more = Beta('coef_age_65_more', 0.0, None, None, 0)
coef_haveGA = Beta('coef_haveGA', 0.0, None, None, 0)
coef_moreThanOneCar = Beta('coef_moreThanOneCar', 0.0, None, None, 0)
coef_moreThanOneBike = Beta('coef_moreThanOneBike', 0.0, None, None, 0)
coef_individualHouse = Beta('coef_individualHouse', 0.0, None, None, 0)
coef_male = Beta('coef_male', 0.0, None, None, 0)
coef_haveChildren = Beta('coef_haveChildren', 0.0, None, None, 0)
coef_highEducation = Beta('coef_highEducation', 0.0, None, None, 0)
thresholds = [None, 4, 6, 8, 10, None]
formulaIncome = piecewise_formula(variable=ScaledIncome, thresholds=thresholds)
Latent variable: structural equation.
CARLOVERS = (
coef_intercept
+ coef_age_65_more * age_65_more
+ formulaIncome
+ coef_moreThanOneCar * moreThanOneCar
+ coef_moreThanOneBike * moreThanOneBike
+ coef_individualHouse * individualHouse
+ coef_male * male
+ coef_haveChildren * haveChildren
+ coef_haveGA * haveGA
+ coef_highEducation * highEducation
)
Measurement equations
Intercepts.
INTER_Envir01 = Beta('INTER_Envir01', 0, None, None, 1)
INTER_Envir02 = Beta('INTER_Envir02', 0, None, None, 0)
INTER_Envir03 = Beta('INTER_Envir03', 0, None, None, 0)
INTER_Mobil11 = Beta('INTER_Mobil11', 0, None, None, 0)
INTER_Mobil14 = Beta('INTER_Mobil14', 0, None, None, 0)
INTER_Mobil16 = Beta('INTER_Mobil16', 0, None, None, 0)
INTER_Mobil17 = Beta('INTER_Mobil17', 0, None, None, 0)
Coefficients.
B_Envir01_F1 = Beta('B_Envir01_F1', -1, None, None, 1)
B_Envir02_F1 = Beta('B_Envir02_F1', -1, None, None, 0)
B_Envir03_F1 = Beta('B_Envir03_F1', 1, None, None, 0)
B_Mobil11_F1 = Beta('B_Mobil11_F1', 1, None, None, 0)
B_Mobil14_F1 = Beta('B_Mobil14_F1', 1, None, None, 0)
B_Mobil16_F1 = Beta('B_Mobil16_F1', 1, None, None, 0)
B_Mobil17_F1 = Beta('B_Mobil17_F1', 1, None, None, 0)
Linear models.
MODEL_Envir01 = INTER_Envir01 + B_Envir01_F1 * CARLOVERS
MODEL_Envir02 = INTER_Envir02 + B_Envir02_F1 * CARLOVERS
MODEL_Envir03 = INTER_Envir03 + B_Envir03_F1 * CARLOVERS
MODEL_Mobil11 = INTER_Mobil11 + B_Mobil11_F1 * CARLOVERS
MODEL_Mobil14 = INTER_Mobil14 + B_Mobil14_F1 * CARLOVERS
MODEL_Mobil16 = INTER_Mobil16 + B_Mobil16_F1 * CARLOVERS
MODEL_Mobil17 = INTER_Mobil17 + B_Mobil17_F1 * CARLOVERS
Scale parameters.
SIGMA_STAR_Envir01 = Beta('SIGMA_STAR_Envir01', 10, 1.0e-5, None, 1)
SIGMA_STAR_Envir02 = Beta('SIGMA_STAR_Envir02', 10, 1.0e-5, None, 0)
SIGMA_STAR_Envir03 = Beta('SIGMA_STAR_Envir03', 10, 1.0e-5, None, 0)
SIGMA_STAR_Mobil11 = Beta('SIGMA_STAR_Mobil11', 10, 1.0e-5, None, 0)
SIGMA_STAR_Mobil14 = Beta('SIGMA_STAR_Mobil14', 10, 1.0e-5, None, 0)
SIGMA_STAR_Mobil16 = Beta('SIGMA_STAR_Mobil16', 10, 1.0e-5, None, 0)
SIGMA_STAR_Mobil17 = Beta('SIGMA_STAR_Mobil17', 10, 1.0e-5, None, 0)
Symmetric thresholds.
delta_1 = Beta('delta_1', 0.1, 1.0e-5, None, 0)
delta_2 = Beta('delta_2', 0.2, 1.0e-5, None, 0)
tau_1 = -delta_1 - delta_2
tau_2 = -delta_1
tau_3 = delta_1
tau_4 = delta_1 + delta_2
Ordered probit models.
Envir01_tau_1 = (tau_1 - MODEL_Envir01) / SIGMA_STAR_Envir01
Envir01_tau_2 = (tau_2 - MODEL_Envir01) / SIGMA_STAR_Envir01
Envir01_tau_3 = (tau_3 - MODEL_Envir01) / SIGMA_STAR_Envir01
Envir01_tau_4 = (tau_4 - MODEL_Envir01) / SIGMA_STAR_Envir01
IndEnvir01 = {
1: bioNormalCdf(Envir01_tau_1),
2: bioNormalCdf(Envir01_tau_2) - bioNormalCdf(Envir01_tau_1),
3: bioNormalCdf(Envir01_tau_3) - bioNormalCdf(Envir01_tau_2),
4: bioNormalCdf(Envir01_tau_4) - bioNormalCdf(Envir01_tau_3),
5: 1 - bioNormalCdf(Envir01_tau_4),
6: 1.0,
-1: 1.0,
-2: 1.0,
}
P_Envir01 = Elem(IndEnvir01, Envir01)
Envir02_tau_1 = (tau_1 - MODEL_Envir02) / SIGMA_STAR_Envir02
Envir02_tau_2 = (tau_2 - MODEL_Envir02) / SIGMA_STAR_Envir02
Envir02_tau_3 = (tau_3 - MODEL_Envir02) / SIGMA_STAR_Envir02
Envir02_tau_4 = (tau_4 - MODEL_Envir02) / SIGMA_STAR_Envir02
IndEnvir02 = {
1: bioNormalCdf(Envir02_tau_1),
2: bioNormalCdf(Envir02_tau_2) - bioNormalCdf(Envir02_tau_1),
3: bioNormalCdf(Envir02_tau_3) - bioNormalCdf(Envir02_tau_2),
4: bioNormalCdf(Envir02_tau_4) - bioNormalCdf(Envir02_tau_3),
5: 1 - bioNormalCdf(Envir02_tau_4),
6: 1.0,
-1: 1.0,
-2: 1.0,
}
P_Envir02 = Elem(IndEnvir02, Envir02)
Envir03_tau_1 = (tau_1 - MODEL_Envir03) / SIGMA_STAR_Envir03
Envir03_tau_2 = (tau_2 - MODEL_Envir03) / SIGMA_STAR_Envir03
Envir03_tau_3 = (tau_3 - MODEL_Envir03) / SIGMA_STAR_Envir03
Envir03_tau_4 = (tau_4 - MODEL_Envir03) / SIGMA_STAR_Envir03
IndEnvir03 = {
1: bioNormalCdf(Envir03_tau_1),
2: bioNormalCdf(Envir03_tau_2) - bioNormalCdf(Envir03_tau_1),
3: bioNormalCdf(Envir03_tau_3) - bioNormalCdf(Envir03_tau_2),
4: bioNormalCdf(Envir03_tau_4) - bioNormalCdf(Envir03_tau_3),
5: 1 - bioNormalCdf(Envir03_tau_4),
6: 1.0,
-1: 1.0,
-2: 1.0,
}
P_Envir03 = Elem(IndEnvir03, Envir03)
Mobil11_tau_1 = (tau_1 - MODEL_Mobil11) / SIGMA_STAR_Mobil11
Mobil11_tau_2 = (tau_2 - MODEL_Mobil11) / SIGMA_STAR_Mobil11
Mobil11_tau_3 = (tau_3 - MODEL_Mobil11) / SIGMA_STAR_Mobil11
Mobil11_tau_4 = (tau_4 - MODEL_Mobil11) / SIGMA_STAR_Mobil11
IndMobil11 = {
1: bioNormalCdf(Mobil11_tau_1),
2: bioNormalCdf(Mobil11_tau_2) - bioNormalCdf(Mobil11_tau_1),
3: bioNormalCdf(Mobil11_tau_3) - bioNormalCdf(Mobil11_tau_2),
4: bioNormalCdf(Mobil11_tau_4) - bioNormalCdf(Mobil11_tau_3),
5: 1 - bioNormalCdf(Mobil11_tau_4),
6: 1.0,
-1: 1.0,
-2: 1.0,
}
P_Mobil11 = Elem(IndMobil11, Mobil11)
Mobil14_tau_1 = (tau_1 - MODEL_Mobil14) / SIGMA_STAR_Mobil14
Mobil14_tau_2 = (tau_2 - MODEL_Mobil14) / SIGMA_STAR_Mobil14
Mobil14_tau_3 = (tau_3 - MODEL_Mobil14) / SIGMA_STAR_Mobil14
Mobil14_tau_4 = (tau_4 - MODEL_Mobil14) / SIGMA_STAR_Mobil14
IndMobil14 = {
1: bioNormalCdf(Mobil14_tau_1),
2: bioNormalCdf(Mobil14_tau_2) - bioNormalCdf(Mobil14_tau_1),
3: bioNormalCdf(Mobil14_tau_3) - bioNormalCdf(Mobil14_tau_2),
4: bioNormalCdf(Mobil14_tau_4) - bioNormalCdf(Mobil14_tau_3),
5: 1 - bioNormalCdf(Mobil14_tau_4),
6: 1.0,
-1: 1.0,
-2: 1.0,
}
P_Mobil14 = Elem(IndMobil14, Mobil14)
Mobil16_tau_1 = (tau_1 - MODEL_Mobil16) / SIGMA_STAR_Mobil16
Mobil16_tau_2 = (tau_2 - MODEL_Mobil16) / SIGMA_STAR_Mobil16
Mobil16_tau_3 = (tau_3 - MODEL_Mobil16) / SIGMA_STAR_Mobil16
Mobil16_tau_4 = (tau_4 - MODEL_Mobil16) / SIGMA_STAR_Mobil16
IndMobil16 = {
1: bioNormalCdf(Mobil16_tau_1),
2: bioNormalCdf(Mobil16_tau_2) - bioNormalCdf(Mobil16_tau_1),
3: bioNormalCdf(Mobil16_tau_3) - bioNormalCdf(Mobil16_tau_2),
4: bioNormalCdf(Mobil16_tau_4) - bioNormalCdf(Mobil16_tau_3),
5: 1 - bioNormalCdf(Mobil16_tau_4),
6: 1.0,
-1: 1.0,
-2: 1.0,
}
P_Mobil16 = Elem(IndMobil16, Mobil16)
Mobil17_tau_1 = (tau_1 - MODEL_Mobil17) / SIGMA_STAR_Mobil17
Mobil17_tau_2 = (tau_2 - MODEL_Mobil17) / SIGMA_STAR_Mobil17
Mobil17_tau_3 = (tau_3 - MODEL_Mobil17) / SIGMA_STAR_Mobil17
Mobil17_tau_4 = (tau_4 - MODEL_Mobil17) / SIGMA_STAR_Mobil17
IndMobil17 = {
1: bioNormalCdf(Mobil17_tau_1),
2: bioNormalCdf(Mobil17_tau_2) - bioNormalCdf(Mobil17_tau_1),
3: bioNormalCdf(Mobil17_tau_3) - bioNormalCdf(Mobil17_tau_2),
4: bioNormalCdf(Mobil17_tau_4) - bioNormalCdf(Mobil17_tau_3),
5: 1 - bioNormalCdf(Mobil17_tau_4),
6: 1.0,
-1: 1.0,
-2: 1.0,
}
P_Mobil17 = Elem(IndMobil17, Mobil17)
loglike = (
log(P_Envir01)
+ log(P_Envir02)
+ log(P_Envir03)
+ log(P_Mobil11)
+ log(P_Mobil14)
+ log(P_Mobil16)
+ log(P_Mobil17)
)
Read the data
database = read_data()
Create the Biogeme object
the_biogeme = bio.BIOGEME(database, loglike)
the_biogeme.modelName = 'b02one_latent_ordered'
Biogeme parameters read from biogeme.toml.
Estimate the parameters
results = the_biogeme.estimate()
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 __b02one_latent_ordered.iter
Cannot read file __b02one_latent_ordered.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. Function Relgrad Radius Rho
0 4e+04 0.34 10 1.3 ++
1 4e+04 0.34 5 1.3 -
2 4e+04 0.34 2.5 0.035 -
3 3.3e+04 0.17 25 1.1 ++
4 3.3e+04 0.17 12 1.1 -
5 3.3e+04 0.17 6.2 1.1 -
6 3.3e+04 0.17 3.1 1.1 -
7 3.3e+04 0.17 1.6 -0.26 -
8 2.5e+04 0.12 16 1.2 ++
9 2.5e+04 0.12 7.8 1.2 -
10 2.5e+04 0.12 3.9 1.2 -
11 2.5e+04 0.12 2 1.2 -
12 2.5e+04 0.12 0.98 -0.82 -
13 2.1e+04 0.1 9.8 1.3 ++
14 2.1e+04 0.1 4.9 1.3 -
15 2.1e+04 0.1 2.4 1.3 -
16 2e+04 0.044 2.4 0.49 +
17 2e+04 0.044 1.2 -0.6 -
18 1.8e+04 0.012 1.2 0.87 +
19 1.8e+04 0.014 1.2 0.9 +
20 1.8e+04 0.0024 12 1.2 ++
21 1.8e+04 0.0024 6.1 1.2 -
22 1.8e+04 0.0024 3.1 -14 -
23 1.8e+04 0.0024 1.5 -0.75 -
24 1.8e+04 0.0069 15 0.92 ++
25 1.8e+04 0.0069 7.6 -1e+02 -
26 1.8e+04 0.0069 3.8 -20 -
27 1.8e+04 0.0069 1.9 -1.3 -
28 1.8e+04 0.0012 19 0.96 ++
29 1.8e+04 0.0012 9.5 0.96 -
30 1.8e+04 0.0012 4.8 0.96 -
31 1.8e+04 0.0012 2.4 0.96 -
32 1.8e+04 0.0012 1.2 -50 -
33 1.8e+04 0.0012 0.6 -5.6 -
34 1.8e+04 0.0012 0.3 -0.16 -
35 1.8e+04 0.0015 0.3 0.83 +
36 1.8e+04 0.00092 3 0.99 ++
37 1.8e+04 0.00092 1.5 -61 -
38 1.8e+04 0.00092 0.75 -8.8 -
39 1.8e+04 0.00092 0.37 -0.62 -
40 1.8e+04 0.0019 0.37 0.49 +
41 1.8e+04 0.00026 3.7 1 ++
42 1.8e+04 0.00026 1.9 -6.9 -
43 1.8e+04 0.0013 1.9 0.59 +
44 1.8e+04 0.00024 19 0.98 ++
45 1.8e+04 0.00012 19 0.85 +
46 1.8e+04 3.4e-06 19 1 +
Results saved in file b02one_latent_ordered.html
Results saved in file b02one_latent_ordered.pickle
print(f'Estimated betas: {len(results.data.betaValues)}')
print(f'final log likelihood: {results.data.logLike:.3f}')
print(f'Output file: {results.data.htmlFileName}')
results.write_latex()
print(f'LaTeX file: {results.data.latexFileName}')
Estimated betas: 34
final log likelihood: -17794.883
Output file: b02one_latent_ordered.html
Results saved in file b02one_latent_ordered.tex
LaTeX file: b02one_latent_ordered.tex
results.get_estimated_parameters()
Total running time of the script: (0 minutes 6.222 seconds)