{
"cells": [
{
"cell_type": "markdown",
"id": "18a303b0-cf44-4c68-9d49-ec524a677a0b",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Class 7\n",
"\n",
"## Mar. 4, 2025\n",
"\n",
"_Where we left off last week_\n",
"\n",
"- understand the role of empirical distributions (sample and statistic) and simulation\n",
"\n",
"_New topics for today_\n",
"\n",
"- Models\n",
"\n",
"- Comparing two samples\n",
"\n",
"- Causality\n"
]
},
{
"cell_type": "markdown",
"id": "ace25f10",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Sampling from a Population\n",
"\n",
"The law of averages also holds when the random sample is drawn from individuals in a large population.\n",
"\n",
"```\n",
" persdur Duration - Personal activities\n",
"\n",
" VALUE LABEL\n",
" 0 No time spent doing this activity\n",
" 9996 Valid skip\n",
" 9997 Don't know\n",
" 9998 Refusal\n",
" 9999 Not stated\n",
"\n",
"\n",
" luc_rst Population centre indicator\n",
"\n",
" VALUE LABEL\n",
" 1 Larger urban population centres (CMA/CA)\n",
" 2 Rural areas and small population centres (non CMA/CA)\n",
" 3 Prince Edward Island\n",
" 6 Valid skip\n",
" 7 Don't know\n",
" 8 Refusal\n",
" 9 Not stated\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "6391b618",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(17390, 3)\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/var/folders/0j/ybsv4ncn5w50v40vdh5jjlww0000gn/T/ipykernel_4201/752036752.py:1: DeprecationWarning: \n",
"Pyarrow will become a required dependency of pandas in the next major release of pandas (pandas 3.0),\n",
"(to allow more performant data types, such as the Arrow string type, and better interoperability with other libraries)\n",
"but was not found to be installed on your system.\n",
"If this would cause problems for you,\n",
"please provide us feedback at https://github.com/pandas-dev/pandas/issues/54466\n",
" \n",
" import pandas as pd\n"
]
},
{
"data": {
"text/html": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"persdur_zoom.sample(\n",
" 100, \n",
" replace=True\n",
" ).plot.hist(bins=10, edgecolor=\"white\", color=\"darkgrey\");"
]
},
{
"cell_type": "markdown",
"id": "7ab338c2",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- The values that occur with the least frequency are less likely to occur in small random samples.\n",
"- As the size of the random sample increases the sample will resemble the population, with high probability."
]
},
{
"cell_type": "markdown",
"id": "23dbe931",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Interactive exploration of empirical distribution of a sample"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b401d770",
"metadata": {},
"outputs": [],
"source": [
"from ipywidgets import interact\n",
"import ipywidgets as widgets\n",
"\n",
"def emp_hist_plot(n):\n",
" persdur_zoom.sample(n, replace = True).plot.hist(bins=10, edgecolor=\"white\", color=\"darkgrey\");\n",
"\n",
"interact(emp_hist_plot, n = widgets.IntSlider(min = 10, max= 500, step=50, value=10));"
]
},
{
"cell_type": "markdown",
"id": "786fc970",
"metadata": {},
"source": [
"The histogram of the \"population\"."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "a467477a",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAGdCAYAAADzOWwgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAnoklEQVR4nO3de3BUZZ7/8U+HpBMu6YSLSaAIkBVFIrcFB9LehptEyVogWCsjQhQcBza4hIyg1FowwuxGcQBxRHFXIbjKIuwMzgjDJQaIoyQggchFBx1lDbO5MYtJQ0aSkJzfH/7SZZOASdvJ6fC8X1VdZZ9++uTbp9rwru7THYdlWZYAAAAMFmL3AAAAAHYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYL9TuAdqD+vp6FRcXKzIyUg6Hw+5xAABAM1iWpfPnz6tXr14KCbn6a0AEUTMUFxcrPj7e7jEAAIAfzpw5o969e191DUHUDJGRkZK+PaAul8vmaQAAQHN4PB7Fx8d7/x2/GoKoGRreJnO5XAQRAADtTHNOd+GkagAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gCgKWZdk9Qou1x5kBALiSULsHgORwOJSfny+Px2P3KM3icrmUlJRk9xgAAAQMQRQkPB6PKioq7B4DAAAj8ZYZAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADBe0ATRs88+K4fDofT0dO+2ixcvKi0tTd27d1eXLl00depUlZWV+dyvqKhIKSkp6tSpk2JiYrRw4UJdunTJZ83+/fs1fPhwhYeHq3///srKymqDRwQAANqLoAiijz76SK+++qqGDBnis33BggV69913tXXrVuXm5qq4uFhTpkzx3l5XV6eUlBTV1NTowIED2rhxo7KysrRkyRLvmtOnTyslJUVjxoxRYWGh0tPT9eijj2r37t1t9vgAAEBwsz2ILly4oOnTp+s//uM/1LVrV+/2yspKvf7661q1apXGjh2rESNGaMOGDTpw4IDy8/MlSXv27NEnn3yiN998U8OGDdM999yj5cuXa+3ataqpqZEkrVu3TgkJCVq5cqUGDhyoefPm6f7779fq1attebwAACD42B5EaWlpSklJ0fjx4322FxQUqLa21mf7TTfdpD59+igvL0+SlJeXp8GDBys2Nta7Jjk5WR6PRydPnvSuuXzfycnJ3n00pbq6Wh6Px+cCAACuXaF2/vDNmzfryJEj+uijjxrdVlpaKqfTqejoaJ/tsbGxKi0t9a75bgw13N5w29XWeDweffPNN+rYsWOjn52ZmalnnnnG78cFAADaF9teITpz5ozmz5+vt956SxEREXaN0aTFixersrLSezlz5ozdIwEAgFZkWxAVFBSovLxcw4cPV2hoqEJDQ5Wbm6sXX3xRoaGhio2NVU1NjSoqKnzuV1ZWpri4OElSXFxco0+dNVz/vjUul6vJV4ckKTw8XC6Xy+cCAACuXbYF0bhx43T8+HEVFhZ6L7fccoumT5/u/e+wsDDl5OR473Pq1CkVFRXJ7XZLktxut44fP67y8nLvmuzsbLlcLiUmJnrXfHcfDWsa9gEAAGDbOUSRkZEaNGiQz7bOnTure/fu3u2zZ89WRkaGunXrJpfLpccff1xut1tJSUmSpAkTJigxMVEzZszQihUrVFpaqqefflppaWkKDw+XJM2ZM0cvvfSSFi1apFmzZmnv3r3asmWLduzY0bYPGAAABC1bT6r+PqtXr1ZISIimTp2q6upqJScn6+WXX/be3qFDB23fvl1z586V2+1W586dlZqaqmXLlnnXJCQkaMeOHVqwYIHWrFmj3r1767XXXlNycrIdDwkAAAQhh2VZlt1DBDuPx6OoqChVVla22vlEe/bsaXS+VLCKjo7WhAkT7B4DAICrasm/37Z/DxEAAIDdCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxrM1iF555RUNGTJELpdLLpdLbrdbO3fu9N5+8eJFpaWlqXv37urSpYumTp2qsrIyn30UFRUpJSVFnTp1UkxMjBYuXKhLly75rNm/f7+GDx+u8PBw9e/fX1lZWW3x8AAAQDthaxD17t1bzz77rAoKCnT48GGNHTtWkyZN0smTJyVJCxYs0LvvvqutW7cqNzdXxcXFmjJlivf+dXV1SklJUU1NjQ4cOKCNGzcqKytLS5Ys8a45ffq0UlJSNGbMGBUWFio9PV2PPvqodu/e3eaPFwAABCeHZVmW3UN8V7du3fT888/r/vvv13XXXadNmzbp/vvvlyT96U9/0sCBA5WXl6ekpCTt3LlT//AP/6Di4mLFxsZKktatW6cnn3xSZ8+eldPp1JNPPqkdO3boxIkT3p8xbdo0VVRUaNeuXc2ayePxKCoqSpWVlXK5XIF/0JL27NmjioqKVtl3oEVHR2vChAl2jwEAwFW15N/voDmHqK6uTps3b1ZVVZXcbrcKCgpUW1ur8ePHe9fcdNNN6tOnj/Ly8iRJeXl5Gjx4sDeGJCk5OVkej8f7KlNeXp7PPhrWNOyjKdXV1fJ4PD4XAABw7bI9iI4fP64uXbooPDxcc+bM0bZt25SYmKjS0lI5nU5FR0f7rI+NjVVpaakkqbS01CeGGm5vuO1qazwej7755psmZ8rMzFRUVJT3Eh8fH4iHCgAAgpTtQTRgwAAVFhbq4MGDmjt3rlJTU/XJJ5/YOtPixYtVWVnpvZw5c8bWeQAAQOsKtXsAp9Op/v37S5JGjBihjz76SGvWrNEDDzygmpoaVVRU+LxKVFZWpri4OElSXFycDh065LO/hk+hfXfN5Z9MKysrk8vlUseOHZucKTw8XOHh4QF5fAAAIPjZ/grR5err61VdXa0RI0YoLCxMOTk53ttOnTqloqIiud1uSZLb7dbx48dVXl7uXZOdnS2Xy6XExETvmu/uo2FNwz4AAABsfYVo8eLFuueee9SnTx+dP39emzZt0v79+7V7925FRUVp9uzZysjIULdu3eRyufT444/L7XYrKSlJkjRhwgQlJiZqxowZWrFihUpLS/X0008rLS3N+wrPnDlz9NJLL2nRokWaNWuW9u7dqy1btmjHjh12PnQAABBEbA2i8vJyzZw5UyUlJYqKitKQIUO0e/du3XXXXZKk1atXKyQkRFOnTlV1dbWSk5P18ssve+/foUMHbd++XXPnzpXb7Vbnzp2VmpqqZcuWedckJCRox44dWrBggdasWaPevXvrtddeU3Jycps/XgAAEJyC7nuIghHfQ+SL7yECALQH7fJ7iAAAAOxCEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIznVxB9+eWXgZ4DAADANn4FUf/+/TVmzBi9+eabunjxYqBnAgAAaFN+BdGRI0c0ZMgQZWRkKC4uTj/72c906NChQM8GAADQJvwKomHDhmnNmjUqLi7W+vXrVVJSottvv12DBg3SqlWrdPbs2UDPCQAA0Gp+0EnVoaGhmjJlirZu3arnnntOf/7zn/XEE08oPj5eM2fOVElJSaDmBAAAaDU/KIgOHz6sf/qnf1LPnj21atUqPfHEE/riiy+UnZ2t4uJiTZo0KVBzAgAAtJpQf+60atUqbdiwQadOndLEiRP1xhtvaOLEiQoJ+bavEhISlJWVpX79+gVyVgAAgFbhVxC98sormjVrlh5++GH17NmzyTUxMTF6/fXXf9BwAAAAbcGvIPr888+/d43T6VRqaqo/uwcAAGhTfp1DtGHDBm3durXR9q1bt2rjxo0/eCgAAIC25FcQZWZmqkePHo22x8TE6N/+7d9+8FAAAABtya8gKioqUkJCQqPtffv2VVFR0Q8eCmgNlmXZPUKLtceZAaA98uscopiYGB07dqzRp8g+/vhjde/ePRBzAQHncDiUn58vj8dj9yjN4nK5lJSUZPcYAGAEv4LoJz/5if75n/9ZkZGRuvPOOyVJubm5mj9/vqZNmxbQAYFA8ng8qqiosHsMAECQ8SuIli9frv/5n//RuHHjFBr67S7q6+s1c+ZMziECAADtjl9B5HQ69fbbb2v58uX6+OOP1bFjRw0ePFh9+/YN9HwAAACtzq8ganDjjTfqxhtvDNQsAAAAtvAriOrq6pSVlaWcnByVl5ervr7e5/a9e/cGZDgAAIC24FcQzZ8/X1lZWUpJSdGgQYPkcDgCPRcAAECb8SuINm/erC1btmjixImBngcAAKDN+fXFjE6nU/379w/0LAAAALbwK4h+/vOfa82aNXyLLgAAuCb49ZbZBx98oH379mnnzp26+eabFRYW5nP7b3/724AMBwAA0Bb8CqLo6Gjdd999gZ4FAADAFn4F0YYNGwI9BwAAgG38OodIki5duqT33ntPr776qs6fPy9JKi4u1oULFwI2HAAAQFvw6xWir776SnfffbeKiopUXV2tu+66S5GRkXruuedUXV2tdevWBXpOAACAVuPXK0Tz58/XLbfcoq+//lodO3b0br/vvvuUk5MTsOEAAADagl+vEP3xj3/UgQMH5HQ6fbb369dP//u//xuQwQAAANqKX68Q1dfXq66urtH2v/zlL4qMjPzBQwEAALQlv4JowoQJeuGFF7zXHQ6HLly4oKVLl/LnPAAAQLvj11tmK1euVHJyshITE3Xx4kU9+OCD+vzzz9WjRw/913/9V6BnBAAAaFV+BVHv3r318ccfa/PmzTp27JguXLig2bNna/r06T4nWQMAALQHfgWRJIWGhuqhhx4K5CwAAAC28CuI3njjjavePnPmTL+GAQAAsINfQTR//nyf67W1tfrb3/4mp9OpTp06EUQAAKBd8etTZl9//bXP5cKFCzp16pRuv/12TqoGAADtjt9/y+xyN9xwg5599tlGrx4BAAAEu4AFkfTtidbFxcWB3CUAAECr8+scot///vc+1y3LUklJiV566SXddtttARkMAACgrfgVRJMnT/a57nA4dN1112ns2LFauXJlIOYCAABoM34FUX19faDnAAAAsE1AzyECAABoj/x6hSgjI6PZa1etWuXPjwAAAGgzfgXR0aNHdfToUdXW1mrAgAGSpM8++0wdOnTQ8OHDvescDkdgpgQAAGhFfgXRvffeq8jISG3cuFFdu3aV9O2XNT7yyCO644479POf/zygQwIAALQmv84hWrlypTIzM70xJEldu3bVL3/5Sz5lBgAA2h2/gsjj8ejs2bONtp89e1bnz5//wUMhuEVERMiyLLvHAAAgYPx6y+y+++7TI488opUrV2rkyJGSpIMHD2rhwoWaMmVKQAdE8AkLC5PD4VB+fr48Ho/d4zRLXFychgwZYvcYAIAg5VcQrVu3Tk888YQefPBB1dbWfruj0FDNnj1bzz//fEAHRPDyeDyqqKiwe4xmiYyMtHsEAEAQ8yuIOnXqpJdfflnPP/+8vvjiC0nS9ddfr86dOwd0OAAAgLbwg76YsaSkRCUlJbrhhhvUuXNnzisBAADtkl9B9H//938aN26cbrzxRk2cOFElJSWSpNmzZ/ORewAA0O74FUQLFixQWFiYioqK1KlTJ+/2Bx54QLt27QrYcAAAAG3Br3OI9uzZo927d6t3794+22+44QZ99dVXARkMAACgrfj1ClFVVZXPK0MNzp07p/Dw8GbvJzMzUz/60Y8UGRmpmJgYTZ48WadOnfJZc/HiRaWlpal79+7q0qWLpk6dqrKyMp81RUVFSklJUadOnRQTE6OFCxfq0qVLPmv279+v4cOHKzw8XP3791dWVlbzHzAAALim+RVEd9xxh9544w3vdYfDofr6eq1YsUJjxoxp9n5yc3OVlpam/Px8ZWdnq7a2VhMmTFBVVZV3zYIFC/Tuu+9q69atys3NVXFxsc93HdXV1SklJUU1NTU6cOCANm7cqKysLC1ZssS75vTp00pJSdGYMWNUWFio9PR0Pfroo9q9e7c/Dx8AAFxj/HrLbMWKFRo3bpwOHz6smpoaLVq0SCdPntS5c+f04YcfNns/l59vlJWVpZiYGBUUFOjOO+9UZWWlXn/9dW3atEljx46VJG3YsEEDBw5Ufn6+kpKStGfPHn3yySd67733FBsbq2HDhmn58uV68skn9Ytf/EJOp1Pr1q1TQkKC98+KDBw4UB988IFWr16t5ORkfw4BAAC4hvj1CtGgQYP02Wef6fbbb9ekSZNUVVWlKVOm6OjRo7r++uv9HqayslKS1K1bN0lSQUGBamtrNX78eO+am266SX369FFeXp4kKS8vT4MHD1ZsbKx3TXJysjwej06ePOld8919NKxp2Mflqqur5fF4fC4AAODa1eJXiGpra3X33Xdr3bp1+pd/+ZeADVJfX6/09HTddtttGjRokCSptLRUTqdT0dHRPmtjY2NVWlrqXfPdGGq4veG2q63xeDz65ptv1LFjR5/bMjMz9cwzzwTssQEAgODW4leIwsLCdOzYsYAPkpaWphMnTmjz5s0B33dLLV68WJWVld7LmTNn7B4JAAC0Ir/eMnvooYf0+uuvB2yIefPmafv27dq3b5/PR/nj4uJUU1PT6O9llZWVKS4uzrvm8k+dNVz/vjUul6vRq0OSFB4eLpfL5XMBAADXLr9Oqr506ZLWr1+v9957TyNGjGj0N8xWrVrVrP1YlqXHH39c27Zt0/79+5WQkOBz+4gRIxQWFqacnBxNnTpVknTq1CkVFRXJ7XZLktxut/71X/9V5eXliomJkSRlZ2fL5XIpMTHRu+YPf/iDz76zs7O9+wAAAGZrURB9+eWX6tevn06cOKHhw4dLkj777DOfNQ6Ho9n7S0tL06ZNm/S73/1OkZGR3nN+oqKi1LFjR0VFRWn27NnKyMhQt27d5HK59Pjjj8vtdispKUmSNGHCBCUmJmrGjBlasWKFSktL9fTTTystLc37nUhz5szRSy+9pEWLFmnWrFnau3evtmzZoh07drTk4QMAgGtUi4LohhtuUElJifbt2yfp2z/V8eKLLzY6Ybm5XnnlFUnS6NGjfbZv2LBBDz/8sCRp9erVCgkJ0dSpU1VdXa3k5GS9/PLL3rUdOnTQ9u3bNXfuXLndbnXu3FmpqalatmyZd01CQoJ27NihBQsWaM2aNerdu7dee+01PnIPAAAktTCILv9r9jt37vT5EsWWunx/TYmIiNDatWu1du3aK67p27dvo7fELjd69GgdPXq0xTMCAIBrn18nVTdoTtAAAAAEuxYFkcPhaHSOUEvOGQIAAAhGLX7L7OGHH/aerHzx4kXNmTOn0afMfvvb3wZuQgAAgFbWoiBKTU31uf7QQw8FdBgAAAA7tCiINmzY0FpzAAAA2OYHnVQNAABwLSCIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiIAgFRERIcuy7B6jxdrjzAAQavcAAJoWFhYmh8Oh/Px8eTweu8dpFpfLpaSkJLvHAIAWI4iAIOfxeFRRUWH3GABwTeMtMwAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxrM1iN5//33de++96tWrlxwOh9555x2f2y3L0pIlS9SzZ0917NhR48eP1+eff+6z5ty5c5o+fbpcLpeio6M1e/ZsXbhwwWfNsWPHdMcddygiIkLx8fFasWJFaz80AADQjtgaRFVVVRo6dKjWrl3b5O0rVqzQiy++qHXr1ungwYPq3LmzkpOTdfHiRe+a6dOn6+TJk8rOztb27dv1/vvv67HHHvPe7vF4NGHCBPXt21cFBQV6/vnn9Ytf/EL//u//3uqPDwAAtA+hdv7we+65R/fcc0+Tt1mWpRdeeEFPP/20Jk2aJEl64403FBsbq3feeUfTpk3Tp59+ql27dumjjz7SLbfcIkn69a9/rYkTJ+pXv/qVevXqpbfeeks1NTVav369nE6nbr75ZhUWFmrVqlU+4QQAAMwVtOcQnT59WqWlpRo/frx3W1RUlEaNGqW8vDxJUl5enqKjo70xJEnjx49XSEiIDh486F1z5513yul0etckJyfr1KlT+vrrr5v82dXV1fJ4PD4XAABw7QraICotLZUkxcbG+myPjY313lZaWqqYmBif20NDQ9WtWzefNU3t47s/43KZmZmKioryXuLj43/4AwIAAEEraIPITosXL1ZlZaX3cubMGbtHAgAArShogyguLk6SVFZW5rO9rKzMe1tcXJzKy8t9br906ZLOnTvns6apfXz3Z1wuPDxcLpfL5wIAAK5dQRtECQkJiouLU05Ojnebx+PRwYMH5Xa7JUlut1sVFRUqKCjwrtm7d6/q6+s1atQo75r3339ftbW13jXZ2dkaMGCAunbt2kaPBgAABDNbg+jChQsqLCxUYWGhpG9PpC4sLFRRUZEcDofS09P1y1/+Ur///e91/PhxzZw5U7169dLkyZMlSQMHDtTdd9+tn/70pzp06JA+/PBDzZs3T9OmTVOvXr0kSQ8++KCcTqdmz56tkydP6u2339aaNWuUkZFh06MGAADBxtaP3R8+fFhjxozxXm+IlNTUVGVlZWnRokWqqqrSY489poqKCt1+++3atWuXIiIivPd56623NG/ePI0bN04hISGaOnWqXnzxRe/tUVFR2rNnj9LS0jRixAj16NFDS5Ys4SP3AADAy9YgGj16tCzLuuLtDodDy5Yt07Jly664plu3btq0adNVf86QIUP0xz/+0e85AQDAtS1ozyECAABoKwQRgICJiIi46qu+wao9zgwgsGx9ywzAtSUsLEwOh0P5+fnt5hveXS6XkpKS7B4DgM0IIgAB5/F4VFFRYfcYANBsvGUGAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABQDtkWZbdI/ilvc6Na1+o3QMAAFrO4XAoPz9fHo/H7lGazeVyKSkpye4xgCYRRADQTnk8HlVUVNg9BnBN4C0zAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAJgtIiICL4sEADfQwTAbGFhYe3uSw7j4uI0ZMgQu8cArikEEQCofX3JYWRkpN0jANcc3jIDAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAuArLsuweocXa48x2C7V7AAAAgpnD4VB+fr48Ho/dozSLy+VSUlKS3WO0OwQRAKBNREREyLIsORwOu0dpMY/Ho4qKCrvHaJb2epztnpkgAgC0ibCwsHb3aktcXJyGDBli9xgt0h6PczC8qkUQAQDaVHt6tSUyMtLuEfzWno5zMOCkagAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPKOCaO3aterXr58iIiI0atQoHTp0yO6RAABAEDAmiN5++21lZGRo6dKlOnLkiIYOHark5GSVl5fbPRoAALCZMUG0atUq/fSnP9UjjzyixMRErVu3Tp06ddL69evtHg0AANjMiC9mrKmpUUFBgRYvXuzdFhISovHjxysvL6/R+urqalVXV3uvV1ZWSlKrfuNnhw4dFBYW1mr7DyTLsuTxeJi5lTFz22DmttMe52bmttGhQ4dW+Te2YZ/N+WO3RgTRX//6V9XV1Sk2NtZne2xsrP70pz81Wp+Zmalnnnmm0fb4+PhWmxEAALSO8+fPKyoq6qprjAiillq8eLEyMjK81+vr63Xu3Dl179693f2xvEDweDyKj4/XmTNn5HK57B6n3eI4BgbHMTA4joHBcQyM1jqOlmXp/Pnz6tWr1/euNSKIevTooQ4dOqisrMxne1lZmeLi4hqtDw8PV3h4uM+26Ojo1hyxXXC5XPwPHwAcx8DgOAYGxzEwOI6B0RrH8fteGWpgxEnVTqdTI0aMUE5OjndbfX29cnJy5Ha7bZwMAAAEAyNeIZKkjIwMpaam6pZbbtHIkSP1wgsvqKqqSo888ojdowEAAJsZE0QPPPCAzp49qyVLlqi0tFTDhg3Trl27Gp1ojcbCw8O1dOnSRm8jomU4joHBcQwMjmNgcBwDIxiOo8NqzmfRAAAArmFGnEMEAABwNQQRAAAwHkEEAACMRxABAADjEUSQJK1du1b9+vVTRESERo0apUOHDl1xbVZWlhwOh88lIiKiDacNPu+//77uvfde9erVSw6HQ++888733mf//v0aPny4wsPD1b9/f2VlZbX6nMGupcdx//79jZ6LDodDpaWlbTNwkMrMzNSPfvQjRUZGKiYmRpMnT9apU6e+935bt27VTTfdpIiICA0ePFh/+MMf2mDa4OXPceT3Y2OvvPKKhgwZ4v3SRbfbrZ07d171PnY8Fwki6O2331ZGRoaWLl2qI0eOaOjQoUpOTlZ5efkV7+NyuVRSUuK9fPXVV204cfCpqqrS0KFDtXbt2matP336tFJSUjRmzBgVFhYqPT1djz76qHbv3t3Kkwa3lh7HBqdOnfJ5PsbExLTShO1Dbm6u0tLSlJ+fr+zsbNXW1mrChAmqqqq64n0OHDign/zkJ5o9e7aOHj2qyZMna/LkyTpx4kQbTh5c/DmOEr8fL9e7d289++yzKigo0OHDhzV27FhNmjRJJ0+ebHK9bc9FC8YbOXKklZaW5r1eV1dn9erVy8rMzGxy/YYNG6yoqKg2mq79kWRt27btqmsWLVpk3XzzzT7bHnjgASs5ObkVJ2tfmnMc9+3bZ0myvv766zaZqb0qLy+3JFm5ublXXPOP//iPVkpKis+2UaNGWT/72c9ae7x2oznHkd+PzdO1a1frtddea/I2u56LvEJkuJqaGhUUFGj8+PHebSEhIRo/frzy8vKueL8LFy6ob9++io+Pv2rpo2l5eXk+x1ySkpOTr3rMcWXDhg1Tz549ddddd+nDDz+0e5ygU1lZKUnq1q3bFdfwnPx+zTmOEr8fr6aurk6bN29WVVXVFf90ll3PRYLIcH/9619VV1fX6Bu7Y2Njr3gexoABA7R+/Xr97ne/05tvvqn6+nrdeuut+stf/tIWI18TSktLmzzmHo9H33zzjU1TtT89e/bUunXr9Jvf/Ea/+c1vFB8fr9GjR+vIkSN2jxY06uvrlZ6erttuu02DBg264rorPSdNPx+rQXOPI78fm3b8+HF16dJF4eHhmjNnjrZt26bExMQm19r1XDTmT3cgcNxut0/Z33rrrRo4cKBeffVVLV++3MbJYJoBAwZowIAB3uu33nqrvvjiC61evVr/+Z//aeNkwSMtLU0nTpzQBx98YPco7VpzjyO/H5s2YMAAFRYWqrKyUv/93/+t1NRU5ebmXjGK7MArRIbr0aOHOnTooLKyMp/tZWVliouLa9Y+wsLC9Pd///f685//3BojXpPi4uKaPOYul0sdO3a0aaprw8iRI3ku/n/z5s3T9u3btW/fPvXu3fuqa6/0nGzu74FrWUuO4+X4/fgtp9Op/v37a8SIEcrMzNTQoUO1Zs2aJtfa9VwkiAzndDo1YsQI5eTkeLfV19crJyfniu/vXq6urk7Hjx9Xz549W2vMa47b7fY55pKUnZ3d7GOOKyssLDT+uWhZlubNm6dt27Zp7969SkhI+N778JxszJ/jeDl+Pzatvr5e1dXVTd5m23OxVU/ZRruwefNmKzw83MrKyrI++eQT67HHHrOio6Ot0tJSy7Isa8aMGdZTTz3lXf/MM89Yu3fvtr744guroKDAmjZtmhUREWGdPHnSrodgu/Pnz1tHjx61jh49akmyVq1aZR09etT66quvLMuyrKeeesqaMWOGd/2XX35pderUyVq4cKH16aefWmvXrrU6dOhg7dq1y66HEBRaehxXr15tvfPOO9bnn39uHT9+3Jo/f74VEhJivffee3Y9hKAwd+5cKyoqytq/f79VUlLivfztb3/zrrn8/+sPP/zQCg0NtX71q19Zn376qbV06VIrLCzMOn78uB0PISj4cxz5/djYU089ZeXm5lqnT5+2jh07Zj311FOWw+Gw9uzZY1lW8DwXCSJYlmVZv/71r60+ffpYTqfTGjlypJWfn++97cc//rGVmprqvZ6enu5dGxsba02cONE6cuSIDVMHj4aPf19+aThuqamp1o9//ONG9xk2bJjldDqtv/u7v7M2bNjQ5nMHm5Yex+eee866/vrrrYiICKtbt27W6NGjrb1799ozfBBp6hhK8nmOXf7/tWVZ1pYtW6wbb7zRcjqd1s0332zt2LGjbQcPMv4cR34/NjZr1iyrb9++ltPptK677jpr3Lhx3hiyrOB5Ljosy7Ja9zUoAACA4MY5RAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOP9P1CetU2+SpbcAAAAAElFTkSuQmCC",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"persdur_zoom.plot.hist(bins=10, edgecolor=\"white\", color=\"darkgrey\");"
]
},