{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "74a2907c",
   "metadata": {},
   "source": [
    "\n",
    "File: 02-aggregation_solution.py\n",
    "\n",
    "\n",
    "Michel Bierlaire\n",
    "\n",
    "Sat Aug 02 2025, 18:41:06\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8abd211d",
   "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 Variable, exp\n",
    "from biogeme.results_processing import EstimationResults\n",
    "\n",
    "from netherlands_model import (\n",
    "    beta_cost,\n",
    "    beta_time_rail,\n",
    "    logprob,\n",
    "    prob_car,\n",
    "    prob_rail,\n",
    "    rail_cost_euro,\n",
    "    rail_time,\n",
    "    v_car,\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1bcbea8",
   "metadata": {},
   "source": [
    "The objective of this laboratory is to calculate aggregate indicators from a choice model."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba0bcbb6",
   "metadata": {},
   "source": [
    "# Calculate the weights of the sample in order to perform sample enumeration."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00c379ee",
   "metadata": {},
   "source": [
    "We consider four segments of the population, characterized by two age categories  (41 or older, and 40 or younger),\n",
    "and gender."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56044e63",
   "metadata": {},
   "source": [
    "The census data report the following number $N_g$ of individuals in each segment $g$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b22d439",
   "metadata": {},
   "outputs": [],
   "source": [
    "census = {\n",
    "    'male_41_more': 4092390,\n",
    "    'male_40_less': 4151092,\n",
    "    'female_41_more': 3984028,\n",
    "    'female_40_less': 4428289,\n",
    "}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3921e300",
   "metadata": {},
   "source": [
    "The total size $N$ of the population is the sum over all segments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34e0b750",
   "metadata": {},
   "outputs": [],
   "source": [
    "population_size = sum(census.values())\n",
    "print(f'Size of the population: {population_size:_}')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48d7c494",
   "metadata": {},
   "source": [
    "We read the data file, and consider only the revealed preferences data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe054a19",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('netherlands.dat', sep='\\t')\n",
    "df = df[df['sp'] == 0]\n",
    "display(df)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef254c1f",
   "metadata": {},
   "source": [
    "Identify each segment in the database."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aedd89e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "filters = {\n",
    "    'male_41_more': (df.age == 0) & (df.gender == 0),\n",
    "    'male_40_less': (df.age == 1) & (df.gender == 0),\n",
    "    'female_41_more': (df.age == 0) & (df.gender == 1),\n",
    "    'female_40_less': (df.age == 1) & (df.gender == 1),\n",
    "}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d307f323",
   "metadata": {},
   "source": [
    "We count the number $S_g$ of individuals in each segment $g$ in the sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c2a4f62",
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_segments = {\n",
    "    segment_name: segment_rows.sum() for segment_name, segment_rows in filters.items()\n",
    "}\n",
    "print(sample_segments)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f2f8b2b",
   "metadata": {},
   "source": [
    "We check the total sample size $S$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8056d10",
   "metadata": {},
   "outputs": [],
   "source": [
    "total_sample = sum(sample_segments.values())\n",
    "print(f'Sample size: {total_sample}')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62d1d9f3",
   "metadata": {},
   "source": [
    "The weight $w_g$ associated with segment $g$ is defined as\n",
    "$$\n",
    "w_g = \\frac{N_g}{N}\\frac{S}{S_g}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4f813c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "weights = {\n",
    "    segment_name: census[segment_name] * total_sample / (segment_size * population_size)\n",
    "    for segment_name, segment_size in sample_segments.items()\n",
    "}\n",
    "print(weights)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e2e12f2",
   "metadata": {},
   "source": [
    "We insert the weight as a new column in the database."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91ff427d",
   "metadata": {},
   "outputs": [],
   "source": [
    "for segment_name, segment_rows in filters.items():\n",
    "    df.loc[segment_rows, 'weight'] = weights[segment_name]\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd43fda5",
   "metadata": {},
   "source": [
    "We check that the sum of the weights is the sample size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c012f4bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "sum_weights = df['weight'].sum()\n",
    "print(f'Sum of the weights: {sum_weights}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1db69249",
   "metadata": {},
   "source": [
    "Here is the database. The new column with weights is the last."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c180c0b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "display(df)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20e39584",
   "metadata": {},
   "source": [
    "# Perform sample enumeration"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b42337b",
   "metadata": {},
   "source": [
    "We obtain the model from the specification file, and estimate its parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3d13872",
   "metadata": {},
   "outputs": [],
   "source": [
    "database = Database('netherlands', df)\n",
    "biogeme = BIOGEME(database, logprob)\n",
    "biogeme.model_name = 'netherlands_model'\n",
    "results: EstimationResults = biogeme.estimate()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "740b2fea",
   "metadata": {},
   "source": [
    "We calculate the choice probability of each alternative for each observation. We include the weight in the simulation\n",
    "results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce638192",
   "metadata": {},
   "outputs": [],
   "source": [
    "weight = Variable('weight')\n",
    "simulate = {\n",
    "    'weight': weight,\n",
    "    'Prob. rail': prob_rail,\n",
    "    'Prob. car': prob_car,\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e1c822d",
   "metadata": {},
   "source": [
    "We perform the simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c4a98cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "biosim = BIOGEME(database, simulate, bootstrap_samples=100)\n",
    "simulated_values = biosim.simulate(results.get_beta_values())\n",
    "display(simulated_values)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15f5b7da",
   "metadata": {},
   "source": [
    "Market shares are calculated using the weighted mean of the\n",
    "individual probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58c463e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulated_values['Weighted rail'] = (\n",
    "    simulated_values['weight'] * simulated_values['Prob. rail']\n",
    ")\n",
    "simulated_values['Weighted car'] = (\n",
    "    simulated_values['weight'] * simulated_values['Prob. car']\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cc8d080",
   "metadata": {},
   "source": [
    "For rail:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e99757e",
   "metadata": {},
   "outputs": [],
   "source": [
    "market_share_rail = simulated_values['Weighted rail'].mean()\n",
    "print(f'Market share for rail: {100*market_share_rail:.1f}%')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "896712c1",
   "metadata": {},
   "source": [
    "For car:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdfa138f",
   "metadata": {},
   "outputs": [],
   "source": [
    "market_share_car = simulated_values['Weighted car'].mean()\n",
    "print(f'Market share for car: {100*market_share_car:.1f}%')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa6222c6",
   "metadata": {},
   "source": [
    "# Confidence intervals"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c0f9254",
   "metadata": {},
   "source": [
    "In order to calculate confidence intervals, we need to re-estimate the parameters using bootstrapping."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07ece9d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "results_bootstrapping = biogeme.estimate(run_bootstrap=True)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "730d1019",
   "metadata": {},
   "source": [
    "We obtain a sample of values for the parameters. We use them to calculate empirically (that is, using simulation),\n",
    "the 90% confidence intervals on the simulated quantities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "567d4c1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "b = results_bootstrapping.get_betas_for_sensitivity_analysis()\n",
    "left, right = biosim.confidence_intervals(b, 0.9)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ee2e9b2",
   "metadata": {},
   "source": [
    "The `left`and `right` data frames have the same structure as the `simulated_values`data frame, and\n",
    "contain the lower bound (left) and the upper bound (right) of the confidence interval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98194e02",
   "metadata": {},
   "outputs": [],
   "source": [
    "display(left)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bd2737f",
   "metadata": {},
   "outputs": [],
   "source": [
    "display(right)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c56762e6",
   "metadata": {},
   "source": [
    "We calculate the weighted sum for the lower bound of the interval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2af61029",
   "metadata": {},
   "outputs": [],
   "source": [
    "left['Weighted rail'] = left['weight'] * left['Prob. rail']\n",
    "left['Weighted car'] = left['weight'] * left['Prob. car']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a327e0c8",
   "metadata": {},
   "source": [
    "And for the upper bound:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "903810b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "right['Weighted rail'] = right['weight'] * right['Prob. rail']\n",
    "right['Weighted car'] = right['weight'] * right['Prob. car']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70ab5287",
   "metadata": {},
   "source": [
    "We calculate the lower bounds on the market shares."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8525470",
   "metadata": {},
   "outputs": [],
   "source": [
    "left_market_share_rail = left['Weighted rail'].mean()\n",
    "left_market_share_car = left['Weighted car'].mean()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1526690a",
   "metadata": {},
   "source": [
    "And the upper bounds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df4d2bae",
   "metadata": {},
   "outputs": [],
   "source": [
    "right_market_share_rail = right['Weighted rail'].mean()\n",
    "right_market_share_car = right['Weighted car'].mean()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89436a6f",
   "metadata": {},
   "source": [
    "Market share for car"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5a7d177",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\n",
    "    f'Market share for car: {100*market_share_car:.1f}% '\n",
    "    f'CI: ['\n",
    "    f'{100*left_market_share_car:.1f}%-'\n",
    "    f'{100*right_market_share_car:.1f}'\n",
    "    f']'\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dead8ade",
   "metadata": {},
   "source": [
    "Market share for rail."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c87b2373",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\n",
    "    f'Market share for rail: {100*market_share_rail:.1f}% '\n",
    "    f'CI: ['\n",
    "    f'{100*left_market_share_rail:.1f}%-'\n",
    "    f'{100*right_market_share_rail:.1f}'\n",
    "    f']'\n",
    ")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d72ea309",
   "metadata": {},
   "source": [
    "# Scenario: rail cost decreased by 10%"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9008052",
   "metadata": {},
   "source": [
    "We write the choice probability under the proposed scenario."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53700ba7",
   "metadata": {},
   "outputs": [],
   "source": [
    "new_rail_cost_euro = rail_cost_euro * 0.9\n",
    "v_rail_scenario = beta_cost * new_rail_cost_euro + beta_time_rail * rail_time\n",
    "prob_car_scenario = 1 / (1 + exp(v_rail_scenario - v_car))\n",
    "prob_rail_scenario = 1 - prob_car_scenario\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5da02b5c",
   "metadata": {},
   "source": [
    "We perform the simulation, exactly as before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0eac0065",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulate_scenario = {\n",
    "    'weight': weight,\n",
    "    'Prob. rail': prob_rail_scenario,\n",
    "    'Prob. car': prob_car_scenario,\n",
    "}\n",
    "biosim_scenario = BIOGEME(database, simulate_scenario)\n",
    "simulated_scenario = biosim_scenario.simulate(results.get_beta_values())\n",
    "left_scenario, right_scenario = biosim_scenario.confidence_intervals(b, 0.9)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e37e7d3c",
   "metadata": {},
   "source": [
    "We calculate the weighted probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f468131a",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulated_scenario['Weighted rail'] = (\n",
    "    simulated_scenario['weight'] * simulated_scenario['Prob. rail']\n",
    ")\n",
    "simulated_scenario['Weighted car'] = (\n",
    "    simulated_scenario['weight'] * simulated_scenario['Prob. car']\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fff58e0",
   "metadata": {},
   "source": [
    "And the corresponding confidence intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "547102ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "left_scenario['Weighted rail'] = left_scenario['weight'] * left_scenario['Prob. rail']\n",
    "left_scenario['Weighted car'] = left_scenario['weight'] * left_scenario['Prob. car']\n",
    "right_scenario['Weighted rail'] = (\n",
    "    right_scenario['weight'] * right_scenario['Prob. rail']\n",
    ")\n",
    "right_scenario['Weighted car'] = right_scenario['weight'] * right_scenario['Prob. car']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c28591e6",
   "metadata": {},
   "source": [
    "Finally, the market shares and the confidence intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83974107",
   "metadata": {},
   "outputs": [],
   "source": [
    "market_share_car_scenario = simulated_scenario['Weighted car'].mean()\n",
    "market_share_rail_scenario = simulated_scenario['Weighted rail'].mean()\n",
    "left_market_share_rail_scenario = left_scenario['Weighted rail'].mean()\n",
    "left_market_share_car_scenario = left_scenario['Weighted car'].mean()\n",
    "right_market_share_rail_scenario = right_scenario['Weighted rail'].mean()\n",
    "right_market_share_car_scenario = right_scenario['Weighted car'].mean()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6748e907",
   "metadata": {},
   "source": [
    "Market share for car."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "410ed944",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\n",
    "    f'Market share for car (scenario): {100*market_share_car_scenario:.1f}% '\n",
    "    f'CI: ['\n",
    "    f'{100*left_market_share_car_scenario:.1f}%-'\n",
    "    f'{100*right_market_share_car_scenario:.1f}'\n",
    "    f']'\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c30615e0",
   "metadata": {},
   "source": [
    "Market share for rail"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5571435e",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\n",
    "    f'Market share for rail (scenario): {100*market_share_rail_scenario:.1f}% '\n",
    "    f'CI: ['\n",
    "    f'{100*left_market_share_rail_scenario:.1f}%-'\n",
    "    f'{100*right_market_share_rail_scenario:.1f}'\n",
    "    f']'\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47f7d0aa",
   "metadata": {},
   "source": [
    "When we decrease the price of rail by 10%, the market shares change from 62.6%/37.4% to 59.4%/40.6%."
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}