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 piecewiseFormula
import biogeme.biogeme as bio
from biogeme.expressions import Beta, log, Elem, bioNormalCdf
from optima import (
database,
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 = piecewiseFormula(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', 1, 1.0e-5, None, 1)
SIGMA_STAR_Envir02 = Beta('SIGMA_STAR_Envir02', 1, 1.0e-5, None, 0)
SIGMA_STAR_Envir03 = Beta('SIGMA_STAR_Envir03', 1, 1.0e-5, None, 0)
SIGMA_STAR_Mobil11 = Beta('SIGMA_STAR_Mobil11', 1, 1.0e-5, None, 0)
SIGMA_STAR_Mobil14 = Beta('SIGMA_STAR_Mobil14', 1, 1.0e-5, None, 0)
SIGMA_STAR_Mobil16 = Beta('SIGMA_STAR_Mobil16', 1, 1.0e-5, None, 0)
SIGMA_STAR_Mobil17 = Beta('SIGMA_STAR_Mobil17', 1, 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)
)
Create the Biogeme object
the_biogeme = bio.BIOGEME(database, loglike)
the_biogeme.modelName = 'b02one_latent_ordered'
File biogeme.toml has been parsed.
Estimate the parameters
results = the_biogeme.estimate()
*** Initial values of the parameters are obtained from the file __b02one_latent_ordered.iter
Cannot read file __b02one_latent_ordered.iter. Statement is ignored.
Optimization algorithm: hybrid Newton/BFGS with simple bounds [simple_bounds]
** Optimization: Newton with trust region for simple bounds
Iter. Function Relgrad Radius Rho
0 2.7e+04 1.1 0.5 -1.8e+304 -
1 2.4e+04 1.2 0.5 0.41 +
2 2.4e+04 1.2 0.25 -1.2e+304 -
3 1.9e+04 0.33 0.25 0.76 +
4 1.9e+04 0.33 0.12 -1.6 -
5 1.8e+04 0.22 0.12 0.51 +
6 1.8e+04 0.089 1.2 1.3 ++
7 1.8e+04 0.089 0.62 -6.6 -
8 1.8e+04 0.089 0.31 -0.2 -
9 1.8e+04 0.017 3.1 0.91 ++
10 1.8e+04 0.017 1.6 -1.6e+306 -
11 1.8e+04 0.017 0.78 -19 -
12 1.8e+04 0.017 0.39 -1.5 -
13 1.8e+04 0.031 0.39 0.63 +
14 1.8e+04 0.017 0.39 0.19 +
15 1.8e+04 0.013 0.39 0.39 +
16 1.8e+04 0.0023 3.9 1.1 ++
17 1.8e+04 0.00071 39 1.1 ++
18 1.8e+04 6e-06 39 1 ++
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.writeLaTeX()
print(f'LaTeX file: {results.data.latexFileName}')
Estimated betas: 34
final log likelihood: -17794.883
Output file: None
Results saved in file b02one_latent_ordered.tex
LaTeX file: b02one_latent_ordered.tex
results.getEstimatedParameters()
Total running time of the script: (0 minutes 3.966 seconds)