{ "cells": [ { "cell_type": "markdown", "id": "b3c50f98", "metadata": {}, "source": [ "\n", "File: 04-asv_solution.py\n", "\n", "\n", "Michel Bierlaire\n", "\n", "Tue Aug 05 2025, 15:57:43\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "b01d13b0", "metadata": {}, "outputs": [], "source": [ "\n", "import biogeme.biogeme_logging as blog\n", "from IPython.core.display_functions import display\n", "from biogeme import models\n", "from biogeme.biogeme import BIOGEME\n", "from biogeme.expressions import Beta, Draws, Expression, MonteCarlo, log\n", "from biogeme.results_processing import (\n", " EstimationResults,\n", " get_pandas_estimated_parameters,\n", ")\n" ] }, { "cell_type": "markdown", "id": "45b8af60", "metadata": {}, "source": [ "The variables of the model are available from the file `airline_variables.py`." ] }, { "cell_type": "code", "execution_count": null, "id": "f10d5830", "metadata": {}, "outputs": [], "source": [ "from airline_variables import (\n", " Fare_1,\n", " Fare_2,\n", " Fare_3,\n", " Legroom_1,\n", " Legroom_2,\n", " Legroom_3,\n", " Opt1_SchedDelayEarly,\n", " Opt1_SchedDelayLate,\n", " Opt2_SchedDelayEarly,\n", " Opt2_SchedDelayLate,\n", " Opt3_SchedDelayEarly,\n", " Opt3_SchedDelayLate,\n", " TripTimeHours_1,\n", " TripTimeHours_2,\n", " TripTimeHours_3,\n", " chosenAlternative,\n", " database,\n", ")\n" ] }, { "cell_type": "markdown", "id": "76dd2538", "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": "01564a60", "metadata": {}, "outputs": [], "source": [ "logger = blog.get_screen_logger(level=blog.INFO)\n", "logger.info('Example b01logit.py')\n", "\n" ] }, { "cell_type": "markdown", "id": "efaf5619", "metadata": {}, "source": [ "**Tip:**
It is advised to start working with a low number of draws, until\n", "the script is working well.\n", "Then, increase the number of draws to 10_000, say, and set `recycle` to False (if `recycle` is set to True, and a\n", "pickle file with the results exists, the results are imported from the file, and the estimation is not performed.)\n", "Then, execute the script overnight.
" ] }, { "cell_type": "code", "execution_count": null, "id": "b1eada1e", "metadata": {}, "outputs": [], "source": [ "number_of_draws = 10\n", "recycle = False\n", "\n" ] }, { "cell_type": "markdown", "id": "d5a7cafa", "metadata": {}, "source": [ "# The model" ] }, { "cell_type": "markdown", "id": "4fd0dadf", "metadata": {}, "source": [ "Parameters" ] }, { "cell_type": "code", "execution_count": null, "id": "c5553eb2", "metadata": {}, "outputs": [], "source": [ "constant2 = Beta('constant2', 0, None, None, 0)\n", "constant3 = Beta('constant3', 0, None, None, 0)\n", "fare = Beta('fare', 0, None, None, 0)\n", "legroom = Beta('legroom', 0, None, None, 0)\n", "schedule_delay_early = Beta('schedule_delay_early', 0, None, None, 0)\n", "schedule_delay_late = Beta('schedule_delay_late', 0, None, None, 0)\n", "total_tt1 = Beta('total_tt1', 0, None, None, 0)\n", "total_tt2 = Beta('total_tt2', 0, None, None, 0)\n", "total_tt3 = Beta('total_tt3', 0, None, None, 0)\n", "\n" ] }, { "cell_type": "markdown", "id": "526e6b40", "metadata": {}, "source": [ "Error components" ] }, { "cell_type": "code", "execution_count": null, "id": "d8967055", "metadata": {}, "outputs": [], "source": [ "sigma_1 = Beta('sigma_1', 1, None, None, 0)\n", "ec_1 = sigma_1 * Draws('ec_1', 'NORMAL')\n", "sigma_2 = Beta('sigma_2', 1, None, None, 0)\n", "ec_2 = sigma_2 * Draws('ec_2', 'NORMAL')\n", "sigma_3 = Beta('sigma_3', 1, None, None, 0)\n", "ec_3 = sigma_3 * Draws('ec_3', 'NORMAL')\n", "\n" ] }, { "cell_type": "markdown", "id": "bf4bc81e", "metadata": {}, "source": [ "Utility functions. The normalization identifies the id of the alternative such that the error component is not\n", "included, which is equivalent to normalize its scale to zero." ] }, { "cell_type": "code", "execution_count": null, "id": "11869a92", "metadata": {}, "outputs": [], "source": [ "def utility(normalization: int | None) -> dict[int, Expression]:\n", " \"\"\"Utility functions.\n", "\n", " :param normalization: The normalization identifies the id of the\n", " alternative such that the error component is not included,\n", " which is equivalent to normalize its scale to zero. If None, no normalization is considered.\n", " :return: dict of utility functions\n", " \"\"\"\n", " opt1 = (\n", " fare * Fare_1\n", " + legroom * Legroom_1\n", " + schedule_delay_early * Opt1_SchedDelayEarly\n", " + schedule_delay_late * Opt1_SchedDelayLate\n", " + total_tt1 * TripTimeHours_1\n", " )\n", " if normalization != 1:\n", " opt1 += ec_1\n", " opt2 = (\n", " constant2\n", " + fare * Fare_2\n", " + legroom * Legroom_2\n", " + schedule_delay_early * Opt2_SchedDelayEarly\n", " + schedule_delay_late * Opt2_SchedDelayLate\n", " + total_tt2 * TripTimeHours_2\n", " )\n", " if normalization != 2:\n", " opt2 += ec_2\n", " opt3 = (\n", " constant3\n", " + fare * Fare_3\n", " + legroom * Legroom_3\n", " + schedule_delay_early * Opt3_SchedDelayEarly\n", " + schedule_delay_late * Opt3_SchedDelayLate\n", " + total_tt3 * TripTimeHours_3\n", " )\n", " if normalization != 3:\n", " opt3 += ec_3\n", " return {1: opt1, 2: opt2, 3: opt3}\n", "\n" ] }, { "cell_type": "markdown", "id": "bba5fb85", "metadata": {}, "source": [ "# Estimation" ] }, { "cell_type": "code", "execution_count": null, "id": "89641a41", "metadata": {}, "outputs": [], "source": [ "def estimate(\n", " normalization: int | None, init_betas: dict[str, float] | None\n", ") -> EstimationResults:\n", " \"\"\"Estimate a normalized model\n", "\n", " :param normalization: The normalization identifies the id of the\n", " alternative such that the error component is not included,\n", " which is equivalent to normalize its scale to zero. If None,\n", " the model is not normalized.\n", " :param init_betas: if not None, starting values for the estimation.\n", " :return: estimation results\n", " \"\"\"\n", " prob = models.logit(utility(normalization=normalization), None, chosenAlternative)\n", " logprob: Expression = log(MonteCarlo(prob))\n", " if init_betas is not None:\n", " logprob.change_init_values(init_betas)\n", " biogeme = BIOGEME(database, logprob, number_of_draws=number_of_draws)\n", " if normalization is None:\n", " biogeme.model_name = '04-asv'\n", " else:\n", " biogeme.model_name = f'04-asv_{normalization}'\n", " results: EstimationResults = biogeme.estimate(recycle=recycle)\n", " return results\n", "\n" ] }, { "cell_type": "markdown", "id": "2b785449", "metadata": {}, "source": [ "## Estimation without normalization" ] }, { "cell_type": "code", "execution_count": null, "id": "0ca8426c", "metadata": {}, "outputs": [], "source": [ "unnormalized_results: EstimationResults = estimate(normalization=None, init_betas=None)\n" ] }, { "cell_type": "markdown", "id": "4edc4ee7", "metadata": {}, "source": [ "General statistics" ] }, { "cell_type": "code", "execution_count": null, "id": "3c2a03d2", "metadata": {}, "outputs": [], "source": [ "print(unnormalized_results.short_summary())\n" ] }, { "cell_type": "markdown", "id": "6a329e48", "metadata": {}, "source": [ "Number of draws are reported" ] }, { "cell_type": "code", "execution_count": null, "id": "eee62068", "metadata": {}, "outputs": [], "source": [ "reported_number_of_draws = unnormalized_results.number_of_draws\n", "print(f'{reported_number_of_draws}')\n" ] }, { "cell_type": "markdown", "id": "bfa1acdd", "metadata": {}, "source": [ "Final log likelihood" ] }, { "cell_type": "code", "execution_count": null, "id": "c189aa9e", "metadata": {}, "outputs": [], "source": [ "LL = unnormalized_results.final_loglikelihood\n", "print(f'{LL:.6g}')\n" ] }, { "cell_type": "markdown", "id": "6bfaaa71", "metadata": {}, "source": [ "Estimated parameters" ] }, { "cell_type": "code", "execution_count": null, "id": "a2ec504b", "metadata": {}, "outputs": [], "source": [ "betas = get_pandas_estimated_parameters(estimation_results=unnormalized_results)\n", "display(betas)\n", "\n" ] }, { "cell_type": "markdown", "id": "17e1fb1f", "metadata": {}, "source": [ "If the number of draws is sufficiently high, we see that $\\sigma_2^2$ is lesser than $\\sigma_1^2$ and $\\sigma_3^3$.\n", "Therefore, the correct normalization consists in setting $\\sigma_2 = 0$." ] }, { "cell_type": "markdown", "id": "a10190a2", "metadata": {}, "source": [ "We extract the sigmas" ] }, { "cell_type": "code", "execution_count": null, "id": "3d873dd5", "metadata": {}, "outputs": [], "source": [ "sigma1_value = unnormalized_results.get_parameter_value('sigma_1')\n", "sigma2_value = unnormalized_results.get_parameter_value('sigma_2')\n", "sigma3_value = unnormalized_results.get_parameter_value('sigma_3')\n", "\n" ] }, { "cell_type": "markdown", "id": "b155d15e", "metadata": {}, "source": [ "Relevant quantities from the variance-covariance matrix: \\begin{align}s_3 &= \\sigma_1^2 + \\sigma_2^2 \\\\ s_2 &=\n", "\\sigma_1^2 + \\sigma_3^2 \\\\ s_1 &= \\sigma_2^2 + \\sigma_3^2\\end{align}" ] }, { "cell_type": "code", "execution_count": null, "id": "5b47bbba", "metadata": {}, "outputs": [], "source": [ "s1 = sigma2_value**2 + sigma3_value**2\n", "s2 = sigma1_value**2 + sigma3_value**2\n", "s3 = sigma1_value**2 + sigma2_value**2\n", "print(f's1 = {s1:.3g}')\n", "print(f's2 = {s2:.3g}')\n", "print(f's3 = {s3:.3g}')\n", "\n" ] }, { "cell_type": "markdown", "id": "ade24c8b", "metadata": {}, "source": [ "## Results with normalization" ] }, { "cell_type": "markdown", "id": "50cc7f5e", "metadata": {}, "source": [ "Normalize $\\sigma_1=0$." ] }, { "cell_type": "code", "execution_count": null, "id": "9a1fb123", "metadata": {}, "outputs": [], "source": [ "normalized_results_1: EstimationResults = estimate(\n", " normalization=1, init_betas=unnormalized_results.get_beta_values()\n", ")\n" ] }, { "cell_type": "markdown", "id": "1cdff4c4", "metadata": {}, "source": [ "General statistics" ] }, { "cell_type": "code", "execution_count": null, "id": "7edc4ac5", "metadata": {}, "outputs": [], "source": [ "print(normalized_results_1.short_summary())\n" ] }, { "cell_type": "markdown", "id": "c66dfcc2", "metadata": {}, "source": [ "Estimated parameters" ] }, { "cell_type": "code", "execution_count": null, "id": "ca28e0b1", "metadata": {}, "outputs": [], "source": [ "display(get_pandas_estimated_parameters(estimation_results=normalized_results_1))\n" ] }, { "cell_type": "markdown", "id": "72490821", "metadata": {}, "source": [ "Normalize $\\sigma_2=0$." ] }, { "cell_type": "code", "execution_count": null, "id": "deabba66", "metadata": {}, "outputs": [], "source": [ "normalized_results_2: EstimationResults = estimate(\n", " normalization=2, init_betas=unnormalized_results.get_beta_values()\n", ")\n" ] }, { "cell_type": "markdown", "id": "010d43ca", "metadata": {}, "source": [ "General statistics" ] }, { "cell_type": "code", "execution_count": null, "id": "e91078a3", "metadata": {}, "outputs": [], "source": [ "print(normalized_results_2.short_summary())\n" ] }, { "cell_type": "markdown", "id": "0b0d04e6", "metadata": {}, "source": [ "Estimated parameters" ] }, { "cell_type": "code", "execution_count": null, "id": "5628bf2b", "metadata": {}, "outputs": [], "source": [ "display(get_pandas_estimated_parameters(estimation_results=normalized_results_2))\n", "\n" ] }, { "cell_type": "markdown", "id": "648ea84a", "metadata": {}, "source": [ "Normalize $\\sigma_3=0$." ] }, { "cell_type": "code", "execution_count": null, "id": "66ca3b3d", "metadata": {}, "outputs": [], "source": [ "normalized_results_3: EstimationResults = estimate(\n", " normalization=3, init_betas=unnormalized_results.get_beta_values()\n", ")\n" ] }, { "cell_type": "markdown", "id": "126c5f6b", "metadata": {}, "source": [ "General statistics" ] }, { "cell_type": "code", "execution_count": null, "id": "7fd75942", "metadata": {}, "outputs": [], "source": [ "print(normalized_results_3.short_summary())\n" ] }, { "cell_type": "markdown", "id": "db707d2e", "metadata": {}, "source": [ "Estimated parameters" ] }, { "cell_type": "code", "execution_count": null, "id": "baeaed59", "metadata": {}, "outputs": [], "source": [ "display(get_pandas_estimated_parameters(estimation_results=normalized_results_3))\n", "\n", "\n", "normalized_results = [normalized_results_1, normalized_results_2, normalized_results_3]\n", "\n" ] }, { "cell_type": "markdown", "id": "98ae7679", "metadata": {}, "source": [ "We extract the final log likelihood and the three scale parameters" ] }, { "cell_type": "code", "execution_count": null, "id": "21317788", "metadata": {}, "outputs": [], "source": [ "def extract(the_results: EstimationResults) -> tuple[float, float, float, float]:\n", " \"\"\"Extracts the final log likelihood and the three scale parameters.\n", "\n", " :param the_results: estimation results.\n", " :return: log likelihood and the three scale parameters.\n", " \"\"\"\n", " the_final_ll = the_results.final_loglikelihood\n", " try:\n", " s_1 = the_results.get_parameter_value('sigma_1')\n", " except ValueError:\n", " s_1 = 0\n", " try:\n", " s_2 = the_results.get_parameter_value('sigma_2')\n", " except ValueError:\n", " s_2 = 0\n", "\n", " try:\n", " s_3 = the_results.get_parameter_value('sigma_3')\n", " except ValueError:\n", " s_3 = 0\n", " return the_final_ll, s_1, s_2, s_3\n", "\n" ] }, { "cell_type": "markdown", "id": "2332b992", "metadata": {}, "source": [ "We gather all the results together" ] }, { "cell_type": "code", "execution_count": null, "id": "afbf854c", "metadata": {}, "outputs": [], "source": [ "all_results = [extract(r) for r in normalized_results]\n", "for a_normalization in [1, 2, 3]:\n", " print(f'Results for sigma_{a_normalization} = 0')\n", " print(f'\\tFinal log likelihood: {all_results[a_normalization-1][0]:.6g}')\n", " print(f'\\tsigma_1: {all_results[a_normalization-1][1]:.3g}')\n", " print(f'\\tsigma_2: {all_results[a_normalization-1][2]:.3g}')\n", " print(f'\\tsigma_3: {all_results[a_normalization-1][3]:.3g}')\n", "\n" ] }, { "cell_type": "markdown", "id": "177d5e37", "metadata": {}, "source": [ "The following function calculates, for a given normalization, the relevant quantities from the variance-covariance\n", "matrix: \\begin{align}s_3 &= \\sigma_1^2 + \\sigma_2^2 \\\\ s_2 &= \\sigma_1^2 + \\sigma_3^2 \\\\ s_1 &=\n", "\\sigma_2^2 + \\sigma_3^2\\end{align}" ] }, { "cell_type": "code", "execution_count": null, "id": "e91d9719", "metadata": {}, "outputs": [], "source": [ "def relevant_quantities(normalization: int) -> tuple[float, float, float]:\n", " \"\"\"Calculates, for a given normalization, the relevant quantities\n", " from the variance-covariance matrix\n", "\n", " :param normalization: The normalization identifies the id of the\n", " alternative such that the error component is not included,\n", " which is equivalent to normalize its scale to zero. If None,\n", " the model is not normalized.^\n", " :return: the three relevant quantities.\n", " \"\"\"\n", " s_1 = (\n", " all_results[normalization - 1][2] ** 2 + all_results[normalization - 1][3] ** 2\n", " )\n", " s_2 = (\n", " all_results[normalization - 1][1] ** 2 + all_results[normalization - 1][3] ** 2\n", " )\n", " s_3 = (\n", " all_results[normalization - 1][1] ** 2 + all_results[normalization - 1][2] ** 2\n", " )\n", " return s_1, s_2, s_3\n", "\n" ] }, { "cell_type": "markdown", "id": "86fc97a3", "metadata": {}, "source": [ "We now compare the relevant quantities of the variance-covariance matrix" ] }, { "cell_type": "markdown", "id": "dc60c829", "metadata": {}, "source": [ "Correct normalization: $\\sigma_2 = 0$." ] }, { "cell_type": "code", "execution_count": null, "id": "fb2b5f4f", "metadata": {}, "outputs": [], "source": [ "s1_normalization_2, s2_normalization_2, s3_normalization_2 = relevant_quantities(\n", " normalization=2\n", ")\n", "print(f's1 = {s1_normalization_2:.3g}')\n", "print(f's2 = {s2_normalization_2:.3g}')\n", "print(f's3 = {s3_normalization_2:.3g}')\n" ] }, { "cell_type": "markdown", "id": "0233569f", "metadata": {}, "source": [ "Same quantities from the unnormalized model" ] }, { "cell_type": "code", "execution_count": null, "id": "eb63ba35", "metadata": {}, "outputs": [], "source": [ "print(f's1 = {s1:.3g}')\n", "print(f's2 = {s2:.3g}')\n", "print(f's3 = {s3:.3g}')\n" ] }, { "cell_type": "markdown", "id": "f8bb63b6", "metadata": {}, "source": [ "Difference" ] }, { "cell_type": "code", "execution_count": null, "id": "f6c25bcc", "metadata": {}, "outputs": [], "source": [ "print(f'delta s1 = {s1_normalization_2 - s1:.3g}')\n", "print(f'delta s2 = {s2_normalization_2 - s2:.3g}')\n", "print(f'delta s3 = {s3_normalization_2 - s3:.3g}')\n" ] }, { "cell_type": "markdown", "id": "1d678987", "metadata": {}, "source": [ "The differences are due to simulation errors." ] }, { "cell_type": "markdown", "id": "88d8a067", "metadata": {}, "source": [ "Incorrect normalization: $\\sigma_1 = 0$" ] }, { "cell_type": "code", "execution_count": null, "id": "7909b906", "metadata": {}, "outputs": [], "source": [ "s1_normalization_1, s2_normalization_1, s3_normalization_1 = relevant_quantities(\n", " normalization=1\n", ")\n", "print(f's1 = {s1_normalization_1:.3g}')\n", "print(f's2 = {s2_normalization_1:.3g}')\n", "print(f's3 = {s3_normalization_1:.3g}')\n", "\n" ] }, { "cell_type": "markdown", "id": "1c23db36", "metadata": {}, "source": [ "Difference with unnormalized model" ] }, { "cell_type": "code", "execution_count": null, "id": "f665939f", "metadata": {}, "outputs": [], "source": [ "print(f'delta s1 = {s1_normalization_1 - s1:.3g}')\n", "print(f'delta s2 = {s2_normalization_1 - s2:.3g}')\n", "print(f'delta s3 = {s3_normalization_1 - s3:.3g}')\n", "\n" ] }, { "cell_type": "markdown", "id": "6bf2d8fd", "metadata": {}, "source": [ "The level of magnitude of the differences is higher than for the correct normalization, and cannot be attributed\n", "to simulation error." ] }, { "cell_type": "markdown", "id": "ac25debf", "metadata": {}, "source": [ "Incorrect normalization: $\\sigma_3=0$" ] }, { "cell_type": "code", "execution_count": null, "id": "baca7d66", "metadata": {}, "outputs": [], "source": [ "s1_normalization_3, s2_normalization_3, s3_normalization_3 = relevant_quantities(\n", " normalization=3\n", ")\n", "print(f's1 = {s1_normalization_3:.3g}')\n", "print(f's2 = {s2_normalization_3:.3g}')\n", "print(f's3 = {s3_normalization_3:.3g}')\n", "\n" ] }, { "cell_type": "markdown", "id": "4d9ed4b2", "metadata": {}, "source": [ "Difference with unnormalized model" ] }, { "cell_type": "code", "execution_count": null, "id": "7e73ae47", "metadata": {}, "outputs": [], "source": [ "print(f'delta s1 = {s1_normalization_3 - s1:.3g}')\n", "print(f'delta s2 = {s2_normalization_3 - s2:.3g}')\n", "print(f'delta s3 = {s3_normalization_3 - s3:.3g}')\n" ] }, { "cell_type": "markdown", "id": "603d5d99", "metadata": {}, "source": [ "In this case, the differences are lower, and it is difficult to say if they are due to simulation or normalization\n", "error." ] } ], "metadata": {}, "nbformat": 4, "nbformat_minor": 5 }