Mixture with lognormal distribution

Example of a mixture of logit models. The mixing distribution is distributed as a log normal. Compared to Mixture with lognormal distribution, the integration is performed using numerical integration instead of Monte-Carlo approximation.

author:

Michel Bierlaire, EPFL

date:

Mon Apr 10 12:13:23 2023

import biogeme.biogeme_logging as blog
import biogeme.biogeme as bio
import biogeme.distributions as dist
from biogeme import models
from biogeme.expressions import (
    Beta,
    RandomVariable,
    exp,
    log,
    Integrate,
)

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 b17lognormal_mixture_integral.py')

Example b17lognormal_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 Monte-Carlo simulation

B_TIME = Beta('B_TIME', 0, None, None, 0)

It is advised not to use 0 as starting value for the following parameter..

B_TIME_S = Beta('B_TIME_S', 1, -2, 2, 0)

Define a random parameter, log normally distributed, designed to be used for numerical integration.

omega = RandomVariable('omega')
B_TIME_RND = -exp(B_TIME + B_TIME_S * omega)
density = dist.normalpdf(omega)

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 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 * density, 'omega'))

Create the Biogeme object.

the_biogeme = bio.BIOGEME(database, logprob)
the_biogeme.modelName = 'b17lognormal_mixture_integral'

File biogeme.toml has been parsed.

Estimate the parameters

results = the_biogeme.estimate()

*** Initial values of the parameters are obtained from the file __b17lognormal_mixture_integral.iter
Cannot read file __b17lognormal_mixture_integral.iter. Statement is ignored.
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        B_TIME_S     Function    Relgrad   Radius      Rho
    0            0.18            -0.4              -1            0.36               1      5.3e+03      0.018       10        1   ++
    1            0.18           -0.34            -1.3            0.57             1.1      5.2e+03     0.0026    1e+02      1.1   ++
    2            0.17           -0.34            -1.4            0.58             1.2      5.2e+03    0.00017    1e+03        1   ++
    3            0.17           -0.34            -1.4            0.58             1.2      5.2e+03    2.9e-06    1e+03        1   ++
Results saved in file b17lognormal_mixture_integral.html
Results saved in file b17lognormal_mixture_integral.pickle

print(results.short_summary())

Results for model b17lognormal_mixture_integral
Nbr of parameters:              5
Sample size:                    6768
Excluded data:                  3960
Final log likelihood:           -5231.419
Akaike Information Criterion:   10472.84
Bayesian Information Criterion: 10506.94

pandas_results = results.getEstimatedParameters()
pandas_results

	Value	Rob. Std err	Rob. t-test	Rob. p-value
ASC_CAR	0.174559	0.062655	2.786034	5.335733e-03
ASC_TRAIN	-0.345857	0.073252	-4.721472	2.341435e-06
B_COST	-1.380475	0.097826	-14.111584	0.000000e+00
B_TIME	0.575771	0.071228	8.083497	6.661338e-16
B_TIME_S	1.239158	0.132469	9.354327	0.000000e+00