
{
"cell_type": "markdown",
"id": "cafd2d6f",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Questions\n",
"\n",
"- Describe how the empirical histogram changes as the the value `n` gets larger."
]
},
{
"cell_type": "markdown",
"id": "fef3a8a1",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Empirical Distribution of a Statistic\n",
"\n",
"The Law of Averages implies that with high probability, the empirical distribution of a large random sample will resemble the distribution of the population from which the sample was drawn.\n",
"\n",
"The resemblance is visible in two histograms: the empirical histogram of a large random sample is likely to resemble the histogram of the population."
]
},
{
"cell_type": "markdown",
"id": "944a7331",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Statistical Parameter\n",
"\n",
"Numerical quantities associated with a (statistical) population are called *statistical parameters* or *parameters*. For the population of respondents in `persdur_zoom`, we know the value of the parameter \"median time (hours) spent on personal activities\":"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "637eee6b",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.75"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"persdur_zoom.median()"
]
},
{
"cell_type": "markdown",
"id": "2010df46",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"More frequently, we are interested in parameters that are unknown to us.\n",
"\n",
"- In a population of voters, what percent will vote for Candidate A?\n",
"\n",
"- In a population of TikTok users, what is the largest number of followers for a user?\n",
"\n",
"- In a population of Air Canada flights, what is the median departure delay?\n",
"\n",
"- In a population of commercial air flights, what is the average [fatal accident rate](https://www.cnn.com/2025/02/19/business/airplane-crashes-statistics/index.html)?"
]
},
{
"cell_type": "markdown",
"id": "3fad42c7",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Statistic\n",
"\n",
"For this, we will rely on data from a large random sample drawn from the population.\n",
"\n",
"A *statistic* (note the singular!) is any number computed using the data in a sample. The sample median, therefore, is a statistic. \n",
"\n",
"Remember that `persdur_zoom.sample(100)` contains a random sample of 100 respondents from `persdur_zoom`. The observed value of the sample median is:"
]
},
{
"cell_type": "markdown",
"id": "8170f8fc",
"metadata": {},
"source": [
""
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "7ad59370",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.75"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"persdur_zoom.sample(100).median()"
]
},
{
"cell_type": "markdown",
"id": "f836e0e0",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Our sample – one set of 100 people – gave us one observed value of the statistic. This raises an important problem of inference:\n",
"\n",
"**The statistic could have been different.**\n",
"\n",
"A fundamental consideration in using any statistic based on a random sample is that *the sample could have come out differently*, and therefore the statistic could have come out differently too."
]
},
{
"cell_type": "markdown",
"id": "2d8fe60f",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Run the cell above a few times to see how the answer varies. Often it is equal to 0.75, the same value as the population parameter. But sometimes it is different."
]
},
{
"cell_type": "markdown",
"id": "c1c15f98",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Example: Poll Tracker\n",
"\n",
"- Consider [the polling averages of the Canadian parties](https://newsinteractives.cbc.ca/elections/poll-tracker/canada/) as of today.\n",
"\n",
"- Each poll consists of responses from a random sample of survey participants.\n",
"\n",
"- Polls even conducted during the same period may not result in the same results."
]
},
{
"cell_type": "markdown",
"id": "29a7b193",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"\n",
"\n",
"_Source: https://poll-mbc.co.kr/bk/2022_president.html_"
]
},
{
"cell_type": "markdown",
"id": "917e4166",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Uh oh...\n",
"\n",
"- We don't know the population parameter.\n",
"\n",
"- Our statistic from a random sample may be a good estimate, but it may be not because it's going to be different with a different sample."
]
},
{
"cell_type": "markdown",
"id": "272a6b93",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"\n",
"**Just how different could the statistic have been?** \n",
"\n",
"One way to answer this is to simulate the statistic many times and note the values. \n",
"A histogram of those values will tell us about the **distribution of the statistic**."
]
},
{
"cell_type": "markdown",
"id": "9f6da1dd",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Simulating a Statistic\n",
"\n",
"We will simulate the sample median using the steps below. You can replace the sample size of 1000 by any other sample size, and the sample median by any other statistic.\n",
"\n",
"**Step 1: Decide which statistic to simulate.** We have already decided that: we are going to simulate the median of a random sample of size 1000 drawn from the population of time use survey respondents that had a median time spent on personal activities between 0.25 and 3.0 hours.\n",
"\n",
"**Step 2: Define a function that returns one simulated value of the statistic.** Draw a random sample of size 1000 and compute the median of the sample. We did this in the code cell above. Here it is again, encapsulated in a function."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "5eb378fb",
"metadata": {},
"outputs": [],
"source": [
"def random_sample_median():\n",
" return persdur_zoom.sample(100).median()"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "1eec73a6",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"0.7083333333333333"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"random_sample_median()"
]
},
{
"cell_type": "markdown",
"id": "3979a0fb",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"**Step 3: Decide how many simulated values to generate.** Let's do 5,000 repetitions.\n",
"\n",
"**Step 4: Use a `for` loop to generate a list of simulated values.** Start by creating an empty list in which to collect our results. We will then set up a `for` loop for generating all the simulated values. The body of the loop will consist of generating one simulated value of the sample median, and appending it to our collection list.\n",
"\n",
"The simulation takes a noticeable amount of time to run. That is because it is performing 5000 repetitions of the process of drawing a sample of size 1000 and computing its median. That's a lot of sampling and repeating!"
]
},
{
"cell_type": "markdown",
"id": "1c6a6530",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Let's break down this step a bit further:\n",
"\n",
"- set up an empty list called `sim_medians`."
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "f20eea7b",
"metadata": {},
"outputs": [],
"source": [
"sim_medians = [] # empty list"
]
},
{
"cell_type": "markdown",
"id": "306bbb0d",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- Use the `append` function to append values to `sim_medians`."
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "e5b24aeb",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[0.9166666666666666]"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sim_medians.append(random_sample_median())\n",
"sim_medians"
]
},
{
"cell_type": "markdown",
"id": "82fd17f1",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- Run the cell above several times and you will see that values keep getting appended to `sim_medians`.\n",
"\n",
"- Each time the cell is run a random sample of 100 is drawn, the median is calculated then appended to the list."
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "a93270f7",
"metadata": {},
"outputs": [],
"source": [
"sim_medians = []\n",
"\n",
"num_sims = 5000 \n",
"\n",
"for _ in range(num_sims):\n",
" sim_medians.append(random_sample_median())"
]
},
{
"cell_type": "markdown",
"id": "fcc87430",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Empirical Distribution of a Statistic\n",
"\n",
"We can now examine the **empirical frequency distribution** of the median statistic."
]
},
{
"cell_type": "markdown",
"id": "0d479dcb",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Create a `pandas` series using the list `sim_medians` then use the `describe` function to describe the distribution."
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "6524a594",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"count 5000.000000\n",
"mean 0.800458\n",
"std 0.081846\n",
"min 0.500000\n",
"25% 0.750000\n",
"50% 0.750000\n",
"75% 0.833333\n",
"max 1.000000\n",
"dtype: float64"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pd.Series(sim_medians).describe()"
]
},