Estimation

book/_toc.yml

@@ -2,13 +2,23 @@ format: jb-book
 root: intro.md
 
 parts:
-  - caption: Subjects
+  - caption: Contents
     chapters:
     - file: simulation/overview.md
       sections:
       - file: simulation/Simulating_Dice.ipynb
       - file: simulation/Birth_Month_Simulation.ipynb 
-  - caption: Miscallaneous
+    - file: estimation/overview.md
+      sections:
+      - file: estimation/Birth_Month_Estimation.ipynb
+      - file: estimation/911_Calls_Estimation.ipynb
+      - file: estimation/Maximum_Likelihood_Estimate.ipynb
+      - file: estimation/Bias.ipynb
+      - file: estimation/Mean_Squared_Error.ipynb
+      - file: estimation/Confidence_Intervals.ipynb
+      - file: estimation/avg_sd_correlation.ipynb
+      - file: estimation/ecdf_kde.ipynb
+  - caption: Miscellaneous
     chapters:
     - file: references.md
     #- file: changelog.md

book/estimation/911_Calls_Estimation.ipynb

@@ -0,0 +1,286 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "9f0496fc",
+   "metadata": {},
+   "source": [
+    "# 911 Calls\n",
+    "\n",
+    "In this section we will be using data on 911 calls in Fernhaven to estimate the probability of having no 911 calls in a minute. \n",
+    "\n",
+    "<!-- This distribution is modeled using a Poisson distribution with probability mass function \n",
+    "\n",
+    "$$\n",
+    "    p(k)=\\frac{\\mu^k}{k!}e^{-\\mu}\n",
+    "$$ \n",
+    "where $\\mu$ is the expected number of calls in one minute. -->\n",
+    "\n",
+    "We will be covering two methods of estimation: proportion and method of moments/maximum likelihood estimate (MLE)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "35f13c30",
+   "metadata": {},
+   "source": [
+    "```{note}\n",
+    "In this case the method of moments and maximum likelihood estimator coincide. However, this is **not** always the case.\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1901c494",
+   "metadata": {},
+   "source": [
+    "The first step is to perform *data wrangling* on the dataset to get several datasets of the number of calls per minute with $n=60$ (i.e. an hour long). In order to assume the same distribution across these datasets we will use the same time of day for each dataset."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7856f2b4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import matplotlib.pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c0d2f978",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Read in the data and extract calls between 6 and 7\n",
+    "df = pd.read_csv(\"data/911_Calls.csv\")\n",
+    "\n",
+    "# Convert 'Date' and 'Call time' to a datetime object\n",
+    "df[\"datetime\"] = pd.to_datetime(df[\"Date\"] + \" \" + df[\"Call time\"])\n",
+    "\n",
+    "# Extract date, hour, and minute from the datetime object\n",
+    "df[\"date\"] = df[\"datetime\"].dt.date.astype(str)\n",
+    "df[\"hour\"] = df[\"datetime\"].dt.hour\n",
+    "df[\"minute\"] = df[\"datetime\"].dt.minute\n",
+    "\n",
+    "# Filter calls that are between 6:00 and 6:59\n",
+    "df_hour = df[df[\"hour\"] == 6]\n",
+    "\n",
+    "# Group by date and minute, counting the number of calls\n",
+    "counts = df_hour.groupby([\"date\", \"minute\"]).size().reset_index(name=\"count\")\n",
+    "\n",
+    "hourly_counts = {}\n",
+    "for date, grp in counts.groupby(\"date\"):\n",
+    "    series = grp.set_index(\"minute\")[\"count\"]\n",
+    "    full = series.reindex(range(60), fill_value=0)\n",
+    "    hourly_counts[date] = full.tolist()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "74f8f04d",
+   "metadata": {},
+   "source": [
+    "Now that we have our datasets we can compute our estimates. \n",
+    "\n",
+    "The first estimate we will compute is the proportion of minutes with no calls in each dataset ($\\hat{p}$). \n",
+    "\n",
+    "$$\n",
+    "    \\hat{p}=\\frac{\\text{number of minutes with 0 calls}}{60}\n",
+    "$$\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ea3c1dad",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "p_hat = [counts.count(0) / 60.0 for counts in hourly_counts.values()]\n",
+    "\n",
+    "p_hat"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b270c19e",
+   "metadata": {},
+   "source": [
+    "The second estimate is the MLE ($e^{-\\hat{\\mu}}$). To do so, we need to calculate an estimate for $\\mu$. For each dataset we calculate $\\hat{\\mu}$ with the formula: \n",
+    "\n",
+    "$$\n",
+    "    \\hat{\\mu}=\\frac{\\text{number of calls}}{60}=\\frac{\\text{sum of calls per minute}}{60}\n",
+    "$$\n",
+    "\n",
+    "Once we have calculated an estimate for $\\mu$, we can use the probability mass function of the Poisson distribution where $k=0$ to calculate $e^{-\\hat{\\mu}}$. \n",
+    "<!-- \n",
+    "$$\n",
+    "    \\hat{p}=\\frac{\\mu^0}{0!}e^{-\\mu}\n",
+    "$$ \n",
+    "\n",
+    "$$\n",
+    "    \\hat{p}=e^{-\\mu}\n",
+    "$$ -->"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e940e960",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mle = [np.exp(-1 * (sum(counts) / 60.0)) for counts in hourly_counts.values()]\n",
+    "\n",
+    "mle"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "946682b2",
+   "metadata": {},
+   "source": [
+    "We can now make histograms of the obtained estimates for both cases to estimate the distribution and indicate the deviation from the true parameter. The true parameter in this case will be obtained using the full data-set. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "91407ce1",
+   "metadata": {},
+   "source": [
+    "```{note}\n",
+    "Although not exact, using the full dataset will be a very very good estimate of the true parameter. For this reason we will use it as the true value.\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7b91c8d1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# compute the true parameter over the full dataset\n",
+    "call_counts_full_dataset = (\n",
+    "    df.groupby([\"date\", \"hour\", \"minute\"]).size().reset_index(name=\"count\")\n",
+    ")\n",
+    "full_counts = []\n",
+    "for (_, _), group in call_counts_full_dataset.groupby([\"date\", \"hour\"]):\n",
+    "    series = group.set_index(\"minute\")[\"count\"]\n",
+    "    full_counts.extend(series.reindex(range(60), fill_value=0).tolist())\n",
+    "\n",
+    "true_parameter = full_counts.count(0) / len(full_counts)\n",
+    "print(f\"Estimated p (full dataset) = {true_parameter}\")\n",
+    "\n",
+    "# plot estimates\n",
+    "plt.figure(figsize=(15, 6))\n",
+    "plt.subplot(1, 2, 1)\n",
+    "plt.hist(p_hat)\n",
+    "plt.title(\"Estimation $\\hat{p}$\")\n",
+    "plt.ylabel(\"frequency\")\n",
+    "\n",
+    "plt.axvline(x=true_parameter, color=\"r\", linestyle=\"--\")\n",
+    "\n",
+    "plt.subplot(1, 2, 2)\n",
+    "plt.hist(mle)\n",
+    "plt.title(\"Estimation $e^{-\\mu}$\")\n",
+    "plt.ylabel(\"frequency\")\n",
+    "\n",
+    "plt.axvline(x=true_parameter, color=\"r\", linestyle=\"--\")\n",
+    "\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fa8fe489",
+   "metadata": {},
+   "source": [
+    "Given these two estimators, we need to determine the best one. $\\hat{p}$ is an unbiased estimator, but $e^{-\\hat{\\mu}}$ is positively biased. From the histograms, we also observe that $\\hat{p}$ has a larger variance than $e^{-\\hat{\\mu}}$. This translates to being typically far away from the true value versus being typically close to a value above the true value. We typically want to select the estimator with the lowest mean squared error (MSE).\n",
+    "\n",
+    "Using the true value from the dataset we can estimate the MSE value for both estimators. Our true value will be the calculated by considering all phone calls taken between 6 and 7. Using the MSE and bias we can calculate the variance of both estimators due to the fact that\n",
+    "\n",
+    "$$\n",
+    "\\text{MSE} - \\text{bias}^2 = \\text{variance}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7f707b53",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# calculate the true parameter over the full dataset\n",
+    "df_full = pd.read_csv(\"data/911_Calls.csv\")\n",
+    "df_full[\"datetime\"] = pd.to_datetime(df_full[\"Date\"] + \" \" + df_full[\"Call time\"])\n",
+    "df_full[\"date\"] = df_full[\"datetime\"].dt.date.astype(str)\n",
+    "df_full[\"hour\"] = df_full[\"datetime\"].dt.hour\n",
+    "df_full[\"minute\"] = df_full[\"datetime\"].dt.minute\n",
+    "df_hour = df_full[df_full[\"hour\"] == 6]\n",
+    "counts = df_hour.groupby([\"date\", \"minute\"]).size().reset_index(name=\"count\")\n",
+    "\n",
+    "hourly_counts = {}\n",
+    "for date, grp in counts.groupby(\"date\"):\n",
+    "    series = grp.set_index(\"minute\")[\"count\"]\n",
+    "    full = series.reindex(range(60), fill_value=0)\n",
+    "    hourly_counts[date] = full.tolist()\n",
+    "\n",
+    "# calculate the number of times there were 0 calls in a minute\n",
+    "total_zero_minutes = np.sum([counts.count(0) for counts in hourly_counts.values()])\n",
+    "\n",
+    "# calculate the total number of minutes\n",
+    "total_minutes = len(hourly_counts) * 60\n",
+    "\n",
+    "# calculate the actual proportion of zero calls\n",
+    "actual_prop = total_zero_minutes / total_minutes\n",
+    "\n",
+    "# calculate MSE and bias for both estimators\n",
+    "mse_prop = np.sum((actual_prop - p_hat) ** 2) / len(p_hat)\n",
+    "bias_prop = np.mean(p_hat) - actual_prop\n",
+    "\n",
+    "mse_mle = np.sum((actual_prop - mle) ** 2) / len(mle)\n",
+    "bias_mle = np.mean(mle) - actual_prop\n",
+    "\n",
+    "print(f\"Mean Squared Error for proportion: {mse_prop}\")\n",
+    "print(f\"Mean Squared Error for MLE: {mse_mle}\")\n",
+    "\n",
+    "print(f\"Bias for proportion: {bias_prop}\")\n",
+    "print(f\"Bias for MLE: {bias_mle}\")\n",
+    "\n",
+    "\n",
+    "print(f\"MSE - Bias^2 for proportion: {mse_prop - bias_prop ** 2}\")\n",
+    "print(f\"MSE - Bias^2 for MLE: {mse_mle - bias_mle ** 2}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "66508de3",
+   "metadata": {},
+   "source": [
+    "```{note}\n",
+    "You may notice that the bias calculated for proportion is 0. This holds because the estimator is unbiased!\n",
+    "```"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "prob_stat_book",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.13.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}