<jupyter_start><jupyter_text>An introduction to Biogeme  Biogeme Basics: Logit Model<jupyter_code>import pandas  as pd
import numpy as np
import biogeme.database as db
import biogeme.biogeme as bio
import biogeme.models as models
import biogeme.expressions as exp
import biogeme.tools as tools
import seaborn as sns
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>**Import Swissmetro data**<jupyter_code>pandas = pd.read_csv("../../Data/8-DiscreteChoiceModels/swissmetro.dat",sep='\t')
database = db.Database("data/swissmetro", pandas)<jupyter_output><empty_output><jupyter_text>Let's see what this dataset has * dataset consists of survey data collected on the trains between St. Gallen and Geneva, Switzerland, during March 1998* It is necessary to obtain data from surveys of hypothetical markets/situations, which include the innovation, to assess the impact. * Survey data were collected on rail-based travels, interviewing 470 respondents. Due to data problems, only 441 are used here. A similar method for relevant car trips. A total of 1070 persons filled in the survey completely and were willing to participate in the second SP survey, which was generated using the same approach used for the rail interviews. 750 usable SP surveys were returned, from the license-plate based survey.* Nine stated choice situations were generated for each the respondents, offering three alternatives: rail, Swissmetro and carBierlaire, M., <NAME>. and <NAME>. (2001), The acceptance of modal innovation: The case of Swissmetro, in ‘Proceedings of the Swiss Transport Research Conference’, Ascona, Switzerland. ![](img/swissmetro_var1.png) ![](img/swissmetro_var2.png)<jupyter_code>plt.figure(figsize=(10,5))
chart = sns.countplot(pandas['PURPOSE'])
chart.set_xticklabels(chart.get_xticklabels(), rotation=30, horizontalalignment='right');
chart.set_xticklabels(['Commuter', 'Shopping', 'Business', 'Leisure', 'Return from work','Return from shopping', 'Return from business','Return from leisure','Other']);<jupyter_output>C:\code\anaconda3\lib\site-packages\seaborn\_decorators.py:36: FutureWarning: Pass the following variable as a keyword arg: x. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.
  warnings.warn(<jupyter_text>**Use collumn names as variables**<jupyter_code>globals().update(database.variables)<jupyter_output><empty_output><jupyter_text>**Exclude some unwanted entries**<jupyter_code>exclude = (( PURPOSE != 1 ) * ( PURPOSE != 3 ) + ( CHOICE == 0 )) > 0

database.remove(exclude)<jupyter_output><empty_output><jupyter_text>**Define some dummy variables**<jupyter_code>SM_COST = SM_CO * ( GA == 0 )
TRAIN_COST = TRAIN_CO * ( GA == 0 )

CAR_AV_SP = exp.DefineVariable ('CAR_AV_SP', CAR_AV * ( SP !=0 ), database)
TRAIN_AV_SP = exp.DefineVariable ('TRAIN_AV_SP', TRAIN_AV * ( SP != 0 ), database)<jupyter_output><empty_output><jupyter_text>**Rescale some data**<jupyter_code>TRAIN_TT_SCALED   = exp.DefineVariable('TRAIN_TT_SCALED',   TRAIN_TT / 100.0, database)
TRAIN_COST_SCALED = exp.DefineVariable('TRAIN_COST_SCALED', TRAIN_COST / 100, database)
SM_TT_SCALED      = exp.DefineVariable('SM_TT_SCALED',      SM_TT / 100.0   , database)
SM_COST_SCALED    = exp.DefineVariable('SM_COST_SCALED',    SM_COST / 100   , database)
CAR_TT_SCALED     = exp.DefineVariable('CAR_TT_SCALED',     CAR_TT / 100    , database)
CAR_CO_SCALED     = exp.DefineVariable('CAR_CO_SCALED',     CAR_CO / 100    , database)
pandas = database.data
plt.figure(figsize=(10,5))
chart = sns.countplot(pandas['PURPOSE'])
chart.set_xticklabels(chart.get_xticklabels(), rotation=30, horizontalalignment='right');
chart.set_xticklabels(['Commuter', 'Business']);
plt.figure(1, figsize=(10,5))
chart = sns.countplot(pandas['TRAIN_AV'])
chart.set_xticklabels(chart.get_xticklabels(), rotation=30, horizontalalignment='right');
chart.set_xticklabels(['Yes']);
chart.set(title='Train Availability', xlabel="Available", ylabel = "Count");

plt.figure(2, figsize=(10,5))
chart = sns.countplot(pandas['CAR_AV'])
chart.set_xticklabels(chart.get_xticklabels(), rotation=30, horizontalalignment='right');
chart.set_xticklabels(['No', 'Yes']);
chart.set(title='Car Availability', xlabel="Available", ylabel = "Count");

plt.figure(3, figsize=(10,5))
chart = sns.countplot(pandas['SM_AV'])
chart.set_xticklabels(chart.get_xticklabels(), rotation=30, horizontalalignment='right');
chart.set_xticklabels(['Yes']);
chart.set(title='Swissmetro Availability', xlabel="Available", ylabel = "Count");
pd.set_option('display.max_rows', 500)
pd.set_option('display.max_columns', 500)
pd.set_option('display.width', 1000)

pandas.describe()<jupyter_output><empty_output><jupyter_text>**Define the utility functions** \begin{align}V_1 & = \beta_{Train} + \beta_{time}X_{Train_{TT}} + \beta_{cost}X_{Train_{cost}}\\V_2 & = \beta_{SM} +  \beta_{time}X_{SM_{TT}} + \beta_{cost}X_{SM_{cost}}\\V_3 & = \beta_{Car} + \beta_{time}X_{Car_{TT}} + \beta_{cost}X_{Car_{cost}}\\\end{align} **Create parameters to be estimated**`Beta`1. name of parameter2. default value for the parameter3. lower bound4. upper bound5. flag indicating if parameter is to be estimated<jupyter_code>ASC_CAR   = exp.Beta('ASC_CAR',0,None ,None ,0)
ASC_TRAIN = exp.Beta('ASC_TRAIN',0,None ,None ,0)
ASC_SM    = exp.Beta('ASC_SM',0,None ,None ,1)
B_TIME    = exp.Beta('B_TIME',0,None ,None ,0)
B_COST    = exp.Beta('B_COST',0,None ,None ,0)<jupyter_output><empty_output><jupyter_text>**Define the utility functions**<jupyter_code>V1 = ASC_TRAIN + \
     B_TIME * TRAIN_TT_SCALED + \
     B_COST * TRAIN_COST_SCALED
V2 = ASC_SM + \
     B_TIME * SM_TT_SCALED + \
     B_COST * SM_COST_SCALED
V3 = ASC_CAR + \
     B_TIME * CAR_TT_SCALED + \
     B_COST * CAR_CO_SCALED<jupyter_output><empty_output><jupyter_text>**Associate utility functions with alternatives and associate availability of alternatives**Create a python dictionary with all utility functionsCreate a python dictionary with availability of choices<jupyter_code>V = {1: V1,
     2: V2,
     3: V3}

av = {1: TRAIN_AV_SP,
      2: SM_AV,
      3: CAR_AV_SP}<jupyter_output><empty_output><jupyter_text>**Define the model**<jupyter_code>logprob = models.loglogit(V, av, CHOICE)<jupyter_output><empty_output><jupyter_text>**Define the Biogeme object*** Give the database with all variables* Give the log likelihood model<jupyter_code>biogeme  = bio.BIOGEME(database, logprob)

biogeme.modelName = "swissmetro_logit_basic"<jupyter_output><empty_output><jupyter_text>**Estimate the model**1. A `.html` can be generated with a report of the results and can be opened with a browser2. A `.pickle` file can also be generaetd with a snapshot with the results. This file can then be used in other scripts<jupyter_code>biogeme.generateHtml = True
biogeme.generatePickle = False

results = biogeme.estimate()

print(f"HTML file:    {results.data.htmlFileName}")
print(f"Pickle file:  {results.data.pickleFileName }")<jupyter_output>HTML file:    swissmetro_logit_basic.html
Pickle file:  None<jupyter_text>**Print results**<jupyter_code>betas = results.getBetaValues()
for k,v in betas.items():
    print(f"{k:10}=\t{v:.3g}")<jupyter_output>ASC_CAR   =	-0.155
ASC_TRAIN =	-0.701
B_COST    =	-1.08
B_TIME    =	-1.28<jupyter_text>**Get the variance-covariance matrix**<jupyter_code>results.getVarCovar()<jupyter_output><empty_output><jupyter_text>**Model results** * Number of estimated parameters($K$)* Sample size ($N$)* Number of excluded observations* Log likelihood of the sample with the default values for the parameters ($\mathcal{L}^i)$)* Log likelihood for the final estimated model ($\mathcal{L}^*)$)* Likelihood ratio:\begin{align}-2 (\mathcal{L}^i-\mathcal{L}^*)\end{align}* Rho-square for the init model\begin{align}\rho^2 = 1- \frac{\mathcal{L}^*}{\mathcal{L}^i}\end{align}* Rho-square adjusted for the init model\begin{align}\rho^2 = 1- \frac{\mathcal{L}^* - K}{\mathcal{L}^i}\end{align}* Akaike Information Criterion\begin{align}2 K - 2 \mathcal{L}^*\end{align}* Bayesian Information Criterion\begin{align}K N- 2 \mathcal{L}^*\end{align}* Final gradient norm* Iterations* Optimization Time<jupyter_code>gs = results.getGeneralStatistics()

for k,v in gs.items():
    print("{}= {}".format(k.ljust(45),v[0]))<jupyter_output>Number of estimated parameters               = 4
Sample size                                  = 6768
Excluded observations                        = 3960
Init log likelihood                          = -6964.662979192295
Final log likelihood                         = -5331.252006916162
Likelihood ratio test for the init. model    = 3266.8219445522664
Rho-square for the init. model               = 0.23452835796306726
Rho-square-bar for the init. model           = 0.23395403010083615
Akaike Information Criterion                 = 10670.504013832324
Bayesian Information Criterion               = 10697.783857437467
Final gradient norm                          = 0.0006288284529603027
Nbr of threads                               = 8<jupyter_text>**Likelihood ratio test** \begin{align}-2 (\mathcal{L}^i-\mathcal{L}^*)\end{align}\begin{align}Asymptotically \chi^2 Distributed\end{align} H0 -> Both models are equal<jupyter_code>tools.likelihood_ratio_test((gs['Final log likelihood'][0],gs['Number of estimated parameters'][0]),
                            (gs['Init log likelihood'][0], 0), 
                            significance_level=0.95)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.128761665093845, 'fraction_numerical': 0.03302925448254168, 'mean_word_length': 5.149580915538363, 'pattern_counts': {'":': 0, '<': 67, '<?xml version=': 0, '>': 69, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '506680', 'license': 'MIT', 'n_tokens_mistral': 3436, 'n_tokens_neox': 3151, 'n_words': 750, 'path': 'Code/8.1-MultinomialLogitAndProbitModels/.ipynb_checkpoints/01-biogeme-basics-checkpoint.ipynb', 'repo_name': 'U-Shift/Transport-Demand-Modelling', 'size': '9543'}