
{
"cell_type": "markdown",
"id": "03518bb4",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"- The distribution of the median can be visualized using a histogram. \n",
"\n",
"- The histogram can tell us how frequent certain values of the median occur in random samples."
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "5127a82f",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGdCAYAAAAMm0nCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAkoUlEQVR4nO3df3DU1b3/8dcmZDdgs0kRkk2ukQas/FAEjTVsRUTNJEDq1ZE7lYJALcLYm/ReSC8g1SKiNTO0otai1KkY74iCXhEt0JAYiqmaSBtNhahcgfQGL2ywYLJAJYHkfP+4Xz5lS9BszGb3xOdjZmfcz+fs5r1n4uxz9kdwGWOMAAAALBIX7QEAAADCRcAAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsE6/aA8QKR0dHTpw4ICSkpLkcrmiPQ4AAOgCY4yOHj2qjIwMxcWd+3WWPhswBw4cUGZmZrTHAAAA3bB//35dcMEF5zzfZwMmKSlJ0v9tgNfrjfI0AACgK4LBoDIzM53n8XPpswFz+m0jr9dLwAAAYJkv+vgHH+IFAADWIWAAAIB1CBgAAGAdAgYAAFiHgAEAANYhYAAAgHUIGAAAYB0CBgAAWIeAAQAA1iFgAACAdQgYAABgHQIGAABYh4ABAADWIWAAAIB1CBgghhhjoj2CNdgr4KutX7QHAPB3LpdLNTU1CgaD0R4lpnm9Xo0bNy7aYwCIIgIGiDHBYFDNzc3RHgMAYhpvIQEAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6YQVMSUmJvvWtbykpKUmpqam6+eabtXv37pA1EydOlMvlCrnceeedIWsaGxtVUFCgAQMGKDU1VQsXLtSpU6dC1mzfvl1XXHGFPB6PLrroIpWWlnbvEQIAgD4nrIB5/fXXVVhYqJqaGlVUVOjkyZPKy8vT8ePHQ9bNnTtXBw8edC4rVqxwzrW3t6ugoEBtbW1666239Mwzz6i0tFRLly511jQ0NKigoEDXXXed6urqNH/+fN1xxx3aunXrl3y4AACgL+gXzuKysrKQ66WlpUpNTVVtba0mTJjgHB8wYIB8Pl+n91FeXq73339fr732mtLS0jR27Fjdf//9Wrx4sZYtWya3263Vq1crKytLDz30kCRp5MiReuONN/Twww8rPz8/3McIAAD6mC/1GZiWlhZJ0sCBA0OOr127VoMGDdKll16qJUuW6G9/+5tzrrq6WqNHj1ZaWppzLD8/X8FgUPX19c6a3NzckPvMz89XdXX1lxkXAAD0EWG9AnOmjo4OzZ8/X1dffbUuvfRS5/j06dM1ZMgQZWRk6L333tPixYu1e/dubdiwQZIUCARC4kWScz0QCHzummAwqM8++0z9+/c/a57W1la1trY614PBYHcfGgAAiHHdDpjCwkLt2rVLb7zxRsjxefPmOf89evRopaen64YbbtDevXs1bNiw7k/6BUpKSnTfffdF7P4BAEDs6NZbSEVFRdq0aZN+//vf64ILLvjctTk5OZKkPXv2SJJ8Pp+amppC1py+fvpzM+da4/V6O331RZKWLFmilpYW57J///7wHxgAALBCWAFjjFFRUZFefvllbdu2TVlZWV94m7q6OklSenq6JMnv92vnzp06dOiQs6aiokJer1ejRo1y1lRWVobcT0VFhfx+/zl/jsfjkdfrDbkAAIC+KayAKSws1LPPPqvnnntOSUlJCgQCCgQC+uyzzyRJe/fu1f3336/a2lr95S9/0auvvqpZs2ZpwoQJuuyyyyRJeXl5GjVqlGbOnKk///nP2rp1q+655x4VFhbK4/FIku68807t27dPixYt0ocffqjHH39cL7zwghYsWNDDDx8AANgorIB54okn1NLSookTJyo9Pd25rF+/XpLkdrv12muvKS8vTyNGjNCPf/xjTZ06Vb/97W+d+4iPj9emTZsUHx8vv9+v2267TbNmzdLy5cudNVlZWdq8ebMqKio0ZswYPfTQQ/rNb37DV6gBAIAkyWWMMdEeIhKCwaCSk5PV0tLC20mwSnl5uZqbm6M9RkxLSUlRXl5etMcAEAFdff7m30ICAADWIWAAAIB1CBgAAGAdAgYAAFiHgAEAANYhYAAAgHUIGAAAYB0CBgAAWIeAAQAA1iFgAACAdQgYAABgHQIGAABYh4ABAADWIWAAAIB1CBgAAGAdAgYAAFiHgAEAANYhYAAAgHUIGAAAYB0CBgAAWIeAAQAA1iFgAACAdQgYAABgHQIGAABYh4ABAADWIWAAAIB1CBgAAGAdAgYAAFiHgAEAANYhYAAAgHUIGAAAYB0CBgAAWIeAAQAA1iFgAACAdQgYAABgHQIGAABYh4ABAADWIWAAAIB1CBgAAGAdAgYAAFiHgAEAANYhYAAAgHUIGAAAYB0CBgAAWIeAAQAA1iFgAACAdQgYAABgHQIGAABYh4ABAADWIWAAAIB1CBgAAGAdAgYAAFiHgAEAANYhYAAAgHXCCpiSkhJ961vfUlJSklJTU3XzzTdr9+7dIWtOnDihwsJCnX/++fra176mqVOnqqmpKWRNY2OjCgoKNGDAAKWmpmrhwoU6depUyJrt27friiuukMfj0UUXXaTS0tLuPUIAANDnhBUwr7/+ugoLC1VTU6OKigqdPHlSeXl5On78uLNmwYIF+u1vf6sXX3xRr7/+ug4cOKBbbrnFOd/e3q6CggK1tbXprbfe0jPPPKPS0lItXbrUWdPQ0KCCggJdd911qqur0/z583XHHXdo69atPfCQAQCA7VzGGNPdG3/yySdKTU3V66+/rgkTJqilpUWDBw/Wc889p3/5l3+RJH344YcaOXKkqqurNW7cOP3ud7/Td77zHR04cEBpaWmSpNWrV2vx4sX65JNP5Ha7tXjxYm3evFm7du1yfta0adPU3NyssrKyLs0WDAaVnJyslpYWeb3e7j5EoNeVl5erubk52mPEtJSUFOXl5UV7DAAR0NXn7y/1GZiWlhZJ0sCBAyVJtbW1OnnypHJzc501I0aM0IUXXqjq6mpJUnV1tUaPHu3EiyTl5+crGAyqvr7eWXPmfZxec/o+OtPa2qpgMBhyAQAAfVO3A6ajo0Pz58/X1VdfrUsvvVSSFAgE5Ha7lZKSErI2LS1NgUDAWXNmvJw+f/rc560JBoP67LPPOp2npKREycnJziUzM7O7Dw0AAMS4bgdMYWGhdu3apXXr1vXkPN22ZMkStbS0OJf9+/dHeyQAABAh/bpzo6KiIm3atElVVVW64IILnOM+n09tbW1qbm4OeRWmqalJPp/PWbNjx46Q+zv9LaUz1/zjN5eamprk9XrVv3//TmfyeDzyeDzdeTgAAMAyYb0CY4xRUVGRXn75ZW3btk1ZWVkh57Ozs5WQkKDKykrn2O7du9XY2Ci/3y9J8vv92rlzpw4dOuSsqaiokNfr1ahRo5w1Z97H6TWn7wMAAHy1hfUKTGFhoZ577jm98sorSkpKcj6zkpycrP79+ys5OVlz5sxRcXGxBg4cKK/Xqx/96Efy+/0aN26cJCkvL0+jRo3SzJkztWLFCgUCAd1zzz0qLCx0XkG588479atf/UqLFi3SD37wA23btk0vvPCCNm/e3MMPHwAA2CisV2CeeOIJtbS0aOLEiUpPT3cu69evd9Y8/PDD+s53vqOpU6dqwoQJ8vl82rBhg3M+Pj5emzZtUnx8vPx+v2677TbNmjVLy5cvd9ZkZWVp8+bNqqio0JgxY/TQQw/pN7/5jfLz83vgIQMAANt9qb8DE8v4OzCwFX8H5ovxd2CAvqtX/g4MAABANBAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrhB0wVVVVuvHGG5WRkSGXy6WNGzeGnP/+978vl8sVcpk0aVLImiNHjmjGjBnyer1KSUnRnDlzdOzYsZA17733nq655holJiYqMzNTK1asCP/RAQCAPinsgDl+/LjGjBmjVatWnXPNpEmTdPDgQefy/PPPh5yfMWOG6uvrVVFRoU2bNqmqqkrz5s1zzgeDQeXl5WnIkCGqra3Vz3/+cy1btkxPPvlkuOMCAIA+qF+4N5g8ebImT578uWs8Ho98Pl+n5z744AOVlZXpj3/8o6688kpJ0mOPPaYpU6boF7/4hTIyMrR27Vq1tbVpzZo1crvduuSSS1RXV6eVK1eGhA4AAPhqishnYLZv367U1FQNHz5cP/zhD3X48GHnXHV1tVJSUpx4kaTc3FzFxcXp7bffdtZMmDBBbrfbWZOfn6/du3fr008/jcTIAADAImG/AvNFJk2apFtuuUVZWVnau3evfvKTn2jy5Mmqrq5WfHy8AoGAUlNTQ4fo108DBw5UIBCQJAUCAWVlZYWsSUtLc859/etfP+vntra2qrW11bkeDAZ7+qEBAIAY0eMBM23aNOe/R48ercsuu0zDhg3T9u3bdcMNN/T0j3OUlJTovvvui9j9AwCA2BHxr1EPHTpUgwYN0p49eyRJPp9Phw4dCllz6tQpHTlyxPncjM/nU1NTU8ia09fP9dmaJUuWqKWlxbns37+/px8KAACIEREPmI8//liHDx9Wenq6JMnv96u5uVm1tbXOmm3btqmjo0M5OTnOmqqqKp08edJZU1FRoeHDh3f69pH0fx8c9nq9IRcAANA3hR0wx44dU11dnerq6iRJDQ0NqqurU2Njo44dO6aFCxeqpqZGf/nLX1RZWambbrpJF110kfLz8yVJI0eO1KRJkzR37lzt2LFDb775poqKijRt2jRlZGRIkqZPny632605c+aovr5e69ev16OPPqri4uKee+QAAMBaYQfMn/70J11++eW6/PLLJUnFxcW6/PLLtXTpUsXHx+u9997TP//zP+viiy/WnDlzlJ2drT/84Q/yeDzOfaxdu1YjRozQDTfcoClTpmj8+PEhf+MlOTlZ5eXlamhoUHZ2tn784x9r6dKlfIUaAABIklzGGBPtISIhGAwqOTlZLS0tvJ0Eq5SXl6u5uTnaY8S0lJQU5eXlRXsMABHQ1edv/i0kAABgHQIGAABYh4ABAADWIWAAAIB1CBgAAGAdAgYAAFiHgAEAANYhYAAAgHUIGAAAYB0CBgAAWIeAAQAA1iFgAACAdQgYAABgHQIGAABYh4ABAADWIWAAAIB1CBgAAGAdAgYAAFiHgAEAANYhYAAAgHUIGAAAYB0CBgAAWIeAAQAA1iFgAFgnMTFRxphoj2EN9gp9Ub9oDwAA4UpISJDL5VJNTY2CwWC0x4lpXq9X48aNi/YYQI8jYABYKxgMqrm5OdpjAIgC3kICAADWIWAAAIB1CBgAAGAdAgYAAFiHgAEAANYhYAAAgHUIGAAAYB0CBgAAWIeAAQAA1iFgAACAdQgYAABgHQIGAABYh4ABAADWIWAAAIB1CBgAAGAdAgYAAFiHgAEAANYhYAAAgHUIGAAAYB0CBgAAWIeAAQAA1iFgAACAdQgYAABgHQIGAABYh4ABAADWIWAAAIB1CBgAAGAdAgYAAFgn7ICpqqrSjTfeqIyMDLlcLm3cuDHkvDFGS5cuVXp6uvr376/c3Fx99NFHIWuOHDmiGTNmyOv1KiUlRXPmzNGxY8dC1rz33nu65pprlJiYqMzMTK1YsSL8RwcAAPqksAPm+PHjGjNmjFatWtXp+RUrVuiXv/ylVq9erbffflvnnXee8vPzdeLECWfNjBkzVF9fr4qKCm3atElVVVWaN2+ecz4YDCovL09DhgxRbW2tfv7zn2vZsmV68sknu/EQAQBAX9Mv3BtMnjxZkydP7vScMUaPPPKI7rnnHt10002SpP/8z/9UWlqaNm7cqGnTpumDDz5QWVmZ/vjHP+rKK6+UJD322GOaMmWKfvGLXygjI0Nr165VW1ub1qxZI7fbrUsuuUR1dXVauXJlSOgAAICvph79DExDQ4MCgYByc3OdY8nJycrJyVF1dbUkqbq6WikpKU68SFJubq7i4uL09ttvO2smTJggt9vtrMnPz9fu3bv16aefdvqzW1tbFQwGQy4AAKBv6tGACQQCkqS0tLSQ42lpac65QCCg1NTUkPP9+vXTwIEDQ9Z0dh9n/ox/VFJSouTkZOeSmZn55R8QAACISX3mW0hLlixRS0uLc9m/f3+0RwIAABHSowHj8/kkSU1NTSHHm5qanHM+n0+HDh0KOX/q1CkdOXIkZE1n93Hmz/hHHo9HXq835AIAAPqmHg2YrKws+Xw+VVZWOseCwaDefvtt+f1+SZLf71dzc7Nqa2udNdu2bVNHR4dycnKcNVVVVTp58qSzpqKiQsOHD9fXv/71nhwZAABYKOyAOXbsmOrq6lRXVyfp/z64W1dXp8bGRrlcLs2fP18PPPCAXn31Ve3cuVOzZs1SRkaGbr75ZknSyJEjNWnSJM2dO1c7duzQm2++qaKiIk2bNk0ZGRmSpOnTp8vtdmvOnDmqr6/X+vXr9eijj6q4uLjHHjgAALBX2F+j/tOf/qTrrrvOuX46KmbPnq3S0lItWrRIx48f17x589Tc3Kzx48errKxMiYmJzm3Wrl2roqIi3XDDDYqLi9PUqVP1y1/+0jmfnJys8vJyFRYWKjs7W4MGDdLSpUv5CjUAAJAkuYwxJtpDREIwGFRycrJaWlr4PAysUl5erubm5miPEdMyMzPl9/vZqy5ISUlRXl5etMcAuqyrz9995ltIAADgq4OAAQAA1iFgAACAdQgYAABgHQIGAABYh4ABAADWIWAAAIB1CBgAAGAdAgYAAFiHgAEAANYhYAAAgHUIGAAAYB0CBgAAWIeAAQAA1iFgAACAdQgYAABgHQIGAABYh4ABAADWIWAAAIB1CBgAAGAdAgYAAFiHgAEAANYhYAAAgHUIGAAAYB0CBgAAWIeAAQAA1iFgAACAdQgYAABgHQIGAABYh4ABAADWIWAAAIB1CBgAAGAdAgYAAFiHgAEAANYhYAAAgHUIGAAAYB0CBgAAWIeAAQAA1iFgAACAdQgYAABgHQIGAABYh4ABAADWIWAAAIB1CBgAAGAdAgYAAFiHgAEAANYhYAAAgHUIGAAAYB0CBgAAWIeAAQAA1iFgAACAdQgYAABgHQIGAABYh4ABAADW6fGAWbZsmVwuV8hlxIgRzvkTJ06osLBQ559/vr72ta9p6tSpampqCrmPxsZGFRQUaMCAAUpNTdXChQt16tSpnh4VAABYql8k7vSSSy7Ra6+99vcf0u/vP2bBggXavHmzXnzxRSUnJ6uoqEi33HKL3nzzTUlSe3u7CgoK5PP59NZbb+ngwYOaNWuWEhIS9OCDD0ZiXAAAYJmIBEy/fv3k8/nOOt7S0qKnnnpKzz33nK6//npJ0tNPP62RI0eqpqZG48aNU3l5ud5//3299tprSktL09ixY3X//fdr8eLFWrZsmdxudyRGBgAAFonIZ2A++ugjZWRkaOjQoZoxY4YaGxslSbW1tTp58qRyc3OdtSNGjNCFF16o6upqSVJ1dbVGjx6ttLQ0Z01+fr6CwaDq6+vP+TNbW1sVDAZDLgAAoG/q8YDJyclRaWmpysrK9MQTT6ihoUHXXHONjh49qkAgILfbrZSUlJDbpKWlKRAISJICgUBIvJw+f/rcuZSUlCg5Odm5ZGZm9uwDAwAAMaPH30KaPHmy89+XXXaZcnJyNGTIEL3wwgvq379/T/84x5IlS1RcXOxcDwaDRAwAAH1UxL9GnZKSoosvvlh79uyRz+dTW1ubmpubQ9Y0NTU5n5nx+XxnfSvp9PXOPldzmsfjkdfrDbkAAIC+KeIBc+zYMe3du1fp6enKzs5WQkKCKisrnfO7d+9WY2Oj/H6/JMnv92vnzp06dOiQs6aiokJer1ejRo2K9LgAAMACPf4W0n/8x3/oxhtv1JAhQ3TgwAHde++9io+P1/e+9z0lJydrzpw5Ki4u1sCBA+X1evWjH/1Ifr9f48aNkyTl5eVp1KhRmjlzplasWKFAIKB77rlHhYWF8ng8PT0uAACwUI8HzMcff6zvfe97Onz4sAYPHqzx48erpqZGgwcPliQ9/PDDiouL09SpU9Xa2qr8/Hw9/vjjzu3j4+O1adMm/fCHP5Tf79d5552n2bNna/ny5T09KgAAsFSPB8y6des+93xiYqJWrVqlVatWnXPNkCFDtGXLlp4eDQAA9BH8W0gAAMA6BAwA9GGJiYkyxkR7DGuwV/aIyD8lAACIDQkJCXK5XKqpqeEvlH8Br9frfKEEsY+AAYCvgGAweNbf4AJsxltIAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAEhKTEyUMSbaY1gj2nvVL6o/HQCAGJGQkCCXy6WamhoFg8FojxPTvF6vxo0bF9UZCBgAAM4QDAbV3Nwc7THwBXgLCREX7ZcZAQB9D6/AIOJ4SbZrfD6fLrvssmiPAQBWIGDQK3hJ9oslJSVFewQAsAZvIQEAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsA4BAwAArEPAAAAA6xAwAADAOgQMAACwDgEDAACsQ8AAAADrEDAAAMA6BAwAALAOAQMAAKxDwAAAAOsQMAAAwDoEDAAAsE5MB8yqVav0jW98Q4mJicrJydGOHTuiPRIAAIgBMRsw69evV3Fxse6991698847GjNmjPLz83Xo0KFojwYAAKIsZgNm5cqVmjt3rm6//XaNGjVKq1ev1oABA7RmzZpojwYAAKKsX7QH6ExbW5tqa2u1ZMkS51hcXJxyc3NVXV3d6W1aW1vV2trqXG9paZEkBYPBHp/PGCOXy9Xj99uXxcfHKyEhIdpjxDRjjILBIHvVBexV17FXXcdedV18fHxEnl+lvz9vG2M+d11MBsxf//pXtbe3Ky0tLeR4WlqaPvzww05vU1JSovvuu++s45mZmRGZEQAARM7Ro0eVnJx8zvMxGTDdsWTJEhUXFzvXOzo6dOTIEZ1//vk9+mpJMBhUZmam9u/fL6/X22P321exX13HXnUde9V17FXXsVddF8m9Msbo6NGjysjI+Nx1MRkwgwYNUnx8vJqamkKONzU1yefzdXobj8cjj8cTciwlJSVSI8rr9fILHgb2q+vYq65jr7qOveo69qrrIrVXn/fKy2kx+SFet9ut7OxsVVZWOsc6OjpUWVkpv98fxckAAEAsiMlXYCSpuLhYs2fP1pVXXqmrrrpKjzzyiI4fP67bb7892qMBAIAoi9mAufXWW/XJJ59o6dKlCgQCGjt2rMrKys76YG9v83g8uvfee896uwqdY7+6jr3qOvaq69irrmOvui4W9splvuh7SgAAADEmJj8DAwAA8HkIGAAAYB0CBgAAWIeAAQAA1iFgOrFq1Sp94xvfUGJionJycrRjx45zri0tLZXL5Qq5JCYm9uK00RXOXklSc3OzCgsLlZ6eLo/Ho4svvlhbtmzppWmjL5z9mjhx4lm/Wy6XSwUFBb04cfSE+7v1yCOPaPjw4erfv78yMzO1YMECnThxopemja5w9urkyZNavny5hg0bpsTERI0ZM0ZlZWW9OG30VFVV6cYbb1RGRoZcLpc2btz4hbfZvn27rrjiCnk8Hl100UUqLS2N+JyxINy9OnjwoKZPn66LL75YcXFxmj9/fuSHNAixbt0643a7zZo1a0x9fb2ZO3euSUlJMU1NTZ2uf/rpp43X6zUHDx50LoFAoJenjo5w96q1tdVceeWVZsqUKeaNN94wDQ0NZvv27aaurq6XJ4+OcPfr8OHDIb9Xu3btMvHx8ebpp5/u3cGjINy9Wrt2rfF4PGbt2rWmoaHBbN261aSnp5sFCxb08uS9L9y9WrRokcnIyDCbN282e/fuNY8//rhJTEw077zzTi9P3vu2bNli7r77brNhwwYjybz88sufu37fvn1mwIABpri42Lz//vvmscceM/Hx8aasrKx3Bo6icPeqoaHB/Nu//Zt55plnzNixY82///u/R3xGAuYfXHXVVaawsNC53t7ebjIyMkxJSUmn659++mmTnJzcS9PFlnD36oknnjBDhw41bW1tvTViTAl3v/7Rww8/bJKSksyxY8ciNWLMCHevCgsLzfXXXx9yrLi42Fx99dURnTMWhLtX6enp5le/+lXIsVtuucXMmDEjonPGmq48KS9atMhccsklIcduvfVWk5+fH8HJYk9X9upM1157ba8EDG8hnaGtrU21tbXKzc11jsXFxSk3N1fV1dXnvN2xY8c0ZMgQZWZm6qabblJ9fX1vjBtV3dmrV199VX6/X4WFhUpLS9Oll16qBx98UO3t7b01dtR093frTE899ZSmTZum8847L1JjxoTu7NW3v/1t1dbWOm+d7Nu3T1u2bNGUKVN6ZeZo6c5etba2nvU2d//+/fXGG29EdFYbVVdXh+ytJOXn53f5/1lEFgFzhr/+9a9qb28/66/9pqWlKRAIdHqb4cOHa82aNXrllVf07LPPqqOjQ9/+9rf18ccf98bIUdOdvdq3b5/+67/+S+3t7dqyZYt++tOf6qGHHtIDDzzQGyNHVXf260w7duzQrl27dMcdd0RqxJjRnb2aPn26li9frvHjxyshIUHDhg3TxIkT9ZOf/KQ3Ro6a7uxVfn6+Vq5cqY8++kgdHR2qqKjQhg0bdPDgwd4Y2SqBQKDTvQ0Gg/rss8+iNBVOI2C+JL/fr1mzZmns2LG69tprtWHDBg0ePFi//vWvoz1azOno6FBqaqqefPJJZWdn69Zbb9Xdd9+t1atXR3u0mPfUU09p9OjRuuqqq6I9Skzavn27HnzwQT3++ON65513tGHDBm3evFn3339/tEeLOY8++qi++c1vasSIEXK73SoqKtLtt9+uuDieDmCXmP23kKJh0KBBio+PV1NTU8jxpqYm+Xy+Lt1HQkKCLr/8cu3ZsycSI8aM7uxVenq6EhISFB8f7xwbOXKkAoGA2tra5Ha7IzpzNH2Z363jx49r3bp1Wr58eSRHjBnd2auf/vSnmjlzpvMK1ejRo3X8+HHNmzdPd999d599cu7OXg0ePFgbN27UiRMndPjwYWVkZOiuu+7S0KFDe2Nkq/h8vk731uv1qn///lGaCqf1zf+ru8ntdis7O1uVlZXOsY6ODlVWVsrv93fpPtrb27Vz506lp6dHasyY0J29uvrqq7Vnzx51dHQ4x/77v/9b6enpfTpepC/3u/Xiiy+qtbVVt912W6THjAnd2au//e1vZ0XK6VA2ffife/syv1eJiYn6p3/6J506dUovvfSSbrrppkiPax2/3x+yt5JUUVHR5ecDRFjEPyZsmXXr1hmPx2NKS0vN+++/b+bNm2dSUlKcr0bPnDnT3HXXXc76++67z2zdutXs3bvX1NbWmmnTppnExERTX18frYfQa8Ldq8bGRpOUlGSKiorM7t27zaZNm0xqaqp54IEHovUQelW4+3Xa+PHjza233trb40ZVuHt17733mqSkJPP888+bffv2mfLycjNs2DDz3e9+N1oPodeEu1c1NTXmpZdeMnv37jVVVVXm+uuvN1lZWebTTz+N0iPoPUePHjXvvvuueffdd40ks3LlSvPuu++a//mf/zHGGHPXXXeZmTNnOutPf4164cKF5oMPPjCrVq36ynyNOty9MsY467Ozs8306dPNu+++G9HnQgKmE4899pi58MILjdvtNldddZWpqalxzl177bVm9uzZzvX58+c7a9PS0syUKVO+En9P4bRw9soYY9566y2Tk5NjPB6PGTp0qPnZz35mTp061ctTR0+4+/Xhhx8aSaa8vLyXJ42+cPbq5MmTZtmyZWbYsGEmMTHRZGZmmn/913/9SjwpGxPeXm3fvt2MHDnSeDwec/7555uZM2ea//3f/43C1L3v97//vZF01uX0/syePdtce+21Z91m7Nixxu12m6FDh34l/g6TMd3bq87WDxkyJGIzuv7/DwUAALAGn4EBAADWIWAAAIB1CBgAAGAdAgYAAFiHgAEAANYhYAAAgHUIGAAAYB0CBgAAWIeAAQAA1iFgAACAdQgYAABgHQIGAABY5/8Bercd6HLU97gAAAAASUVORK5CYII=",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"# create bins of length 0.1 \n",
"# starting at 0.5 and ending at 1.1\n",
"\n",
"median_bins = np.arange(start = 0.5, stop = 1.1, step = 0.1) # use a numpy function here\n",
"\n",
"# use the bins to plot empirical distribution of medians\n",
"\n",
"plt.hist(sim_medians, bins=median_bins, edgecolor=\"white\", color=\"darkgrey\");"
]
},
{
"cell_type": "markdown",
"id": "96d02581",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"The exact counts of medians in each interval can be examined using `pd.cut` and `value_counts`."
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "0cfab888",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.7, 0.8] 2727\n",
"(0.8, 0.9] 1275\n",
"(0.9, 1.0] 437\n",
"(1.0, 1.1] 285\n",
"(0.6, 0.7] 269\n",
"(0.5, 0.6] 5\n",
"Name: count, dtype: int64"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# create a pandas series so we can use cut function\n",
"sim_means_series = pd.Series(sim_medians)\n",
"\n",
"# frequency of values in each bin\n",
"pd.cut(sim_means_series, median_bins).value_counts()"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "94226cd1",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.7, 0.8] 0.5454\n",
"(0.8, 0.9] 0.2550\n",
"(0.9, 1.0] 0.0874\n",
"(1.0, 1.1] 0.0570\n",
"(0.6, 0.7] 0.0538\n",
"(0.5, 0.6] 0.0010\n",
"Name: count, dtype: float64"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# relative frequency (proportion) of values in each bin\n",
"pd.cut(sim_means_series, median_bins).value_counts() / num_sims"
]
},
{
"cell_type": "markdown",
"id": "68d81e97",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- The histogram shows that the median values between 0.70 and 0.80 have the highest probability of occurring.\n",
"\n",
"- This means that a random sample would rarely yield a median value in the range of 0.5 - 0.6 or 1.0 - 1.1."
]
},
{
"cell_type": "markdown",
"id": "5888f9e1",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## What can we learn from simulation?\n",
"\n",
"- If we could generate all possible random samples of size 100, we would know all possible values of the statistic (the sample median), as well as the probabilities of all those values. We could visualize all the values and probabilities in the probability histogram of the statistic."
]
},
{
"cell_type": "markdown",
"id": "3c898c41",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- But in many situations including this one, the number of all possible samples is large enough to exceed the capacity of the computer, and purely mathematical calculations of the probabilities can be intractably difficult --- that is, if you have access to the whole population."
]
},
{
"cell_type": "markdown",
"id": "1d28b7fc",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- This is where empirical histograms come in."
]
},
{
"cell_type": "markdown",
"id": "24b20368",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- We know that by the Law of Averages, the empirical histogram of the statistic is likely to resemble the probability histogram of the statistic, if the sample size is large and if you repeat the random sampling process numerous times."
]
},
{
"cell_type": "markdown",
"id": "ce9dbb41",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- This means that simulating random processes repeatedly is a way of approximating probability distributions *without figuring out the probabilities mathematically or generating all possible random samples*. "
]
},
{
"cell_type": "markdown",
"id": "f2b39896",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- Thus computer simulations become a powerful tool in data science. They can help data scientists understand the properties of random quantities that would be complicated to analyze in other ways."
]
},
{
"cell_type": "markdown",
"id": "f6d31097-9363-4fb2-9ecc-1aac7b179527",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Models\n",
"\n",
"- A model is a set of assumptions about the data. In many cases models include assumptions about random (stochastic) processes used to generate the data.\n",
"\n",
"- Data scientists are often in a position of formulating and assessing models."
]
},
{
"cell_type": "markdown",
"id": "4cebbb19-5679-4ed2-9ffd-012b985fb312",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Goals of Data Science\n",
"\n",
"- Deeper understanding of the world.\n",
"\n",
"- Make the world a better place to live.\n",
"\n",
"- For example, help expose injustice.\n",
"\n",
"- The following example demonstrates how a model can achieve such goals."
]
},
{
"cell_type": "markdown",
"id": "8c9e9d64",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Jury Selection\n",
"\n",
"- U.S. Constitution grants equal protection under the law"
]
},
{
"cell_type": "markdown",
"id": "ecbd30b8",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- All defendants have the right to due process "
]
},
{
"cell_type": "markdown",
"id": "a2ee959c",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- Robert Swain, a Black man, was convicted in Talladega County, AL"
]
},
{
"cell_type": "markdown",
"id": "a4ce2705",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- He appealed to the U.S. Supreme Court\n",
"\n",
"- Main reason: Unfair jury selection in the County’s trials"
]
},
{
"cell_type": "markdown",
"id": "97cb022e",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- At the time of the trial, only men aged 21 or more were eligible to serve on juries in Talladega County. In that population, 26% of the men were Black. \n",
"\n",
"- But only eight men among the panel of 100 men (that is, 8%) were Black."
]
},
{
"cell_type": "markdown",
"id": "df094ad7",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- The U.S. Supreme Court reviewed the appeal and concluded, “the overall percentage disparity has been small.” But was this assertion reasonable? "
]
},
{
"cell_type": "markdown",
"id": "a6eef555-860c-47d6-93f4-3f82f0c227b0",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- If jury panelists were selected at random from the county’s eligible population, there would be some chance variation. We wouldn’t get exactly 26 Black panelists on every 100-person panel. But would we expect as few as eight?"
]
},
{
"cell_type": "markdown",
"id": "cc99dd90-1b5c-467e-b07d-14024fc83bb7",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## A model of random selection\n",
"\n",
"- A model of the data is that the panel was selected at random and ended up with a small number of Black panelists just due to chance.\n",
"\n",
"- Since the panel was supposed to resemble the population of all eligible jurors, the model of random selection is important to assess. Let’s see if it stands up to scrutiny.\n",
"\n",
"- The `numpy.random` function `multinomial(n, pvals, size)` can be used to simulate sample proportions or counts with two or more categories."
]
},
{
"cell_type": "markdown",
"id": "64577a5a-a0c2-400b-a2f5-f7281d742c86",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Example 1: Rolling a six-sided die 20 times"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "237a8f9b",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[3, 3, 4, 4, 5, 1]])"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import numpy as np\n",
"\n",
"# number of times die is rolled\n",
"sample_size = 20\n",
"\n",
"# number of experiments\n",
"num_simulations = 1\n",
"\n",
"# probability of each side\n",
"true_probabilities = [1/6] * 6\n",
"\n",
"# number of times each side appears\n",
"counts = np.random.multinomial(sample_size, true_probabilities, size=num_simulations)\n",
"\n",
"counts"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "51f82363-2dcb-4f92-8c66-80f0acd9897e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Sample counts: \n",
" [[3 3 4 4 5 1]]\n",
"Sample proportions: \n",
" [[0.15 0.15 0.2 0.2 0.25 0.05]]\n",
"True probabilities: \n",
" [0.16666666666666666, 0.16666666666666666, 0.16666666666666666, 0.16666666666666666, 0.16666666666666666, 0.16666666666666666]\n"
]
}
],
"source": [
"proportions = counts / sample_size\n",
"\n",
"print(\"Sample counts: \\n\", counts) \n",
"print(\"Sample proportions: \\n\", proportions)\n",
"print(\"True probabilities: \\n\", true_probabilities)"
]
},
{
"cell_type": "markdown",
"id": "0f2119c0-adcc-4095-90cf-b5e6f1dace60",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Example 2: Rolling a loaded six-side 100 times more likely to land on 6 - repeated 3 times"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "970611e6-05ef-4836-a1a4-57f076c7e0f9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Sample counts: \n",
" [[10 17 9 14 16 34]\n",
" [12 18 10 14 13 33]\n",
" [16 14 16 14 15 25]]\n",
"Sample proportions: \n",
" [[0.1 0.17 0.09 0.14 0.16 0.34]\n",
" [0.12 0.18 0.1 0.14 0.13 0.33]\n",
" [0.16 0.14 0.16 0.14 0.15 0.25]]\n",
"True probabilities: \n",
" [0.14285714285714285, 0.14285714285714285, 0.14285714285714285, 0.14285714285714285, 0.14285714285714285, 0.2857142857142857]\n"
]
}
],
"source": [
"sample_size = 100\n",
"\n",
"num_simulations = 3\n",
"\n",
"true_probabilities = [1 / 7] * 5 + [2 / 7]\n",
"\n",
"counts = np.random.multinomial(sample_size, true_probabilities, size=num_simulations)\n",
"\n",
"proportions = counts / sample_size\n",
"\n",
"print(\"Sample counts: \\n\", counts) \n",
"print(\"Sample proportions: \\n\", proportions)\n",
"print(\"True probabilities: \\n\", true_probabilities)"
]
},
{
"cell_type": "markdown",
"id": "4c6f2a05-cf84-4a68-8f8c-12eeb4313acc",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"- Let's use this to simulate the jury selection process.\n",
"\n",
"- The size of the jury panel is 100, so `sample_size` is 100. \n",
"\n",
"- The distribution from which we will draw the sample is the distribution in the population of eligible jurors: 26% of them were Black, so 100% - 26% = 74% are white (very simplistic assumption, but let's go with it for now). \n",
"\n",
"- This means `true_pobabilities` is `[0.26, 0.74]`.\n",
"\n",
"- One simulation is below."
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "eca2d986",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[23, 77]])"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sample_size = 100\n",
"\n",
"true_probabilities = [0.26, 0.74]\n",
"\n",
"num_simulations = 1\n",
"\n",
"counts = np.random.multinomial(sample_size, true_probabilities, size=num_simulations)\n",
"\n",
"counts"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "98626e7c-772e-4dda-b191-1c39e8fa735b",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Black
\n",
"
White
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
0.23
\n",
"
0.77
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" Black White\n",
"0 0.23 0.77"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"proportions = counts / sample_size\n",
"\n",
"sim_counts = pd.DataFrame(proportions, columns = [\"Black\", \"White\"])\n",
"\n",
"sim_counts"
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "e974b9ae",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.23"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sim_counts.iloc[0,0]"
]
},
{
"cell_type": "markdown",
"id": "d54a719d",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Quick side note on Indexing a `pandas` `DataFrame`\n",
"\n",
"- `iloc` and `loc` are both functions that can access data in a `DataFrame`.\n",
"\n",
"- `loc` uses labels to access rows and columns.\n",
"\n",
"- `iloc` uses integers to access rows and columns."
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "f8833462",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(1, 2)\n"
]
}
],
"source": [
"print(sim_counts.shape)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"id": "38a9f6d3",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"RangeIndex(start=0, stop=1, step=1)\n"
]
}
],
"source": [
"print(sim_counts.index)"
]
},
{
"cell_type": "code",
"execution_count": 27,
"id": "1426a8f6",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Index(['Black', 'White'], dtype='object')\n"
]
}
],
"source": [
"print(sim_counts.columns)"
]
},
{
"cell_type": "code",
"execution_count": 28,
"id": "29aaea0e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.77\n",
"0.77\n"
]
}
],
"source": [
"# row 0 column 1\n",
"print(sim_counts.iloc[0, 1])\n",
"\n",
"# access row 0 column 1 using label names\n",
"print(sim_counts.loc[0, \"White\"])"
]
},
{
"cell_type": "code",
"execution_count": 29,
"id": "69f40599",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.23\n",
"0.23\n"
]
}
],
"source": [
"print(sim_counts.iloc[0, 0])\n",
"\n",
"print(sim_counts.loc[0, \"Black\"])"
]
},
{
"cell_type": "markdown",
"id": "1e1069d8-4602-40fb-be47-7e21430eefc7",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Back to simulation ... Simulate one value\n",
"\n",
"- Let's write a function to simulate one value."
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "f3d23692-67a0-4205-b466-7f4a5db54f54",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"source": [
"def simulate_one_count():\n",
" sample_size = 100 # jury size\n",
"\n",
" true_probabilities = [0.26, 0.74] #true prob of race\n",
" \n",
" num_simulations = 1 # number of simulations\n",
" # get the random counts\n",
" counts = np.random.multinomial(sample_size, true_probabilities, size=num_simulations) \n",
" \n",
" # store in data frame\n",
" sim_counts = pd.DataFrame(counts, columns = [\"Black\", \"White\"])\n",
" \n",
" return sim_counts.iloc[0,0]\n"
]
},
{
"cell_type": "code",
"execution_count": 31,
"id": "3e371a58-6c50-458f-ad07-e0dd7592e1ce",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"29"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"simulate_one_count()"
]
},
{
"cell_type": "markdown",
"id": "06961dc2-9562-46ef-adc4-a190aa44d76d",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Simulate multiple values\n",
"\n",
"- Our analysis is focused on the variability in the counts. \n",
"\n",
"- Let’s generate 10,000 simulated values of the count and see how they vary.\n",
"\n",
"- We will do this by using a for loop and collecting all the simulated counts in a list called `sim_counts`"
]
},
{
"cell_type": "code",
"execution_count": 32,
"id": "ccea85b6-57ee-4291-86a7-3aa6c56a31d8",
"metadata": {},
"outputs": [],
"source": [
"sim_counts = []\n",
"\n",
"number_sims = 10000\n",
"\n",
"for _ in np.arange(number_sims):\n",
" sim_counts.append(simulate_one_count())"
]
},