{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3416d103",
   "metadata": {},
   "source": [
    "\n",
    "File: 09-lecture_solution.py\n",
    "\n",
    "\n",
    "Michel Bierlaire\n",
    "\n",
    "Thu Aug 07 2025, 08:47:55\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "555f4bf0",
   "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, 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": "1db8b61d",
   "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": "5c24da75",
   "metadata": {},
   "outputs": [],
   "source": [
    "from swissmetro_variables import (\n",
    "    BUSINESS,\n",
    "    CAR_AV_SP,\n",
    "    CAR_CO_SCALED,\n",
    "    CAR_TT_SCALED,\n",
    "    CHOICE,\n",
    "    FIRST,\n",
    "    GA,\n",
    "    LOW_INC,\n",
    "    MALE,\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": "b2874c23",
   "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": "7095c465",
   "metadata": {},
   "outputs": [],
   "source": [
    "logger = blog.get_screen_logger(level=blog.INFO)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b967852",
   "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": "de82e5c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "number_of_draws = 10\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f49b7962",
   "metadata": {},
   "source": [
    "# Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "915f0d5a",
   "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": "7919e164",
   "metadata": {},
   "source": [
    "# Availability conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9df1e24a",
   "metadata": {},
   "outputs": [],
   "source": [
    "av = {1: TRAIN_AV_SP, 2: SM_AV, 3: CAR_AV_SP}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a078ac78",
   "metadata": {},
   "source": [
    "# Logit model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bd41ae0",
   "metadata": {},
   "source": [
    "## Utility functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e0c694e",
   "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": "800e249f",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fb583dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "logprob = loglogit(v, av, CHOICE)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dedc3281",
   "metadata": {},
   "source": [
    "## Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e4fc39f",
   "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": "6f49aa55",
   "metadata": {},
   "source": [
    "## Results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4696429",
   "metadata": {},
   "source": [
    "General statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc3e76e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results_logit.short_summary())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0013fb40",
   "metadata": {},
   "source": [
    "Estimated parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52a91d49",
   "metadata": {},
   "outputs": [],
   "source": [
    "param_logit = get_pandas_estimated_parameters(estimation_results=results_logit)\n",
    "display(param_logit)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a4fe951",
   "metadata": {},
   "source": [
    "# Random parameter: normal distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "104a9117",
   "metadata": {},
   "source": [
    "Read the results from file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c7e7317",
   "metadata": {},
   "outputs": [],
   "source": [
    "results_normal = EstimationResults.from_yaml_file(filename='02normal.yaml')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef230c56",
   "metadata": {},
   "source": [
    "# Random parameter: lognormal distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca910176",
   "metadata": {},
   "source": [
    "Read the results from file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0447e948",
   "metadata": {},
   "outputs": [],
   "source": [
    "results_lognormal = EstimationResults.from_yaml_file(filename='03lognormal.yaml')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ac01890",
   "metadata": {},
   "source": [
    "# Latent classes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b04a0a39",
   "metadata": {},
   "source": [
    "Read the results from file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbebda6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "results_latent = EstimationResults.from_yaml_file(filename='04latentClass.yaml')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07cec85f",
   "metadata": {},
   "source": [
    "# Latent classes with class membership model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0066be2",
   "metadata": {},
   "source": [
    "Read the results from file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33557778",
   "metadata": {},
   "outputs": [],
   "source": [
    "results_latentsocio = EstimationResults.from_yaml_file(filename='05latentClass.yaml')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e45dda41",
   "metadata": {},
   "source": [
    "# Latent classes with random parameter"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "734d8c1c",
   "metadata": {},
   "source": [
    "We consider again two classes in the population. The first class of individuals have considered all variables\n",
    "when making their choice. For them, the specification of the utility function is the same as for the logit model,\n",
    "where the time coefficient is now distributed in the population."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e00bb54a",
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5464aa05",
   "metadata": {},
   "source": [
    "Utility function, class 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d289b51f",
   "metadata": {},
   "outputs": [],
   "source": [
    "v_train_class_1 = (\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_class_1 = (\n",
    "    b_time * SM_TT_SCALED + b_cost * SM_COST_SCALED + b_fr * SM_HE_SCALED\n",
    ")\n",
    "v_car_class_1 = asc_car + b_time * CAR_TT_SCALED + b_cost * CAR_CO_SCALED\n",
    "v_class_1 = {1: v_train_class_1, 2: v_swissmetro_class_1, 3: v_car_class_1}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a87dca6f",
   "metadata": {},
   "source": [
    "The second class of individuals ignored the travel time variable when making the choice. Therefore, this variable\n",
    "is removed from the utility function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be54ab05",
   "metadata": {},
   "outputs": [],
   "source": [
    "v_train_class_2 = asc_train + b_cost * TRAIN_COST_SCALED + b_fr * TRAIN_HE_SCALED\n",
    "v_swissmetro_class_2 = b_cost * SM_COST_SCALED + b_fr * SM_HE_SCALED\n",
    "v_car_class_2 = asc_car + b_cost * CAR_CO_SCALED\n",
    "v_class_2 = {1: v_train_class_2, 2: v_swissmetro_class_2, 3: v_car_class_2}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7287edb5",
   "metadata": {},
   "source": [
    "The following parameters are involved in the class membership model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9abc2355",
   "metadata": {},
   "outputs": [],
   "source": [
    "g_intercept = Beta('g_intercept', 0, None, None, 0)\n",
    "g_male = Beta('g_male', 0, None, None, 0)\n",
    "g_ga = Beta('g_ga', 0, None, None, 0)\n",
    "g_business = Beta('g_business', 0, None, None, 0)\n",
    "g_low_inc = Beta('g_low_inc', 0, None, None, 0)\n",
    "g_first = Beta('g_first', 0, None, None, 0)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c066f81",
   "metadata": {},
   "source": [
    "Class membership model. Note that `omega` can potentially take any real value. We have to transform it into a\n",
    "probability using the transform `1 / (1 + exp(omega))`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a22ef42",
   "metadata": {},
   "outputs": [],
   "source": [
    "omega = (\n",
    "    g_intercept\n",
    "    + g_male * MALE\n",
    "    + g_ga * GA\n",
    "    + g_business * BUSINESS\n",
    "    + g_low_inc * LOW_INC\n",
    "    + g_first * FIRST\n",
    ")\n",
    "prob_class_1 = 1 / (1 + exp(omega))\n",
    "prob_class_2 = 1 - prob_class_1\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "895cb463",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93bd47d1",
   "metadata": {},
   "source": [
    "We first calculate the choice probability for each class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b984973",
   "metadata": {},
   "outputs": [],
   "source": [
    "cond_choice_prob_class_1 = logit(v_class_1, av, CHOICE)\n",
    "choice_prob_class_1 = MonteCarlo(cond_choice_prob_class_1)\n",
    "choice_prob_class_2 = logit(v_class_2, av, CHOICE)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8bb230b",
   "metadata": {},
   "source": [
    "The choice probability is obtained by using the class membership model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c46e5289",
   "metadata": {},
   "outputs": [],
   "source": [
    "prob = prob_class_1 * choice_prob_class_1 + prob_class_2 * choice_prob_class_2\n",
    "logprob = log(prob)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a72f34e0",
   "metadata": {},
   "source": [
    "## Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f6d902a",
   "metadata": {},
   "outputs": [],
   "source": [
    "biogeme = BIOGEME(database, logprob, number_of_draws=number_of_draws)\n",
    "biogeme.model_name = '06mixedLatentClass'\n",
    "results_latentrandom: EstimationResults = biogeme.estimate()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21c0e5c4",
   "metadata": {},
   "source": [
    "## Results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb836772",
   "metadata": {},
   "source": [
    "General statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f393aeb8",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results_latentrandom.short_summary())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90306458",
   "metadata": {},
   "source": [
    "Estimated parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f511ca4",
   "metadata": {},
   "outputs": [],
   "source": [
    "param_latentrandom = get_pandas_estimated_parameters(\n",
    "    estimation_results=results_latentrandom\n",
    ")\n",
    "display(param_latentrandom)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db637001",
   "metadata": {},
   "source": [
    "# Comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3176223a",
   "metadata": {},
   "source": [
    "We build a summary data frame. We first gather the parameter estimates for each model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c1a4655",
   "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",
    "        'Latent class': results_latent,\n",
    "        'Latent class with class mbship': results_latentsocio,\n",
    "        'Latent with random param.': results_latentrandom,\n",
    "    }\n",
    ")\n",
    "display(summary[0])"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}