{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "571a13c5",
   "metadata": {},
   "source": [
    "\n",
    "File: 08-lecture_solution.py\n",
    "\n",
    "\n",
    "Michel Bierlaire\n",
    "\n",
    "Thu Aug 07 2025, 08:37:00\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f6e931b",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import itertools\n",
    "\n",
    "import biogeme.biogeme_logging as blog\n",
    "import pandas as pd\n",
    "from IPython.core.display_functions import display\n",
    "from biogeme.biogeme import BIOGEME\n",
    "from biogeme.expressions import Beta, Expression, 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",
    "from pandas import Series\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27bcb0bc",
   "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": "6077ce77",
   "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": "9c123acb",
   "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": "30b9e50b",
   "metadata": {},
   "outputs": [],
   "source": [
    "logger = blog.get_screen_logger(level=blog.INFO)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3d1847a",
   "metadata": {},
   "source": [
    "# Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcbfb62a",
   "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": "e66704a8",
   "metadata": {},
   "source": [
    "# Availability conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "721ff501",
   "metadata": {},
   "outputs": [],
   "source": [
    "av = {1: TRAIN_AV_SP, 2: SM_AV, 3: CAR_AV_SP}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7222e5f2",
   "metadata": {},
   "source": [
    "# Logit model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e54bfb96",
   "metadata": {},
   "source": [
    "## Utility functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b01b0c9d",
   "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": "13dff512",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f690a7db",
   "metadata": {},
   "outputs": [],
   "source": [
    "logprob = loglogit(V, av, CHOICE)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f79b8d3b",
   "metadata": {},
   "source": [
    "## Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3abb79e",
   "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": "69dd5cad",
   "metadata": {},
   "source": [
    "## Results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "968d90da",
   "metadata": {},
   "source": [
    "General statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14962e11",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results_logit.short_summary())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54e89c21",
   "metadata": {},
   "source": [
    "Estimated parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec186b57",
   "metadata": {},
   "outputs": [],
   "source": [
    "param_logit = get_pandas_estimated_parameters(estimation_results=results_logit)\n",
    "display(param_logit)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e01b48b",
   "metadata": {},
   "source": [
    "# Random parameter: normal distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9589c900",
   "metadata": {},
   "source": [
    "Read the results from file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "082baba3",
   "metadata": {},
   "outputs": [],
   "source": [
    "results_normal = EstimationResults.from_yaml_file(filename='02normal.yaml')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c10cb0dc",
   "metadata": {},
   "source": [
    "# Random parameter: lognormal distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b91a15ae",
   "metadata": {},
   "source": [
    "Read the results from file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b038d5c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "results_lognormal = EstimationResults.from_yaml_file(filename='03lognormal.yaml')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20ad5f0b",
   "metadata": {},
   "source": [
    "# Latent classes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbc9430c",
   "metadata": {},
   "source": [
    "Read the results from file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "441960e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "results_latent = EstimationResults.from_yaml_file(filename='04latentClass.yaml')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "024bf0fe",
   "metadata": {},
   "source": [
    "# Latent classes with class membership model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3ef3e08",
   "metadata": {},
   "source": [
    "We consider again two classes in the population. The first class of individuals have considered all variables when\n",
    "making their choice. For them, the specification of the utility function is the same as for the logit model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1cec84c",
   "metadata": {},
   "outputs": [],
   "source": [
    "v_train_class_1 = (\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_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": "7c1f88fc",
   "metadata": {},
   "source": [
    "The second class of individuals ignored the travel time variable when making the choice. Therefore, this variable is\n",
    "removed from the utility function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "982ed14b",
   "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": "033e0496",
   "metadata": {},
   "source": [
    "The following parameters are involved in the class membership model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5fc50247",
   "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": "5f740782",
   "metadata": {},
   "source": [
    "The following function returns the expressions for the class membership probabilities. If `value` is set to `True`,\n",
    "the values instead of the expressions are returned."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad194650",
   "metadata": {},
   "source": [
    "Note that `w` can potentially take any real value. We have to transform it into a probability using the\n",
    "transform `1 / (1 + exp(w))`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6ac1356",
   "metadata": {},
   "outputs": [],
   "source": [
    "def omega(\n",
    "    male: float | Expression,\n",
    "    ga: float | Expression,\n",
    "    business: float | Expression,\n",
    "    low_inc: float | Expression,\n",
    "    first: float | Expression,\n",
    "    estimates: Series | None = None,\n",
    ") -> tuple[float, float] | tuple[Expression, Expression]:\n",
    "    \"\"\"The following function returns the expressions for the class\n",
    "     membership probabilities. If `value` is set to `True`, the values\n",
    "    instead of the expressions are returned.\n",
    "\n",
    "    Note that `w` can potentially take any real value. We have to\n",
    "    transform it into a probability using the transform `1 / (1 +\n",
    "    exp(w))`\n",
    "\n",
    "    :param male: value or expression of the variable MALE\n",
    "    :param ga: value or expression of the variable GA\n",
    "    :param business: value or expression of the variable BUSINESS\n",
    "    :param low_inc: value or expression of the variable LOW_INC\n",
    "    :param first: value or expression of the variable FIRST\n",
    "    :param estimates: estimated value of the parameters. If None, the expression is built.\n",
    "    :return: values or expressions for the class membership probability\n",
    "        for each of the two classes.\n",
    "\n",
    "    \"\"\"\n",
    "    # We calculate the values of the parameters. They are provided either as parameters, or are available in the\n",
    "    # estimates.\n",
    "    g_intercept_w = g_intercept if estimates is None else estimates['g_intercept']\n",
    "    g_male_w = g_male if estimates is None else estimates['g_male']\n",
    "    g_ga_w = g_ga if estimates is None else estimates['g_ga']\n",
    "    g_business_w = g_business if estimates is None else estimates['g_business']\n",
    "    g_low_inc_w = g_low_inc if estimates is None else estimates['g_low_inc']\n",
    "    g_first_w = g_first if estimates is None else estimates['g_first']\n",
    "    w = (\n",
    "        g_intercept_w\n",
    "        + g_male_w * male\n",
    "        + g_ga_w * ga\n",
    "        + g_business_w * business\n",
    "        + g_low_inc_w * low_inc\n",
    "        + g_first_w * first\n",
    "    )\n",
    "    omega_1 = 1 / (1 + exp(w))\n",
    "    omega_2 = 1 - omega_1\n",
    "    if estimates is None:\n",
    "        return omega_1, omega_2\n",
    "    return omega_1.get_value(), omega_2.get_value()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0e16cf4",
   "metadata": {},
   "source": [
    "We use the function to obtain the expressions of the class membership probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "937684ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "prob_class_1, prob_class_2 = omega(MALE, GA, BUSINESS, LOW_INC, FIRST)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ab6e4c9",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c7cb1c9",
   "metadata": {},
   "source": [
    "We first calculate the choice probability for each class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2db0ffc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "choice_prob_class_1 = logit(v_class_1, av, CHOICE)\n",
    "choice_prob_class_2 = logit(v_class_2, av, CHOICE)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f558c271",
   "metadata": {},
   "source": [
    "The choice probability is obtained by using the class membership model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00daa011",
   "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": "54fb08a7",
   "metadata": {},
   "source": [
    "## Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef64e3f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "biogeme = BIOGEME(database, logprob)\n",
    "biogeme.model_name = '05latentClass'\n",
    "results_latentsocio: EstimationResults = biogeme.estimate()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24f61aec",
   "metadata": {},
   "source": [
    "## Results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c667a31b",
   "metadata": {},
   "source": [
    "General statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d32343c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results_latentsocio.short_summary())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06cc3646",
   "metadata": {},
   "source": [
    "Estimated parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6de197c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "param_latentsocio = get_pandas_estimated_parameters(\n",
    "    estimation_results=results_latentsocio\n",
    ")\n",
    "display(param_latentsocio)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ca160f4",
   "metadata": {},
   "source": [
    "# Membership probability"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4167f939",
   "metadata": {},
   "source": [
    "We use the `itertools.product` function to enumerate all the combinations of values of the binary variables. We\n",
    "also use the `omega` function defined above, that returns the class membership probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d6c0cd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "rows = []\n",
    "variables = [' MALE', ' GA', ' BUSINESS', ' LOW_INC', ' FIRST']\n",
    "estimates = param_latentsocio.set_index('Name')['Value']\n",
    "for x in itertools.product([0, 1], [0, 1], [0, 1], [0, 1], [0, 1]):\n",
    "    prob = omega(*x, estimates=estimates)\n",
    "    prob_dict = {f'Class {i+1}': v for i, v in enumerate(prob)}\n",
    "    vars_dict = dict(zip(variables, x))\n",
    "    row = {**prob_dict, **vars_dict}\n",
    "    rows.append(row)\n",
    "simulation = pd.DataFrame(rows)\n",
    "for i in range(2):\n",
    "    key = f'Class {i+1}'\n",
    "    simulation[key] = simulation[key].apply(lambda the_prob: f'{100 * the_prob:.1f}%')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "869f42bd",
   "metadata": {},
   "source": [
    "Here are the values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de9866f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "display(simulation)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cafc7ef",
   "metadata": {},
   "source": [
    "# Comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e655f690",
   "metadata": {},
   "source": [
    "We build a summary data frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7a619f0",
   "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",
    "    }\n",
    ")\n",
    "display(summary[0])"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}