
{
"cell_type": "code",
"execution_count": 33,
"id": "cbc747ac-ea7c-4d32-8b37-15f6db0d9dc0",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAYAAABB4NqyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA1y0lEQVR4nO3de1wVdf7H8fcB5GaKl7OCKHgpFE3E9bqYRZsUppXIZq7rhcxfl0eSGGWtpdKu7dLlpw8t3cj9peX+1nTd0K2tNCK0LUlTtHIXDbvsMQX0bKXCghLM7w9/nu0kwuF44ByY1/PxmMd6Zj4z5zPMSm9nvjNjMQzDEAAAgIn4ebsBAACAlkYAAgAApkMAAgAApkMAAgAApkMAAgAApkMAAgAApkMAAgAAphPg7QZ8UV1dnY4dO6YOHTrIYrF4ux0AAOACwzB0+vRpRUZGys+v4XM8BKB6HDt2TFFRUd5uAwAAuOHIkSPq2bNngzUEoHp06NBB0rkfYMeOHb3cDQAAcMWpU6cUFRXl+O94QwhA9Th/2atjx44EIAAAWhlXhq8wCBoAAJgOAQgAAJgOAQgAAJgOAQgAAJgOAQgAAJgOAQgAAJgOAQgAAJgOAQgAAJgOAQgAAJgOAQgAAJgOAQgAAJgOAQgAAJgOAQgAAJgOAQgAAJhOgLcbAGBONptNdrvdpVqr1aro6Ohm7giAmRCAALQ4m82m2NhYVVVVuVQfEhKigwcPEoIAeAwBCECLs9vtqqqqUmpqqqxWa6O1ubm5stvtBCAAHkMAAuA1VqtVkZGR3m4DgAkRgAC0Ga6OK2JMEQACEIA2oSnjihhTBIAABKBNcHVcEWOKAEgEIABtDOOKALiCByECAADTIQABAADTIQABAADTIQABAADTIQABAADTIQABAADT4TZ4AK1CcXHxJS0HgO8jAAHwaRUVFbJYLJo+fbq3WwHQhhCAAHiUK+/jasrZmurqahmG0egTnktKSlRQUODydgGYGwEIgMc05X1cTdXYE55deQkqAJxHAALgMa6+j4uzNQC8jQAEwOM4WwPA13EbPAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB0CEAAAMB2vB6BVq1apd+/eCg4O1qhRo7R79+4G6zdt2qTY2FgFBwcrLi5Ob7zxhtPyiooKpaenq2fPngoJCdHAgQOVk5PTnLsAAABaGa8GoI0bNyozM1NZWVkqKipSfHy8kpOTdfz48Xrrd+7cqalTp2r27Nnat2+fUlJSlJKSogMHDjhqMjMztXXrVv3v//6viouLNW/ePKWnp+vVV19tqd0CAAA+zqsBaNmyZbrzzjs1a9Ysx5ma0NBQrVmzpt76FStWaNy4cZo/f74GDBigJUuWaOjQoVq5cqWjZufOnUpLS9O1116r3r1766677lJ8fHyDZ5bOnDmjU6dOOU0AAKDt8loAOnv2rPbu3aukpKT/NOPnp6SkJBUWFta7TmFhoVO9JCUnJzvVjx49Wq+++qqOHj0qwzBUUFCgTz/9VDfccMNFe8nOzlZYWJhjioqKusS9AwAAvsxrAchut6u2tlbh4eFO88PDw1VWVlbvOmVlZY3WP/vssxo4cKB69uypwMBAjRs3TqtWrdI111xz0V4WLFigkydPOqYjR45cwp4BAABfF+DtBjzt2Wef1QcffKBXX31VvXr10rvvvqs5c+YoMjLygrNH5wUFBSkoKKiFOwUAAN7itQBktVrl7++v8vJyp/nl5eWKiIiod52IiIgG66uqqvTII49o8+bNmjBhgiRp8ODB2r9/v/77v//7ogEIQONsNpvsdnuDNcXFxS3UDQBcGq8FoMDAQA0bNkz5+flKSUmRJNXV1Sk/P1/p6en1rpOQkKD8/HzNmzfPMS8vL08JCQmSpJqaGtXU1MjPz/nKnr+/v+rq6pplPwAzsNlsio2NVVVVlbdbAQCP8OolsMzMTKWlpWn48OEaOXKkli9frsrKSs2aNUuSNHPmTPXo0UPZ2dmSpIyMDCUmJmrp0qWaMGGCNmzYoD179mj16tWSpI4dOyoxMVHz589XSEiIevXqpR07dmjdunVatmyZ1/YTaO3sdruqqqqUmpoqq9V60bqSkhIVFBS0YGcA4B6vBqApU6boxIkTWrx4scrKyjRkyBBt3brVMdDZZrM5nc0ZPXq01q9fr4ULF+qRRx5RTEyMtmzZokGDBjlqNmzYoAULFmjatGn6+uuv1atXL/3mN7/RPffc0+L7B7Q1VqtVkZGRF13e2CUyAPAVXh8EnZ6eftFLXtu3b79g3uTJkzV58uSLbi8iIkJr1671VHsAAKAN8vqrMAAAAFoaAQgAAJgOAQgAAJiO18cAAYA3uPLMIqvVqujo6BboBkBLIwABMJWKigpZLBZNnz690dqQkBAdPHiQEAS0QQQgAKZSXV0twzAafaaR3W5Xbm6u7HY7AQhogwhAAEypsWcaAWjbGAQNAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMhwAEAABMJ8DbDQCALysuLnapzmq1Kjo6upm7AeApBCAAqEdFRYUsFoumT5/uUn1ISIgOHjxICAJaCQIQANSjurpahmEoNTVVVqu1wVq73a7c3FzZ7XYCENBKEIAAoAFWq1WRkZHebgOAhzEIGgAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmA4BCAAAmE6AtxsA4F02m012u73BmuLi4hbqBgBaBgEIMDGbzabY2FhVVVV5uxUAaFEEIMDE7Ha7qqqqlJqaKqvVetG6kpISFRQUtGBnANC8CEAAZLVaFRkZedHljV0iA4DWhkHQAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdLwegFatWqXevXsrODhYo0aN0u7duxus37Rpk2JjYxUcHKy4uDi98cYbF9QUFxfrlltuUVhYmNq3b68RI0bIZrM11y4AAIBWxqsBaOPGjcrMzFRWVpaKiooUHx+v5ORkHT9+vN76nTt3aurUqZo9e7b27dunlJQUpaSk6MCBA46azz77TGPGjFFsbKy2b9+ujz/+WIsWLVJwcHBL7RYAAPBxXg1Ay5Yt05133qlZs2Zp4MCBysnJUWhoqNasWVNv/YoVKzRu3DjNnz9fAwYM0JIlSzR06FCtXLnSUfPoo49q/Pjxeuqpp/TjH/9Yl19+uW655RZ169atpXYLAAD4OK8FoLNnz2rv3r1KSkr6TzN+fkpKSlJhYWG96xQWFjrVS1JycrKjvq6uTq+//rr69eun5ORkdevWTaNGjdKWLVsa7OXMmTM6deqU0wQAANourwUgu92u2tpahYeHO80PDw9XWVlZveuUlZU1WH/8+HFVVFToiSee0Lhx4/TWW29p0qRJSk1N1Y4dOy7aS3Z2tsLCwhxTVFTUJe4dAADwZV4fBO1JdXV1kqSJEyfq/vvv15AhQ/TLX/5SN910k3Jyci663oIFC3Ty5EnHdOTIkZZqGQAAeIHX3gZvtVrl7++v8vJyp/nl5eWKiIiod52IiIgG661WqwICAjRw4ECnmgEDBui99967aC9BQUEKCgpyZzcAAEAr5LUzQIGBgRo2bJjy8/Md8+rq6pSfn6+EhIR610lISHCql6S8vDxHfWBgoEaMGKFDhw451Xz66afq1auXh/cAAAC0Vl47AyRJmZmZSktL0/DhwzVy5EgtX75clZWVmjVrliRp5syZ6tGjh7KzsyVJGRkZSkxM1NKlSzVhwgRt2LBBe/bs0erVqx3bnD9/vqZMmaJrrrlGP/3pT7V161a99tpr2r59uzd2EQAA+CCvBqApU6boxIkTWrx4scrKyjRkyBBt3brVMdDZZrPJz+8/J6lGjx6t9evXa+HChXrkkUcUExOjLVu2aNCgQY6aSZMmKScnR9nZ2Zo7d6769++vV155RWPGjGnx/QMAAL7JqwFIktLT05Wenl7vsvrO2kyePFmTJ09ucJt33HGH7rjjDk+0BwAA2qA2dRcYAACAKwhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdAhAAADAdLz+KgwAaCuKi4sbrbFarYqOjm6BbgA0xK0A9Pnnn6tv376e7gUAWqWKigpZLBZNnz690dqQkBAdPHiQEAR4mVsB6IorrlBiYqJmz56tW2+9VcHBwZ7uCwBajerqahmGodTUVFmt1ovW2e125ebmym63E4AAL3NrDFBRUZEGDx6szMxMRURE6O6779bu3bs93RsAtCpWq1WRkZEXnRoKRwBallsBaMiQIVqxYoWOHTumNWvWqLS0VGPGjNGgQYO0bNkynThxwtN9AgAAeMwl3QUWEBCg1NRUbdq0SU8++aQOHz6sBx98UFFRUZo5c6ZKS0s91ScAAIDHXFIA2rNnj+699151795dy5Yt04MPPqjPPvtMeXl5OnbsmCZOnOipPgEAADzGrUHQy5Yt09q1a3Xo0CGNHz9e69at0/jx4+Xndy5P9enTRy+++KJ69+7tyV4BAAA8wq0A9Nxzz+mOO+7Q7bffru7du9db061bN73wwguX1BwAAEBzcCsAlZSUNFoTGBiotLQ0dzYPAADQrNwaA7R27Vpt2rTpgvmbNm3SSy+9dMlNAQAANCe3AlB2dna9z7Po1q2bfvvb315yUwAAAM3JrQBks9nUp0+fC+b36tVLNpvtkpsCAABoTm4FoG7duunjjz++YP5HH32krl27XnJTAAAAzcmtADR16lTNnTtXBQUFqq2tVW1trd555x1lZGTo5z//uad7BAAA8Ci37gJbsmSJvvzyS40dO1YBAec2UVdXp5kzZzIGCAAA+Dy3AlBgYKA2btyoJUuW6KOPPlJISIji4uLUq1cvT/cHAADgcW4FoPP69eunfv36eaoXAACAFuFWAKqtrdWLL76o/Px8HT9+XHV1dU7L33nnHY80BwAA0BzcCkAZGRl68cUXNWHCBA0aNEgWi8XTfQEAADQbtwLQhg0b9Kc//Unjx4/3dD8AAADNzq3b4AMDA3XFFVd4uhcAAIAW4dYZoAceeEArVqzQypUrufwF+CCbzSa73d5oXXFxcQt0AwC+x60A9N5776mgoEBvvvmmrrzySrVr185peW5urkeaA9B0NptNsbGxqqqq8nYrAOCz3ApAnTp10qRJkzzdCwAPsNvtqqqqUmpqar0vLf6+kpISFRQUtFBnAOA73ApAa9eu9XQfADzMarUqMjKywRpXLpMBQFvk1iBoSfruu+/09ttv6/nnn9fp06clSceOHVNFRYXHmgMAAGgObp0B+uc//6lx48bJZrPpzJkzuv7669WhQwc9+eSTOnPmjHJycjzdJwAAgMe4dQYoIyNDw4cP1zfffKOQkBDH/EmTJik/P99jzQEAADQHt84A/e1vf9POnTsVGBjoNL937946evSoRxoDAABoLm6dAaqrq1Ntbe0F87/66it16NDhkpsCAABoTm4FoBtuuEHLly93fLZYLKqoqFBWVhavxwAAAD7PrUtgS5cuVXJysgYOHKjq6mr94he/UElJiaxWq15++WVP9wgAAOBRbgWgnj176qOPPtKGDRv08ccfq6KiQrNnz9a0adOcBkUDAC7kyitIrFaroqOjW6AbwJzcCkCSFBAQoOnTp3uyFwBo0yoqKmSxWFz63RkSEqKDBw8SgoBm4lYAWrduXYPLZ86c6VYzANCWVVdXyzCMRl9TYrfblZubK7vdTgACmolbASgjI8Ppc01Njf79738rMDBQoaGhBCAAaIArrykB0Lzcugvsm2++cZoqKip06NAhjRkzhkHQAADA57n9LrAfiomJ0RNPPHHB2SEAAABf47EAJJ0bGH3s2DFPbhIAAMDj3BoD9Oqrrzp9NgxDpaWlWrlypa666iqPNAYAANBc3ApAKSkpTp8tFot+9KMf6brrrtPSpUs90RcAAECzcSsA1dXVeboPAACAFuPRMUAAAACtgVtngDIzM12uXbZsmTtfAQAA0GzcCkD79u3Tvn37VFNTo/79+0uSPv30U/n7+2vo0KGOOovF4pkuAQAAPMitAHTzzTerQ4cOeumll9S5c2dJ5x6OOGvWLF199dV64IEHPNokAACAJ7k1Bmjp0qXKzs52hB9J6ty5sx5//HHuAgMAAD7PrQB06tQpnThx4oL5J06c0OnTpy+5KQAAgObkVgCaNGmSZs2apdzcXH311Vf66quv9Morr2j27NlKTU31dI8AAAAe5dYYoJycHD344IP6xS9+oZqamnMbCgjQ7Nmz9fTTT3u0QQAAAE9zKwCFhobqd7/7nZ5++ml99tlnkqTLL79c7du392hzAAAAzeGSHoRYWlqq0tJSxcTEqH379jIMw1N9AQAANBu3AtC//vUvjR07Vv369dP48eNVWloqSZo9eza3wAMAAJ/nVgC6//771a5dO9lsNoWGhjrmT5kyRVu3bvVYcwAAAM3BrTFAb731lrZt26aePXs6zY+JidE///lPjzQGAADQXNw6A1RZWel05ue8r7/+WkFBQZfcFAAAQHNyKwBdffXVWrduneOzxWJRXV2dnnrqKf30pz9t8vZWrVql3r17Kzg4WKNGjdLu3bsbrN+0aZNiY2MVHBysuLg4vfHGGxetveeee2SxWLR8+fIm9wUAANomtwLQU089pdWrV+vGG2/U2bNn9dBDD2nQoEF699139eSTTzZpWxs3blRmZqaysrJUVFSk+Ph4JScn6/jx4/XW79y5U1OnTtXs2bO1b98+paSkKCUlRQcOHLigdvPmzfrggw8UGRnpzm4CAIA2yq0ANGjQIH366acaM2aMJk6cqMrKSqWmpmrfvn26/PLLm7StZcuW6c4779SsWbM0cOBA5eTkKDQ0VGvWrKm3fsWKFRo3bpzmz5+vAQMGaMmSJRo6dKhWrlzpVHf06FHdd999+uMf/6h27dq5s5sAAKCNavIg6JqaGo0bN045OTl69NFHL+nLz549q71792rBggWOeX5+fkpKSlJhYWG96xQWFiozM9NpXnJysrZs2eL4XFdXpxkzZmj+/Pm68sorG+3jzJkzOnPmjOPzqVOnmrgnAACgNWnyGaB27drp448/9siX2+121dbWKjw83Gl+eHi4ysrK6l2nrKys0fonn3xSAQEBmjt3rkt9ZGdnKywszDFFRUU1cU8AAEBr4tYlsOnTp+uFF17wdC8esXfvXq1YsUIvvviiLBaLS+ssWLBAJ0+edExHjhxp5i4BAIA3ufUcoO+++05r1qzR22+/rWHDhl3wDrBly5a5tB2r1Sp/f3+Vl5c7zS8vL1dERES960RERDRY/7e//U3Hjx9XdHS0Y3ltba0eeOABLV++XF9++eUF2wwKCuL2fQAATKRJAejzzz9X7969deDAAQ0dOlSS9OmnnzrVuHrWRZICAwM1bNgw5efnKyUlRdK58Tv5+flKT0+vd52EhATl5+dr3rx5jnl5eXlKSEiQJM2YMUNJSUlO6yQnJ2vGjBmaNWuWy70BAIC2q0kBKCYmRqWlpSooKJB07tUXzzzzzAVjcpoiMzNTaWlpGj58uEaOHKnly5ersrLSEVZmzpypHj16KDs7W5KUkZGhxMRELV26VBMmTNCGDRu0Z88erV69WpLUtWtXde3a1ek72rVrp4iICPXv39/tPgEAQNvRpAD0w7e9v/nmm6qsrLykBqZMmaITJ05o8eLFKisr05AhQ7R161ZHqLLZbPLz+89QpdGjR2v9+vVauHChHnnkEcXExGjLli0aNGjQJfUBAADMw60xQOf9MBC5Kz09/aKXvLZv337BvMmTJ2vy5Mkub7++cT8AAMC8mnQXmMViuWCMT1PG/AAAAPiCJl8Cu/322x13TFVXV+uee+654C6w3Nxcz3UIAADgYU0KQGlpaU6fp0+f7tFmAAAAWkKTAtDatWubqw8AAIAWc0mDoAEAzae4uLjRGqvV6vTgVwCuIQABgI+pqKiQxWJxaZhBSEiIDh48SAgCmogABAA+prq6WoZhKDU1VVar9aJ1drtdubm5stvtBCCgiQhAAOCjrFarIiMjvd0G0Ca59TZ4AACA1owABAAATIdLYEArYrPZZLfbG6xx5c4hADA7AhDQSthsNsXGxqqqqsrbrQBAq0cAAloJu92uqqqqRu8MKikpUUFBQQt2BgCtDwEIaGUauzOosUtkAAAGQQMAABMiAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMhAAEAANMJ8HYDAIBLU1xc7FKd1WpVdHR0M3cDtA4EIABopSoqKmSxWDR9+nSX6kNCQnTw4EFCECACEAC0WtXV1TIMQ6mpqbJarQ3W2u125ebmym63E4AAEYAAoNWzWq2KjIz0dhtAq8IgaAAAYDoEIAAAYDoEIAAAYDoEIAAAYDoMggZ8gM1mk91ub7DG1We9AAAaRwACvMxmsyk2NlZVVVXebgUATIMABHiZ3W5XVVVVo89yKSkpUUFBQQt2BgBtFwEI8BGNPculsUtkAADXMQgaAACYjk8EoFWrVql3794KDg7WqFGjtHv37gbrN23apNjYWAUHBysuLk5vvPGGY1lNTY0efvhhxcXFqX379oqMjNTMmTN17Nix5t4NAADQSng9AG3cuFGZmZnKyspSUVGR4uPjlZycrOPHj9dbv3PnTk2dOlWzZ8/Wvn37lJKSopSUFB04cECS9O9//1tFRUVatGiRioqKlJubq0OHDumWW25pyd0CAAA+zOsBaNmyZbrzzjs1a9YsDRw4UDk5OQoNDdWaNWvqrV+xYoXGjRun+fPna8CAAVqyZImGDh2qlStXSpLCwsKUl5en2267Tf3799dPfvITrVy5Unv37pXNZqt3m2fOnNGpU6ecJgAA0HZ5NQCdPXtWe/fuVVJSkmOen5+fkpKSVFhYWO86hYWFTvWSlJycfNF6STp58qQsFos6depU7/Ls7GyFhYU5pqioqKbvDAAAaDW8GoDsdrtqa2sVHh7uND88PFxlZWX1rlNWVtak+urqaj388MOaOnWqOnbsWG/NggULdPLkScd05MgRN/YGAAC0Fm36NviamhrddtttMgxDzz333EXrgoKCFBQU1IKdAQAAb/JqALJarfL391d5ebnT/PLyckVERNS7TkREhEv158PPP//5T73zzjsXPfsDAADMx6uXwAIDAzVs2DDl5+c75tXV1Sk/P18JCQn1rpOQkOBUL0l5eXlO9efDT0lJid5++2117dq1eXYAAAC0Sl6/BJaZmam0tDQNHz5cI0eO1PLly1VZWalZs2ZJkmbOnKkePXooOztbkpSRkaHExEQtXbpUEyZM0IYNG7Rnzx6tXr1a0rnwc+utt6qoqEh//etfVVtb6xgf1KVLFwUGBnpnRwEAgM/wegCaMmWKTpw4ocWLF6usrExDhgzR1q1bHQOdbTab/Pz+c6Jq9OjRWr9+vRYuXKhHHnlEMTEx2rJliwYNGiRJOnr0qF599VVJ0pAhQ5y+q6CgQNdee22L7BcAAPBdXg9AkpSenq709PR6l23fvv2CeZMnT9bkyZPrre/du7cMw/BkewAAoI3x+oMQAQAAWhoBCAAAmA4BCAAAmA4BCAAAmI5PDIIGALSM4uLiRmusVquio6NboBvAewhAAGACFRUVslgsmj59eqO1ISEhOnjwICEIbRoBCGhGNptNdru9wRpX/kUOXKrq6moZhqHU1FRZrdaL1tntduXm5sputxOA0KYRgIBmYrPZFBsbq6qqKm+3AjhYrVZFRkZ6uw3A6whAQDOx2+2qqqpq9F/cJSUlKigoaMHOAAAEIKCZNfYv7sYukQEAPI/b4AEAgOkQgAAAgOkQgAAAgOkQgAAAgOkQgAAAgOkQgAAAgOkQgAAAgOkQgAAAgOkQgAAAgOnwJGigiVx5wanES04BwJcRgIAm4AWnANA2EICAJnD1BacSLzkFAF9GAALc0NgLTiVecgoAvoxB0AAAwHQIQAAAwHQIQAAAwHQIQAAAwHQYBA0AuIArz7GyWq2Kjo5ugW4AzyMAAQAcKioqZLFYNH369EZrQ0JCdPDgQUIQWiUCEADAobq6WoZhNPqsK7vdrtzcXNntdgIQWiUCEADgAq486wpozRgEDQAATIcABAAATIcABAAATIcABAAATIcABAAATIcABAAATIcABAAATIfnAAEA3MYrM9BaEYAAAE3GKzPQ2hGAAABNxisz0NoRgAAAbuOVGWitCEDA/7PZbLLb7Q3WuDLeAQDg+whAgM6Fn9jYWFVVVXm7FQBACyAAATo3TqGqqqrR8QwlJSUqKChowc4AAM2BAAR8T2PjGRq7RAYAaB14ECIAADAdAhAAADAdAhAAADAdxgABAJqdq4+Q4LUZaCkEIABAs2nKKzMkXpuBlkMAAgA0G1dfmSHx2gy0LAIQAKDZ8coM+BoGQQMAANMhAAEAANMhAAEAANMhAAEAANNhEDQAwKe48swgnheES0UAQptns9kafYmpqw9pA9B8mvLMIJ4XhEtFAEKr5UqwKS0t1a233qrq6uoW6gqAu1x9ZhDPC4InEIDQKtlsNsXGxqqqqsql+sZ+oZaUlKigoMBT7QG4BDwzCC2BAIRWyW63q6qqyuVg09gv1MbOJAEA2hYCUFtSaZM+f1E6XSJ9d1oK6CB1iJH63i61b5uniQk2AAB3EIDagvId0sGl0tG/Spb/f7KBUStZ/M/9+ZPHpB43SQMelLpd47U2AQDwFQSg1swwzgWfffMlS4Ak41zwcSz/3p+PvSkdfU368X9LsZmSxdLi7QKAJ3G7PC4FAag1O7jsXPiRJOO7hmvPL9/34Ln/HfBA8/UFAM2I2+XhCQSg1qp8x3/CTFPte1DqOoLLYQBaJW6Xhyf4RABatWqVnn76aZWVlSk+Pl7PPvusRo4cedH6TZs2adGiRfryyy8VExOjJ598UuPHj3csNwxDWVlZ+v3vf69vv/1WV111lZ577jnFxMS0xO60jINLz132auzMT30sAVLxUp8MQK4820fiwYUAXL9dnktlqI/XA9DGjRuVmZmpnJwcjRo1SsuXL1dycrIOHTqkbt26XVC/c+dOTZ06VdnZ2brpppu0fv16paSkqKioSIMGDZIkPfXUU3rmmWf00ksvqU+fPlq0aJGSk5P1j3/8Q8HBwS29i55XaTs34FmGe+sb350bD1R5RGof5dHWLkVTn+0DAA1pyqWyoKAgvfLKK+revXuDdQSltsPrAWjZsmW68847NWvWLElSTk6OXn/9da1Zs0a//OUvL6hfsWKFxo0bp/nzz419WbJkifLy8rRy5Url5OTIMAwtX75cCxcu1MSJEyVJ69atU3h4uLZs2aKf//znLbdzzeXzF8/d7fX9Qc5NZfGTPl8rxS32WFsNcfV1FK4820fiwYUAGufqpTKbzaZt27bppptuanSbrgYlyfNhydUz5IQ013g1AJ09e1Z79+7VggULHPP8/PyUlJSkwsLCetcpLCxUZmam07zk5GRt2bJFkvTFF1+orKxMSUlJjuVhYWEaNWqUCgsL6w1AZ86c0ZkzZxyfT548KUk6deqU2/vWkLKyMpWVlTVa5+fnp7q6ugvmR5e9p86VdbqU+7gMGfqm+H3Zvnm3Sd/tTl15eblmzJjh9DNuSE1Njc6ePdtgzXffnbv0V1pa2mDtiRMnvFLnze+mR3r0pe/2hR4b+51SWVkpwzA0evRodezYscHt7d2716WgJJ0LS3/4wx8UHh7eYJ2nf4+6+r2ufndz1EVERCgiIqLRuqY6/99tw3DhConhRUePHjUkGTt37nSaP3/+fGPkyJH1rtOuXTtj/fr1TvNWrVpldOvWzTAMw3j//fcNScaxY8ecaiZPnmzcdttt9W4zKyvL0LnrSUxMTExMTEytfDpy5EijGcTrl8B8wYIFC5zOKtXV1enrr79W165dZeF5OZLOpeqoqCgdOXKkwX8hoWVxXHwXx8Y3cVx8lyeOjWEYOn36tEuD470agKxWq/z9/VVeXu40v7y8/KKnxiIiIhqsP/+/5eXlTtdoy8vLNWTIkHq3GRQUpKCgIKd5nTp1asqumEbHjh35peGDOC6+i2PjmzguvutSj01YWJhLdX5uf4MHBAYGatiwYcrPz3fMq6urU35+vhISEupdJyEhwalekvLy8hz1ffr0UUREhFPNqVOntGvXrotuEwAAmIvXL4FlZmYqLS1Nw4cP18iRI7V8+XJVVlY67gqbOXOmevTooezsbElSRkaGEhMTtXTpUk2YMEEbNmzQnj17tHr1akmSxWLRvHnz9PjjjysmJsZxG3xkZKRSUlK8tZsAAMCHeD0ATZkyRSdOnNDixYtVVlamIUOGaOvWrY7R6zabTX5+/zlRNXr0aK1fv14LFy7UI488opiYGG3ZssXxDCBJeuihh1RZWam77rpL3377rcaMGaOtW7e2jWcAeUlQUJCysrIuuFQI7+K4+C6OjW/iuPiulj42FsNw5V4xAACAtsOrY4AAAAC8gQAEAABMhwAEAABMhwAEAABMhwAEh3fffVc333yzIiMjZbFYHO9XO88wDC1evFjdu3dXSEiIkpKSVFJS4p1mTSQ7O1sjRoxQhw4d1K1bN6WkpOjQoUNONdXV1ZozZ466du2qyy67TD/72c8ueGAoPO+5557T4MGDHQ9uS0hI0JtvvulYznHxDU888YTjESnncWy847HHHpPFYnGaYmNjHctb8rgQgOBQWVmp+Ph4rVq1qt7lTz31lJ555hnl5ORo165dat++vZKTk1VdXd3CnZrLjh07NGfOHH3wwQfKy8tTTU2NbrjhBlVWVjpq7r//fr322mvatGmTduzYoWPHjik1NdWLXZtDz5499cQTT2jv3r3as2ePrrvuOk2cOFF///vfJXFcfMGHH36o559/XoMHD3aaz7HxniuvvFKlpaWO6b333nMsa9Hj0ujbwmBKkozNmzc7PtfV1RkRERHG008/7Zj37bffGkFBQcbLL7/shQ7N6/jx44YkY8eOHYZhnDsO7dq1MzZt2uSoKS4uNiQZhYWF3mrTtDp37mz8z//8D8fFB5w+fdqIiYkx8vLyjMTERCMjI8MwDP7OeFNWVpYRHx9f77KWPi6cAYJLvvjiC5WVlSkpKckxLywsTKNGjVJhYaEXOzOfkydPSpK6dOkiSdq7d69qamqcjk1sbKyio6M5Ni2otrZWGzZsUGVlpRISEjguPmDOnDmaMGGC0zGQ+DvjbSUlJYqMjFTfvn01bdo02Ww2SS1/XLz+JGi0DmVlZZLkeEL3eeHh4Y5laH51dXWaN2+errrqKsfTz8vKyhQYGHjBC3w5Ni3jk08+UUJCgqqrq3XZZZdp8+bNGjhwoPbv389x8aINGzaoqKhIH3744QXL+DvjPaNGjdKLL76o/v37q7S0VL/61a909dVX68CBAy1+XAhAQCsyZ84cHThwwOmaObyrf//+2r9/v06ePKk///nPSktL044dO7zdlqkdOXJEGRkZysvL4xVIPubGG290/Hnw4MEaNWqUevXqpT/96U8KCQlp0V64BAaXRERESNIFo/HLy8sdy9C80tPT9de//lUFBQXq2bOnY35ERITOnj2rb7/91qmeY9MyAgMDdcUVV2jYsGHKzs5WfHy8VqxYwXHxor179+r48eMaOnSoAgICFBAQoB07duiZZ55RQECAwsPDOTY+olOnTurXr58OHz7c4n9nCEBwSZ8+fRQREaH8/HzHvFOnTmnXrl1KSEjwYmdtn2EYSk9P1+bNm/XOO++oT58+TsuHDRumdu3aOR2bQ4cOyWazcWy8oK6uTmfOnOG4eNHYsWP1ySefaP/+/Y5p+PDhmjZtmuPPHBvfUFFRoc8++0zdu3dv8b8zXAKDQ0VFhQ4fPuz4/MUXX2j//v3q0qWLoqOjNW/ePD3++OOKiYlRnz59tGjRIkVGRiolJcV7TZvAnDlztH79ev3lL39Rhw4dHNfCw8LCFBISorCwMM2ePVuZmZnq0qWLOnbsqPvuu08JCQn6yU9+4uXu27YFCxboxhtvVHR0tE6fPq3169dr+/bt2rZtG8fFizp06OAYI3de+/bt1bVrV8d8jo13PPjgg7r55pvVq1cvHTt2TFlZWfL399fUqVNb/u+Mx+8rQ6tVUFBgSLpgSktLMwzj3K3wixYtMsLDw42goCBj7NixxqFDh7zbtAnUd0wkGWvXrnXUVFVVGffee6/RuXNnIzQ01Jg0aZJRWlrqvaZN4o477jB69eplBAYGGj/60Y+MsWPHGm+99ZZjOcfFd3z/NnjD4Nh4y5QpU4zu3bsbgYGBRo8ePYwpU6YYhw8fdixvyeNiMQzD8HysAgAA8F2MAQIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAAIAAKZDAALgcV9++aUsFov279/v7VaaxGKxaMuWLd5uw2fcfvvtvOoGbRYBCGgFysrKdN9996lv374KCgpSVFSUbr75ZqeXBrYUV0JCVFSUSktLL3gfEwD4Cl6GCvi4L7/8UldddZU6deqkp59+WnFxcaqpqdG2bds0Z84cHTx40NstXsDf318REREe325NTY3atWvn8e0CMB/OAAE+7t5775XFYtHu3bv1s5/9TP369dOVV16pzMxMffDBB446m82miRMn6rLLLlPHjh112223qby83LG8vssZ8+bN07XXXuv4fO2112ru3Ll66KGH1KVLF0VEROixxx5zLO/du7ckadKkSbJYLI7PP/TDS2Dbt2+XxWJRfn6+hg8frtDQUI0ePVqHDh266H6f38bGjRuVmJio4OBg/fGPf9S//vUvTZ06VT169FBoaKji4uL08ssvO63b2H5IUklJia655hoFBwdr4MCBysvLu6CHTz75RNddd51CQkLUtWtX3XXXXaqoqLjgZ/rb3/5W4eHh6tSpk37961/ru+++0/z589WlSxf17NlTa9euveh+StKf//xnxcXFOb4nKSlJlZWVkqQPP/xQ119/vaxWq8LCwpSYmKiioiKn9S0Wi55//nnddNNNCg0N1YABA1RYWKjDhw/r2muvVfv27TV69Gh99tlnjnUee+wxDRkyRM8//7yioqIUGhqq2267TSdPnrxon3V1dcrOzlafPn0UEhKi+Ph4/fnPf25w3wBfRQACfNjXX3+trVu3as6cOWrfvv0Fyzt16iTp3H+YJk6cqK+//lo7duxQXl6ePv/8c02ZMqXJ3/nSSy+pffv22rVrl5566in9+te/doSDDz/8UJK0du1alZaWOj676tFHH9XSpUu1Z88eBQQE6I477mh0nV/+8pfKyMhQcXGxkpOTVV1drWHDhun111/XgQMHdNddd2nGjBnavXu3y/tRV1en1NRUBQYGateuXcrJydHDDz/stH5lZaWSk5PVuXNnffjhh9q0aZPefvttpaenO9W98847OnbsmN59910tW7ZMWVlZuummm9S5c2ft2rVL99xzj+6++2599dVX9e5faWmppk6dqjvuuEPFxcXavn27UlNTdf491adPn1ZaWpree+89ffDBB4qJidH48eN1+vRpp+0sWbJEM2fO1P79+xUbG6tf/OIXuvvuu7VgwQLt2bNHhmFc0Pvhw4f1pz/9Sa+99pq2bt2qffv26d57773oscjOzta6deuUk5Ojv//977r//vs1ffp07dixo4EjCPioZnnHPACP2LVrlyHJyM3NbbDurbfeMvz9/Q2bzeaY9/e//92QZOzevdswDMNIS0szJk6c6LReRkaGkZiY6PicmJhojBkzxqlmxIgRxsMPP+z4LMnYvHlzg/188cUXhiRj3759hmEYRkFBgSHJePvttx01r7/+uiHJqKqqanAby5cvb/C7DMMwJkyYYDzwwAMu78e2bduMgIAA4+jRo47lb775ptO+rV692ujcubNRUVHh1LOfn59RVlZmGMa5n2mvXr2M2tpaR03//v2Nq6++2vH5u+++M9q3b2+8/PLL9fa+d+9eQ5Lx5ZdfNrqfhmEYtbW1RocOHYzXXnvNMU+SsXDhQsfnwsJCQ5LxwgsvOOa9/PLLRnBwsONzVlaW4e/vb3z11VdOPwM/Pz+jtLTUsX/n/z9TXV1thIaGGjt37nTqZ/bs2cbUqVNd6h3wJZwBAnyY8f9nARpTXFysqKgoRUVFOeYNHDhQnTp1UnFxcZO+c/DgwU6fu3fvruPHjzdpG65su3v37pLU6LaHDx/u9Lm2tlZLlixRXFycunTpossuu0zbtm2TzWa76Hed/77z33X+5xUZGelYnpCQ4FRfXFys+Ph4pzNvV111lerq6pwu3V155ZXy8/vPr9Lw8HDFxcU5Pvv7+6tr164X3c/4+HiNHTtWcXFxmjx5sn7/+9/rm2++cSwvLy/XnXfeqZiYGIWFhaljx46qqKhocH/Dw8MlyamP8PBwVVdX69SpU4550dHR6tGjh9PP4If7d97hw4f173//W9dff70uu+wyx7Ru3TqnS2tAa8EgaMCHxcTEyGKxeGSgs5+f3wWBqqam5oK6Hw4ytlgsqquru+Tv/+G2LRaLJDW67R9e+nv66ae1YsUKLV++XHFxcWrfvr3mzZuns2fPXvS7zn+fp/ajse9pynf7+/srLy9PO3fu1FtvvaVnn31Wjz76qHbt2qU+ffooLS1N//rXv7RixQr16tVLQUFBSkhIaHB/z/9s3fl5X8z5sU+vv/66U2iSpKCgILe2CXgTZ4AAH9alSxclJydr1apVjkGx3/ftt99KkgYMGKAjR47oyJEjjmX/+Mc/9O2332rgwIGSpB/96EcqLS11Wt+d5/S0a9dOtbW1TV7PU95//31NnDhR06dPV3x8vPr27atPP/20Sds4//P6/s/j+wPKz9d89NFHTj/3999/X35+furfv/+l7cQPWCwWXXXVVfrVr36lffv2KTAwUJs3b3Z859y5czV+/HhdeeWVCgoKkt1u98j32mw2HTt2zPH5gw8+uOj+DRw4UEFBQbLZbLriiiucpu+feQRaCwIQ4ONWrVql2tpajRw5Uq+88opKSkpUXFysZ555xnHZJikpSXFxcZo2bZqKioq0e/duzZw5U4mJiY5LSNddd5327NmjdevWqaSkRFlZWTpw4ECT++ndu7fy8/NVVlbmdKmmpcTExDjOmBQXF+vuu+92utvNFUlJSerXr5/S0tL00Ucf6W9/+5seffRRp5pp06YpODhYaWlpOnDggAoKCnTfffdpxowZjktMnrBr1y799re/1Z49e2Sz2ZSbm6sTJ05owIABjv39wx/+oOLiYu3atUvTpk1TSEiIR777/P6d/xnMnTtXt912W72PMOjQoYMefPBB3X///XrppZf02WefqaioSM8++6xeeuklj/QDtCQCEODj+vbtq6KiIv30pz/VAw88oEGDBun6669Xfn6+nnvuOUnnziD85S9/UefOnXXNNdcoKSlJffv21caNGx3bSU5O1qJFi/TQQw9pxIgROn36tGbOnNnkfpYuXaq8vDxFRUXpxz/+scf201ULFy7U0KFDlZycrGuvvVYRERFNflqxn5+fNm/erKqqKo0cOVL/9V//pd/85jdONaGhodq2bZu+/vprjRgxQrfeeqvGjh2rlStXenBvpI4dO+rdd9/V+PHj1a9fPy1cuFBLly7VjTfeKEl64YUX9M0332jo0KGaMWOG5s6dq27dunnku6+44gqlpqZq/PjxuuGGGzR48GD97ne/u2j9kiVLtGjRImVnZ2vAgAEaN26cXn/9dfXp08cj/QAtyWK4OsoSANBmPPbYY9qyZUure10J4CmcAQIAAKZDAAIAAKbDJTAAAGA6nAECAACmQwACAACmQwACAACmQwACAACmQwACAACmQwACAACmQwACAACmQwACAACm83+VzQYKfNdtrQAAAABJRU5ErkJggg==",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"plt.hist(sim_counts, bins = np.arange(5.5, 50, 1), edgecolor = 'black', \n",
" color = 'grey', density = True);\n",
"plt.xlabel('Count in random sample')\n",
"plt.ylabel('Frequency')\n",
"plt.scatter(8, 0, color = 'orange', s = 150);"
]
},
{
"cell_type": "markdown",
"id": "830c8bf5-d743-41db-a2df-24019c281c07",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Conclusion of the data analysis\n",
"\n",
"- The histogram shows that if we select a panel of size 100 at random from the eligible population, we are very unlikely to get counts of Black panelists that are as low as the eight that were observed on the panel in the trial.\n",
"\n",
"- This is evidence that the model of random selection of the jurors in the panel is not consistent with the data from the panel. While it is possible that the panel could have been generated by chance, our simulation demonstrates that it is hugely unlikely.\n",
"\n",
"- Therefore the most *reasonable* conclusion is that the assumption of random selection is unjustified for this jury panel."
]
},