{ "cells": [ { "cell_type": "markdown", "id": "460b4c98", "metadata": {}, "source": [ "\n", "File: 02-maximum_likelihood.py\n", "\n", "\n", "Michel Bierlaire\n", "\n", "Mon Aug 04 2025, 10:28:28\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "5b918a32", "metadata": {}, "outputs": [], "source": [ "\n", "import os\n", "import pickle\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "import tqdm\n", "from IPython.core.display_functions import display\n", "from biogeme.biogeme import BIOGEME\n", "from biogeme.database import Database\n", "from biogeme.models import loglogit\n", "from biogeme.results_processing import (\n", " EstimationResults,\n", " get_pandas_estimated_parameters,\n", ")\n", "from pandas import DataFrame, Series\n", "\n", "from spec_optima import Choice, av, v\n", "from stratified_sampling import Population, SamplingStrategy, SegmentSample\n" ] }, { "cell_type": "markdown", "id": "763c3f7e", "metadata": {}, "source": [ "The objective of this laboratory is to investigate maximum likelihood estimation under various sampling strategies." ] }, { "cell_type": "markdown", "id": "eb92a53c", "metadata": {}, "source": [ "The specification of the true model is available in `spec_optima.py`." ] }, { "cell_type": "markdown", "id": "56ccece4", "metadata": {}, "source": [ "The procedure to perform stratified sampling is implemented in the file `stratidied_sampling.py`." ] }, { "cell_type": "markdown", "id": "93f7022a", "metadata": {}, "source": [ "We implement the ingredients of the following experiment." ] }, { "cell_type": "markdown", "id": "9190b7e2", "metadata": {}, "source": [ "One instance of the experiment consists in:\n", "- sampling from the population using a given sampling strategy\n", "- estimate the unknown parameters using the generated sample.\n", "\n", "This is repeated several times. The mean and the standard deviation of the estimated value of each parameter across\n", "all instances are then calculated. As we know the true value of each parameter (that is, the value that has been\n", "used to generate the synthetic data), we can perform a formal t-test to test the hypothesis of the true value\n", "of the parameter." ] }, { "cell_type": "markdown", "id": "30baaf53", "metadata": {}, "source": [ "We provide below the full experiment for the simple random sampling strategy, as well as the functions implementing\n", "three other sampling strategies. And we ask you to run the same\n", "experiment for one exogenous and two endogenous sampling strategies, and to analyze and comment the results." ] }, { "cell_type": "markdown", "id": "3559e941", "metadata": {}, "source": [ "Full population" ] }, { "cell_type": "code", "execution_count": null, "id": "3e518b6b", "metadata": {}, "outputs": [], "source": [ "population_data_file = 'synthetic_population_with_choice.zip'\n", "population_data = pd.read_csv(population_data_file)\n", "population = Population(data=population_data, choice_variable='Choice')\n" ] }, { "cell_type": "markdown", "id": "0a82e110", "metadata": {}, "source": [ "We set the seed so that the results are reproducible" ] }, { "cell_type": "code", "execution_count": null, "id": "33fe6e3e", "metadata": {}, "outputs": [], "source": [ "np.random.seed(seed=90267)\n" ] }, { "cell_type": "markdown", "id": "e7f99575", "metadata": {}, "source": [ "We define the sample size for all experiments" ] }, { "cell_type": "code", "execution_count": null, "id": "74551f7b", "metadata": {}, "outputs": [], "source": [ "sample_size = 200\n" ] }, { "cell_type": "markdown", "id": "d8a43436", "metadata": {}, "source": [ "# Sampling procedures" ] }, { "cell_type": "markdown", "id": "ceea4534", "metadata": {}, "source": [ "## Simple random sampling" ] }, { "cell_type": "code", "execution_count": null, "id": "477c09f4", "metadata": {}, "outputs": [], "source": [ "\n", "srs_sampling = SamplingStrategy(\n", " name='srs', population=population, sample_size=sample_size, segments=None\n", ")\n", "\n" ] }, { "cell_type": "markdown", "id": "ac236017", "metadata": {}, "source": [ "## Exogenously stratified Sample" ] }, { "cell_type": "markdown", "id": "0dc19f3c", "metadata": {}, "source": [ "The strata are defined based on trip purpose and car availability. As the choice is not involved in the definition\n", "of strata, the stratified sampling is *exogenous*. Target shares in the sample: 25% for each stratum." ] }, { "cell_type": "markdown", "id": "88fe7ebc", "metadata": {}, "source": [ "We need to create mask generators implementing those segments for Pandas.\n", "purposes = {1: 'work', 2: 'other'}\n", "car_availability = {1: 'car_available', 3: 'car_unavailable'}" ] }, { "cell_type": "code", "execution_count": null, "id": "9a6698b5", "metadata": {}, "outputs": [], "source": [ "def mask_work_car_available(dataframe: DataFrame):\n", " return (dataframe['TripPurpose'] == 1) & (dataframe['CarAvail'] == 1)\n", "\n", "\n", "segment_work_car_available = SegmentSample(\n", " name=f'work_car_available',\n", " mask_generator=mask_work_car_available,\n", " target_share=0.25,\n", " target_share_alternatives=None,\n", ")\n", "\n", "\n", "def mask_work_car_unavailable(dataframe: DataFrame):\n", " return (dataframe['TripPurpose'] == 1) & (dataframe['CarAvail'] == 3)\n", "\n", "\n", "segment_work_car_unavailable = SegmentSample(\n", " name=f'work_car_unavailable',\n", " mask_generator=mask_work_car_unavailable,\n", " target_share=0.25,\n", " target_share_alternatives=None,\n", ")\n", "\n", "\n", "def mask_other_car_available(dataframe: DataFrame):\n", " return (dataframe['TripPurpose'] == 2) & (dataframe['CarAvail'] == 1)\n", "\n", "\n", "segment_other_car_available = SegmentSample(\n", " name=f'other_car_available',\n", " mask_generator=mask_other_car_available,\n", " target_share=0.25,\n", " target_share_alternatives=None,\n", ")\n", "\n", "\n", "def mask_other_car_unavailable(dataframe: DataFrame):\n", " return (dataframe['TripPurpose'] == 2) & (dataframe['CarAvail'] == 3)\n", "\n", "\n", "segment_other_car_unavailable = SegmentSample(\n", " name=f'other_car_unavailable',\n", " mask_generator=mask_other_car_unavailable,\n", " target_share=0.25,\n", " target_share_alternatives=None,\n", ")\n", "\n", "xss_segments = [\n", " segment_work_car_available,\n", " segment_work_car_unavailable,\n", " segment_other_car_available,\n", " segment_other_car_unavailable,\n", "]\n", "\n", "xss_sampling = SamplingStrategy(\n", " name='xss', population=population, sample_size=sample_size, segments=xss_segments\n", ")\n", "\n", "print(xss_sampling)\n" ] }, { "cell_type": "markdown", "id": "bbb53515", "metadata": {}, "source": [ "## Endogenously stratified Sample 1" ] }, { "cell_type": "markdown", "id": "41bc20da", "metadata": {}, "source": [ "The strata are defined based on car availability and the chosen alternative. As the choice is involved in the\n", "definition of strata, the stratified sampling is *endogenous*. Target shares in the sample: there are 5 non empty\n", "strata, taking 20% each" ] }, { "cell_type": "markdown", "id": "aa40aab8", "metadata": {}, "source": [ "Note that we define the masks for the exogenous strata. And within each of them, we specify the share of each\n", "alternative." ] }, { "cell_type": "code", "execution_count": null, "id": "dfa1749c", "metadata": {}, "outputs": [], "source": [ "\n", "\n", "def mask_car_available(dataframe: DataFrame):\n", " return dataframe['CarAvail'] == 1\n", "\n", "\n", "car_available_segment = SegmentSample(\n", " name=f'car_available',\n", " mask_generator=mask_car_available,\n", " target_share=0.6,\n", " target_share_alternatives={0: 1 / 3, 1: 1 / 3, 2: 1 / 3},\n", ")\n", "\n", "\n", "def mask_car_unavailable(dataframe: DataFrame):\n", " return dataframe['CarAvail'] != 1\n", "\n", "\n", "car_unavailable_segment = SegmentSample(\n", " name=f'car_unavailable',\n", " mask_generator=mask_car_unavailable,\n", " target_share=0.4,\n", " target_share_alternatives={0: 1 / 2, 1: 0, 2: 1 / 2},\n", ")\n", "\n", "\n", "ess_segments = [car_available_segment, car_unavailable_segment]\n", "\n", "\n", "ess_sampling = SamplingStrategy(\n", " name='ess', population=population, sample_size=sample_size, segments=ess_segments\n", ")\n", "print(ess_sampling)\n" ] }, { "cell_type": "markdown", "id": "1eb0c8b2", "metadata": {}, "source": [ "## Endogenously stratified Sample 2" ] }, { "cell_type": "markdown", "id": "ed879bc3", "metadata": {}, "source": [ "The strata are defined based on the trip purpose, the language and the chosen alternative. As the choice is\n", "involved in the definition of strata, the stratified sampling is *endogenous*. Target shares in the sample: there\n", "are 10 non empty strata, taking 10% each" ] }, { "cell_type": "code", "execution_count": null, "id": "80bef746", "metadata": {}, "outputs": [], "source": [ "\n", "\n", "def mask_french_car_available(dataframe: DataFrame):\n", " return (dataframe['LangCode'] == 1) & (dataframe['CarAvail'] == 1)\n", "\n", "\n", "french_car_available_segment = SegmentSample(\n", " name=f'french_car_available',\n", " mask_generator=mask_french_car_available,\n", " target_share=0.3,\n", " target_share_alternatives={0: 1 / 3, 1: 1 / 3, 2: 1 / 3},\n", ")\n", "\n", "\n", "def mask_german_car_available(dataframe: DataFrame):\n", " return (dataframe['LangCode'] == 2) & (dataframe['CarAvail'] == 1)\n", "\n", "\n", "german_car_available_segment = SegmentSample(\n", " name=f'german_car_available',\n", " mask_generator=mask_german_car_available,\n", " target_share=0.3,\n", " target_share_alternatives={0: 1 / 3, 1: 1 / 3, 2: 1 / 3},\n", ")\n", "\n", "\n", "def mask_french_car_unavailable(dataframe: DataFrame):\n", " return (dataframe['LangCode'] == 1) & (dataframe['CarAvail'] == 3)\n", "\n", "\n", "french_car_unavailable_segment = SegmentSample(\n", " name=f'french_car_unavailable',\n", " mask_generator=mask_french_car_unavailable,\n", " target_share=0.2,\n", " target_share_alternatives={0: 1 / 2, 1: 0, 2: 1 / 2},\n", ")\n", "\n", "\n", "def mask_german_car_unavailable(dataframe: DataFrame):\n", " return (dataframe['LangCode'] == 2) & (dataframe['CarAvail'] == 3)\n", "\n", "\n", "german_car_unavailable_segment = SegmentSample(\n", " name=f'german_car_unavailable',\n", " mask_generator=mask_german_car_unavailable,\n", " target_share=0.2,\n", " target_share_alternatives={0: 1 / 2, 1: 0, 2: 1 / 2},\n", ")\n", "\n", "ess2_segments = [\n", " french_car_available_segment,\n", " french_car_unavailable_segment,\n", " german_car_available_segment,\n", " german_car_unavailable_segment,\n", "]\n", "\n", "\n", "ess2_sampling = SamplingStrategy(\n", " name='ess2', population=population, sample_size=sample_size, segments=ess2_segments\n", ")\n", "\n" ] }, { "cell_type": "markdown", "id": "9c8cb51c", "metadata": {}, "source": [ "# Choice model" ] }, { "cell_type": "markdown", "id": "0af9aba0", "metadata": {}, "source": [ "The specification of the choice model is available in the file spec_optima.py. Note that it is the\n", "exact same specification that has been used to generate the synthetic choices." ] }, { "cell_type": "markdown", "id": "86795a12", "metadata": {}, "source": [ "The following function extracts a sample from the population and estimates the choice model using this sample.\n", "For each sample, we estimate the parameters with ESML." ] }, { "cell_type": "code", "execution_count": null, "id": "741675dc", "metadata": {}, "outputs": [], "source": [ "def sample_and_estimate(sampling_strategy: SamplingStrategy) -> Series:\n", " \"\"\"Estimate the choice model\n", "\n", " The specification of the choice model is available in the file\n", " spec_optima.py. Note that it is the exact same\n", " specification that has been used to generate the synthetic\n", " choices.\n", "\n", " This function extracts a sample from the population and\n", " estimates the choice model using this sample.\n", "\n", " :param sampling_strategy: sampling_strategy to apply\n", " :return: estimated parameters\n", " \"\"\"\n", " # Extract the sample using the provided function.\n", "\n", " the_sample = sampling_strategy.sample()\n", " if the_sample.sample.isna().any().any():\n", " nan_columns = the_sample.sample.columns[the_sample.sample.isna().any()].tolist()\n", " print(f'NaN values in column(s) {nan_columns}')\n", "\n", " # Prepare the data for Biogeme\n", " database = Database(the_sample.name, the_sample.sample)\n", "\n", " # Estimate the parameters\n", " logprob = loglogit(v, av, Choice)\n", " biogeme = BIOGEME(database, logprob, generate_html=False, generate_yaml=False)\n", " biogeme.model_name = the_sample.name\n", " results: EstimationResults = biogeme.quick_estimate()\n", "\n", " # Extract the estimated parameters\n", " pandas_results = get_pandas_estimated_parameters(estimation_results=results)\n", " return pandas_results['Value']\n", "\n" ] }, { "cell_type": "markdown", "id": "c78c419a", "metadata": {}, "source": [ "# Code the experiment" ] }, { "cell_type": "markdown", "id": "0f73e7a3", "metadata": {}, "source": [ "We first load the true values of the parameters" ] }, { "cell_type": "code", "execution_count": null, "id": "79ebcc06", "metadata": {}, "outputs": [], "source": [ "true_values_of_the_parameters = {\n", " 'asc_car': -1.8610422074360626,\n", " 'asc_car_german': -1.0188194471148966,\n", " 'asc_pt': -1.6550877823146726,\n", " 'asc_pt_car_unavail': 3.5259706461037417,\n", " 'asc_pt_german': 0.21100250245252905,\n", " 'beta_cost_car': -1.5653923336388653,\n", " 'beta_cost_pt': -0.45824318794797375,\n", " 'beta_dist': -1.3658487759164655,\n", " 'beta_dist_high_school': -1.0242967065023554,\n", " 'beta_dist_higher_education': -0.011682106903893624,\n", " 'beta_dist_university': 0.43619104548817356,\n", " 'beta_time': -1.6340227136243364,\n", " 'beta_time_others': -0.588784620423224,\n", " 'beta_time_part_time': -0.2778165351485947,\n", " 'beta_waiting': -0.22948530603741885,\n", " 'beta_waiting_not_work': -0.1501487521368328,\n", "}\n" ] }, { "cell_type": "markdown", "id": "1e3cf95d", "metadata": {}, "source": [ "We select some parameters that we will analyze more closely" ] }, { "cell_type": "code", "execution_count": null, "id": "5baf63bb", "metadata": {}, "outputs": [], "source": [ "selected_parameters = [\n", " 'beta_cost_car',\n", " 'beta_cost_pt',\n", " 'asc_car',\n", " 'asc_pt_car_unavail',\n", "]\n", "\n" ] }, { "cell_type": "markdown", "id": "058d3410", "metadata": {}, "source": [ "The following procedure implements the experiment. It repeats sampling and estimation, stores all estimation results.\n", "When it is done, the mean and the standard deviation is calculated for each parameter, as well as the *t*-statistic." ] }, { "cell_type": "code", "execution_count": null, "id": "31014609", "metadata": {}, "outputs": [], "source": [ "def run_experiment(\n", " sampling_strategy: SamplingStrategy,\n", " repetitions: int,\n", ") -> tuple[DataFrame, DataFrame]:\n", " \"\"\"Function that implements the experiment.\n", "\n", " It repeats sampling and estimation, stores all estimation\n", " results. When it is done, the mean and the standard deviation is\n", " calculated for each parameter, as well as the *t*-statistic.\n", "\n", " :param sampling_strategy: sampling strategy to use\n", " :param repetitions: number of times the experiment is performed\n", " :return: two data frames: one with the complete estimation results\n", " (the estimate of each parameter for each repetition), and one\n", " with a comparison between the true value and the average of\n", " the estimated values.\n", " \"\"\"\n", " file_name = f'{sampling_strategy.name}_{repetitions}.pickle'\n", " if os.path.exists(file_name):\n", " with open(file_name, 'rb') as f:\n", " results, comparisons = pickle.load(f)\n", " print(f'Loaded cached results from {file_name}')\n", " return results, comparisons\n", "\n", " list_of_results = []\n", "\n", " # Loop to generate rows\n", " for _ in tqdm.tqdm(range(repetitions)):\n", " row = sample_and_estimate(sampling_strategy)\n", " list_of_results.append(row.to_frame().T)\n", "\n", " # Concatenate all the collected results into a data frame\n", " results = pd.concat(list_of_results, ignore_index=True)\n", "\n", " comparisons = pd.DataFrame.from_dict(\n", " true_values_of_the_parameters, orient='index', columns=['True']\n", " )\n", " # Set the columns of results to match the parameter names (comparisons.index)\n", " results.columns = comparisons.index\n", "\n", " comparisons['Estimated'] = results.mean()\n", " comparisons['StdDev'] = results.std()\n", " comparisons['t-test'] = (\n", " comparisons['Estimated'] - comparisons['True']\n", " ) / comparisons['StdDev']\n", "\n", " with open(file_name, 'wb') as f:\n", " pickle.dump((results, comparisons), f)\n", " print(f'Saved results to {file_name}')\n", " return results, comparisons\n", "\n" ] }, { "cell_type": "markdown", "id": "78196bf0", "metadata": {}, "source": [ "We also plot the distribution of the estimates, and compare it with the true value and the mean." ] }, { "cell_type": "code", "execution_count": null, "id": "db1eb44e", "metadata": {}, "outputs": [], "source": [ "def report(results: DataFrame, comparisons: DataFrame, parameters: list[str]) -> None:\n", " \"\"\"Prepare the plots for a selected list of parameters.\n", "\n", " We plot the distribution of the estimates, and compare it with the\n", " true value and the mean.\n", "\n", " :param results: results of the experiment\n", " :param comparisons: data frame containing the comparison between\n", " the true value and the empirical average of the estimates of\n", " each coefficient. Generated by the function\n", " :func:run_experiment.\n", " :param parameters: list of parameters for which we plot the distribution\n", " \"\"\"\n", " print('*** Results ***')\n", " display(results)\n", " print('*** Comparisons ***')\n", " display(comparisons)\n", " for param in parameters:\n", " xrangespan = 5\n", " plt.figure()\n", " plt.title(param)\n", " results[param].plot.hist(\n", " bins=20,\n", " range=[\n", " comparisons.loc[param, 'Estimated']\n", " - xrangespan * comparisons.loc[param, 'StdDev'],\n", " comparisons.loc[param, 'Estimated']\n", " + xrangespan * comparisons.loc[param, 'StdDev'],\n", " ],\n", " )\n", " plt.axvline(\n", " true_values_of_the_parameters[param],\n", " color='r',\n", " linestyle='dashed',\n", " label='True Value',\n", " )\n", " plt.axvline(\n", " comparisons.loc[param, 'Estimated'],\n", " color='g',\n", " linestyle='dotted',\n", " label='Estimated Value',\n", " )\n", " plt.legend()\n", "\n" ] }, { "cell_type": "markdown", "id": "76476b40", "metadata": {}, "source": [ "# We have prepared all the functions. We now run the experiment" ] }, { "cell_type": "markdown", "id": "36904b8f", "metadata": {}, "source": [ "First, we define the number of repetitions." ] }, { "cell_type": "code", "execution_count": null, "id": "652dc0c8", "metadata": {}, "outputs": [], "source": [ "the_repetitions = 200\n" ] }, { "cell_type": "markdown", "id": "eed7b6e3", "metadata": {}, "source": [ "## Simple random sampling" ] }, { "cell_type": "markdown", "id": "a80f304c", "metadata": {}, "source": [ "We run the experiment." ] }, { "cell_type": "code", "execution_count": null, "id": "a60d36e9", "metadata": {}, "outputs": [], "source": [ "results_srs, comparisons_srs = run_experiment(\n", " sampling_strategy=srs_sampling,\n", " repetitions=the_repetitions,\n", ")\n" ] }, { "cell_type": "markdown", "id": "33253a8a", "metadata": {}, "source": [ "The complete results of the estimations are available in the following table." ] }, { "cell_type": "code", "execution_count": null, "id": "d5cf4396", "metadata": {}, "outputs": [], "source": [ "display(results_srs)\n" ] }, { "cell_type": "markdown", "id": "ddff057a", "metadata": {}, "source": [ "The average estimated value of the parameters is therefore:" ] }, { "cell_type": "code", "execution_count": null, "id": "44a3bf3d", "metadata": {}, "outputs": [], "source": [ "display(results_srs.mean())\n" ] }, { "cell_type": "markdown", "id": "fc3a6d2e", "metadata": {}, "source": [ "In order to compare those values with the true value of the parameters (that is, those used to generate the data), we\n", "use the comparison table" ] }, { "cell_type": "code", "execution_count": null, "id": "48cb9f5c", "metadata": {}, "outputs": [], "source": [ "display(comparisons_srs)\n" ] }, { "cell_type": "markdown", "id": "251b1629", "metadata": {}, "source": [ "The results above show the good adequacy between the mean of the estimated parameters and the true value.\n", "In particular, the hypothesis that the true value of the parameter is what it is cannot be rejected for any parameter.\n", "And when the difference is large (for instance, look at `BETA_DIST_high_school`) it corresponds to a parameter\n", "estimated with high standard error, that is with low precision." ] }, { "cell_type": "markdown", "id": "4200afea", "metadata": {}, "source": [ "Note that the share of individuals in the population with no car available is about 5%. The parameters associated\n", "with this group are estimated will less precision than the others." ] }, { "cell_type": "markdown", "id": "91f06871", "metadata": {}, "source": [ "Increasing the sample size and the number of repetitions would increase the precision. We have kept those values\n", "relatively low for the sake of computational time." ] }, { "cell_type": "markdown", "id": "8549c0bd", "metadata": {}, "source": [ "We plot the results for selected parameters." ] }, { "cell_type": "code", "execution_count": null, "id": "1984411e", "metadata": {}, "outputs": [], "source": [ "report(results=results_srs, comparisons=comparisons_srs, parameters=selected_parameters)\n" ] }, { "cell_type": "markdown", "id": "2bd09c72", "metadata": {}, "source": [ "We observe that the true value of the parameters are fairly well recovered." ] }, { "cell_type": "markdown", "id": "69d51465", "metadata": {}, "source": [ "Now, run the same experiment for the following sampling strategies:\n", "- exogenous sampling defined by the function `xss`,\n", "- endogenous sampling defined by the function `ess1`,\n", "- endogenous sampling defined by the function `ess2`." ] }, { "cell_type": "markdown", "id": "c800338f", "metadata": {}, "source": [ "Analyze the results in light of the theory seen in class ." ] } ], "metadata": {}, "nbformat": 4, "nbformat_minor": 5 }