{
  "cells": [
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#%%\n",
        "\"\"\"File 05normalMixtureIntegral.py\n",
        "\n",
        ":author: Michel Bierlaire, EPFL\n",
        ":date: Sat Sep  7 18:23:01 2019\n",
        "\n",
        " Example of a mixture of logit models, using numerical integration.\n",
        " Three alternatives: Train, Car and Swissmetro\n",
        " SP data\n",
        "\"\"\"\n",
        "\n",
        "import pandas as pd\n",
        "import biogeme.database as db\n",
        "import biogeme.biogeme as bio\n",
        "import biogeme.distributions as dist\n",
        "import biogeme.models as models\n",
        "import biogeme.messaging as msg\n",
        "from biogeme.expressions import Beta, DefineVariable, RandomVariable, log, Integrate\n",
        "\n",
        "# Read the data\n",
        "df = pd.read_csv('swissmetro.dat', '\\t')\n",
        "database = db.Database('swissmetro', df)\n",
        "\n",
        "# The Pandas data structure is available as database.data. Use all the\n",
        "# Pandas functions to investigate the database. For example:\n",
        "#print(database.data.describe())\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",
        "# Removing some observations can be done directly using pandas.\n",
        "#remove = (((database.data.PURPOSE != 1) &\n",
        "#           (database.data.PURPOSE != 3)) |\n",
        "#          (database.data.CHOICE == 0))\n",
        "#database.data.drop(database.data[remove].index,inplace=True)\n",
        "\n",
        "# Here we use the \"biogeme\" way for backward compatibility\n",
        "exclude = ((PURPOSE != 1) * (PURPOSE != 3) + (CHOICE == 0)) > 0\n",
        "database.remove(exclude)\n",
        "\n",
        "# Parameters to be estimated\n",
        "ASC_CAR = Beta('ASC_CAR', 0, None, None, 0)\n",
        "ASC_TRAIN = Beta('ASC_TRAIN', 0, None, None, 0)\n",
        "ASC_SM = Beta('ASC_SM', 0, None, None, 1)\n",
        "B_COST = Beta('B_COST', 0, None, None, 0)\n",
        "\n",
        "# Define a random parameter, normally distributed, designed to be used\n",
        "# for numerical integration\n",
        "B_TIME = Beta('B_TIME', 0, None, None, 0)\n",
        "# It is advised not to use 0 as starting value for the following parameter.\n",
        "B_TIME_S = Beta('B_TIME_S', 1, None, None, 0)\n",
        "omega = RandomVariable('omega')\n",
        "density = dist.normalpdf(omega)\n",
        "B_TIME_RND = B_TIME + B_TIME_S * omega\n",
        "\n",
        "# Definition of new variables\n",
        "SM_COST = SM_CO * (GA == 0)\n",
        "TRAIN_COST = TRAIN_CO * (GA == 0)\n",
        "\n",
        "# Definition of new variables: adding columns to the database\n",
        "CAR_AV_SP = DefineVariable('CAR_AV_SP', CAR_AV * (SP != 0), database)\n",
        "TRAIN_AV_SP = DefineVariable('TRAIN_AV_SP', TRAIN_AV * (SP != 0), database)\n",
        "TRAIN_TT_SCALED = DefineVariable('TRAIN_TT_SCALED', TRAIN_TT / 100.0, database)\n",
        "TRAIN_COST_SCALED = DefineVariable('TRAIN_COST_SCALED', TRAIN_COST / 100, database)\n",
        "SM_TT_SCALED = DefineVariable('SM_TT_SCALED', SM_TT / 100.0, database)\n",
        "SM_COST_SCALED = DefineVariable('SM_COST_SCALED', SM_COST / 100, database)\n",
        "CAR_TT_SCALED = DefineVariable('CAR_TT_SCALED', CAR_TT / 100, database)\n",
        "CAR_CO_SCALED = DefineVariable('CAR_CO_SCALED', CAR_CO / 100, database)\n",
        "\n",
        "# Definition of the utility functions\n",
        "V1 = ASC_TRAIN + \\\n",
        "     B_TIME_RND * TRAIN_TT_SCALED + \\\n",
        "     B_COST * TRAIN_COST_SCALED\n",
        "V2 = ASC_SM + \\\n",
        "     B_TIME_RND * SM_TT_SCALED + \\\n",
        "     B_COST * SM_COST_SCALED\n",
        "V3 = ASC_CAR + \\\n",
        "     B_TIME_RND * CAR_TT_SCALED + \\\n",
        "     B_COST * CAR_CO_SCALED\n",
        "\n",
        "# Associate utility functions with the numbering of alternatives\n",
        "V = {1: V1,\n",
        "     2: V2,\n",
        "     3: V3}\n",
        "\n",
        "# Associate the availability conditions with the alternatives\n",
        "av = {1: TRAIN_AV_SP,\n",
        "      2: SM_AV,\n",
        "      3: CAR_AV_SP}\n",
        "\n",
        "# Conditional to omega, we have a logit model (called the kernel)\n",
        "condprob = models.logit(V, av, CHOICE)\n",
        "# We integrate over omega using numerical integration\n",
        "logprob = 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, logprob)\n",
        "biogeme.modelName = '05normalMixtureIntegral'\n",
        "\n",
        "# Estimate the parameters\n",
        "results = biogeme.estimate()\n",
        "pandasResults = results.getEstimatedParameters()\n",
        "print(pandasResults)\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": 1
}