
{
"cell_type": "markdown",
"id": "8c372dea-0534-444e-b3ba-0ebecb785b89",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"- The simulation also could have been done using `np.random.multinomial`.\n",
"\n",
"- This is an example of a 'vectorized' computation, and are usually faster than non-vectorized computations."
]
},
{
"cell_type": "code",
"execution_count": 34,
"id": "f665c7b0-8d5f-4e85-94a7-4127c52615e0",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"10000"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sample_size = 100\n",
"\n",
"true_probabilities = [0.26, 0.74]\n",
"\n",
"num_simulations = 10000\n",
"\n",
"counts = np.random.multinomial(sample_size, true_probabilities, size=num_simulations)\n",
"\n",
"len(counts)"
]
},
{
"cell_type": "code",
"execution_count": 35,
"id": "be364d12",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[26, 74],\n",
" [24, 76],\n",
" [24, 76],\n",
" ...,\n",
" [34, 66],\n",
" [22, 78],\n",
" [27, 73]])"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"counts"
]
},
{
"cell_type": "markdown",
"id": "338e88fa-54e3-4150-8ba0-c35aab5f3600",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Comparing two samples \n",
"\n",
"### Are mental health visit rates higher in Toronto Neighbourhoods with higher material deprivation?\n",
"\n",
"> Material deprivation is closely connected to poverty and it refers to inability for individuals and communities to access and attain basic material needs. The indicators included in this dimension measure income, quality of housing, educational attainment, and family structure characteristics. \n",
"\n",
"See [2011 Ontario Marginalization Index Documentation](https://www.publichealthontario.ca/-/media/Documents/O/2017/on-marg-technical.pdf?la=en≻lang=en&hash=EED54DF437EDEDA2DFE1A00A4B14A50A) and [Toronto Health Profiles website](http://www.torontohealthprofiles.ca/index.php?varTab=HPDtbl)"
]
},
{
"cell_type": "markdown",
"id": "b99a6b64",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Data wrangling\n",
"\n",
"- The next few slides involve reading the data into pandas and getting it ready for analysis.\n",
"- This won't be covered in detail in class, but we have already covered this process in previous classes."
]
},
{
"cell_type": "markdown",
"id": "29dc46a3-8122-450f-b237-71bb36b92dd4",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- The data is stored in `1_marg_neighb_toronto_2006_OnMarg.xls` - a Microsoft Excel file format with file extension `.xls`.\n",
"\n",
"- Use `pandas` function `read_excel` with `sheet_name` parameter."
]
},
{
"cell_type": "markdown",
"id": "2b7c8916-4b64-4249-96b5-b2cd552620b2",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### Neighbourhood deprivation scores"
]
},
{
"cell_type": "code",
"execution_count": 36,
"id": "c9df7fb0",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Requirement already satisfied: xlrd in /Users/moon/anaconda3/envs/ggr274/lib/python3.10/site-packages (2.0.1)\n",
"\n",
"\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m A new release of pip is available: \u001b[0m\u001b[31;49m23.3.2\u001b[0m\u001b[39;49m -> \u001b[0m\u001b[32;49m25.0.1\u001b[0m\n",
"\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m To update, run: \u001b[0m\u001b[32;49mpip install --upgrade pip\u001b[0m\n",
"Note: you may need to restart the kernel to use updated packages.\n"
]
}
],
"source": [
"# this python library is not pre-installed\n",
"# and is required to use `pandas` `read_excel`\n",
"%pip install xlrd"
]
},
{
"cell_type": "code",
"execution_count": 37,
"id": "4a8ca612-4deb-42ea-9735-cd5baaa9c8ea",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Unnamed: 0
\n",
"
Ontario Marginalization Index. Toronto Neighbourhoods 2006, \\nQuintiles: Material Deprivation, Residential Instability, Dependency, Ethnic Concentration
"
],
"text/plain": [
" n_id name mh_visit_rates_mf deprivation depriv_binary\n",
"131 136 West Hill 8.9 0.8138 High\n",
"132 137 Woburn 7.2 0.5257 High\n",
"133 138 Eglinton East 8.3 1.0344 High\n",
"134 139 Scarborough Village 9.2 1.3915 High\n",
"135 140 Guildwood 8.2 -0.4834 Low"
]
},
"execution_count": 49,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mh_visit_depriv = mh_visit_depriv.copy()\n",
"\n",
"# create a new column in mhvisitdepric called depriv_HL\n",
"# classify neighbourhoods above the median as High\n",
"\n",
"mh_visit_depriv.loc[mh_visit_depriv[\"deprivation\"] > median_depriv, \"depriv_binary\"] = \"High\"\n",
"\n",
"# classify neighbourhoods at or below the median as Low\n",
"mh_visit_depriv.loc[mh_visit_depriv[\"deprivation\"] <= median_depriv, \"depriv_binary\"] = \"Low\"\n",
"\n",
"mh_visit_depriv.tail(n=5) # print last 3 rows of mhvisitdepriv"
]
},
{
"cell_type": "markdown",
"id": "416f67d7-72b4-44cc-8682-55adc6d50eea",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Do Neighbourhoods with high deprivation have more mental health visits compared to Neighbourhoods with low deprivation?\n",
"\n",
"- We can compare the mean mental health visit rates between the two groups."
]
},
{
"cell_type": "code",
"execution_count": 50,
"id": "29eb5419-a1c2-41ca-aa88-217c17a63e76",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [
{
"data": {
"text/plain": [
"depriv_binary\n",
"High 8.480882\n",
"Low 7.733824\n",
"Name: mh_visit_rates_mf, dtype: float64"
]
},
"execution_count": 50,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mean_table = mh_visit_depriv.groupby(\"depriv_binary\")[\"mh_visit_rates_mf\"].mean()\n",
"mean_table"
]
},
{
"cell_type": "code",
"execution_count": 51,
"id": "5f2e46d8-d86f-427a-b962-cabc24798599",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.7470588235294136"
]
},
"execution_count": 51,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"observed_mean_difference = mean_table.iloc[0] - mean_table.iloc[1]\n",
"observed_mean_difference"
]
},
{
"cell_type": "markdown",
"id": "815883ca",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- Does the difference represent a _real_ diffrence between the two populations or did we have a peculiar group in 2012?"
]
},
{
"cell_type": "markdown",
"id": "1cf9fe2a-cecd-49d5-9300-9590f4959bdd",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## The Logic of Hypothesis Testing\n",
"\n",
"### 1. Hypotheses\n",
"\n",
"Two claims:\n",
"\n",
"1. There is no difference in the mean mental health visit rates between high and low deprivation neighbourhoods. This is called the **null** hypothesis.\n",
"\n",
"2. There is a difference in the mean mental health visit rates between high and low deprivation neighbourhoods. This is called the **alternative** hypothesis."
]
},
{
"cell_type": "markdown",
"id": "396a61b0-0ea1-4fba-b654-3765b3b4bbc1",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### 2. Test statistic\n",
"\n",
"The test statistic is a number, **calculated from the data**, that captures what we're interested in.\n",
"\n",
"What would be a useful test statistic for this study?"
]
},
{
"cell_type": "markdown",
"id": "b2e72994",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### 3. Simulate what the null hypothesis predicts will happen\n",
"\n",
"- If the null hypothesis is true then the mean values of high and low deprivation neighbourhoods will be the same regardless of how they are named or labelled. \n",
"\n",
"- That means we can randomly assign (or shuffle) the neigbourhood names to high and low deprivation and the mean difference should be close to 0."
]
},
{
"cell_type": "markdown",
"id": "5e099fb7",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"+ Imagine we have 68 playing cards labelled `High` and 68 cards labelled `Low`."
]
},
{
"cell_type": "markdown",
"id": "00fcfb17",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- Shuffle the cards ..."
]
},
{
"cell_type": "markdown",
"id": "f5a9c833",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- Assign the cards to the 136 neighbourhoods then calculate the mean difference between high and low. This is one simulated value of the test statistic. "
]
},
{
"cell_type": "markdown",
"id": "5727b614",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- Shuffle the cards again ..."
]
},
{
"cell_type": "markdown",
"id": "49af6fa4",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- Assign the cards to the 136 neighbourhoods then calculate the mean difference between high and low. This is another simulated value of the test statistic. "
]
},
{
"cell_type": "markdown",
"id": "88bfd64b-4672-4fda-aec2-98e8d0b835da",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- Continue shuffling, assigning to neigbourhoods, and computing the mean difference."
]
},
{
"cell_type": "markdown",
"id": "44c9b5ce-be09-41e1-973d-26f2da32b541",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Shuffling"
]
},
{
"cell_type": "markdown",
"id": "03a939ea-d800-4c95-b7f0-0a5097ccd498",
"metadata": {},
"source": [
"- The _observed_ difference in mean rate of mental health visits between high and low deprivation is 0.75.\n",
"\n",
"- Could this difference be due to chance?\n",
"\n",
"- Let's repeat this study **assuming that the difference is due to chance**.\n",
"\n",
"- Suppose that the (true) mean mental health visit rates in high deprivation neigbourhoods is equal to the (true) mean mental health visit rates in low deprivation neigbourhoods. Then the **labels** of `\"depriv_binary\"` `High` and `Low` are ..."
]
},
{
"cell_type": "markdown",
"id": "b1a14244",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Quiz\n",
"\n",
"If the (true) mean mental health visit rates in high deprivation neigbourhoods is equal to the (true) mean mental health visit rates in low deprivation neigbourhoods. Then the **labels** of `\"depriv_binary\"` `High` and `Low` on neighbourhoods are ...\n",
"\n",
"A. interchangable (`High` can be changed for `Low` without effecting the mean mental health visit rates)\n",
"\n",
"B. not interchangable (Changing any lablel from `High` to `Low` affects the mean mental health visit rates)\n"
]
},
{
"cell_type": "markdown",
"id": "ba4c9cab-b103-46ac-bb42-bff179ec26e7",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"So, if \n",
"\n",
"> (The true) mean mental health visit rates in high deprivation neigbourhoods are equal to the (true) mean mental health visit rates in low deprivation neigbourhoods\n",
"\n",
"Then\n",
"\n",
"> The difference in mean mental health rates between `High` and `Low` deprivation neighbourhoods after switching any two `\"deperiv_binary\"` labels was a possible observation.\n",
"\n",
"Thus,\n",
"\n",
"> Re-calculating the difference in mean mental health rates between `High` and `Low` deprivation neighbourhoods for **all possible** combinations of the **labels** will provide the distribution of the possible differences."
]
},
{
"cell_type": "code",
"execution_count": 52,
"id": "6d21e588",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"5.94910575592826e+39"
]
},
"execution_count": 52,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# the number of possible combinations...\n",
"import math\n",
"math.factorial(136) / (\n",
" math.factorial(68) * math.factorial(68)\n",
")"
]
},
{
"cell_type": "markdown",
"id": "dcf90b52",
"metadata": {},
"source": [
"_It will take ages (of the universe) even with the fastest supercomputer._"
]
},
{
"cell_type": "markdown",
"id": "423918e0",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Approximate via simulation\n",
"\n",
"So, if \n",
"\n",
"> (The true) mean mental health visit rates in high deprivation neigbourhoods are equal to the (true) mean mental health visit rates in low deprivation neigbourhoods\n",
"\n",
"Then\n",
"\n",
"> The difference in mean mental health rates between `High` and `Low` deprivation neighbourhoods after switching any two `\"deperiv_binary\"` labels was a possible observation.\n",
"\n",
"Thus,\n",
"\n",
"> Re-calculating the difference in mean mental health rates between `High` and `Low` deprivation neighbourhoods for **a set of randomly shuffled labels** will provide the _approximaate_ distribution of the possible differences.\n"
]
},
{
"cell_type": "markdown",
"id": "60311806",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Random shuffling\n",
"\n",
"- We can randlomly shuffle using the `sample` function in `pandas`.\n",
"- The parameter `frac` in the `pandas` `sample` function refers to the fraction of rows to return. `frac=1` means all the rows are returned."
]
},
{
"cell_type": "code",
"execution_count": 53,
"id": "6848f37f-6c35-4994-aa0a-4c5043acb247",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 High\n",
"1 High\n",
"2 High\n",
"3 High\n",
"4 High\n",
"5 High\n",
"6 Low\n",
"7 High\n",
"8 Low\n",
"9 Low\n",
"Name: depriv_binary, dtype: object"
]
},
"execution_count": 53,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mh_visit_depriv.iloc[0:10, 4]"
]
},
{
"cell_type": "code",
"execution_count": 54,
"id": "bc77ed95",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"8 Low\n",
"5 High\n",
"0 High\n",
"2 High\n",
"1 High\n",
"9 Low\n",
"7 High\n",
"3 High\n",
"6 Low\n",
"4 High\n",
"Name: depriv_binary, dtype: object"
]
},
"execution_count": 54,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.random.seed(7) # for reproducability\n",
"\n",
"mh_visit_depriv.iloc[0:10, 4].sample(frac=1, replace=False)"
]
},
{
"cell_type": "markdown",
"id": "de341713",
"metadata": {},
"source": [
"- `reset_index` resets the index. \n",
"- the argument `drop=True` in `reset_index` indicates not to save the index after running `sample` as a column."
]
},
{
"cell_type": "code",
"execution_count": 55,
"id": "c4adb67d-31e7-41c8-9059-b7cf4da5fe98",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 Low\n",
"1 High\n",
"2 High\n",
"3 High\n",
"4 High\n",
"5 Low\n",
"6 High\n",
"7 High\n",
"8 Low\n",
"9 High\n",
"Name: depriv_binary, dtype: object"
]
},
"execution_count": 55,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.random.seed(7)\n",
"\n",
"mh_visit_depriv.iloc[0:10,4].sample(frac=1, replace=False).reset_index(drop=True)"
]
},
{
"cell_type": "markdown",
"id": "3db92599-56a9-4b80-a2d0-df97f548b63d",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"We can investigate what happens if the labels for `High` and `Low` are randomly shuffled for the first 10 rows. "
]
},
{
"cell_type": "code",
"execution_count": 56,
"id": "0f2511b5-0b93-4998-bafc-178d676dea54",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
mh_visit_rates_mf
\n",
"
depriv_binary
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
7.4
\n",
"
High
\n",
"
\n",
"
\n",
"
1
\n",
"
7.8
\n",
"
High
\n",
"
\n",
"
\n",
"
2
\n",
"
7.8
\n",
"
High
\n",
"
\n",
"
\n",
"
3
\n",
"
8.9
\n",
"
High
\n",
"
\n",
"
\n",
"
4
\n",
"
8.5
\n",
"
Low
\n",
"
\n",
"
\n",
"
5
\n",
"
7.6
\n",
"
High
\n",
"
\n",
"
\n",
"
6
\n",
"
7.8
\n",
"
High
\n",
"
\n",
"
\n",
"
7
\n",
"
8.9
\n",
"
Low
\n",
"
\n",
"
\n",
"
8
\n",
"
8.2
\n",
"
High
\n",
"
\n",
"
\n",
"
9
\n",
"
8.0
\n",
"
Low
\n",
"
\n",
"
\n",
"
10
\n",
"
7.8
\n",
"
Low
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" mh_visit_rates_mf depriv_binary\n",
"0 7.4 High\n",
"1 7.8 High\n",
"2 7.8 High\n",
"3 8.9 High\n",
"4 8.5 Low\n",
"5 7.6 High\n",
"6 7.8 High\n",
"7 8.9 Low\n",
"8 8.2 High\n",
"9 8.0 Low\n",
"10 7.8 Low"
]
},
"execution_count": 56,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# mental health visit rates for first 10 rows\n",
"\n",
"visits = mh_visit_depriv.loc[0:10, \"mh_visit_rates_mf\"]\n",
"\n",
"# randomly shuffle depriv_HL \n",
"# previous index is dropped because we don't need it\n",
"\n",
"shuffled_depriv = mh_visit_depriv.loc[0:10, \"depriv_binary\"].sample(\n",
" frac=1, replace=False).reset_index(drop=True)\n",
"\n",
"# put two dataframes in a list as input to\n",
"# pd.concat\n",
"\n",
"L = [visits, shuffled_depriv]\n",
"\n",
"# combine two columns\n",
"pd.concat(L, axis = 1)\n"
]
},
{
"cell_type": "markdown",
"id": "acbf5c74-0237-4ada-96cd-28724d635ac6",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Set up the simulation in python\n",
"\n",
"**Step 1:** Shuffle the column `\"depriv_binary\"`\n",
"\n",
"- To do this we will use `pandas` `sample` with the following parameters for the sample function:\n",
"\n",
" + `frac=1` (sample 100% of values)\n",
" + `replace=True` (we want the same number of `High` `Low` as original sample)\n",
" \n",
"- `reset_index(drop=True)` (use the dafault index - row labels in `pandas` so that we can assign shuffled labels to )"
]
},
{
"cell_type": "code",
"execution_count": 57,
"id": "9a2421bd-4015-453a-b9b1-4fbf1e8aef8e",
"metadata": {},
"outputs": [],
"source": [
"depriv_binary_shuffle = mh_visit_depriv[\"depriv_binary\"].sample(\n",
" frac=1, replace=False).reset_index(drop=True)"
]
},
{
"cell_type": "markdown",
"id": "ed59b00b-538b-4803-8170-5a924ae6bc32",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"**Step 2:** Assign the shuffled labels to the neighbourhoods and compute the mean rate of mental health visits"
]
},
{
"cell_type": "code",
"execution_count": 58,
"id": "dd0d0037-9c0d-434b-ae0d-7018289ef8a5",
"metadata": {},
"outputs": [],
"source": [
"visit_rate_low_shuffle = mh_visit_depriv.loc[depriv_binary_shuffle == \"High\", \"mh_visit_rates_mf\"].mean()\n",
"visit_rate_high_shuffle = mh_visit_depriv.loc[depriv_binary_shuffle == \"Low\", \"mh_visit_rates_mf\"].mean()"
]
},
{
"cell_type": "markdown",
"id": "9c881906-05e4-4e81-abd2-9bf49552cffb",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"**Step 3:** Compute the mean difference between the groups"
]
},
{
"cell_type": "code",
"execution_count": 59,
"id": "1591d57d-717a-4e67-8e1c-6c4cea8ff10a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-0.08823529411764675"
]
},
"execution_count": 59,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"visit_rate_high_shuffle - visit_rate_low_shuffle "
]
},
{
"cell_type": "markdown",
"id": "1b6e1bb3-5f8f-47fa-95f9-94b34a7ac4e8",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"- Steps 1 - 3 is an algorithm for computing one simulated value of the mean difference when the null hypothesis is true.\n",
"\n",
"- Let's create a function to do Steps 1 - 3 that returns a simulated value of the mean difference (a simulated value of the test statistic)."
]
},
{
"cell_type": "code",
"execution_count": 60,
"id": "6c49132b-0559-4632-b6ae-887f77eb0d07",
"metadata": {},
"outputs": [],
"source": [
"def random_shuffle_mean():\n",
" # step 1\n",
" depriv_binary_shuffle = mh_visit_depriv[\"depriv_binary\"].sample(\n",
" frac=1, replace=False).reset_index(drop=True)\n",
" \n",
" #step 2\n",
" visit_rate_low_shuffle = mh_visit_depriv.loc[\n",
" depriv_binary_shuffle == \"High\", \"mh_visit_rates_mf\"].mean()\n",
" visit_rate_high_shuffle = mh_visit_depriv.loc[\n",
" depriv_binary_shuffle == \"Low\", \"mh_visit_rates_mf\"].mean()\n",
" \n",
" #step 3\n",
" shuffled_diff = visit_rate_high_shuffle - visit_rate_low_shuffle \n",
" \n",
" return shuffled_diff "
]
},
{
"cell_type": "code",
"execution_count": 61,
"id": "4612d41d-46c9-4633-a8eb-257cada992c5",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-0.00588235294117645"
]
},
"execution_count": 61,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"random_shuffle_mean()"
]
},
{
"cell_type": "markdown",
"id": "0b0acd61-a822-45c7-a1fc-2a217ea55ef7",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Distribution of simulated values of the mean difference assuming the null hypothesis is true\n",
"\n",
"- Statistical tests based on shuffles or permutations of the data are called permutation tests. "
]
},
{
"cell_type": "code",
"execution_count": 62,
"id": "8dc984a6-12e7-44c6-8689-6e9f9670dae6",
"metadata": {},
"outputs": [],
"source": [
"shuffled_diffs = []\n",
"\n",
"for _ in range(5000):\n",
" shuffled_diffs.append(random_shuffle_mean())"
]
},