{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e73700fb",
   "metadata": {},
   "source": [
    "\n",
    "File: 06-lecture_solution.py\n",
    "\n",
    "\n",
    "Michel Bierlaire\n",
    "\n",
    "Wed Aug 06 2025, 11:49:08\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "957a824c",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import biogeme.biogeme_logging as blog\n",
    "import numpy as np\n",
    "from IPython.core.display_functions import display\n",
    "from biogeme.biogeme import BIOGEME\n",
    "from biogeme.expressions import Beta, Draws, MonteCarlo, exp, 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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "783376d2",
   "metadata": {},
   "source": [
    "Variables used for the specification of the Swissmetro model are defined in the file `swissmetro_variables.py`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4020f09e",
   "metadata": {},
   "outputs": [],
   "source": [
    "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": "cd27cfe0",
   "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": "8f0d2170",
   "metadata": {},
   "outputs": [],
   "source": [
    "logger = blog.get_screen_logger(level=blog.INFO)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11ea4d3b",
   "metadata": {},
   "source": [
    "**Tip:**<div class=\"alert alert-block alert-info\">It is advised to start working with a low number of draws, until\n",
    "the script is working well. Then, increase the number of draws to 10000, say.\n",
    "Then, execute the script overnight.</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22f69100",
   "metadata": {},
   "outputs": [],
   "source": [
    "number_of_draws = 10\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88d854b0",
   "metadata": {},
   "source": [
    "# Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "713d939a",
   "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": "5b6bd448",
   "metadata": {},
   "source": [
    "# Availability conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fcaf30f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "av = {1: TRAIN_AV_SP, 2: SM_AV, 3: CAR_AV_SP}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fff4cdf",
   "metadata": {},
   "source": [
    "# Logit model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8907bb1",
   "metadata": {},
   "source": [
    "## Utility functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acbdedcf",
   "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": "3314cec4",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10760eba",
   "metadata": {},
   "outputs": [],
   "source": [
    "logprob = loglogit(v, av, CHOICE)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b60401a0",
   "metadata": {},
   "source": [
    "## Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f2c8267",
   "metadata": {},
   "outputs": [],
   "source": [
    "biogeme = BIOGEME(database, logprob)\n",
    "biogeme.model_name = '01logit'\n",
    "results_logit: EstimationResults = biogeme.estimate(recycle=True)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3a9f160",
   "metadata": {},
   "source": [
    "## Results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9d7a81b",
   "metadata": {},
   "source": [
    "General statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "827d3e91",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results_logit.short_summary())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "071c197d",
   "metadata": {},
   "source": [
    "Estimated parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6618e724",
   "metadata": {},
   "outputs": [],
   "source": [
    "param_logit = get_pandas_estimated_parameters(estimation_results=results_logit)\n",
    "display(param_logit)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbc75fc1",
   "metadata": {},
   "source": [
    "# Random parameter: normal distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2d0ef37",
   "metadata": {},
   "source": [
    "Read the results from file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1541605d",
   "metadata": {},
   "outputs": [],
   "source": [
    "results_normal = EstimationResults.from_yaml_file(filename='02normal.yaml')\n",
    "param_normal = get_pandas_estimated_parameters(estimation_results=results_normal)\n",
    "display(param_normal)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb9193de",
   "metadata": {},
   "source": [
    "# Random parameter: lognormal distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2963bcec",
   "metadata": {},
   "outputs": [],
   "source": [
    "b_time_s = Beta('b_time_s', 1, None, None, 0)\n",
    "b_time_rnd = -exp(b_time + b_time_s * Draws('b_time_rnd', 'NORMAL'))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fdaf4c5",
   "metadata": {},
   "source": [
    "## Utility functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8862d278",
   "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": "631bc715",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5c22498",
   "metadata": {},
   "outputs": [],
   "source": [
    "prob = logit(v, av, CHOICE)\n",
    "logprob = log(MonteCarlo(prob))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c99ea4ab",
   "metadata": {},
   "source": [
    "## Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9f13076",
   "metadata": {},
   "outputs": [],
   "source": [
    "biogeme = BIOGEME(database, logprob, number_of_draws=number_of_draws)\n",
    "biogeme.model_name = '03lognormal'\n",
    "results_lognormal: EstimationResults = biogeme.estimate()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "521d919e",
   "metadata": {},
   "source": [
    "## Results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c744cf09",
   "metadata": {},
   "source": [
    "General statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "565937ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results_lognormal.short_summary())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5661e37e",
   "metadata": {},
   "source": [
    "Estimated parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e95a92ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "param_lognormal = get_pandas_estimated_parameters(estimation_results=results_lognormal)\n",
    "display(param_lognormal)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be986704",
   "metadata": {},
   "source": [
    "# Comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfb1dd9a",
   "metadata": {},
   "source": [
    "We build a summary data frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c37f334",
   "metadata": {},
   "outputs": [],
   "source": [
    "summary = compile_estimation_results(\n",
    "    {\n",
    "        'Logit': results_logit,\n",
    "        'Random param. (normal)': results_normal,\n",
    "        'Random param. (lognormal)': results_lognormal,\n",
    "    }\n",
    ")\n",
    "display(summary[0])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "121aee3e",
   "metadata": {},
   "source": [
    "The values of `b_time`and `b_time_s` cannot be directly compared. Indeed, in the case of the log normal distribution,\n",
    "they do not capture the mean and the standard deviation of the underlying distribution. We need to calculate them\n",
    "explicitly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cb64ee7",
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_normal = results_normal.get_parameter_value('b_time')\n",
    "stddev_normal = results_normal.get_parameter_value('b_time_s')\n",
    "\n",
    "location = results_lognormal.get_parameter_value('b_time')\n",
    "scale = results_lognormal.get_parameter_value('b_time_s')\n",
    "mean_lognormal = -np.exp(location + 0.5 * scale**2)\n",
    "variance_lognormal = np.exp(2 * location + scale**2) * (np.exp(scale**2) - 1)\n",
    "stddev_lognormal = np.sqrt(variance_lognormal)\n",
    "print(f'Mean log normal: {mean_lognormal:.3g}')\n",
    "print(f'Std. dev. log normal: {stddev_lognormal:.3g}')\n",
    "print(f'Mean normal: {mean_normal:.3g}')\n",
    "print(f'Std. dev. normal: {stddev_normal:.3g}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "455b139c",
   "metadata": {},
   "source": [
    "The mean of the lognormal is more negative than the mean of the normal. Also, the standard deviation is larger."
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}