.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/programmers/plot_models.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_programmers_plot_models.py: biogeme.models ============== Examples of use of several functions. This is designed for programmers who need examples of use of the functions of the module. The examples are designed to illustrate the syntax. They do not correspond to any meaningful model. Michel Bierlaire Sun Jun 29 2025, 11:21:51 .. GENERATED FROM PYTHON SOURCE LINES 15-53 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np import pandas as pd from IPython.core.display_functions import display from biogeme.calculator import evaluate_simple_expression_per_row from biogeme.database import Database from biogeme.expressions import Beta, Variable from biogeme.models import ( boxcox, cnl, cnlmu, get_mev_for_cross_nested, get_mev_for_cross_nested_mu, get_mev_for_nested, get_mev_for_nested_mu, logit, loglogit, logmev_endogenous_sampling, lognested, lognested_mev_mu, mev_endogenous_sampling, nested, nested_mev_mu, piecewise_formula, piecewise_function, piecewise_variables, ) from biogeme.nests import ( NestsForCrossNestedLogit, NestsForNestedLogit, OneNestForCrossNestedLogit, OneNestForNestedLogit, ) from biogeme.second_derivatives import SecondDerivativesMode from biogeme.version import get_text .. GENERATED FROM PYTHON SOURCE LINES 54-55 Version of Biogeme. .. GENERATED FROM PYTHON SOURCE LINES 55-58 .. code-block:: Python print(get_text()) .. rst-class:: sphx-glr-script-out .. code-block:: none biogeme 3.3.1 [2025-09-03] Home page: http://biogeme.epfl.ch Submit questions to https://groups.google.com/d/forum/biogeme Michel Bierlaire, Transport and Mobility Laboratory, Ecole Polytechnique Fédérale de Lausanne (EPFL) .. GENERATED FROM PYTHON SOURCE LINES 59-61 Definition of a database ++++++++++++++++++++++++ .. GENERATED FROM PYTHON SOURCE LINES 61-75 .. code-block:: Python df = pd.DataFrame( { 'Person': [1, 1, 1, 2, 2], 'Exclude': [0, 0, 1, 0, 1], 'Variable1': [1, 2, 3, 4, 5], 'Variable2': [10, 20, 30, 40, 50], 'Choice': [2, 2, 3, 1, 2], 'Av1': [0, 1, 1, 1, 1], 'Av2': [1, 1, 0, 1, 1], 'Av3': [0, 1, 1, 1, 0], } ) display(df) .. rst-class:: sphx-glr-script-out .. code-block:: none Person Exclude Variable1 Variable2 Choice Av1 Av2 Av3 0 1 0 1 10 2 0 1 0 1 1 0 2 20 2 1 1 1 2 1 1 3 30 3 1 0 1 3 2 0 4 40 1 1 1 1 4 2 1 5 50 2 1 1 0 .. GENERATED FROM PYTHON SOURCE LINES 76-78 .. code-block:: Python my_data = Database('test', df) .. GENERATED FROM PYTHON SOURCE LINES 79-81 Piecewise linear specification ++++++++++++++++++++++++++++++ .. GENERATED FROM PYTHON SOURCE LINES 83-126 A piecewise linear specification (sometimes called 'spline') is a continuous but not differentiable function of the variable. It is defined based on thresholds. Between two thresholds, the function is linear. And the slope is changing after each threshold. Consider a variable :math:`t` and an interval :math:`[a, a+b]`. We define a new variable .. math:: x_{[a,b]}(t) = \max(0,\min(t-a,b)) = \left\{ \begin{array}{ll} 0 & \text{if } t < a, \\ t-a & \text{if } a \leq t < a+b, \\ b & \text{otherwise}. \end{array} \right. For each interval :math:`]-\infty,a]`, we have .. math:: x_{]-\infty,a]}(t) = \min(t,a) = \left\{ \begin{array}{ll} t & \text{if } t < a, \\ a & \text{otherwise}. \end{array} \right.. For each interval :math:`[a,+\infty[`, we have .. math:: x_{]-\infty,a]}(t) = \max(0,t-a) = \left\{ \begin{array}{ll} 0& \text{if } t < a, \\ t-a & \text{otherwise}. \end{array} \right.. If we consider a series of threshold .. math:: \alpha_1 < \alpha_2 < \ldots <\alpha_K, the piecewise linear transform of variable :math:`t` is .. math:: \sum_{k=1}^{K-1} \beta_k x_{[\alpha_k,\alpha_{k+1}]}, where :math:`\beta_k` is the slope of the linear function in interval :math:`[\alpha_k,\alpha_{k+1}]`. .. GENERATED FROM PYTHON SOURCE LINES 128-131 The next statement generates the variables, given the thresholds. A `None` is equivalent to :math:`\infty`, and can only appear first (and it means :math:`-\infty`) or last (and it means :math:`+\infty`). .. GENERATED FROM PYTHON SOURCE LINES 131-136 .. code-block:: Python x = Variable('x') thresholds = [None, 90, 180, 270, None] variables = piecewise_variables(x, thresholds) print(variables) .. rst-class:: sphx-glr-script-out .. code-block:: none [BinaryMin(, ), BinaryMax(, BinaryMin(( - ), )), BinaryMax(, BinaryMin(( - ), )), BinaryMax(, ( - ))] .. GENERATED FROM PYTHON SOURCE LINES 137-139 The next statement automatically generates the formula, including the Beta parameters, that are initialized to zero. .. GENERATED FROM PYTHON SOURCE LINES 139-142 .. code-block:: Python formula = piecewise_formula('x', thresholds) print(formula) .. rst-class:: sphx-glr-script-out .. code-block:: none MultipleSum((Beta('beta_x_minus_inf_90', 0, None, None, 0) * BinaryMin(x, `90.0`)), (Beta('beta_x_90_180', 0, None, None, 0) * BinaryMax(`0.0`, BinaryMin((x - `90.0`), `90.0`))), (Beta('beta_x_180_270', 0, None, None, 0) * BinaryMax(`0.0`, BinaryMin((x - `180.0`), `90.0`))), (Beta('beta_x_270_inf', 0, None, None, 0) * BinaryMax(`0.0`, (x - `270.0`)))) .. GENERATED FROM PYTHON SOURCE LINES 143-146 It is also possible to initialize the Beta parameters with other values. Note also that the first argument can be either the name of the variable (as in the previous call) or the variable itself. .. GENERATED FROM PYTHON SOURCE LINES 146-150 .. code-block:: Python betas = [-0.016806308, -0.010491137, -0.002012234, -0.020051303] formula = piecewise_formula(x, thresholds, betas) print(formula) .. rst-class:: sphx-glr-script-out .. code-block:: none MultipleSum((`-0.016806308` * BinaryMin(x, `90.0`)), (`-0.010491137` * BinaryMax(`0.0`, BinaryMin((x - `90.0`), `90.0`))), (`-0.002012234` * BinaryMax(`0.0`, BinaryMin((x - `180.0`), `90.0`))), (`-0.020051303` * BinaryMax(`0.0`, (x - `270.0`)))) .. GENERATED FROM PYTHON SOURCE LINES 151-152 We provide a plot of a piecewise linear specification. .. GENERATED FROM PYTHON SOURCE LINES 154-163 .. code-block:: Python X = np.arange(0, 300, 0.1) Y = [ piecewise_function( x, thresholds, [-0.016806308, -0.010491137, -0.002012234, -0.020051303] ) for x in X ] plt.plot(X, Y) .. image-sg:: /auto_examples/programmers/images/sphx_glr_plot_models_001.png :alt: plot models :srcset: /auto_examples/programmers/images/sphx_glr_plot_models_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 164-166 Logit +++++ .. GENERATED FROM PYTHON SOURCE LINES 168-171 .. code-block:: Python v = {1: Variable('Variable1'), 2: 0.1, 3: -0.1} av = {1: Variable('Av1'), 2: Variable('Av2'), 3: Variable('Av3')} .. GENERATED FROM PYTHON SOURCE LINES 172-174 Calculation of the (log of the) logit for the three alternatives, based on their availability. .. GENERATED FROM PYTHON SOURCE LINES 176-177 Alternative 1 .. GENERATED FROM PYTHON SOURCE LINES 177-187 .. code-block:: Python p1 = logit(v, av, 1) prob_1_value = evaluate_simple_expression_per_row( expression=p1, database=my_data, numerically_safe=False, second_derivatives_mode=SecondDerivativesMode.NEVER, use_jit=True, ) display(f'Probability of alternative 1: {prob_1_value}') .. rst-class:: sphx-glr-script-out .. code-block:: none Probability of alternative 1: [0. 0.78614804 0.95689275 0.9644926 0.99260846] .. GENERATED FROM PYTHON SOURCE LINES 188-189 Alternative 2 .. GENERATED FROM PYTHON SOURCE LINES 189-199 .. code-block:: Python p2 = logit(v, av, 2) prob_2_value = evaluate_simple_expression_per_row( expression=p2, database=my_data, numerically_safe=False, second_derivatives_mode=SecondDerivativesMode.NEVER, use_jit=True, ) display(f'Probability of alternative 2: {prob_2_value}') .. rst-class:: sphx-glr-script-out .. code-block:: none Probability of alternative 2: [1. 0.11758308 0. 0.01952317 0.00739154] .. GENERATED FROM PYTHON SOURCE LINES 200-201 Alternative 3 .. GENERATED FROM PYTHON SOURCE LINES 201-211 .. code-block:: Python p3 = logit(v, av, 3) prob_3_value = evaluate_simple_expression_per_row( expression=p3, database=my_data, numerically_safe=False, second_derivatives_mode=SecondDerivativesMode.NEVER, use_jit=True, ) display(f'Probability of alternative 3: {prob_3_value}') .. rst-class:: sphx-glr-script-out .. code-block:: none Probability of alternative 3: [0. 0.09626888 0.04310725 0.01598422 0. ] .. GENERATED FROM PYTHON SOURCE LINES 212-216 Calculation of the log of the logit for the three alternatives. If `av` is set to None, it means that all alternatives are always available. %% Alternative 1 .. GENERATED FROM PYTHON SOURCE LINES 216-226 .. code-block:: Python p1 = loglogit(util=v, av=None, i=1) log_prob_1_value = evaluate_simple_expression_per_row( expression=p1, database=my_data, numerically_safe=False, second_derivatives_mode=SecondDerivativesMode.NEVER, use_jit=True, ) display(f'Log probability of alternative 1: {log_prob_1_value}') .. rst-class:: sphx-glr-script-out .. code-block:: none Log probability of alternative 1: [-0.55356365 -0.24061016 -0.09537602 -0.03615312 -0.01345244] .. GENERATED FROM PYTHON SOURCE LINES 227-228 Alternative 2 .. GENERATED FROM PYTHON SOURCE LINES 228-238 .. code-block:: Python p2 = loglogit(util=v, av=None, i=2) log_prob_2_value = evaluate_simple_expression_per_row( expression=p2, database=my_data, numerically_safe=False, second_derivatives_mode=SecondDerivativesMode.NEVER, use_jit=True, ) display(f'Log probability of alternative 2: {log_prob_2_value}') .. rst-class:: sphx-glr-script-out .. code-block:: none Log probability of alternative 2: [-1.45356365 -2.14061016 -2.99537602 -3.93615312 -4.91345244] .. GENERATED FROM PYTHON SOURCE LINES 239-240 Alternative 3 .. GENERATED FROM PYTHON SOURCE LINES 240-250 .. code-block:: Python p3 = loglogit(util=v, av=None, i=3) log_prob_3_value = evaluate_simple_expression_per_row( expression=p3, database=my_data, numerically_safe=False, second_derivatives_mode=SecondDerivativesMode.NEVER, use_jit=True, ) display(f'Log probability of alternative 3: {log_prob_3_value}') .. rst-class:: sphx-glr-script-out .. code-block:: none Log probability of alternative 3: [-1.65356365 -2.34061016 -3.19537602 -4.13615312 -5.11345244] .. GENERATED FROM PYTHON SOURCE LINES 251-253 Box-Cox transform +++++++++++++++++ .. GENERATED FROM PYTHON SOURCE LINES 255-263 The Box-Cox transform of a variable :math:`x` is defined as .. math:: B(x,\ell) =\frac{x^{\ell}-1}{\ell}, where :math:`\ell` is a parameter that can be estimated from data. It has the property that .. math:: \lim_{\ell \to 0} B(x,\ell)=\log(x). .. GENERATED FROM PYTHON SOURCE LINES 263-266 .. code-block:: Python x = Variable('Variable1') display(boxcox(x, 4)) .. rst-class:: sphx-glr-script-out .. code-block:: none {{0:{{0:((Variable1**4.0 - `1.0`) / `4.0`), 1:(((log(Variable1) + (`4.0` * log(Variable1)**2.0)) + ((`16.0` * log(Variable1)**3.0) / `6.0`)) + ((`64.0` * log(Variable1)**4.0) / `24.0`))}[((`1e-05` > `4.0`) * ((-`1e-05`) < `4.0`))], 1:`0.0`}[(Variable1 == `0.0`)] .. GENERATED FROM PYTHON SOURCE LINES 267-270 .. code-block:: Python x = Variable('Variable1') display(boxcox(x, 0)) .. rst-class:: sphx-glr-script-out .. code-block:: none {{0:{{0:((Variable1**0.0 - `1.0`) / `0.0`), 1:(((log(Variable1) + (`0.0` * log(Variable1)**2.0)) + ((`0.0` * log(Variable1)**3.0) / `6.0`)) + ((`0.0` * log(Variable1)**4.0) / `24.0`))}[((`1e-05` > `0.0`) * ((-`1e-05`) < `0.0`))], 1:`0.0`}[(Variable1 == `0.0`)] .. GENERATED FROM PYTHON SOURCE LINES 271-275 .. code-block:: Python ell = Variable('Variable2') e = boxcox(x, ell) display(e) .. rst-class:: sphx-glr-script-out .. code-block:: none {{0:{{0:(((Variable1 ** Variable2) - `1.0`) / Variable2), 1:(((log(Variable1) + (Variable2 * log(Variable1)**2.0)) + ((Variable2**2.0 * log(Variable1)**3.0) / `6.0`)) + ((Variable2**3.0 * log(Variable1)**4.0) / `24.0`))}[((Variable2 < `1e-05`) * (Variable2 > (-`1e-05`)))], 1:`0.0`}[(Variable1 == `0.0`)] .. GENERATED FROM PYTHON SOURCE LINES 276-285 .. code-block:: Python boxcox_value = evaluate_simple_expression_per_row( expression=e, database=my_data, numerically_safe=False, second_derivatives_mode=SecondDerivativesMode.NEVER, use_jit=True, ) display(f'Box-Cox transform of Variable2: {boxcox_value}') .. rst-class:: sphx-glr-script-out .. code-block:: none Box-Cox transform of Variable2: [0.00000000e+00 5.24287500e+04 6.86303774e+12 3.02231455e+22 1.77635684e+33] .. GENERATED FROM PYTHON SOURCE LINES 286-288 We numerically illustrate that, when :math:`\lambda` goes to 0, the BoxCox transform of :math:`x` converges to the log of :math:`x`. .. GENERATED FROM PYTHON SOURCE LINES 290-296 .. code-block:: Python for ell in range(1, 16): x = 3 bc = boxcox(x, 10**-ell).get_value() print(f'ell=l0^(-{ell}): {bc:.6g} - {np.log(x):.6g} = {bc - np.log(x):.6g}') .. rst-class:: sphx-glr-script-out .. code-block:: none ell=l0^(-1): 1.16123 - 1.09861 = 0.0626195 ell=l0^(-2): 1.10467 - 1.09861 = 0.00605691 ell=l0^(-3): 1.09922 - 1.09861 = 0.000603696 ell=l0^(-4): 1.09867 - 1.09861 = 6.03497e-05 ell=l0^(-5): 1.09862 - 1.09861 = 6.03476e-06 ell=l0^(-6): 1.09861 - 1.09861 = 1.20695e-06 ell=l0^(-7): 1.09861 - 1.09861 = 1.20695e-07 ell=l0^(-8): 1.09861 - 1.09861 = 1.20695e-08 ell=l0^(-9): 1.09861 - 1.09861 = 1.20695e-09 ell=l0^(-10): 1.09861 - 1.09861 = 1.20695e-10 ell=l0^(-11): 1.09861 - 1.09861 = 1.20695e-11 ell=l0^(-12): 1.09861 - 1.09861 = 1.20703e-12 ell=l0^(-13): 1.09861 - 1.09861 = 1.20792e-13 ell=l0^(-14): 1.09861 - 1.09861 = 1.19904e-14 ell=l0^(-15): 1.09861 - 1.09861 = 1.11022e-15 .. GENERATED FROM PYTHON SOURCE LINES 297-299 MEV models ++++++++++ .. GENERATED FROM PYTHON SOURCE LINES 301-309 MEV models are defined as .. math:: \frac{e^{V_i + \ln G_i(e^{V_1},\ldots,e^{V_J})}}{\sum_j e^{V_j + \ln G_j(e^{V_1},\ldots,e^{V_J})}}, where :math:`G` is a generating function, and .. math:: G_i=\frac{\partial G}{\partial y_i}(e^{V_1},\ldots,e^{V_J}). .. GENERATED FROM PYTHON SOURCE LINES 311-323 **Nested logit model**: the :math:`G` function for the nested logit model is defined such that .. math:: G_i=\frac{\partial G}{\partial y_i}(e^{V_1},\ldots,e^{V_J}) = \mu e^{(\mu_m-1)V_i} \left(\sum_{i=1}^{J_m} e^{\mu_m V_i}\right)^{\frac{\mu}{\mu_m}-1}, where the choice set is partitioned into :math:`J_m` nests, each associated with a parameter :math:`\mu_m`, and :math:`\mu` is the scale parameter. The condition is :math:`0 \leq \mu \leq \mu_m` must be verified. In general, :math:`\mu` is normalized to 1.0. .. GENERATED FROM PYTHON SOURCE LINES 325-329 This is an example with 5 alternatives. Nest A contains alternatives 1, 2 and 4, and is associated with a scale parameter :math:`\mu_A=1.2`. Nest B contains alternatives 3 and 5, and is associated with a scale parameter :math:`\mu_B=2.3`. .. GENERATED FROM PYTHON SOURCE LINES 331-334 .. code-block:: Python v = {1: Variable('Variable1'), 2: 0.1, 3: -0.1, 4: -0.2, 5: 0.2} av = {1: 1, 2: 0, 3: 1, 4: 1, 5: 1} .. GENERATED FROM PYTHON SOURCE LINES 335-336 Definition of the nests. .. GENERATED FROM PYTHON SOURCE LINES 336-345 .. code-block:: Python nest_a = OneNestForNestedLogit( nest_param=1.2, list_of_alternatives=[1, 2, 4], name='nest_a' ) nest_b = OneNestForNestedLogit( nest_param=2.3, list_of_alternatives=[3, 5], name='name_b' ) nests = NestsForNestedLogit(choice_set=list(v), tuple_of_nests=(nest_a, nest_b)) .. GENERATED FROM PYTHON SOURCE LINES 346-356 .. code-block:: Python p1 = nested(v, availability=av, nests=nests, choice=1) p1_value = evaluate_simple_expression_per_row( expression=p1, database=my_data, numerically_safe=False, second_derivatives_mode=SecondDerivativesMode.NEVER, use_jit=True, ) display(f'Nested logit probability of alternative 1: {p1_value}') .. rst-class:: sphx-glr-script-out .. code-block:: none Nested logit probability of alternative 1: [0.55789028 0.78684631 0.9138115 0.96786857 0.98836318] .. GENERATED FROM PYTHON SOURCE LINES 357-358 If all the alternatives are available, define the availability dictionary as None. .. GENERATED FROM PYTHON SOURCE LINES 358-369 .. code-block:: Python p1_value = evaluate_simple_expression_per_row( expression=p1, database=my_data, numerically_safe=False, second_derivatives_mode=SecondDerivativesMode.NEVER, use_jit=True, ) display( f'Nested logit probability of alternative 1, all alternatives are available: {p1_value}' ) .. rst-class:: sphx-glr-script-out .. code-block:: none Nested logit probability of alternative 1, all alternatives are available: [0.55789028 0.78684631 0.9138115 0.96786857 0.98836318] .. GENERATED FROM PYTHON SOURCE LINES 370-371 The syntax is similar to obtain the log of the probability. .. GENERATED FROM PYTHON SOURCE LINES 371-381 .. code-block:: Python p2 = lognested(v, availability=av, nests=nests, choice=1) p2_value = evaluate_simple_expression_per_row( expression=p2, database=my_data, numerically_safe=False, second_derivatives_mode=SecondDerivativesMode.NEVER, use_jit=True, ) display(f'Nested logit log probability of alternative 1: {p2_value}') .. rst-class:: sphx-glr-script-out .. code-block:: none Nested logit log probability of alternative 1: [-0.58359296 -0.23972233 -0.09013096 -0.03265898 -0.01170506] .. GENERATED FROM PYTHON SOURCE LINES 382-394 .. code-block:: Python p2 = lognested(v, availability=None, nests=nests, choice=1) p2_value = evaluate_simple_expression_per_row( expression=p2, database=my_data, numerically_safe=False, second_derivatives_mode=SecondDerivativesMode.NEVER, use_jit=True, ) display( f'Nested logit log probability of alternative 1, all alternatives are available: {p2_value}' ) .. rst-class:: sphx-glr-script-out .. code-block:: none Nested logit log probability of alternative 1, all alternatives are available: [-0.7677835 -0.31935224 -0.11821542 -0.04163899 -0.01446801] .. GENERATED FROM PYTHON SOURCE LINES 395-399 If the value of the parameter :math:`\mu` is not 1, there is another function to call. Note that, for the sake of computational efficiency, it is not verified by the code if the condition :math:`0 \leq \mu \leq \mu_m` is valid. .. GENERATED FROM PYTHON SOURCE LINES 399-409 .. code-block:: Python p1 = nested_mev_mu(v, availability=av, nests=nests, choice=1, mu=1.1) p1_value = evaluate_simple_expression_per_row( expression=p1, database=my_data, numerically_safe=False, second_derivatives_mode=SecondDerivativesMode.NEVER, use_jit=True, ) display(f'Nested logit probability of alternative 1, mu=1.1: {p1_value}') .. rst-class:: sphx-glr-script-out .. code-block:: none Nested logit probability of alternative 1, mu=1.1: [0.57151598 0.80643395 0.92814743 0.9755475 0.9919295 ] .. GENERATED FROM PYTHON SOURCE LINES 410-420 .. code-block:: Python p1 = lognested_mev_mu(v, availability=av, nests=nests, choice=1, mu=1.1) p1_value = evaluate_simple_expression_per_row( expression=p1, database=my_data, numerically_safe=False, second_derivatives_mode=SecondDerivativesMode.NEVER, use_jit=True, ) display(f'Nested logit log probability of alternative 1, mu=1.1: {p1_value}') .. rst-class:: sphx-glr-script-out .. code-block:: none Nested logit log probability of alternative 1, mu=1.1: [-0.55946283 -0.21513329 -0.07456469 -0.02475643 -0.00810324] .. GENERATED FROM PYTHON SOURCE LINES 421-422 The validity of the nested structure can be verified. .. GENERATED FROM PYTHON SOURCE LINES 422-427 .. code-block:: Python nest_c = OneNestForNestedLogit(nest_param=2.3, list_of_alternatives=[3], name='name_c') nests = NestsForNestedLogit(choice_set=list(v), tuple_of_nests=(nest_a, nest_b)) is_valid, msg = nests.check_partition() .. GENERATED FROM PYTHON SOURCE LINES 428-430 .. code-block:: Python display(is_valid) .. rst-class:: sphx-glr-script-out .. code-block:: none True .. GENERATED FROM PYTHON SOURCE LINES 431-433 .. code-block:: Python print(msg) .. rst-class:: sphx-glr-script-out .. code-block:: none ; .. GENERATED FROM PYTHON SOURCE LINES 434-435 If an alternative belongs to two nests .. GENERATED FROM PYTHON SOURCE LINES 437-442 .. code-block:: Python nest_a = OneNestForNestedLogit(nest_param=1.2, list_of_alternatives=[1, 2, 3, 4]) nest_b = OneNestForNestedLogit(nest_param=2.3, list_of_alternatives=[3, 5]) nests = NestsForNestedLogit(choice_set=list(v), tuple_of_nests=(nest_a, nest_b)) is_valid, msg = nests.check_partition() .. GENERATED FROM PYTHON SOURCE LINES 443-445 .. code-block:: Python display(is_valid) .. rst-class:: sphx-glr-script-out .. code-block:: none False .. GENERATED FROM PYTHON SOURCE LINES 446-448 .. code-block:: Python print(msg) .. rst-class:: sphx-glr-script-out .. code-block:: none ; The following alternatives appear both in nests nest_1 and nest_2: {3}. .. GENERATED FROM PYTHON SOURCE LINES 449-465 **Cross-nested logit model**: the :math:`G` function for the cross nested logit model is defined such that .. math:: G_i=\frac{\partial G}{\partial y_i}(e^{V_1},\ldots,e^{V_J}) = \mu \sum_{m=1}^{M} \alpha_{im}^{\frac{\mu_m}{\mu}} e^{(\mu_m-1) V_i}\left( \sum_{j=1}^{J} \alpha_{jm}^{\frac{\mu_m}{\mu}} e^{\mu_m V_j} \right)^{\frac{\mu}{\mu_m}-1}, where each nest :math:`m` is associated with a parameter :math:`\mu_m` and, for each alternative :math:`i`, a parameter :math:`\alpha_{im} \geq 0` that captures the degree of membership of alternative :math:`i` to nest :math:`m`. :math:`\mu` is the scale parameter. For each alternative :math:`i`, there must be at least one nest :math:`m` such that :math:`\alpha_{im}>0`. The condition :math:`0 \leq \mu \leq \mu_m` must be also verified. In general, :math:`\mu` is normalized to 1.0. .. GENERATED FROM PYTHON SOURCE LINES 467-474 This is an example with 5 alternatives and two nests. - Alt. 1 belongs to nest A. - Alt. 2 belongs to nest A. - Alt. 3 belongs to both nest A (50%) and nest B (50%). - Alt. 4 belongs to nest B. - Alt. 5 belongs to nest B. .. GENERATED FROM PYTHON SOURCE LINES 476-488 .. code-block:: Python v = {1: Variable('Variable1'), 2: 0.1, 3: -0.1, 4: -0.2, 5: 0.2} av = {1: 1, 2: 0, 3: 1, 4: 1, 5: 1} alpha_a = {1: 1, 2: 1, 3: 0.5, 4: 0, 5: 0} alpha_b = {1: 0, 2: 0, 3: 0.5, 4: 1, 5: 1} nest_a = OneNestForCrossNestedLogit( nest_param=1.2, dict_of_alpha=alpha_a, name='Nest a' ) nest_b = OneNestForCrossNestedLogit( nest_param=2.3, dict_of_alpha=alpha_b, name='Nest b' ) nests = NestsForCrossNestedLogit(choice_set=list(v), tuple_of_nests=(nest_a, nest_b)) .. GENERATED FROM PYTHON SOURCE LINES 489-499 .. code-block:: Python p1 = cnl(v, availability=av, nests=nests, choice=1) p1_value = evaluate_simple_expression_per_row( expression=p1, database=my_data, numerically_safe=False, second_derivatives_mode=SecondDerivativesMode.NEVER, use_jit=True, ) display(f'Cross-Nested logit probability of alternative 1: {p1_value}') .. rst-class:: sphx-glr-script-out .. code-block:: none Cross-Nested logit probability of alternative 1: [0.60161076 0.81080413 0.92317655 0.97098919 0.98933903] .. GENERATED FROM PYTHON SOURCE LINES 500-501 If all the alternatives are available, define the availability dictionary as None. .. GENERATED FROM PYTHON SOURCE LINES 501-513 .. code-block:: Python p1 = cnl(v, availability=None, nests=nests, choice=1) p1_value = evaluate_simple_expression_per_row( expression=p1, database=my_data, numerically_safe=False, second_derivatives_mode=SecondDerivativesMode.NEVER, use_jit=True, ) display( f'Cross-Nested logit probability of alternative 1, all alternatives are available: {p1_value}' ) .. rst-class:: sphx-glr-script-out .. code-block:: none Cross-Nested logit probability of alternative 1, all alternatives are available: [0.49345928 0.74695583 0.89735316 0.9622809 0.98660661] .. GENERATED FROM PYTHON SOURCE LINES 514-518 If the value of the parameter :math:`\mu` is not 1, there is another function to call. Note that, for the sake of computational efficiency, it is not verified by the code if the condition :math:`0 \leq \mu \leq \mu_m` is verified. .. GENERATED FROM PYTHON SOURCE LINES 518-528 .. code-block:: Python p1 = cnlmu(v, availability=av, nests=nests, choice=1, mu=1.1) p1_value = evaluate_simple_expression_per_row( expression=p1, database=my_data, numerically_safe=False, second_derivatives_mode=SecondDerivativesMode.NEVER, use_jit=True, ) display(f'Cross-Nested logit probability of alternative 1, mu = 1.1: {p1_value}') .. rst-class:: sphx-glr-script-out .. code-block:: none Cross-Nested logit probability of alternative 1, mu = 1.1: [0.6110354 0.828654 0.93675704 0.97837752 0.99280536] .. GENERATED FROM PYTHON SOURCE LINES 529-534 If the sample is endogenous, a correction must be included in the model, as proposed by `Bierlaire, Bolduc and McFadden (2008) `_. In this case, the generating function must first be defined, and the MEV model with correction is then called. .. GENERATED FROM PYTHON SOURCE LINES 534-537 .. code-block:: Python log_gi = get_mev_for_cross_nested(v, availability=av, nests=nests) display(log_gi) .. rst-class:: sphx-glr-script-out .. code-block:: none {1: logzero(MultipleSum([((PowerConstant(, 1.2) * exp(( * ))) * PowerConstant(MultipleSum([(( * PowerConstant(, 1.2)) * exp(( * ))), (( * PowerConstant(, 1.2)) * exp()), (( * PowerConstant(, 1.2)) * exp()), (( * PowerConstant(, 1.2)) * exp()), (( * PowerConstant(, 1.2)) * exp())]), -0.16666666666666663)), ((PowerConstant(, 2.3) * exp(( * ))) * PowerConstant(MultipleSum([(( * PowerConstant(, 2.3)) * exp(( * ))), (( * PowerConstant(, 2.3)) * exp()), (( * PowerConstant(, 2.3)) * exp()), (( * PowerConstant(, 2.3)) * exp()), (( * PowerConstant(, 2.3)) * exp())]), -0.5652173913043478))])), 2: logzero(MultipleSum([((PowerConstant(, 1.2) * exp()) * PowerConstant(MultipleSum([(( * PowerConstant(, 1.2)) * exp(( * ))), (( * PowerConstant(, 1.2)) * exp()), (( * PowerConstant(, 1.2)) * exp()), (( * PowerConstant(, 1.2)) * exp()), (( * PowerConstant(, 1.2)) * exp())]), -0.16666666666666663)), ((PowerConstant(, 2.3) * exp()) * PowerConstant(MultipleSum([(( * PowerConstant(, 2.3)) * exp(( * ))), (( * PowerConstant(, 2.3)) * exp()), (( * PowerConstant(, 2.3)) * exp()), (( * PowerConstant(, 2.3)) * exp()), (( * PowerConstant(, 2.3)) * exp())]), -0.5652173913043478))])), 3: logzero(MultipleSum([((PowerConstant(, 1.2) * exp()) * PowerConstant(MultipleSum([(( * PowerConstant(, 1.2)) * exp(( * ))), (( * PowerConstant(, 1.2)) * exp()), (( * PowerConstant(, 1.2)) * exp()), (( * PowerConstant(, 1.2)) * exp()), (( * PowerConstant(, 1.2)) * exp())]), -0.16666666666666663)), ((PowerConstant(, 2.3) * exp()) * PowerConstant(MultipleSum([(( * PowerConstant(, 2.3)) * exp(( * ))), (( * PowerConstant(, 2.3)) * exp()), (( * PowerConstant(, 2.3)) * exp()), (( * PowerConstant(, 2.3)) * exp()), (( * PowerConstant(, 2.3)) * exp())]), -0.5652173913043478))])), 4: logzero(MultipleSum([((PowerConstant(, 1.2) * exp()) * PowerConstant(MultipleSum([(( * PowerConstant(, 1.2)) * exp(( * ))), (( * PowerConstant(, 1.2)) * exp()), (( * PowerConstant(, 1.2)) * exp()), (( * PowerConstant(, 1.2)) * exp()), (( * PowerConstant(, 1.2)) * exp())]), -0.16666666666666663)), ((PowerConstant(, 2.3) * exp()) * PowerConstant(MultipleSum([(( * PowerConstant(, 2.3)) * exp(( * ))), (( * PowerConstant(, 2.3)) * exp()), (( * PowerConstant(, 2.3)) * exp()), (( * PowerConstant(, 2.3)) * exp()), (( * PowerConstant(, 2.3)) * exp())]), -0.5652173913043478))])), 5: logzero(MultipleSum([((PowerConstant(, 1.2) * exp()) * PowerConstant(MultipleSum([(( * PowerConstant(, 1.2)) * exp(( * ))), (( * PowerConstant(, 1.2)) * exp()), (( * PowerConstant(, 1.2)) * exp()), (( * PowerConstant(, 1.2)) * exp()), (( * PowerConstant(, 1.2)) * exp())]), -0.16666666666666663)), ((PowerConstant(, 2.3) * exp()) * PowerConstant(MultipleSum([(( * PowerConstant(, 2.3)) * exp(( * ))), (( * PowerConstant(, 2.3)) * exp()), (( * PowerConstant(, 2.3)) * exp()), (( * PowerConstant(, 2.3)) * exp()), (( * PowerConstant(, 2.3)) * exp())]), -0.5652173913043478))]))} .. GENERATED FROM PYTHON SOURCE LINES 538-539 Assume the following correction factors .. GENERATED FROM PYTHON SOURCE LINES 539-550 .. code-block:: Python correction = {1: -0.1, 2: 0.1, 3: 0.2, 4: -0.2, 5: 0} p1 = mev_endogenous_sampling(v, log_gi, av, correction, choice=1) p1_value = evaluate_simple_expression_per_row( expression=p1, database=my_data, numerically_safe=False, second_derivatives_mode=SecondDerivativesMode.NEVER, use_jit=True, ) display(f'MEV model probability with correction: {p1_value}') .. rst-class:: sphx-glr-script-out .. code-block:: none MEV model probability with correction: [0.57460723 0.79415806 0.91584339 0.96822294 0.98835331] .. GENERATED FROM PYTHON SOURCE LINES 551-562 .. code-block:: Python correction = {1: -0.1, 2: 0.1, 3: 0.2, 4: -0.2, 5: 0} p1 = logmev_endogenous_sampling(v, log_gi, av, correction, choice=1) p1_value = evaluate_simple_expression_per_row( expression=p1, database=my_data, numerically_safe=False, second_derivatives_mode=SecondDerivativesMode.NEVER, use_jit=True, ) display(f'MEV model log probability with correction: {p1_value}') .. rst-class:: sphx-glr-script-out .. code-block:: none MEV model log probability with correction: [-0.55406855 -0.23047277 -0.0879099 -0.03229291 -0.01171504] .. GENERATED FROM PYTHON SOURCE LINES 563-564 The MEV generating function for the following models are available. .. GENERATED FROM PYTHON SOURCE LINES 566-567 Nested logit model .. GENERATED FROM PYTHON SOURCE LINES 567-576 .. code-block:: Python v = {1: Variable('Variable1'), 2: 0.1, 3: -0.1, 4: -0.2, 5: 0.2} nest_a = OneNestForNestedLogit( nest_param=Beta('muA', 1.2, 1.0, None, 0), list_of_alternatives=[1, 2, 4] ) nest_b = OneNestForNestedLogit( nest_param=Beta('muB', 2.3, 1.0, None, 0), list_of_alternatives=[3, 5] ) nests = NestsForNestedLogit(choice_set=list(v), tuple_of_nests=(nest_a, nest_b)) .. GENERATED FROM PYTHON SOURCE LINES 577-580 .. code-block:: Python log_gi = get_mev_for_nested(v, availability=None, nests=nests) display(log_gi) .. rst-class:: sphx-glr-script-out .. code-block:: none {1: ((( - ) * ) + ((( / ) - ) * log(MultipleSum([exp(( * )), exp(( * )), exp(( * ))])))), 2: ((( - ) * ) + ((( / ) - ) * log(MultipleSum([exp(( * )), exp(( * )), exp(( * ))])))), 4: ((( - ) * ) + ((( / ) - ) * log(MultipleSum([exp(( * )), exp(( * )), exp(( * ))])))), 3: ((( - ) * ) + ((( / ) - ) * log(MultipleSum([exp(( * )), exp(( * ))])))), 5: ((( - ) * ) + ((( / ) - ) * log(MultipleSum([exp(( * )), exp(( * ))]))))} .. GENERATED FROM PYTHON SOURCE LINES 581-582 And with the :math:`\mu` parameter. .. GENERATED FROM PYTHON SOURCE LINES 584-587 .. code-block:: Python log_gi = get_mev_for_nested_mu(v, availability=None, nests=nests, mu=1.1) display(log_gi) .. rst-class:: sphx-glr-script-out .. code-block:: none {1: ((log() + (( - ) * )) + ((( / ) - ) * log(MultipleSum([exp(( * )), exp(( * )), exp(( * ))])))), 2: ((log() + (( - ) * )) + ((( / ) - ) * log(MultipleSum([exp(( * )), exp(( * )), exp(( * ))])))), 4: ((log() + (( - ) * )) + ((( / ) - ) * log(MultipleSum([exp(( * )), exp(( * )), exp(( * ))])))), 3: ((log() + (( - ) * )) + ((( / ) - ) * log(MultipleSum([exp(( * )), exp(( * ))])))), 5: ((log() + (( - ) * )) + ((( / ) - ) * log(MultipleSum([exp(( * )), exp(( * ))]))))} .. GENERATED FROM PYTHON SOURCE LINES 588-589 Cross nested logit model .. GENERATED FROM PYTHON SOURCE LINES 591-602 .. code-block:: Python v = {1: Variable('Variable1'), 2: 0.1, 3: -0.1, 4: -0.2, 5: 0.2} av = {1: 1, 2: 0, 3: 1, 4: 1, 5: 1} alpha_a = {1: 1, 2: 1, 3: 0.5, 4: 0, 5: 0} alpha_b = {1: 0, 2: 0, 3: 0.5, 4: 1, 5: 1} nest_a = OneNestForCrossNestedLogit( nest_param=Beta('muA', 1.2, 1.0, None, 0), dict_of_alpha=alpha_a ) nest_b = OneNestForCrossNestedLogit( nest_param=Beta('muB', 2.3, 1.0, None, 0), dict_of_alpha=alpha_b ) nests = NestsForCrossNestedLogit(choice_set=list(v), tuple_of_nests=(nest_a, nest_b)) .. GENERATED FROM PYTHON SOURCE LINES 603-606 .. code-block:: Python log_gi = get_mev_for_cross_nested(v, availability=None, nests=nests) display(log_gi) .. rst-class:: sphx-glr-script-out .. code-block:: none {1: logzero(MultipleSum([((( ** ) * exp((( - ) * ))) * (MultipleSum([(( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * )))]) ** (( - ) / ))), ((( ** ) * exp((( - ) * ))) * (MultipleSum([(( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * )))]) ** (( - ) / )))])), 2: logzero(MultipleSum([((( ** ) * exp((( - ) * ))) * (MultipleSum([(( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * )))]) ** (( - ) / ))), ((( ** ) * exp((( - ) * ))) * (MultipleSum([(( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * )))]) ** (( - ) / )))])), 3: logzero(MultipleSum([((( ** ) * exp((( - ) * ))) * (MultipleSum([(( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * )))]) ** (( - ) / ))), ((( ** ) * exp((( - ) * ))) * (MultipleSum([(( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * )))]) ** (( - ) / )))])), 4: logzero(MultipleSum([((( ** ) * exp((( - ) * ))) * (MultipleSum([(( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * )))]) ** (( - ) / ))), ((( ** ) * exp((( - ) * ))) * (MultipleSum([(( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * )))]) ** (( - ) / )))])), 5: logzero(MultipleSum([((( ** ) * exp((( - ) * ))) * (MultipleSum([(( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * )))]) ** (( - ) / ))), ((( ** ) * exp((( - ) * ))) * (MultipleSum([(( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * ))), (( ** ) * exp(( * )))]) ** (( - ) / )))]))} .. GENERATED FROM PYTHON SOURCE LINES 607-608 Cross nested logit model with :math:`\mu` parameter. .. GENERATED FROM PYTHON SOURCE LINES 610-612 .. code-block:: Python log_gi = get_mev_for_cross_nested_mu(v, availability=None, nests=nests, mu=1.1) display(log_gi) .. rst-class:: sphx-glr-script-out .. code-block:: none {1: log(( * MultipleSum([((( ** ( / )) * exp((( - ) * ))) * (MultipleSum([(( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * )))]) ** (( / ) - ))), ((( ** ( / )) * exp((( - ) * ))) * (MultipleSum([(( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * )))]) ** (( / ) - )))]))), 2: log(( * MultipleSum([((( ** ( / )) * exp((( - ) * ))) * (MultipleSum([(( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * )))]) ** (( / ) - ))), ((( ** ( / )) * exp((( - ) * ))) * (MultipleSum([(( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * )))]) ** (( / ) - )))]))), 3: log(( * MultipleSum([((( ** ( / )) * exp((( - ) * ))) * (MultipleSum([(( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * )))]) ** (( / ) - ))), ((( ** ( / )) * exp((( - ) * ))) * (MultipleSum([(( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * )))]) ** (( / ) - )))]))), 4: log(( * MultipleSum([((( ** ( / )) * exp((( - ) * ))) * (MultipleSum([(( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * )))]) ** (( / ) - ))), ((( ** ( / )) * exp((( - ) * ))) * (MultipleSum([(( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * )))]) ** (( / ) - )))]))), 5: log(( * MultipleSum([((( ** ( / )) * exp((( - ) * ))) * (MultipleSum([(( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * )))]) ** (( / ) - ))), ((( ** ( / )) * exp((( - ) * ))) * (MultipleSum([(( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * ))), (( ** ( / )) * exp(( * )))]) ** (( / ) - )))])))} .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 2.203 seconds) .. _sphx_glr_download_auto_examples_programmers_plot_models.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_models.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_models.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_models.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_