
{
"cell_type": "code",
"execution_count": 63,
"id": "f0a79597-b5ed-499f-86d7-3d1ae2dda22a",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGdCAYAAADqsoKGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAwoUlEQVR4nO3df1yU9Z7//ycjAooOiMog56BZmb8yLS0c004pK/6oY0fOlsV6rGW1U1Crnpupe9RMK8pcNY1ybSttb3qs9qhrViZhiiWRYmz+CrPc9OQZoCVAdOXXXN8/+nJ9GsNOg4PAm8f9drtuN+f9fl/X9XqJOk+vmWsmyLIsSwAAAC2co6kLAAAACARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACMFNXUBj8Xq9On36tDp27KigoKCmLgcAAPwMlmXpzJkzio2NlcPh37UXY0PN6dOnFRcX19RlAACABjh16pR++ctf+rWPsaGmY8eOkr7/TXE6nU1cDQAA+DnKy8sVFxdnP4/7w9hQU/eSk9PpJNQAANDCNOStI7xRGAAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACM4Heoyc7O1h133KHY2FgFBQVpy5YtP1pz9OhR/frXv1ZERITCw8N144036uTJk/b8+fPnlZqaqs6dO6tDhw5KSkpSYWGhzzFOnjyp8ePHq3379oqOjtasWbNUU1Pjf4cAAKBV8DvUnD17VgMHDlRGRka9819++aWGDx+uPn36aNeuXfrss880f/58hYWF2WtmzJiht956S2+++aZ2796t06dPa+LEifZ8bW2txo8fr6qqKu3du1fr1q3T2rVrtWDBgga0CAAAWoMgy7KsBu8cFKTNmzfrzjvvtMcmTZqktm3b6j/+4z/q3aesrExdu3bVhg0b9Nvf/laS9Pnnn6tv377KycnR0KFD9e677+r222/X6dOn5XK5JEmrV6/W7NmzVVxcrJCQkL9ZW3l5uSIiIlRWVsYXWgIA0EJcyvN3QN9T4/V69fbbb+uaa65RYmKioqOjFR8f7/MSVV5enqqrq5WQkGCP9enTR927d1dOTo4kKScnRwMGDLADjSQlJiaqvLxchw8frvfclZWVKi8v99kAAEDrEdBQU1RUpIqKCj399NMaM2aMduzYod/85jeaOHGidu/eLUnyeDwKCQlRZGSkz74ul0sej8de88NAUzdfN1ef9PR0RURE2FtcXFwgW/t/amul3fu/32prG+ccAAA0V834eTDgV2okacKECZoxY4YGDRqkOXPm6Pbbb9fq1asDeaofmTt3rsrKyuzt1KlTjXo+AADQvAQ01HTp0kXBwcHq16+fz3jfvn3tu59iYmJUVVWl0tJSnzWFhYWKiYmx11x4N1Td47o1FwoNDZXT6fTZAABA6xHQUBMSEqIbb7xRBQUFPuPHjh1Tjx49JEmDBw9W27ZtlZWVZc8XFBTo5MmTcrvdkiS3262DBw+qqKjIXpOZmSmn0/mjwAQAACBJwf7uUFFRoePHj9uPT5w4ofz8fEVFRal79+6aNWuW7r77bt1yyy267bbbtH37dr311lvatWuXJCkiIkIpKSmaOXOmoqKi5HQ69fDDD8vtdmvo0KGSpNGjR6tfv36aPHmylixZIo/Ho3nz5ik1NVWhoaGB6RwAABjF71Czf/9+3XbbbfbjmTNnSpKmTJmitWvX6je/+Y1Wr16t9PR0PfLII+rdu7f+/Oc/a/jw4fY+y5cvl8PhUFJSkiorK5WYmKgXXnjBnm/Tpo22bdumBx98UG63W+Hh4ZoyZYoWLVp0Kb0CAACDXdLn1DRnjfY5NbW10oeffv/r4ddLbdoE7tgAADR3jfw82Gw+pwYAAKCpEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACP4HWqys7N1xx13KDY2VkFBQdqyZctF1/7+979XUFCQVqxY4TNeUlKi5ORkOZ1ORUZGKiUlRRUVFT5rPvvsM40YMUJhYWGKi4vTkiVL/C0VAAC0In6HmrNnz2rgwIHKyMj4yXWbN2/Wxx9/rNjY2B/NJScn6/Dhw8rMzNS2bduUnZ2tadOm2fPl5eUaPXq0evTooby8PD377LNauHCh1qxZ42+5AACglQj2d4exY8dq7NixP7nmm2++0cMPP6z33ntP48eP95k7evSotm/frn379mnIkCGSpFWrVmncuHFaunSpYmNjtX79elVVVemVV15RSEiI+vfvr/z8fC1btswn/AAAANQJ+HtqvF6vJk+erFmzZql///4/ms/JyVFkZKQdaCQpISFBDodDubm59ppbbrlFISEh9prExEQVFBTou+++q/e8lZWVKi8v99kAAEDrEfBQ88wzzyg4OFiPPPJIvfMej0fR0dE+Y8HBwYqKipLH47HXuFwunzV1j+vWXCg9PV0RERH2FhcXd6mtAACAFiSgoSYvL0/PPfec1q5dq6CgoEAe+m+aO3euysrK7O3UqVOX9fwAAKBpBTTU7NmzR0VFRerevbuCg4MVHBysr7/+Wn/4wx90xRVXSJJiYmJUVFTks19NTY1KSkoUExNjryksLPRZU/e4bs2FQkND5XQ6fTYAANB6BDTUTJ48WZ999pny8/PtLTY2VrNmzdJ7770nSXK73SotLVVeXp69386dO+X1ehUfH2+vyc7OVnV1tb0mMzNTvXv3VqdOnQJZMgAAMITfdz9VVFTo+PHj9uMTJ04oPz9fUVFR6t69uzp37uyzvm3btoqJiVHv3r0lSX379tWYMWM0depUrV69WtXV1UpLS9OkSZPs27/vvfdePf7440pJSdHs2bN16NAhPffcc1q+fPml9AoAAAzmd6jZv3+/brvtNvvxzJkzJUlTpkzR2rVrf9Yx1q9fr7S0NI0aNUoOh0NJSUlauXKlPR8REaEdO3YoNTVVgwcPVpcuXbRgwQJu5wYAABcVZFmW1dRFNIby8nJFRESorKwssO+vqa2VPvz0+18Pv15q0yZwxwYAoLlr5OfBS3n+5rufAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAJfM6/U2dQm25lQLgMsruKkLANDyORwObdq0ScXFxU1aR9euXTVx4sQmrQFA0yHUAAiI4uJieTyepi4DQCvGy08AAMAIhBoAAGAEQg0AADACoQYAABjB71CTnZ2tO+64Q7GxsQoKCtKWLVvsuerqas2ePVsDBgxQeHi4YmNj9bvf/U6nT5/2OUZJSYmSk5PldDoVGRmplJQUVVRU+Kz57LPPNGLECIWFhSkuLk5LlixpWIcAAKBV8DvUnD17VgMHDlRGRsaP5s6dO6cDBw5o/vz5OnDggDZt2qSCggL9+te/9lmXnJysw4cPKzMzU9u2bVN2dramTZtmz5eXl2v06NHq0aOH8vLy9Oyzz2rhwoVas2ZNA1oEAACtgd+3dI8dO1Zjx46tdy4iIkKZmZk+Y88//7xuuukmnTx5Ut27d9fRo0e1fft27du3T0OGDJEkrVq1SuPGjdPSpUsVGxur9evXq6qqSq+88opCQkLUv39/5efna9myZT7hBwAAoE6jv6emrKxMQUFBioyMlCTl5OQoMjLSDjSSlJCQIIfDodzcXHvNLbfcopCQEHtNYmKiCgoK9N1339V7nsrKSpWXl/tsAACg9WjUUHP+/HnNnj1b99xzj5xOpyTJ4/EoOjraZ11wcLCioqLsD+7yeDxyuVw+a+oeX+zDvdLT0xUREWFvcXFxgW4HAAA0Y40Waqqrq3XXXXfJsiy9+OKLjXUa29y5c1VWVmZvp06davRzAgCA5qNRviahLtB8/fXX2rlzp32VRpJiYmJUVFTks76mpkYlJSWKiYmx1xQWFvqsqXtct+ZCoaGhCg0NDWQbAACgBQn4lZq6QPPFF1/o/fffV+fOnX3m3W63SktLlZeXZ4/t3LlTXq9X8fHx9prs7GxVV1fbazIzM9W7d2916tQp0CUDMESHDh2a1bd0N6dagNbA7ys1FRUVOn78uP34xIkTys/PV1RUlLp166bf/va3OnDggLZt26ba2lr7PTBRUVEKCQlR3759NWbMGE2dOlWrV69WdXW10tLSNGnSJMXGxkqS7r33Xj3++ONKSUnR7NmzdejQIT333HNavnx5gNoGYKKwsDC+MRxoxfwONfv379dtt91mP545c6YkacqUKVq4cKG2bt0qSRo0aJDPfh988IFuvfVWSdL69euVlpamUaNGyeFwKCkpSStXrrTXRkREaMeOHUpNTdXgwYPVpUsXLViwgNu5AfwsfGM40Dr5HWpuvfVWWZZ10fmfmqsTFRWlDRs2/OSa6667Tnv27PG3PAAA0Erx3U8AAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAh+h5rs7Gzdcccdio2NVVBQkLZs2eIzb1mWFixYoG7duqldu3ZKSEjQF1984bOmpKREycnJcjqdioyMVEpKiioqKnzWfPbZZxoxYoTCwsIUFxenJUuW+N8dAABoNfwONWfPntXAgQOVkZFR7/ySJUu0cuVKrV69Wrm5uQoPD1diYqLOnz9vr0lOTtbhw4eVmZmpbdu2KTs7W9OmTbPny8vLNXr0aPXo0UN5eXl69tlntXDhQq1Zs6YBLQIAgNYg2N8dxo4dq7Fjx9Y7Z1mWVqxYoXnz5mnChAmSpNdee00ul0tbtmzRpEmTdPToUW3fvl379u3TkCFDJEmrVq3SuHHjtHTpUsXGxmr9+vWqqqrSK6+8opCQEPXv31/5+flatmyZT/gBAACoE9D31Jw4cUIej0cJCQn2WEREhOLj45WTkyNJysnJUWRkpB1oJCkhIUEOh0O5ubn2mltuuUUhISH2msTERBUUFOi7776r99yVlZUqLy/32QAAQOsR0FDj8XgkSS6Xy2fc5XLZcx6PR9HR0T7zwcHBioqK8llT3zF+eI4LpaenKyIiwt7i4uIuvSEAANBiGHP309y5c1VWVmZvp06dauqSAADAZRTQUBMTEyNJKiws9BkvLCy052JiYlRUVOQzX1NTo5KSEp819R3jh+e4UGhoqJxOp88GAABaj4CGmp49eyomJkZZWVn2WHl5uXJzc+V2uyVJbrdbpaWlysvLs9fs3LlTXq9X8fHx9prs7GxVV1fbazIzM9W7d2916tQpkCUDAABD+B1qKioqlJ+fr/z8fEnfvzk4Pz9fJ0+eVFBQkKZPn64nnnhCW7du1cGDB/W73/1OsbGxuvPOOyVJffv21ZgxYzR16lR98skn+uijj5SWlqZJkyYpNjZWknTvvfcqJCREKSkpOnz4sF5//XU999xzmjlzZsAaBwAAZvH7lu79+/frtttusx/XBY0pU6Zo7dq1evTRR3X27FlNmzZNpaWlGj58uLZv366wsDB7n/Xr1ystLU2jRo2Sw+FQUlKSVq5cac9HRERox44dSk1N1eDBg9WlSxctWLCA27kBAMBF+R1qbr31VlmWddH5oKAgLVq0SIsWLbromqioKG3YsOEnz3Pddddpz549/pYHAABaKWPufgIAAK0boQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAVowr9fb1CUAQLPh9xdaAmg+HA6HNm3apOLi4iaroVevXho5cmSTnR8A6hBqgBauuLhYHo+nyc7fpUuXJjs3APwQLz8BAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACAEPNbW1tZo/f7569uypdu3a6aqrrtLixYtlWZa9xrIsLViwQN26dVO7du2UkJCgL774wuc4JSUlSk5OltPpVGRkpFJSUlRRURHocgEAgCECHmqeeeYZvfjii3r++ed19OhRPfPMM1qyZIlWrVplr1myZIlWrlyp1atXKzc3V+Hh4UpMTNT58+ftNcnJyTp8+LAyMzO1bds2ZWdna9q0aYEuFwAAGCI40Afcu3evJkyYoPHjx0uSrrjiCv3pT3/SJ598Iun7qzQrVqzQvHnzNGHCBEnSa6+9JpfLpS1btmjSpEk6evSotm/frn379mnIkCGSpFWrVmncuHFaunSpYmNjA102AABo4QJ+pWbYsGHKysrSsWPHJEn//d//rQ8//FBjx46VJJ04cUIej0cJCQn2PhEREYqPj1dOTo4kKScnR5GRkXagkaSEhAQ5HA7l5ubWe97KykqVl5f7bAAAoPUI+JWaOXPmqLy8XH369FGbNm1UW1urJ598UsnJyZIkj8cjSXK5XD77uVwue87j8Sg6Otq30OBgRUVF2WsulJ6erscffzzQ7QAAgBYi4Fdq3njjDa1fv14bNmzQgQMHtG7dOi1dulTr1q0L9Kl8zJ07V2VlZfZ26tSpRj0fAABoXgJ+pWbWrFmaM2eOJk2aJEkaMGCAvv76a6Wnp2vKlCmKiYmRJBUWFqpbt272foWFhRo0aJAkKSYmRkVFRT7HrampUUlJib3/hUJDQxUaGhrodgAAQAsR8Cs1586dk8Phe9g2bdrI6/VKknr27KmYmBhlZWXZ8+Xl5crNzZXb7ZYkud1ulZaWKi8vz16zc+dOeb1excfHB7pkAABggIBfqbnjjjv05JNPqnv37urfv78+/fRTLVu2TP/4j/8oSQoKCtL06dP1xBNPqFevXurZs6fmz5+v2NhY3XnnnZKkvn37asyYMZo6dapWr16t6upqpaWladKkSdz5BAAA6hXwULNq1SrNnz9fDz30kIqKihQbG6sHHnhACxYssNc8+uijOnv2rKZNm6bS0lINHz5c27dvV1hYmL1m/fr1SktL06hRo+RwOJSUlKSVK1cGulwAAGCIgIeajh07asWKFVqxYsVF1wQFBWnRokVatGjRRddERUVpw4YNgS4PAAAYiu9+AgAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAGkGHDh3k9Xqbugxbc6oFaCzBTV0AAJgoLCxMDodDmzZtUnFxcZPW0rVrV02cOLFJawAuB0INADSi4uJieTyepi4DaBV4+QkAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwQqOEmm+++Ub/8A//oM6dO6tdu3YaMGCA9u/fb89blqUFCxaoW7duateunRISEvTFF1/4HKOkpETJyclyOp2KjIxUSkqKKioqGqNcAABggICHmu+++04333yz2rZtq3fffVdHjhzRv/7rv6pTp072miVLlmjlypVavXq1cnNzFR4ersTERJ0/f95ek5ycrMOHDyszM1Pbtm1Tdna2pk2bFuhyAQCAIYIDfcBnnnlGcXFxevXVV+2xnj172r+2LEsrVqzQvHnzNGHCBEnSa6+9JpfLpS1btmjSpEk6evSotm/frn379mnIkCGSpFWrVmncuHFaunSpYmNjA102AABo4QJ+pWbr1q0aMmSI/v7v/17R0dG6/vrr9dJLL9nzJ06ckMfjUUJCgj0WERGh+Ph45eTkSJJycnIUGRlpBxpJSkhIkMPhUG5ubr3nraysVHl5uc8GAABaj4CHmq+++kovvviievXqpffee08PPvigHnnkEa1bt06S5PF4JEkul8tnP5fLZc95PB5FR0f7zAcHBysqKspec6H09HRFRETYW1xcXKBbAwAAzVjAQ43X69UNN9ygp556Stdff72mTZumqVOnavXq1YE+lY+5c+eqrKzM3k6dOtWo5wMAAM1LwENNt27d1K9fP5+xvn376uTJk5KkmJgYSVJhYaHPmsLCQnsuJiZGRUVFPvM1NTUqKSmx11woNDRUTqfTZwMAAK1HwEPNzTffrIKCAp+xY8eOqUePHpK+f9NwTEyMsrKy7Pny8nLl5ubK7XZLktxut0pLS5WXl2ev2blzp7xer+Lj4wNdMgAAMEDA736aMWOGhg0bpqeeekp33XWXPvnkE61Zs0Zr1qyRJAUFBWn69Ol64okn1KtXL/Xs2VPz589XbGys7rzzTknfX9kZM2aM/bJVdXW10tLSNGnSJO58AgAA9Qp4qLnxxhu1efNmzZ07V4sWLVLPnj21YsUKJScn22seffRRnT17VtOmTVNpaamGDx+u7du3KywszF6zfv16paWladSoUXI4HEpKStLKlSsDXS4AADBEwEONJN1+++26/fbbLzofFBSkRYsWadGiRRddExUVpQ0bNjRGeQAAwEB89xMAADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMEKjh5qnn35aQUFBmj59uj12/vx5paamqnPnzurQoYOSkpJUWFjos9/Jkyc1fvx4tW/fXtHR0Zo1a5Zqamoau1wAANBCNWqo2bdvn/7t3/5N1113nc/4jBkz9NZbb+nNN9/U7t27dfr0aU2cONGer62t1fjx41VVVaW9e/dq3bp1Wrt2rRYsWNCY5QI/i9frbeoSAAD1CG6sA1dUVCg5OVkvvfSSnnjiCXu8rKxML7/8sjZs2KCRI0dKkl599VX17dtXH3/8sYYOHaodO3boyJEjev/99+VyuTRo0CAtXrxYs2fP1sKFCxUSEtJYZQN/k8Ph0KZNm1RcXNykdfTq1cv+OwQAaMRQk5qaqvHjxyshIcEn1OTl5am6uloJCQn2WJ8+fdS9e3fl5ORo6NChysnJ0YABA+Ryuew1iYmJevDBB3X48GFdf/31PzpfZWWlKisr7cfl5eWN1BkgFRcXy+PxNGkNXbp0adLzA0Bz0yihZuPGjTpw4ID27dv3ozmPx6OQkBBFRkb6jLtcLvtJwuPx+ASauvm6ufqkp6fr8ccfD0D1AACgJQr4e2pOnTqlf/7nf9b69esVFhYW6MNf1Ny5c1VWVmZvp06dumznBgAATS/goSYvL09FRUW64YYbFBwcrODgYO3evVsrV65UcHCwXC6XqqqqVFpa6rNfYWGhYmJiJEkxMTE/uhuq7nHdmguFhobK6XT6bAAAoPUIeKgZNWqUDh48qPz8fHsbMmSIkpOT7V+3bdtWWVlZ9j4FBQU6efKk3G63JMntduvgwYMqKiqy12RmZsrpdKpfv36BLhkAABgg4O+p6dixo6699lqfsfDwcHXu3NkeT0lJ0cyZMxUVFSWn06mHH35YbrdbQ4cOlSSNHj1a/fr10+TJk7VkyRJ5PB7NmzdPqampCg0NDXTJAADAAI1299NPWb58uRwOh5KSklRZWanExES98MIL9nybNm20bds2Pfjgg3K73QoPD9eUKVO0aNGipigXAAC0AJcl1OzatcvncVhYmDIyMpSRkXHRfXr06KF33nmnkSsDAACm4LufAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAM16FDB3m93qYuw9acaoFZgpu6AABA4woLC5PD4dCmTZtUXFzcpLV07dpVEydObNIaYC5CDQC0EsXFxfJ4PE1dBtBoePkJAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIwQ8FCTnp6uG2+8UR07dlR0dLTuvPNOFRQU+Kw5f/68UlNT1blzZ3Xo0EFJSUkqLCz0WXPy5EmNHz9e7du3V3R0tGbNmqWamppAlwsAAAwR8FCze/dupaam6uOPP1ZmZqaqq6s1evRonT171l4zY8YMvfXWW3rzzTe1e/dunT592ucLzmprazV+/HhVVVVp7969WrdundauXasFCxYEulwAAGCIgH+h5fbt230er127VtHR0crLy9Mtt9yisrIyvfzyy9qwYYNGjhwpSXr11VfVt29fffzxxxo6dKh27NihI0eO6P3335fL5dKgQYO0ePFizZ49WwsXLlRISEigywYAAC1co7+npqysTJIUFRUlScrLy1N1dbUSEhLsNX369FH37t2Vk5MjScrJydGAAQPkcrnsNYmJiSovL9fhw4frPU9lZaXKy8t9NgAA0Ho0aqjxer2aPn26br75Zl177bWSJI/Ho5CQEEVGRvqsdblc8ng89pofBpq6+bq5+qSnpysiIsLe4uLiAtwNAABozho11KSmpurQoUPauHFjY55GkjR37lyVlZXZ26lTpxr9nAAAoPkI+Htq6qSlpWnbtm3Kzs7WL3/5S3s8JiZGVVVVKi0t9blaU1hYqJiYGHvNJ5984nO8uruj6tZcKDQ0VKGhoQHuAgAAtBQBv1JjWZbS0tK0efNm7dy5Uz179vSZHzx4sNq2bausrCx7rKCgQCdPnpTb7ZYkud1uHTx4UEVFRfaazMxMOZ1O9evXL9AlAwAAAwT8Sk1qaqo2bNig//qv/1LHjh3t98BERESoXbt2ioiIUEpKimbOnKmoqCg5nU49/PDDcrvdGjp0qCRp9OjR6tevnyZPnqwlS5bI4/Fo3rx5Sk1N5WoMAACoV8BDzYsvvihJuvXWW33GX331Vd13332SpOXLl8vhcCgpKUmVlZVKTEzUCy+8YK9t06aNtm3bpgcffFBut1vh4eGaMmWKFi1aFOhyAQCAIQIeaizL+ptrwsLClJGRoYyMjIuu6dGjh955551AlgYAAAzGdz8BAAAjEGoAAIARCDUAAMAIhBoAAGAEQg0AADACoQYAABiBUAMAAIxAqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGrQIXq+3qUsAADRzAf9CS6AxOBwObdq0ScXFxU1aR69evTRy5MgmrQEAUD9CDVqM4uJieTyeJq2hS5cuTXp+AMDF8fITAAAwAqEGAAAYgVADAACMQKgBAABGINQAAAAjEGoAAIARCDUAgMumQ4cOzebDNJtLHQgcPqcGAHDZhIWFNYsP0+zatasmTpzYZOdH4yDUAAAuu+bwYZowDy8/AQAAIxBqAACAEQg1AADACIQaAABgBEINAAAwAqEGAAAYgVCDi+KDqQAALQmfU4OLag4fkCVJvXr10siRI5u0BgBA80eowU9qDh+Q1aVLlyY9PwCgZeDlJwAAYIRmHWoyMjJ0xRVXKCwsTPHx8frkk0+auiQAANBMNdtQ8/rrr2vmzJl67LHHdODAAQ0cOFCJiYkqKipq6tIAAC1cc/q2cIkbMwKl2b6nZtmyZZo6daruv/9+SdLq1av19ttv65VXXtGcOXOauDoAQEvWXL4tXOIbwwOpWYaaqqoq5eXlae7cufaYw+FQQkKCcnJy6t2nsrJSlZWV9uOysjJJUnl5eWCLq62Vzlbo/z+41KZNYI/fzISHhysiIqJJawgJCVF5eTm1NNNamksd1EItDanj3Llz+r//+78mq0OSzp07F/jnqsbUyM+Ddb8XlmX5vW+zDDXffvutamtr5XK5fMZdLpc+//zzevdJT0/X448//qPxuLi4RqkRAAA0njNnzvgdfJtlqGmIuXPnaubMmfZjr9erkpISde7cWUFBQT5ry8vLFRcXp1OnTsnpdF7uUi87+jVfa+uZfs3X2nqm3//HsiydOXNGsbGxfh+3WYaaLl26qE2bNiosLPQZLywsVExMTL37hIaGKjQ01GcsMjLyJ8/jdDpbxR+eOvRrvtbWM/2ar7X1TL/fa+hLk83y7qeQkBANHjxYWVlZ9pjX61VWVpbcbncTVgYAAJqrZnmlRpJmzpypKVOmaMiQIbrpppu0YsUKnT171r4bCgAA4Ieabai5++67VVxcrAULFsjj8WjQoEHavn37j9483BChoaF67LHHfvRylano13ytrWf6NV9r65l+AyPIasg9UwAAAM1Ms3xPDQAAgL8INQAAwAiEGgAAYARCDQAAMEKrCDUlJSVKTk6W0+lUZGSkUlJSVFFR8ZP7eDweTZ48WTExMQoPD9cNN9ygP//5z5ep4kvXkJ4lKScnRyNHjlR4eLicTqduueWWJv9elJ+jof1K33965dixYxUUFKQtW7Y0bqEB4m+/JSUlevjhh9W7d2+1a9dO3bt31yOPPGJ/R1pzlJGRoSuuuEJhYWGKj4/XJ5988pPr33zzTfXp00dhYWEaMGCA3nnnnctUaWD40+9LL72kESNGqFOnTurUqZMSEhL+5u9Pc+Tvz7jOxo0bFRQUpDvvvLNxCwwwf/stLS1VamqqunXrptDQUF1zzTUt6s+1v/2uWLHC/jcqLi5OM2bM0Pnz5/07qdUKjBkzxho4cKD18ccfW3v27LGuvvpq65577vnJff7u7/7OuvHGG63c3Fzryy+/tBYvXmw5HA7rwIEDl6nqS9OQnvfu3Ws5nU4rPT3dOnTokPX5559br7/+unX+/PnLVHXDNaTfOsuWLbPGjh1rSbI2b97cuIUGiL/9Hjx40Jo4caK1detW6/jx41ZWVpbVq1cvKykp6TJW/fNt3LjRCgkJsV555RXr8OHD1tSpU63IyEirsLCw3vUfffSR1aZNG2vJkiXWkSNHrHnz5llt27a1Dh48eJkrbxh/+7333nutjIwM69NPP7WOHj1q3XfffVZERIT1l7/85TJX3nD+9lznxIkT1i9+8QtrxIgR1oQJEy5PsQHgb7+VlZXWkCFDrHHjxlkffvihdeLECWvXrl1Wfn7+Za68Yfztd/369VZoaKi1fv1668SJE9Z7771ndevWzZoxY4Zf5zU+1Bw5csSSZO3bt88ee/fdd62goCDrm2++ueh+4eHh1muvveYzFhUVZb300kuNVmugNLTn+Ph4a968eZejxIBqaL+WZVmffvqp9Ytf/ML661//2mJCzaX0+0NvvPGGFRISYlVXVzdGmZfkpptuslJTU+3HtbW1VmxsrJWenl7v+rvuussaP368z1h8fLz1wAMPNGqdgeJvvxeqqamxOnbsaK1bt66xSgy4hvRcU1NjDRs2zPr3f/93a8qUKS0q1Pjb74svvmhdeeWVVlVV1eUqMaD87Tc1NdUaOXKkz9jMmTOtm2++2a/zGv/yU05OjiIjIzVkyBB7LCEhQQ6HQ7m5uRfdb9iwYXr99ddVUlIir9erjRs36vz587r11lsvQ9WXpiE9FxUVKTc3V9HR0Ro2bJhcLpd+9atf6cMPP7xcZTdYQ3/G586d07333quMjIyLfqdYc9TQfi9UVlYmp9Op4ODm9RmcVVVVysvLU0JCgj3mcDiUkJCgnJycevfJycnxWS9JiYmJF13fnDSk3wudO3dO1dXVioqKaqwyA6qhPS9atEjR0dFKSUm5HGUGTEP63bp1q9xut1JTU+VyuXTttdfqqaeeUm1t7eUqu8Ea0u+wYcOUl5dnv0T11Vdf6Z133tG4ceP8Onfz+tesEXg8HkVHR/uMBQcHKyoqSh6P56L7vfHGG7r77rvVuXNnBQcHq3379tq8ebOuvvrqxi75kjWk56+++kqStHDhQi1dulSDBg3Sa6+9plGjRunQoUPq1atXo9fdUA39Gc+YMUPDhg3ThAkTGrvEgGpovz/07bffavHixZo2bVpjlHhJvv32W9XW1v7o08NdLpc+//zzevfxeDz1rv+5vx9NqSH9Xmj27NmKjY39UbBrrhrS84cffqiXX35Z+fn5l6HCwGpIv1999ZV27typ5ORkvfPOOzp+/LgeeughVVdX67HHHrscZTdYQ/q999579e2332r48OGyLEs1NTX6/e9/r3/5l3/x69wt9krNnDlzFBQU9JPbz/0HoT7z589XaWmp3n//fe3fv18zZ87UXXfdpYMHDwawC/80Zs9er1eS9MADD+j+++/X9ddfr+XLl6t379565ZVXAtnGz9aY/W7dulU7d+7UihUrAlv0JWjsP9N1ysvLNX78ePXr108LFy689MLRpJ5++mlt3LhRmzdvVlhYWFOX0yjOnDmjyZMn66WXXlKXLl2aupzLwuv1Kjo6WmvWrNHgwYN19913649//KNWr17d1KU1il27dumpp57SCy+8oAMHDmjTpk16++23tXjxYr+O02Kv1PzhD3/Qfffd95NrrrzySsXExKioqMhnvKamRiUlJRd9yeHLL7/U888/r0OHDql///6SpIEDB2rPnj3KyMhosj9Ujdlzt27dJEn9+vXzGe/bt69OnjzZ8KIvQWP2u3PnTn355ZeKjIz0GU9KStKIESO0a9euS6i8YRqz3zpnzpzRmDFj1LFjR23evFlt27a91LIDrkuXLmrTpo0KCwt9xgsLCy/aX0xMjF/rm5OG9Ftn6dKlevrpp/X+++/ruuuua8wyA8rfnr/88kv9z//8j+644w57rO4/YsHBwSooKNBVV13VuEVfgob8jLt166a2bduqTZs29ljfvn3l8XhUVVWlkJCQRq35UjSk3/nz52vy5Mn6p3/6J0nSgAEDdPbsWU2bNk1//OMf5XD8vGswLTbUdO3aVV27dv2b69xut0pLS5WXl6fBgwdL+v4Jzev1Kj4+vt59zp07J0k/+k1s06aN/RepKTRmz1dccYViY2NVUFDgM37s2DGNHTv20otvgMbsd86cOfZfnjoDBgzQ8uXLff7hvJwas1/p+ys0iYmJCg0N1datW5vt/+pDQkI0ePBgZWVl2bfser1eZWVlKS0trd593G63srKyNH36dHssMzNTbrf7MlR8aRrSryQtWbJETz75pN577z2f91e1BP723KdPnx9dJZ83b57OnDmj5557TnFxcZej7AZryM/45ptv1oYNG+T1eu3nomPHjqlbt27NOtBIDev33Llz9T7nSt9/7MbP5u87mluiMWPGWNdff72Vm5trffjhh1avXr18bn/9y1/+YvXu3dvKzc21LMuyqqqqrKuvvtoaMWKElZubax0/ftxaunSpFRQUZL399ttN1YZf/O3Zsixr+fLlltPptN58803riy++sObNm2eFhYVZx48fb4oW/NKQfi+kFnL3k2X5329ZWZkVHx9vDRgwwDp+/Lj117/+1d5qamqaqo2L2rhxoxUaGmqtXbvWOnLkiDVt2jQrMjLS8ng8lmVZ1uTJk605c+bY6z/66CMrODjYWrp0qXX06FHrsccea3G3dPvT79NPP22FhIRY//mf/+nzszxz5kxTteA3f3u+UEu7+8nffk+ePGl17NjRSktLswoKCqxt27ZZ0dHR1hNPPNFULfjF334fe+wxq2PHjtaf/vQn66uvvrJ27NhhXXXVVdZdd93l13lbRaj53//9X+uee+6xOnToYDmdTuv+++/3+ct/4sQJS5L1wQcf2GPHjh2zJk6caEVHR1vt27e3rrvuuh/d4t2cNaRny7Ks9PR065e//KXVvn17y+12W3v27LnMlTdMQ/v9oZYUavzt94MPPrAk1budOHGiaZr4G1atWmV1797dCgkJsW666Sbr448/tud+9atfWVOmTPFZ/8Ybb1jXXHONFRISYvXv37/F/Aekjj/99ujRo96f5WOPPXb5C78E/v6Mf6ilhRrL8r/fvXv3WvHx8VZoaKh15ZVXWk8++WSz/E/IxfjTb3V1tbVw4ULrqquussLCwqy4uDjroYcesr777ju/zhlkWf5c1wEAAGieWuzdTwAAAD9EqAEAAEYg1AAAACMQagAAgBEINQAAwAiEGgAAYARCDQAAMAKhBgAAGIFQAwAAjECoAQAARiDUAAAAIxBqAACAEf4/VoKFP6JR5XwAAAAASUVORK5CYII=",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"plt.hist(shuffled_diffs, color=\"grey\", edgecolor=\"white\")\n",
"\n",
"plt.vlines(x=observed_mean_difference, ymin=0, ymax=1600, color=\"pink\")\n",
"\n",
"plt.vlines(x=-1 * observed_mean_difference, ymin=0, ymax=1600, color=\"pink\");"
]
},