{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3b5d9c0d",
   "metadata": {},
   "source": [
    "\n",
    "File: 02-outlier_solution.py\n",
    "\n",
    "\n",
    "Michel Bierlaire\n",
    "\n",
    "Sat Aug 02 2025, 18:19:59\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a114b74",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import pandas as pd\n",
    "from IPython.core.display_functions import display\n",
    "from biogeme.biogeme import BIOGEME\n",
    "from biogeme.database import Database\n",
    "from biogeme.expressions import Beta, Variable\n",
    "from biogeme.models import logit, loglogit\n",
    "from biogeme.results_processing import (\n",
    "    EstimationResults,\n",
    "    compile_estimation_results,\n",
    "    get_pandas_estimated_parameters,\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f02845fc",
   "metadata": {},
   "source": [
    "The objective of this laboratory is to illustrate the outlier analysis. We use the Optima case study, for\n",
    "transportation mode choice in Switzerland."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d809bea7",
   "metadata": {},
   "source": [
    "Read the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6298be11",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('optima.dat', sep='\\t')\n",
    "display(df)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29445304",
   "metadata": {},
   "source": [
    "Prepare the data for Biogeme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48076e15",
   "metadata": {},
   "outputs": [],
   "source": [
    "database = Database('optima', df)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ebf96b9",
   "metadata": {},
   "source": [
    "Identification of the relevant variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d62bfc4b",
   "metadata": {},
   "outputs": [],
   "source": [
    "Choice = Variable('Choice')\n",
    "Weight = Variable('Weight')\n",
    "age = Variable('age')\n",
    "Gender = Variable('Gender')\n",
    "TimePT = Variable('TimePT')\n",
    "TimeCar = Variable('TimeCar')\n",
    "TripPurpose = Variable('TripPurpose')\n",
    "MarginalCostPT = Variable('MarginalCostPT')\n",
    "CostCarCHF = Variable('CostCarCHF')\n",
    "distance_km = Variable('distance_km')\n",
    "Education = Variable('Education')\n",
    "LangCode = Variable('LangCode')\n",
    "OccupStat = Variable('OccupStat')\n",
    "NbTransf = Variable('NbTransf')\n",
    "Income = Variable('Income')\n",
    "WaitingTimePT = Variable('WaitingTimePT')\n",
    "CarAvail = Variable('CarAvail')\n",
    "Subscription = Variable('Subscription')\n",
    "GenAbST = Variable('GenAbST')\n",
    "OwnHouse = Variable('OwnHouse')\n",
    "NbBicy = Variable('NbBicy')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "471f27a7",
   "metadata": {},
   "source": [
    "Removing some incorrectly coded data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ad231e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "exclude = (\n",
    "    (age == -1)\n",
    "    + (Choice == -1)\n",
    "    + (CostCarCHF < 0)\n",
    "    + (Income == -1)\n",
    "    + (CarAvail == 3) * (Choice == 1)\n",
    ") > 0\n",
    "database.remove(exclude)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d673f7e7",
   "metadata": {},
   "source": [
    "We first estimate a simple model:\n",
    "\\begin{align*}\n",
    "V_\\text{PT} &= \\text{Cte}_\\text{PT} + \\beta_{t, \\text{PT}} \\text{Time}_\\text{PT}\n",
    "+ \\beta_{c, \\text{PT}} \\text{Cost}_\\text{PT}, \\\\\n",
    "V_\\text{Car} &= \\text{Cte}_\\text{Car} + \\beta_{t, \\text{Car}} \\text{Time}_\\text{Car}\n",
    "+ \\beta_{c, \\text{Car}} \\text{Cost}_\\text{Car}, \\\\\n",
    "V_\\text{SM} &= \\beta_d \\text{distance}.\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe7d9096",
   "metadata": {},
   "source": [
    "Parameters to be estimated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de9fcf6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "asc_pt = Beta('asc_pt', 0, None, None, 0)\n",
    "asc_car = Beta('asc_car', 0, None, None, 0)\n",
    "beta_time_pt = Beta('beta_time_pt', 0.0, None, None, 0)\n",
    "beta_time_car = Beta('beta_time_car', 0.0, None, None, 0)\n",
    "beta_cost = Beta('beta_cost', 0, None, None, 0)\n",
    "beta_distance = Beta('beta_distance', 0, None, None, 0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65b49ec3",
   "metadata": {},
   "source": [
    "Utility functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e65e9451",
   "metadata": {},
   "outputs": [],
   "source": [
    "v_pt = asc_pt + beta_time_pt * TimePT + beta_cost * MarginalCostPT\n",
    "v_car = asc_car + beta_time_car * TimeCar + beta_cost * CostCarCHF\n",
    "v_sm = beta_distance * distance_km\n",
    "\n",
    "v = {0: v_pt, 1: v_car, 2: v_sm}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50482b3b",
   "metadata": {},
   "source": [
    "Availability conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d687e0f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "av = {0: 1, 1: CarAvail != 3, 2: 1}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b92aa256",
   "metadata": {},
   "source": [
    "Estimation of the parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2258f7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "logprob = loglogit(v, av, Choice)\n",
    "biogeme = BIOGEME(database, logprob)\n",
    "biogeme.model_name = 'logit_optima_base'\n",
    "results: EstimationResults = biogeme.estimate()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "945d041e",
   "metadata": {},
   "source": [
    "General statistics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a85eaef3",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results.short_summary())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5da6371",
   "metadata": {},
   "source": [
    "Estimated parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48bca2fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "display(get_pandas_estimated_parameters(estimation_results=results))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68f43d95",
   "metadata": {},
   "source": [
    "We now simulate the estimated model to obtain the contribution of each observation to the likelihood function. It is\n",
    "the probability predicted by the model to choose the actually chosen alternative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6649d7f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "prob_chosen = logit(v, av, Choice)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d1f26da",
   "metadata": {},
   "source": [
    "We define a dictionary with the formulas to be simulated. Here, there is only one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f61e841d",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulate = {'Prob. chosen': prob_chosen}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4792cd8a",
   "metadata": {},
   "source": [
    "We perform the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89ab8113",
   "metadata": {},
   "outputs": [],
   "source": [
    "biosim = BIOGEME(database, simulate)\n",
    "betas_values = results.get_beta_values()\n",
    "sim_results = biosim.simulate(the_beta_values=betas_values)\n",
    "display(sim_results)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be64570b",
   "metadata": {},
   "source": [
    "We now merge the result of the simulation with the original database.\n",
    "We also sort the rows based on the 'Prob. chosen'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0f7c10a",
   "metadata": {},
   "outputs": [],
   "source": [
    "merged_table = pd.concat([sim_results, database.dataframe], axis='columns').sort_values(\n",
    "    by=['Prob. chosen']\n",
    ")\n",
    "display(merged_table)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4589cdfa",
   "metadata": {},
   "source": [
    "As we are interested in the outliers, we drop from the table all observations such that the probability predicted by\n",
    "the model is above 10%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91d3da26",
   "metadata": {},
   "outputs": [],
   "source": [
    "merged_table.drop(merged_table[merged_table['Prob. chosen'] >= 0.1].index, inplace=True)\n",
    "display(merged_table)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d53fd49",
   "metadata": {},
   "source": [
    "Let's investigate only the variables in the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97ec3785",
   "metadata": {},
   "outputs": [],
   "source": [
    "variables = [\n",
    "    'Choice',\n",
    "    'TimePT',\n",
    "    'MarginalCostPT',\n",
    "    'TimeCar',\n",
    "    'CostCarCHF',\n",
    "    'distance_km',\n",
    "]\n",
    "display(merged_table[variables])\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8c66172",
   "metadata": {},
   "source": [
    "We realize that alternative 2 (slow modes) has been selected for trips with a long distance. In order to\n",
    "investigate why, we extract all observations such that the alternative \"slow modes\" was chosen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "600438f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "outliers_sm = merged_table[merged_table['Choice'] == 2]\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26baac23",
   "metadata": {},
   "source": [
    "We now compare the distribution of the variable `NbBicy` in this table, with the distribution of the same variable\n",
    "in the full sample."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95645b99",
   "metadata": {},
   "source": [
    "Frequency of NbBicy in the outliers' database."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca8f3f20",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(outliers_sm['NbBicy'].value_counts(normalize=True).sort_index().cumsum())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00dfa5ef",
   "metadata": {},
   "source": [
    "Frequency of NbBicy in the complete database."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01c09f88",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(database.dataframe['NbBicy'].value_counts(normalize=True).sort_index().cumsum())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b3e257c",
   "metadata": {},
   "source": [
    "We realize that, among the outliers that have chosen \"slow modes\", the proportion of individuals with 5 bicycles or\n",
    "more is very high, compared with the proportion in the full sample. It seems that the ownership of bicycles is a\n",
    "sign of an attraction to use slow modes, even for long distance."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecc68b17",
   "metadata": {},
   "source": [
    "We investigate another variable: `age`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8568586d",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(outliers_sm['age'].value_counts(normalize=True).sort_index().cumsum())\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad8464e0",
   "metadata": {},
   "source": [
    "We realize that most outliers (~70%) are 45 or more. We can assume that individuals in that age range like to practice\n",
    "slow modes to keep in shape."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c404aed",
   "metadata": {},
   "source": [
    "In order to check these hypotheses, let's introduce the variables `NbBicy` and `age` in the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50c31513",
   "metadata": {},
   "outputs": [],
   "source": [
    "beta_bicy = Beta('beta_bicy', 0, None, None, 0)\n",
    "beta_45_or_more = Beta('beta_45_or_more', 0, None, None, 0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0598b5bb",
   "metadata": {},
   "source": [
    "New specification of the utility function for slow modes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53fe92b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "v_sm_revised = (\n",
    "    beta_distance * distance_km + beta_bicy * NbBicy + beta_45_or_more * (age >= 45)\n",
    ")\n",
    "v_revised = {0: v_pt, 1: v_car, 2: v_sm_revised}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbedc8ae",
   "metadata": {},
   "source": [
    "Estimation of the revised model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1751e1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "logprob_revised = loglogit(v_revised, av, Choice)\n",
    "biogeme_revised = BIOGEME(database, logprob_revised)\n",
    "biogeme_revised.model_name = 'logit_optima_revised'\n",
    "results_revised: EstimationResults = biogeme_revised.estimate()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41443580",
   "metadata": {},
   "source": [
    "General statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4d20ec6",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results_revised.short_summary())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "127f3f67",
   "metadata": {},
   "source": [
    "Estimated parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "928e603f",
   "metadata": {},
   "outputs": [],
   "source": [
    "display(get_pandas_estimated_parameters(estimation_results=results_revised))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88bc99b8",
   "metadata": {},
   "source": [
    "We compare the two models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3615198c",
   "metadata": {},
   "outputs": [],
   "source": [
    "comparison, _ = compile_estimation_results(\n",
    "    {'Base': results, 'Bicycles': results_revised}\n",
    ")\n",
    "display(comparison)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bf19132",
   "metadata": {},
   "source": [
    "We observe a significant improvement in the fit."
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}