{ "cells": [ { "cell_type": "code", "metadata": {}, "source": [ "#%%\n", "\"\"\"File 03choiceOnly.py\n", "\n", "Choice model with the latent variable.\n", "Mixture of logit.\n", "No measurement equation for the indicators.\n", "\n", ":author: Michel Bierlaire, EPFL\n", ":date: Thu Sep 6 15:14:39 2018\n", "\n", "\"\"\"\n", "\n", "import pandas as pd\n", "import biogeme.database as db\n", "import biogeme.biogeme as bio\n", "from biogeme import models\n", "import biogeme.distributions as dist\n", "import biogeme.optimization as opt\n", "import biogeme.messaging as msg\n", "from biogeme.expressions import (\n", " Beta,\n", " DefineVariable,\n", " RandomVariable,\n", " Integrate,\n", " exp,\n", " log,\n", ")\n", "\n", "# Read the data\n", "df = pd.read_csv('optima.dat', sep='\\t')\n", "database = db.Database('optima', df)\n", "\n", "\n", "# The following statement allows you to use the names of the variable\n", "# as Python variable.\n", "globals().update(database.variables)\n", "\n", "# Exclude observations such that the chosen alternative is -1\n", "database.remove(Choice == -1.0)\n", "\n", "### Variables\n", "\n", "# Piecewise linear definition of income\n", "ScaledIncome = DefineVariable(\n", " 'ScaledIncome', CalculatedIncome / 1000, database\n", ")\n", "\n", "thresholds = [None, 4, 6, 8, 10, None]\n", "formulaIncome = models.piecewiseFormula(\n", " ScaledIncome, thresholds, [0.0, 0.0, 0.0, 0.0, 0.0]\n", ")\n", "\n", "# Definition of other variables\n", "age_65_more = DefineVariable('age_65_more', age >= 65, database)\n", "moreThanOneCar = DefineVariable('moreThanOneCar', NbCar > 1, database)\n", "moreThanOneBike = DefineVariable('moreThanOneBike', NbBicy > 1, database)\n", "individualHouse = DefineVariable('individualHouse', HouseType == 1, database)\n", "male = DefineVariable('male', Gender == 1, database)\n", "haveChildren = DefineVariable(\n", " 'haveChildren', ((FamilSitu == 3) + (FamilSitu == 4)) > 0, database\n", ")\n", "haveGA = DefineVariable('haveGA', GenAbST == 1, database)\n", "highEducation = DefineVariable('highEducation', Education >= 6, database)\n", "\n", "# Parameters to be estimated\n", "coef_intercept = Beta('coef_intercept', 0.0, None, None, 1)\n", "coef_age_65_more = Beta('coef_age_65_more', 0.0, None, None, 0)\n", "coef_haveGA = Beta('coef_haveGA', 0.0, None, None, 0)\n", "coef_moreThanOneCar = Beta('coef_moreThanOneCar', 0.0, None, None, 0)\n", "coef_moreThanOneBike = Beta('coef_moreThanOneBike', 0.0, None, None, 0)\n", "coef_individualHouse = Beta('coef_individualHouse', 0.0, None, None, 0)\n", "coef_male = Beta('coef_male', 0.0, None, None, 0)\n", "coef_haveChildren = Beta('coef_haveChildren', 0.0, None, None, 0)\n", "coef_highEducation = Beta('coef_highEducation', 0.0, None, None, 0)\n", "\n", "### Latent variable: structural equation\n", "\n", "# Define a random parameter, normally distributed, designed to be used\n", "# for numerical integration\n", "omega = RandomVariable('omega')\n", "density = dist.normalpdf(omega)\n", "sigma_s = Beta('sigma_s', 1, None, None, 0)\n", "\n", "CARLOVERS = (\n", " coef_intercept\n", " + coef_age_65_more * age_65_more\n", " + formulaIncome\n", " + coef_moreThanOneCar * moreThanOneCar\n", " + coef_moreThanOneBike * moreThanOneBike\n", " + coef_individualHouse * individualHouse\n", " + coef_male * male\n", " + coef_haveChildren * haveChildren\n", " + coef_haveGA * haveGA\n", " + coef_highEducation * highEducation\n", " + sigma_s * omega\n", ")\n", "\n", "# Choice model\n", "\n", "ASC_CAR = Beta('ASC_CAR', 0.0, None, None, 0)\n", "ASC_PT = Beta('ASC_PT', 0.0, None, None, 1)\n", "ASC_SM = Beta('ASC_SM', 0.0, None, None, 0)\n", "BETA_COST_HWH = Beta('BETA_COST_HWH', 0.0, None, None, 0)\n", "BETA_COST_OTHER = Beta('BETA_COST_OTHER', 0.0, None, None, 0)\n", "BETA_DIST = Beta('BETA_DIST', 0.0, None, None, 0)\n", "BETA_TIME_CAR_REF = Beta('BETA_TIME_CAR_REF', 0.0, None, 0, 0)\n", "BETA_TIME_CAR_CL = Beta('BETA_TIME_CAR_CL', 0.0, None, None, 0)\n", "BETA_TIME_PT_REF = Beta('BETA_TIME_PT_REF', 0.0, None, 0, 0)\n", "BETA_TIME_PT_CL = Beta('BETA_TIME_PT_CL', -1.0, None, None, 0)\n", "BETA_WAITING_TIME = Beta('BETA_WAITING_TIME', 0.0, None, None, 0)\n", "\n", "TimePT_scaled = DefineVariable('TimePT_scaled', TimePT / 200, database)\n", "TimeCar_scaled = DefineVariable('TimeCar_scaled', TimeCar / 200, database)\n", "MarginalCostPT_scaled = DefineVariable(\n", " 'MarginalCostPT_scaled', MarginalCostPT / 10, database\n", ")\n", "CostCarCHF_scaled = DefineVariable(\n", " 'CostCarCHF_scaled', CostCarCHF / 10, database\n", ")\n", "distance_km_scaled = DefineVariable(\n", " 'distance_km_scaled', distance_km / 5, database\n", ")\n", "PurpHWH = DefineVariable('PurpHWH', TripPurpose == 1, database)\n", "PurpOther = DefineVariable('PurpOther', TripPurpose != 1, database)\n", "\n", "### Definition of utility functions:\n", "\n", "BETA_TIME_PT = BETA_TIME_PT_REF * exp(BETA_TIME_PT_CL * CARLOVERS)\n", "\n", "V0 = (\n", " ASC_PT\n", " + BETA_TIME_PT * TimePT_scaled\n", " + BETA_WAITING_TIME * WaitingTimePT\n", " + BETA_COST_HWH * MarginalCostPT_scaled * PurpHWH\n", " + BETA_COST_OTHER * MarginalCostPT_scaled * PurpOther\n", ")\n", "\n", "BETA_TIME_CAR = BETA_TIME_CAR_REF * exp(BETA_TIME_CAR_CL * CARLOVERS)\n", "\n", "V1 = (\n", " ASC_CAR\n", " + BETA_TIME_CAR * TimeCar_scaled\n", " + BETA_COST_HWH * CostCarCHF_scaled * PurpHWH\n", " + BETA_COST_OTHER * CostCarCHF_scaled * PurpOther\n", ")\n", "\n", "V2 = ASC_SM + BETA_DIST * distance_km_scaled\n", "\n", "# Associate utility functions with the numbering of alternatives\n", "V = {0: V0, 1: V1, 2: V2}\n", "\n", "# Conditional to omega, we have a logit model (called the kernel)\n", "condprob = models.logit(V, None, Choice)\n", "# We integrate over omega using numerical integration\n", "loglike = log(Integrate(condprob * density, 'omega'))\n", "\n", "# Define level of verbosity\n", "logger = msg.bioMessage()\n", "# logger.setSilent()\n", "# logger.setWarning()\n", "logger.setGeneral()\n", "# logger.setDetailed()\n", "\n", "# Create the Biogeme object\n", "biogeme = bio.BIOGEME(database, loglike)\n", "biogeme.modelName = '03choiceOnly'\n", "\n", "# Estimate the parameters\n", "results = biogeme.estimate(algorithm=opt.bioNewton)\n", "\n", "print(f'Estimated betas: {len(results.data.betaValues)}')\n", "print(f'Final log likelihood: {results.data.logLike:.3f}')\n", "print(f'Output file: {results.data.htmlFileName}')\n", "results.writeLaTeX()\n", "print(f'LaTeX file: {results.data.latexFileName}')\n" ], "outputs": [], "execution_count": null } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.1" } }, "nbformat": 4, "nbformat_minor": 4 }