
{
"cell_type": "markdown",
"id": "e675b176-7eea-46ef-abc6-4ea6d12b2226",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## p-value\n",
"\n",
"- Assuming that the null hypothesis is true, the p-value gives a measure of the probability of getting data that are *at least* as unusual or extreme as the sample data.\n",
"\n",
"- What does \"at least as unusual\" mean?\n",
"\n",
"- An unusual value would be a mean difference (`High` minus `Low`) greater than what we observed or a mean difference (`Low` minus `High`) less than what we observed.\n",
"\n",
"- Values that are larger than the observed value mean difference (the test statistic) 0.75 (`High` minus `Low`) or smaller than - 0.75 (`Low` minus `High`)."
]
},
{
"cell_type": "markdown",
"id": "27462f52-dc20-4723-8310-39e1bb6d0abd",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"The number of simulated mean differences that are greater than the observed mean difference can be computed using:"
]
},
{
"cell_type": "code",
"execution_count": 64,
"id": "b9931f03-9f89-4dda-9e85-4222ea35e5fd",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0"
]
},
"execution_count": 64,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"right_extreme = shuffled_diffs >= observed_mean_difference\n",
"\n",
"right_extreme.sum()"
]
},
{
"cell_type": "markdown",
"id": "78641883-0bb9-4e65-ba63-1d40e3b6c88d",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"The number of simulated mean differences that are less than the -1 times observed mean difference can be computed using:"
]
},
{
"cell_type": "code",
"execution_count": 65,
"id": "fa253024-283a-49de-a895-29c2fe064240",
"metadata": {
"tags": []
},
"outputs": [
{
"data": {
"text/plain": [
"0"
]
},
"execution_count": 65,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"left_extreme = shuffled_diffs <= -1 * observed_mean_difference\n",
"\n",
"left_extreme.sum()"
]
},
{
"cell_type": "markdown",
"id": "82bae078-3e00-4023-bfd0-683af3c5d820",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"The proportion of simulations that are unusual or extreme are:"
]
},
{
"cell_type": "code",
"execution_count": 66,
"id": "01aa3a37-1b8e-4c1d-a96f-cf1f92976b61",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0"
]
},
"execution_count": 66,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(right_extreme.sum() + left_extreme.sum()) / 5000"
]
},
{
"cell_type": "markdown",
"id": "a76cffb3-ef04-4643-9710-781ed34f0d4d",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"We can also write a function to carry out this computation."
]
},
{
"cell_type": "code",
"execution_count": 67,
"id": "1c63e16e-3a87-4feb-8b03-d5019759b713",
"metadata": {},
"outputs": [],
"source": [
"def pvalue_two_sided(shuff_diffs, obs_diffs):\n",
" \n",
" right_extreme_count = (shuff_diffs >= obs_diffs).sum()\n",
" \n",
" left_extreme_count = (shuff_diffs < -1 * obs_diffs).sum()\n",
" \n",
" all_extreme = right_extreme_count + left_extreme_count\n",
" \n",
" pval = all_extreme / len(shuff_diffs)\n",
" \n",
" return pval "
]
},
{
"cell_type": "code",
"execution_count": 68,
"id": "3c7e5631-76a3-4f32-9418-5e4e785cfbfd",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0"
]
},
"execution_count": 68,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pvalue_two_sided(shuff_diffs=shuffled_diffs, obs_diffs=observed_mean_difference)"
]
},
{
"cell_type": "markdown",
"id": "cecd15b0-2bb6-4023-8475-87267138f739",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"- This is called the **P-value**. \n",
"\n",
"- The entire procedure is often referred to as **significance testing**.\n",
"\n",
"- None of the permuted 5000 samples resulted in a mean difference more extreme than what we observed."
]
},
{
"cell_type": "markdown",
"id": "ff7ec95d-3774-43bf-bbab-ceb0a5b9a413",
"metadata": {
"slideshow": {
"slide_type": "subslide"
},
"tags": []
},
"source": [
"## Conclusions\n",
"\n",
"So, there are two possibilities:\n",
"\n",
"1. The null hypothesis of no difference is true, but we have observed a **very rare** value of the mean difference in our study.\n",
"\n",
"2. The null hypothesis is **false**, and there is a difference between high and low deprivation.\n",
"\n",
"It's scientific convention to assume that a small p-value is evidence in favour of 2 \n",
"(i.e., the null hypothesis is false, and there is evidence of a difference in the means)."
]
},
{
"cell_type": "markdown",
"id": "0035ac2f",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## History of Significance Testing\n",
"\n",
"R.A. Fisher made fundamental contributions to statistics including developing significance testing.\n",
"\n",
"> \"_He was from an early age a supporter of certain eugenic ideas, and it is for this reason that he has been accused of being a racist and an advocate of forced sterilisation (Evans 2020). His promotion of eugenics has recently caused various organisations to remove his name from awards and dedications of buildings._\" ([Bodmer et al, 2021](https://www.nature.com/articles/s41437-020-00394-6))"
]
},
{
"cell_type": "markdown",
"id": "d33f18f0-ca00-4489-aad5-9fad858a9b4b",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Causality\n",
"\n",
"Imagine...\n",
"\n",
"
Now, you go back in time to 10:00 and you don't take an Aspirin.
\n",
"
Your pain goes away after 48 minutes.
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "c3ca77ca",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"The **causal effect** of taking an Aspirin is 18 minutes (48 minutes - 30 minutes).\n"
]
},
{
"cell_type": "markdown",
"id": "a07027d9",
"metadata": {},
"source": [
""
]
},
{
"cell_type": "markdown",
"id": "7859b66e",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Potential outcomes and randomized control trials\n",
"\n",
"- Establishing causality involves comparing these **potential outcomes**. "
]
},
{
"cell_type": "markdown",
"id": "b7a9922e",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- The problem is that we can never observe both taking an Aspirin and not taking as Aspirin (in the same person at the same time under the same conditions)."
]
},
{
"cell_type": "markdown",
"id": "9e350596",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- A close approximation to comparing potential outcomes is to compare two groups of people that are similar on average (age, sex, income, etc.) except one group is allowed to take Aspirin after a headache and the other group takes a fake Aspirin (sugar pill/placebo) after a headache. This is an example of a randomized control trial."
]
},
{
"cell_type": "markdown",
"id": "8d4f8cd5-2b31-46d3-aad4-897bbb3021e2",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- Then the mean difference between time to pain relief should be due to Aspirin and not other factors related to why people may or may not take an Aspirin. "
]
},
{
"cell_type": "markdown",
"id": "5c44bf71-07a9-46af-aa37-6f46db59c025",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Does Material Deprivation Cause increases in mental health visits?\n",
"\n",
"- Neighbourhoods can't be randomly assigned to high or low deprivationdepravity, like people could be randomized to take a real or fake Aspirin.\n",
"\n",
"- There could be many other factors related to why a neighbourhood has a low or high deprivation score. \n",
"\n",
"- This means that when we are comparing neighbourhoods with low deprivation to high deprivation the differences in a certain outcome such as mental health visit rates could be due to factors other than deprivation such as environmental factors that are associated with mental illness, but happen to be more or less prevalent in one of the deprivation groups.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "a13007b5-a193-4a36-b630-906be48578dc",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Just in case you have seen a two-sample t-test ...\n",
"\n",
"This is similar to a two-sample t-test (we won't cover in this course except here) except not as flexible. "
]
},
{
"cell_type": "code",
"execution_count": 69,
"id": "9be7e8d7-eb51-45a5-b774-e2c7eb8348d1",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The p-value from the two-sample t-test is: 0.0\n"
]
}
],
"source": [
"from scipy import stats\n",
"\n",
"mh_visits_low = mh_visit_depriv.loc[\n",
" mh_visit_depriv[\"depriv_binary\"] == \"Low\", \"mh_visit_rates_mf\"]\n",
"\n",
"mh_visits_high = mh_visit_depriv.loc[\n",
" mh_visit_depriv[\"depriv_binary\"] == \"High\", \"mh_visit_rates_mf\"]\n",
"\n",
"statistic, pvalue = stats.ttest_ind(mh_visits_high, mh_visits_low)\n",
"\n",
"print(\"The p-value from the two-sample t-test is: \", round(pvalue, 3))"
]
},
{
"cell_type": "markdown",
"id": "798bf90b",
"metadata": {},
"source": [
"The two-sample t-test involves: \n",
"\n",
"1. Computing $t_{\\text observed} = \\frac{\\bar{x_{\\text Low}} - \\bar{x_{\\text High}}}{SD}$, where $SD$ is an estimate of the standard deviation of the mean difference.\n",
"\n",
"2. Computing the p-value using $t_{\\text observed}$ and the appropriate $t$-distribution.\n",
"\n",
"3. Assessing several statistical assumptions about the data to ensure the accuracy of the p-value.\n"
]
},
{
"cell_type": "markdown",
"id": "5a671c9e",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Review of this class\n",
"\n",
"- Comparing two samples by simulating the distribution of a test statistic (e.g., the difference in two means) assuming the null hypothesis is true (i.e., there is no difference in the means) performs a statistical test:\n",
"\n",
" \n",
" 1. Stated the null and alternative hypotheses;\n",
" \n",
" 2. imulated the test statistic assuming the null hypotheis by:\n",
" \n",
" + Randomly shuffling the group labels (e.g., high/low deprivation);\n",
" \n",
" + Calculating the test statistic in the shuffled data se; and\n",
" \n",
" + Repeating the previous two steps a large number of times (e.g., 5,000).\n",
" \n",
" 3. Computed the p-value as a measure of how consistent the data are with the null hypothesis. \n",
"\n",
" + The p-value is computed by summing the number of simulated values that are more extreme in the positive or negative direction.\n"
]
},
{
"cell_type": "markdown",
"id": "9575cf62",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"\n",
"- If the p-value is small then the data are inconsistent with the null hypothesis and there is evidence in support of the alternative hypothesis.\n",
"\n",
"- A large p-value isn't an indication of evidence in support of the null hypothesis. The procedure only tests whether you have data inconsistent with the null hypothesis.\n",
"\n",
"- A small p-value may not be **causal** evidence that for the alternative hypothesis unless random assignment was used part of the study design."
]
}
],
"metadata": {
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "ggr274",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.13"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
![](\"headache.svg\")
\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "4720a94f",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"![](\"pill.svg\")
\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "366ff202",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"![](\"clock.svg\")
\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "c3ca77ca",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"The **causal effect** of taking an Aspirin is 18 minutes (48 minutes - 30 minutes).\n"
]
},
{
"cell_type": "markdown",
"id": "a07027d9",
"metadata": {},
"source": [
""
]
},
{
"cell_type": "markdown",
"id": "7859b66e",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Potential outcomes and randomized control trials\n",
"\n",
"- Establishing causality involves comparing these **potential outcomes**. "
]
},
{
"cell_type": "markdown",
"id": "b7a9922e",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- The problem is that we can never observe both taking an Aspirin and not taking as Aspirin (in the same person at the same time under the same conditions)."
]
},
{
"cell_type": "markdown",
"id": "9e350596",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- A close approximation to comparing potential outcomes is to compare two groups of people that are similar on average (age, sex, income, etc.) except one group is allowed to take Aspirin after a headache and the other group takes a fake Aspirin (sugar pill/placebo) after a headache. This is an example of a randomized control trial."
]
},
{
"cell_type": "markdown",
"id": "8d4f8cd5-2b31-46d3-aad4-897bbb3021e2",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- Then the mean difference between time to pain relief should be due to Aspirin and not other factors related to why people may or may not take an Aspirin. "
]
},
{
"cell_type": "markdown",
"id": "5c44bf71-07a9-46af-aa37-6f46db59c025",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Does Material Deprivation Cause increases in mental health visits?\n",
"\n",
"- Neighbourhoods can't be randomly assigned to high or low deprivationdepravity, like people could be randomized to take a real or fake Aspirin.\n",
"\n",
"- There could be many other factors related to why a neighbourhood has a low or high deprivation score. \n",
"\n",
"- This means that when we are comparing neighbourhoods with low deprivation to high deprivation the differences in a certain outcome such as mental health visit rates could be due to factors other than deprivation such as environmental factors that are associated with mental illness, but happen to be more or less prevalent in one of the deprivation groups.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "a13007b5-a193-4a36-b630-906be48578dc",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Just in case you have seen a two-sample t-test ...\n",
"\n",
"This is similar to a two-sample t-test (we won't cover in this course except here) except not as flexible. "
]
},
{
"cell_type": "code",
"execution_count": 69,
"id": "9be7e8d7-eb51-45a5-b774-e2c7eb8348d1",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The p-value from the two-sample t-test is: 0.0\n"
]
}
],
"source": [
"from scipy import stats\n",
"\n",
"mh_visits_low = mh_visit_depriv.loc[\n",
" mh_visit_depriv[\"depriv_binary\"] == \"Low\", \"mh_visit_rates_mf\"]\n",
"\n",
"mh_visits_high = mh_visit_depriv.loc[\n",
" mh_visit_depriv[\"depriv_binary\"] == \"High\", \"mh_visit_rates_mf\"]\n",
"\n",
"statistic, pvalue = stats.ttest_ind(mh_visits_high, mh_visits_low)\n",
"\n",
"print(\"The p-value from the two-sample t-test is: \", round(pvalue, 3))"
]
},
{
"cell_type": "markdown",
"id": "798bf90b",
"metadata": {},
"source": [
"The two-sample t-test involves: \n",
"\n",
"1. Computing $t_{\\text observed} = \\frac{\\bar{x_{\\text Low}} - \\bar{x_{\\text High}}}{SD}$, where $SD$ is an estimate of the standard deviation of the mean difference.\n",
"\n",
"2. Computing the p-value using $t_{\\text observed}$ and the appropriate $t$-distribution.\n",
"\n",
"3. Assessing several statistical assumptions about the data to ensure the accuracy of the p-value.\n"
]
},
{
"cell_type": "markdown",
"id": "5a671c9e",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Review of this class\n",
"\n",
"- Comparing two samples by simulating the distribution of a test statistic (e.g., the difference in two means) assuming the null hypothesis is true (i.e., there is no difference in the means) performs a statistical test:\n",
"\n",
" \n",
" 1. Stated the null and alternative hypotheses;\n",
" \n",
" 2. imulated the test statistic assuming the null hypotheis by:\n",
" \n",
" + Randomly shuffling the group labels (e.g., high/low deprivation);\n",
" \n",
" + Calculating the test statistic in the shuffled data se; and\n",
" \n",
" + Repeating the previous two steps a large number of times (e.g., 5,000).\n",
" \n",
" 3. Computed the p-value as a measure of how consistent the data are with the null hypothesis. \n",
"\n",
" + The p-value is computed by summing the number of simulated values that are more extreme in the positive or negative direction.\n"
]
},
{
"cell_type": "markdown",
"id": "9575cf62",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"\n",
"- If the p-value is small then the data are inconsistent with the null hypothesis and there is evidence in support of the alternative hypothesis.\n",
"\n",
"- A large p-value isn't an indication of evidence in support of the null hypothesis. The procedure only tests whether you have data inconsistent with the null hypothesis.\n",
"\n",
"- A small p-value may not be **causal** evidence that for the alternative hypothesis unless random assignment was used part of the study design."
]
}
],
"metadata": {
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "ggr274",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.13"
}
},
"nbformat": 4,
"nbformat_minor": 5
}