{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "366516f7",
   "metadata": {},
   "source": [
    "\n",
    "File: 05-lecture_solution.py\n",
    "\n",
    "\n",
    "Michel Bierlaire\n",
    "\n",
    "Tue Aug 05 2025, 16:08:02\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea0c16b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import biogeme.biogeme_logging as blog\n",
    "from IPython.core.display_functions import display\n",
    "from biogeme.biogeme import BIOGEME\n",
    "from biogeme.expressions import Beta, Draws, MonteCarlo, log\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",
    "from scipy.stats import norm\n",
    "\n",
    "from swissmetro_variables import (\n",
    "    CAR_AV_SP,\n",
    "    CAR_CO_SCALED,\n",
    "    CAR_TT_SCALED,\n",
    "    CHOICE,\n",
    "    SM_AV,\n",
    "    SM_COST_SCALED,\n",
    "    SM_HE_SCALED,\n",
    "    SM_TT_SCALED,\n",
    "    TRAIN_AV_SP,\n",
    "    TRAIN_COST_SCALED,\n",
    "    TRAIN_HE_SCALED,\n",
    "    TRAIN_TT_SCALED,\n",
    "    database,\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64c4865f",
   "metadata": {},
   "source": [
    "Variables used for the specification of the Swissmetro model are defined in the file `swissmetro_variables.py`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74e824e3",
   "metadata": {},
   "source": [
    "As the estimation time may be long, we ask Biogeme to report the details of the iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6f1e0f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "logger = blog.get_screen_logger(level=blog.INFO)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38bf1bab",
   "metadata": {},
   "source": [
    "**Tip:**<div class=\"alert alert-block alert-info\">It is advised to start working with a low number of draws,\n",
    "until the script is working well. Then, increase the number of draws to 10000, say. Then, execute the script\n",
    "overnight.  </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6cc999e",
   "metadata": {},
   "outputs": [],
   "source": [
    "number_of_draws = 10_000\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d47ab68a",
   "metadata": {},
   "source": [
    "# Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1b66944",
   "metadata": {},
   "outputs": [],
   "source": [
    "asc_car = Beta('asc_car', 0, None, None, 0)\n",
    "asc_train = Beta('asc_train', 0, None, None, 0)\n",
    "b_time = Beta('b_time', 0, None, None, 0)\n",
    "b_cost = Beta('b_cost', 0, None, None, 0)\n",
    "b_fr = Beta('b_fr', 0, None, None, 0)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da3e8d31",
   "metadata": {},
   "source": [
    "# Availability conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b701d94d",
   "metadata": {},
   "outputs": [],
   "source": [
    "av = {1: TRAIN_AV_SP, 2: SM_AV, 3: CAR_AV_SP}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d847b180",
   "metadata": {},
   "source": [
    "# Logit model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72b57991",
   "metadata": {},
   "source": [
    "## Utility functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85a32af8",
   "metadata": {},
   "outputs": [],
   "source": [
    "v_train = (\n",
    "    asc_train\n",
    "    + b_time * TRAIN_TT_SCALED\n",
    "    + b_cost * TRAIN_COST_SCALED\n",
    "    + b_fr * TRAIN_HE_SCALED\n",
    ")\n",
    "v_swissmetro = b_time * SM_TT_SCALED + b_cost * SM_COST_SCALED + b_fr * SM_HE_SCALED\n",
    "v_car = asc_car + b_time * CAR_TT_SCALED + b_cost * CAR_CO_SCALED\n",
    "v = {1: v_train, 2: v_swissmetro, 3: v_car}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7da9e3d",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5f10a70",
   "metadata": {},
   "outputs": [],
   "source": [
    "logprob = loglogit(v, av, CHOICE)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f00a63dd",
   "metadata": {},
   "source": [
    "## Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ae7ad44",
   "metadata": {},
   "outputs": [],
   "source": [
    "biogeme = BIOGEME(database, logprob)\n",
    "biogeme.model_name = '01logit'\n",
    "results_logit: EstimationResults = biogeme.estimate(recycle=False)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c385b147",
   "metadata": {},
   "source": [
    "## Results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a83a3a2",
   "metadata": {},
   "source": [
    "General statistics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea2ede2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results_logit.short_summary())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70c9a007",
   "metadata": {},
   "source": [
    "Estimated parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74b2ef9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "param_logit = get_pandas_estimated_parameters(estimation_results=results_logit)\n",
    "display(param_logit)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98e10b89",
   "metadata": {},
   "source": [
    "# Random parameter: normal distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31034fcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "b_time_s = Beta('b_time_s', 1, None, None, 0)\n",
    "b_time_rnd = b_time + b_time_s * Draws('b_time_rnd', 'NORMAL')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b16b6d14",
   "metadata": {},
   "source": [
    "## Utility functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9aaa75d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "v_train = (\n",
    "    asc_train\n",
    "    + b_time_rnd * TRAIN_TT_SCALED\n",
    "    + b_cost * TRAIN_COST_SCALED\n",
    "    + b_fr * TRAIN_HE_SCALED\n",
    ")\n",
    "v_swissmetro = b_time_rnd * SM_TT_SCALED + b_cost * SM_COST_SCALED + b_fr * SM_HE_SCALED\n",
    "v_car = asc_car + b_time_rnd * CAR_TT_SCALED + b_cost * CAR_CO_SCALED\n",
    "v = {1: v_train, 2: v_swissmetro, 3: v_car}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "447a6783",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a3d849d",
   "metadata": {},
   "outputs": [],
   "source": [
    "prob = logit(v, av, CHOICE)\n",
    "logprob = log(MonteCarlo(prob))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d6c54c1",
   "metadata": {},
   "source": [
    "## Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9ecd6d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "biogeme = BIOGEME(database, logprob, number_of_draws=number_of_draws)\n",
    "biogeme.model_name = '02normal'\n",
    "results_normal: EstimationResults = biogeme.estimate()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1157b8b",
   "metadata": {},
   "source": [
    "## Results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8a0a193",
   "metadata": {},
   "source": [
    "General statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4b6c0d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results_normal.short_summary())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f863ca32",
   "metadata": {},
   "source": [
    "Estimated parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c9731d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "param_normal = get_pandas_estimated_parameters(estimation_results=results_normal)\n",
    "display(param_normal)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1ec9b35",
   "metadata": {},
   "source": [
    "# Comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6b93178",
   "metadata": {},
   "source": [
    "We build a summary data frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ad994cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "summary = compile_estimation_results(\n",
    "    {'Logit': results_logit, 'Random param. (normal)': results_normal}\n",
    ")\n",
    "display(summary[0])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b63f148",
   "metadata": {},
   "source": [
    "You may note the following:\n",
    "- the large increase of the final log likelihood with the mixture model,\n",
    "- the increase of the (average) value of time when assume the time coefficient to be distributed,\n",
    "- there is a significant portion of the normal distribution that lies in the positive values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1398f1a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_beta_time = results_normal.get_parameter_value('b_time')\n",
    "stddev_beta_time = results_normal.get_parameter_value('b_time_s')\n",
    "probability_to_be_larger_than_zero = 1 - norm.cdf(\n",
    "    0, loc=mean_beta_time, scale=stddev_beta_time\n",
    ")\n",
    "print(\n",
    "    f'Probability for beta_time to be larger than 0: {100 * probability_to_be_larger_than_zero:.3g}%'\n",
    ")"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}