sec2&3_bessel.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3348985a",
   "metadata": {},
   "source": [
    "## Bessel's Correction (Motivation)\n",
    "\n",
    "How do we estimate the variance of a random variable by operating on observations?\n",
    "\n",
    "First attempt:\n",
    "\n",
    "$$ \\hat{\\sigma^2}_{bias} = \\frac{1}{N}\\sum_i (x_i - \\bar{x})^2 $$\n",
    "\n",
    "where $\\bar{x}$ is the sample mean.  Lets see how this works out for us..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "196c1a42",
   "metadata": {},
   "source": [
    "## Building a ground truth distribution:\n",
    "\n",
    "In practice, we don't usually have access to ground truth.  But to help ourselves study, lets build a toy example with known ground truth variance.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8accfb49",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# a fair 6 sided die\n",
    "outcomes = np.array([1, 2, 3, 4, 5, 6])\n",
    "prob = np.array([1/6, 1/6, 1/6, 1/6, 1/6, 1/6])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f43f75da",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.916666666666666"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute exp_value & variance\n",
    "truth_exp_val = np.dot(outcomes, prob)\n",
    "truth_var = np.dot(outcomes ** 2, prob) - truth_exp_val**2\n",
    "truth_var"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e524fb1",
   "metadata": {},
   "source": [
    "## Lets \"roll\" the die N times and estimate the variance from our observations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "cb136d6f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[6, 6, 5, 3, 5]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from random import choices\n",
    "\n",
    "N = 5\n",
    "\n",
    "observations = choices(outcomes, prob, k=N)\n",
    "observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "321ddedb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.0"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute sample mean, our estimate of expected value\n",
    "est_exp_val = sum(observations) / N\n",
    "est_exp_val"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee830657",
   "metadata": {},
   "source": [
    "In the cell below:\n",
    "\n",
    "`est_var_bias` = $$ \\hat{\\sigma^2}_{bias} = \\frac{1}{N}\\sum_i (x_i - \\bar{x})^2 $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "993a0d70",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.2"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute biased variance estimator\n",
    "est_var_bias = sum((np.array(observations) - est_exp_val)**2) / N\n",
    "est_var_bias"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "463d8be5",
   "metadata": {},
   "source": [
    "## Maybe its just unlucky in these first few examples that our estimate variance was often lower than ground truth ... lets repeat that a whole bunch:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "41fe6342",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3wAAAFPCAYAAAAWZcwIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABoO0lEQVR4nO3dd1hUR9sG8JtdYOlVsICAFVRs2BGNIioWAjHYS4zd2DWxpLwxRU0UW6JRgx17QyyoqKDG2LtGrLEDgoXelt39/uBjdUNXdg/u3r/ryvW+nDNn5jlnFuFh5szoKRQKBYiIiIiIiEjriIQOgIiIiIiIiNSDCR8REREREZGWYsJHRERERESkpZjwERERERERaSkmfERERERERFqKCR8REREREZGWYsJHRERERESkpfSFDqAsvH6dBrm8fG0naGtrhpcvU4UOg/6D/VL+sE80K3bFHwCAyiO/KLIc+6X8YZ+UT+yX8od9Uv6wT9RLJNKDtbVpoee1IuGTyxXlLuEDUC5jIvZLecQ+0Zzs14kASvbM2S/lD/ukfGK/lD/sk/KHfSIcTukkItIxW+S1sHTXdaHDICIiIg3QihE+IiIquQzoQ54hFToMIiIi0gCO8BEREREREWkpJnxERERERERaigkfERERERGRluI7fEREOqaaXhIsXayFDoOoxDIy0pCamgiZLEfoUPKJjxdBLpcLHQa9hX1S/rBP3p1YrA8zMysYGxe+7UJxmPAREemYdnoxqNq6mtBhEJVIRkYaUlJew8rKDgYGhtDT0xM6JBX6+iLk5PAX2fKEfVL+sE/ejUKhgFSajcTEBAB456SPUzqJiIio3EpNTYSVlR0MDSXlLtkjIlInPT09GBpKYGVlh9TUxHeuhwkfEZGOCZG7YsG2K0KHQVQiMlkODAwMhQ6DiEgwBgaG7zWlnVM6iT5AMgBZUvVPjZAYiCBWeyukaTkQQaqBzw9RWeHIHhHpsvf9N5AJH9EHKEsqR8SZh2pvp1NLF5gYcCJASWkiEWcSTkRERKXBhI+IqIxoIhFnEk5ERESlwd8aiIiIiLTQpUsX4OXVFJcuXRA6lDIVGOiHWbNmCh2GUnj4Xnh5NUVsbIzQobw3mUyG/v0DsWlTiPLYqlUr4OXVFCkpKcVe7+XVFKtWrVBLbD/88C3+978Zaqlb2zHhIyLSMbX1XqNhzQpCh0FEZeTIkUPYtm2T0GGouHHjOlatWlGiJIHKj4iIA3j58iUCAnoIHUo+/foNwrFjR3H//j2hQ/ngcEonEZGOaa0Xh6otnIQOg+i9aWoBq6IYi/Qg9JIyR49G4O7dO+jVq5/K8UaNPHD06N8wMDDQeEw3b17HmjXB6NrVD+bm5hpvX5M6d+6KDh06wdDww19NdsuWDejY0RcmJu+239vRo39DLFbPm+a1atWGm1tdbN26EV9//b1a2tBWTPiIiIjog6SpBayK4tvKBUb65XPClEgkgkQiEToMrZWRkQFjY2OIxWK1JTmadOtWNO7fv4eJE7965zrU/Xlr394Ha9YEY/LkaTAyMlJrW9qECR8RkY5ZI68DycZLmNbfQ+hQiHTa8+dx+PPPP3D27GmkpaXCyckZAwYMRseOvsoyOTk5WL9+NSIiDiA+/jmMjIzh7OyCIUOGo1mzlhg7dgSuXLkEIPf9KQCoVKkyduzYi0uXLmD8+FH47bfl8PDIPTd27AikpqZixoz/YdGiebhz5xYqVqyEceMmoVUrL5w+fRLLly/FkyeP4OzsgunT/wdXVzdlPPfu3cXWrRtx5cplvHyZADMzc7Rs6YkxYybA0tIKQO47X2vWBAMAevb8WHnt9u17ULlyFQDA/v17sHPnVjx8+BDGxkZo1coLX3wxHjY2tsryCoUC69atQljYLiQnJ6FuXXdMmjS12OcqlUrRpYsP2rfvgGnTvlU59+LFC/To0RVDhozA4MHDkJychPXr1+DcudOIjY2Bnp4I9es3xKhRY1GrVm3ldXnP8scff8GdO7dw4MA+vHr1En/9dR7h4Xsxe/YPKvf311/HsGdPKO7cuY3k5CTY2dmja1c/DBz4uUpymNcf//vfj1iwYC6io/+BubkFevbsg/79P1OJPSsrEyEha3HkyCE8fx4HCwtLNGrUGGPHToKdnb2yzLp1q3H48CG8eBEPGxtbdOnSHYMHD4O+ftG/9v/99wkYGkpQv37DAs+/fv0Kv/76E86ePQOJRIKuXbtjxIgxKvV6eTXF558Px9ChIwEAcXGx2LBhHS5dOo+4uDgYGRnBw6MpxoyZoHxWQPGf8zzNmrXA0qWLcOHCOXh5tS3yfuiNEiV8Dx8+xKJFi3Dp0iUkJyejSpUqCAgIwODBg1WGry9duoR58+bh5s2bMDMzQ5cuXTBlyhQYGxur1JednY3FixcjLCwMycnJcHNzw6RJk9CqVauyvTsiIiKicujFixcYOfJzGBgYoGfPPrCwsMRffx3HDz98C6lUiq5d/QDkJk+bN4egR4+eqF69BlJSUnHr1j+4ffsWmjVric8+G4KMjAw8fx6LceMmAwCMjU2KbDs5OQnTp09Gp05d0KFDR+zYsRXffDMVM2b8D3/88Rs++SQQYrEYISFr8f33X2PTph0QiXJHMc+fP4uYmGfo1s0PNja2ePDgX+zZswsPHvyLP/9cCz09PXz0kTdiYp7i0KEDGD9+sjIRtLKyBgCsXv0n1q5dCR+fzvj440/w8uVLbN++Gbdu3cSqVSGQSHJHblauXI5161bB07MNWrRohTt3bmHy5LHIyZEWeX8GBgZo27YdTpyIwpQp01USkmPHjkAul6NDh04AgJiYZ/jrr2No394HlStXwevXrxAWtgvjxo3Ahg3bUaGCnUrdq1f/CSMjI/Tv/xnS09MKjSE8fB+MjU3Qu3d/mJgY4+LFC1i5cjnS0tIwZsyEfP0xZcp4tG/fAR06dEJU1BEsW/Y7qleviVatWgPIXUzlyy8n4PLli+jUqQt69eqHtLRUnDp1Ek+fPoGdnT3kcjmmTp2MmzdvwN+/B6pWdcLt29FYv3414uOfFzsN8vr1q6hRo2ahieF3301DlSqOGDVqLK5du4JNm0KQnp6BL7+cXmid0dH/4MaNa+jYsRNsbe0RGxuD3bt3Yty4kdiwYbtylK64z3keF5dqkEgkuH79KhO+Uig24Xv+/Dl69uwJc3NzDBgwAJaWlrhw4QLmz5+Pu3fvYt68eQCA6OhoDB48GDVr1sT06dMRFxeH1atX4+nTp1i+fLlKndOnT0dERAQGDRoEZ2dnhIaGYvjw4QgJCUHjxo3Vc6dERERE5URw8B8QiURYvXqj8h23Tz4JxJQp47FixVL4+naDSCTC6dN/w8/vE4wfP6XAepo1a4ldu7YjKSkRnTt3LVHb8fHP8eOPv8Db2wcAULeuO0aMGIxZs2Zi/fotcHJyAQCYmZlj3rzZuH79Kho2zP39rEePQPTtO0Clvnr13DFz5je4du0KGjZsjJo1a8HVtQ4OHTqANm3aqYzkxMbGYN26VRg9erxKPS1aeGL06CE4cGAfAgIC8fr1a2zatB5t2nyE2bODlBtPr1ixFCEha4q9R2/vjggP34uLF8+jRYs3AwqRkUdQu7YrqlbNfY+5evWa2Lx5lzKhBXLfyevfPxD79oVh8OBhKvXm5ORg6dLgYqcuzpz5szJxBYCAgEDMmzcboaHbMXz4aJUBk/j45/j++5+VI7vdu/sjMLA79u8PUyZ8Bw7sxeXLFzFx4pcIDOyjvHbgwM+hUCgA5C64cuXKRfzxxyrUq+euLFOligOWL1+C/v0/g7OzS6ExP378CPXrNyj0vKNjVcyalft7/6ef9oJEIkFY2E706zcQVao4FHiNp6cX2rf3gb6+CDk5ue/btm7dFqNGfY5jx47C17cbABT7Oc+jr68Pe/uKePjw3yLLkapiJ53njcKtWLECI0aMQO/evTFv3jx07twZ4eHhkEpz/8qyYMECWFlZISQkBH379sWkSZPwv//9D1FRUTh9+rSyvmvXrmH//v348ssvMXXqVPTu3Rvr1q1D5cqVERQUpL47JSIiIioHFAoFjh+PQuvWbSGTyZCYmKj8r0WLVnj58gWePHkMADAzM8PNmzcQH/+8zNo3MzNH+/YdlF+7udWFWCyGu3sDZbIH5CaCQO4oWJ63k5isrCwkJiaiXr36AIDbt28V2/aJE1FQKBT46KP2Kvft6FgVtrYVcPnyRQDAhQtnIZVK8emnvZXJHoB8C9MUpmnT5rCwsERk5GHlsRcvEnD9+lV4e3dUHjM0NFQmezKZDElJiTA2NkHVqs4F3k/Xrt1L9J7a288pPT0NiYmJaNiwMTIzM/Ho0UOVsmZmZvDx6az82sDAAHXq1FN57sePR8HGxhaffNIzX1t5z+fYsaOoVq0GHBwcVZ5t06bNAQCXLxe9PUdSUiLMzS0KPd+jR698XysUCpw7d7qQK1SfQ05ODpKScvvazMwcd+68eb6l+Zybm1sgKSmp2HL0RrEjfGlpucPVtra2KscrVKgAfX19iMVipKam4tSpUxg6dChMTd+s6uPv74/Zs2fjwIEDyumaBw8e/P/pC28+sBKJBIGBgVi4cCHi4+Nhb29fJjdHREREVN4kJr5GamoKQkO3IzR0e6FlnJ1dMGzYKEyfPgWfftodtWq5okWLVujUqQuqVav+zu3b29urJFEikQjGxsaoWLGiSjkzMzMAUNlaITk5CatXB+Po0Qi8fv1KpXxaWmqxbT958gRyuRy9evkXeD4xMREAEBcXBwBwdFRdUdja2rrIpCSPvr4+PvqoPY4di8RXX30NfX19REYegUKhUEn45HI5tm/fjNDQHYiNjYFMJlOes7S0zFfv26OVRfn33/sIDl6GS5fOK3+XzvPf52RvX1GlP4DcpObt7QdiYp7Bycm5yMVhnj59gocPH6B7d58Cz+c926LkjRYWxNGxqsrXVavmfh0bG1voNXnvHYaH70VCQrxK/ampb55DaT7nCoUCekIvjfuBKTbha9asGZYvX45vvvkGEyZMgKWlJc6fP6+chikSiXD79m3k5OTA3d1d5VpDQ0PUqVMH0dHRymPR0dGoVq2aSmIIAA0aNIBCoUB0dDQTPiIiNaqn9xLWdZoKHQaRzpLLc6e2de3qp7JAy9uqVasBIHdrhW3bduPkyRM4d+4Mdu/eiU2b1mPq1G/QrdvHBV5bHJGo4KShsOPAm1/Sv/tuBm7cuIp+/QahZs3aMDExhlyuwJQp45T3VRSFQg6xWIygoMVAARtalCSZKylvbx/s3bsb58+fQatWXoiKOoI6deqqTD9cv341Vq5cjm7dPsawYaNgYWEJPT09/Pbb/ALv5+0Rq8KkpKRg3LgRMDExw9Cho+Dg4AhDQ0PcuXMLy5b9nq/ewp57UclXQeRyOWrXdsXo0eMLPF/YtMs8lpZWZb5v4sKF8xAevhe9e/dD3bru//9HBD3MnPm1yv2V5nOekpICFxeXMo1T2xWb8Hl5eWHChAlYsWIFIiMjlcfHjx+PMWPGAAASEhIAAHZ2dvmut7Ozw5UrV5RfJyQk5PsL0tvXxsfHl+4OiIioVJrrxaOqh6PQYRDpLCsra5iYmEKhUKBZsxbFlrewsETXrn7o2tUPGRkZGDduJFatWvHWL8KaGe5ITk7GxYvnMHToSHz++XDl8bzpp6oKjsnBwREymQyOjk5FjpZVqlQJAPD06WPl/weA169fIyUluUTxeng0g5WVFSIjj6Bmzdq4ceMavvhCdcGUY8ci4eHRFDNm/E/leGpqqnKxmdK6fPkikpKSMGvWPDRq9GY15NjYmHeqD8h9brduRSMnJ6fQRVUcHBzx4MG/JfpMFcTJyRmxsc8KPf/06RNUrPimL548eQIgd1XYwuS9pzdhwmTlO3xZWVkqo3t5iv+c504LjY9/jo8+al/q+9NlJVql09HREc2bN0fHjh1hZWWFY8eO4ffff4eNjQ369u2LzMxMAChww0mJRKI8DwCZmZkFbgCaNx86Kyur1Ddha2tW6ms0wc5Ouzca/VBpQ788f5UOU1P1761kZGQAO5uiV3srC9rQJ4Bm+uV9+yTOUB/ZCj2YWxrDyLDoHwHa0i/aRBf7JD5eBP1C9rkT5cghEgk/t6uw+Ioq365dexw9ehiDBw/Jt5DG69evYW2du6JlUlKiSuJhbm4KJycnxMXFKNs1MTFGWlpqvjjEYpHyf/PO6enpQU+v4Jj19PRUjuddLxLlHpdI9P+/nOr1O3ZsVikHAKamuf9OZWSkqZRt394by5cvwbp1K/HttzNV2pfL5UhJSYGlpSVatmwFfX197Nq1DS1bvlmlMa+t/8b6X/r6uffcvn0HHD58CNWr504N7NSpU4H3+Paxo0cPIyEhHo6OVZXH3zzL/O3mfQbznrOBgTjf85BKpdi9e4dKubz7KKg/8qZ45h1v184bp06dxJ49u9CrVx+VsrlTHPXQoYMPfv75Bxw6tB/duvmplElPT4dYLC7y/cMGDRpiw4Z1AGTQ13/zu3re/YWGbkeLFm+Syd27t0NPTw+enp4q8b993yKRWDn9Mu/Yli3bIJPJVPqwJJ9zAHjw4CGys7PQsGGjUn/ffehEItE7/wwoNuHbv38/vv/+exw8eFA5MtepUycoFArMnTsXXbt2VS6pmp2dne/6rKwslY0RjYyMlAu9/Lcc8G4bNr58mQq5vHTD3upmZ2eOhISyHRan96ct/ZIplSMtrfR/HCl1O5lStT8vbekTQDP98r59Is3OwVp5HUj++LvIffi0qV+0ha72iVwuV44M5DunQLn4+V9YfEUZMWIsLlw4j8GDB+Djjz+Bs7MLkpISER19E3fu3MKOHXsBAH36BKJRIw+4udWBhYUlbt+OxpEjEejRo6ey3dq13XDo0AEsXBgEN7e6MDY2gZdXW8hkuedlsjfPUKFQQKEoOGaFQqFyPO96uTz3uERigkaNPLBhwzpkZ0thZ2ePc+fOKEeu8soBQK1argCAZcuWoEOHTtDX10fr1m1RubIjhg4dieDgZXj69Ck8PdvA2NgYz549xfHjkRg0aAj8/AJgbm6Jvn0HIiRkDSZPnoCWLT1x584tnDlzClZWVvlifdvbK0K2a+eD0NCdWLNmFdzd68PW1l7lOk/PNlizJhg//vg93N0b4N9/7yEi4iCqVHFQaePNs8zfbt5nMO85161bH+bmFvjxx/8hMDB30ZlDh8LzlSuqP/KmO+Yd79SpK/bv34sFC+bin3+uw929IdLT03D69N8YOnQkGjdugo4duyIiIgI//zwTZ8+eQb169ZGTI8XDhw8QGXkYq1ZtyPce3ts8Pdti9epgXLx4CU2aNMt3f0+ePMZXX01Gs2YtcPXqZRw9GgF//x6oWLGKSvxvfw48Pb1w8GA4zMzM4ORUDf/8cx0XLpyDpaWlyvMtyeccAM6cOQNjYxM0btz0nb7vPmRyubzQnwEikV6RA2DFJnybNm1CvXr18k3D9Pb2xq5du3Dr1i3ldMy8qZ1vS0hIUHknz87OrsBpm3nX8v09IiIiKgmJgQidWroIG4OhODfzLKUKFSogOHgd1qwJRlTUEbx69RKWllaoUaMWhg8frSwXGNgbJ0+ewPnzZyGVZqNSpcoYNmwU+vUbpCzj7/8p7ty5hfDwfdi6dRMqVaqstj3Kvv/+ZyxcOA+7dm0HoECzZi0RFPQbAgJU30WsXdsNI0eOwa5d23H27On/XxxlD4yNjfHZZ0Ph6OiE7ds3Y/XqFdDTE6FixUpo0+YjlemIedsX7N69ExcvnkPduu5YsGAJpk6dWOJ4GzXygI2NLV69eon27TvmOz9w4OfIyMjA4cMHcfRoBGrXdsPcuYuwfPnv7/qIYGlphblzF2LJkkUIDl4Gc3MLdOrUBU2bNsfkyWPfqU59fX3Mn/871q5diaNHI3D06GFYWVmjYcPGyiROLBbj118XYPPmEEREHEBU1BEYG5vAwcERAwYMLvDVq7e5udVBjRo1ERV1RCXhy/PTT79ixYolWLbsd0gkhujTZwBGjSr6fiZM+BIikQiHDh1AVlY26tdviEWLlmLy5HEq5UryOQeAqKgjaNfOO98e31Q0PUUxb4R27twZNjY22Lx5s8rx8PBwTJo0CWvWrEH9+vXRsmVLDB06FJMnT1aWyc7ORosWLeDn54cff/wRAPDrr78iJCQEZ8+eVVm4Zfny5Vi4cCFOnDhR4Dt+ReEIH5WUtvRLulSOiDMP1d5Op5YuMDFQ75QJbekTQDP98r598mTuHKyR14GkqhNH+D4wutoncXGPUKmSs9BhFOrt0SQqH9gn7+7AgX1YvHg+du3aBxMT0+IvKKGy6JO7d+9g6NABWL16I2rWrFVGkX04ivq3sLgRvmJ/a6hWrRpu3LiBx49VX8jdv38/xGIxXF1dYW5ujlatWiEsLExl6dmwsDCkp6fD1/fNX318fX0hlUqxffubZYizs7Oxa9cueHh4lDrZIyIiKowMuYm4uv+TFRsJEVH516lTF9ja2mL37l1Ch5LPpk3r0a5dB51M9t5XsVM6hw4dihMnTqBv377o378/LC0tcezYMZw4cQJ9+vRR7s83adIk9OnTBwMHDkTPnj0RFxeHNWvWoG3btvD09FTW17BhQ/j6+iIoKAgJCQlwcnJCaGgoYmJiMGfOHPXdKRER6ZwsLRoNJyJSN7FYjI0bdwgdRoG+//5noUP4YJVoH74tW7bg999/x6ZNm5CYmAgHBwdMmTIFQ4cOVZarV68e1qxZg6CgIMyZMwdmZmbo1auXyhTPPHPnzsWiRYsQFhaGpKQkuLq64s8//0STJk3K9u6IiCifRnoJsKn/bst2ExER0YelRNsyNGjQAMHBwcWWa9q0KbZs2VJsOYlEgmnTpmHatGklaZ6IiMpQY70XqNqg8H2TiIiISHtw/gkRkY5JU+gjJT3/NjpERESkfZjwERHpmG2KWvgj9IbQYRAREZEGMOEjIiIiIiLSUkz4iIiIiIiItBQTPiIiIiIiIi3FhI+IiIiIiEhLlWhbBiIi0h7N9J7D1qOV0GEQERGRBnCEj4hIx7jrvULzOhWFDoOItNCqVSvg5dVU6DAKFR6+F15eTXH37m2hQylWYKAfZs2aKXQYpAWY8BER6ZgkhSFeJWcKHQYR6bAjRw5h27ZNaqt/9+4dCA/fq7b6iT4knNJJRKRjdilqQLL3Jqb19xA6FCLSUUePRuDu3Tvo1aufWurfvXsXzMzM0LWrn1rq14RNm3ZCJOLYDL0/foqIiIiIPiAymQzZ2dlCh6ExGRkZQoegUZmZuTMwDA0Noa/PsRl6f0z4iIiIiARw6dIFDB06EN7enujVyx+7d+8s8B04L6+mWLx4Pg4c2Id+/T6Ft7cnbty4BgC4ffsWJk8ei44d26Jjx7aYPHkc7t69o3L9rFkzERiYf6SrqLaOHTuKAQN6oX37VhgwoBfOnDmV7/qrV69g2LBBKvGXxNixI/DXX8cRFxcLL6+m8PJqqozv0qUL8PJqisjII1i+fAn8/X3RqVPbQuMF3ryXFxsbAwAICOiGe/fu4MqVS8r6x44doXJNVlY2Fi+ej+7dfeDj44UZM77E69evi4x706b1aNOmGeLjn+c7FxQ0Bx07tkVWVm6y9tdfx/DVVxPg7++L9u1boVcvf6xduxIymSzfsxg8uB9u3ryB0aOHwtu7NTZuXAcg/zt8yclJWLJkEQYN6o2OHdugU6ePMGXK+Hz9nfcMjx+PxNq1KxEQ0AXe3p6YMGE0nj59ki/2GzeuYcqU8fD1bYeOHdvg88/7Ye/e3Splrl+/iokTv0CnTh/Bx8cLEyZ8gVu3oot8XlR+8M8GRERERBp2584tfPnleFSoYIchQ0ZALpdj7dqVsLKyLrD8+fNnEBkZgU8+6Qlzc3NUqFAB//57H2PHDoe5uQUGDhwMANi9eyfGjBmOP/9cCxeXau8U25UrFxEVdQSffBIIY2MT7NixBd9+OxU7d+6DpaUVAOD+/XuYPHkMrK1tMGTICOTk5GD16j9hbW1TbP2ffTYEGRkZeP48FuPGTQYAGBubqJRZvfpPGBkZoX//z5Cenlaq+CdN+hILFsyDkZERBg0aAgCwsVGNa8GCX2BhYYnPPx+B2NgYbN++GQsXzsWPP84ptF5v747444/fEBV1BL1791cel8lkOH48Cl5ebSGRGAEAwsP3wdjYBL1794eJiTEuXryAlSuXIy0tDWPGTFCpNzHxNaZOnYROnXzh69sNFStWKrD9mJhn+OuvY2jf3geVK1fB69evEBa2C+PGjcCGDdtRoYKdSvl161ZBJBKjX79BSElJxubNIfjhh28RHLxOWebMmVOYPn0y7OwqolevfrCxscG9e/dw6tRJ+PkFAADOnz+Lr76agHr16mPYsJGQy+XYsycUY8cOR3DwelSrVr3I/iDhMeEjIiIi0rBVq1ZALNbHsmWrYGtbAUBuQtG/f2CB5Z88eYyQkG1wcnJWHpsx40vIZHL88cdKVKpUGQDQsaMv+vcPRHDwH5g1a947xfbo0UNs2LAdVao4AAA8PJpi8OC+OHLkED79tDcAYOXK5dDT08OyZatgZ2cPAGjXrgM++6xPsfU3a9YSu3ZtR1JSIjp37lpgmZycHCxdGgyJRFLq+D/6qD2Cg1fAzMys0PqtrGwwf/5v0NPTAwAoFHLs2LEVqampMDMzK/CaSpUqo25dd0RGqiZ8ly9fxOvXr9ChQ0flsZkzf1YmfwAQEBCIefNmIzR0O4YPHw1DQ0PluRcvEjB9+nfo3t2/yPuqXr0mNm/epfJeX+fOXdG/fyD27QvD4MHDVMrnJuHrlNNCLSwssXhxEP799x6qV68JmUyG+fN/gb19RaxevVHlvhUKBQBALpdj/vxf0Lx5K8ydu1B5vnv3APTv/ynWrg3GDz8UniRT+cCEj4hIx3jqxaJCcy+hwyB6L8mn/kbSyRNChwHrjz6CWQvPUl0jk8lw4cI5tG/vo0z2AMDRsSpatvTE33//le8aD4+mKsmeTCbD+fNn8NFH7ZXJHgBUrlwFbdq0w6lTf0Emk0EsFpf6npo3b6lM9gCgZs1aMDU1RUzMM2Xb586dRrt23spkDwBcXKqhefOWOH3671K3+V9du3Z/p2SvpPz9eyiTPQBo0KAxtm7dhLi4WNSsWavQ67y9fbB06WLExcWhUqXckbjIyMMwMzND8+Zv9jd9O9lLT09DdrYUDRs2RljYLjx69BC1atVWnjcyMoKvb7diY347SZTJZEhNTYGxsQmqVnXG7du38pXv1u1jlXcAGzZsBCB3pLB69Zq4c+cWYmNjMGnSV/mS3Lxnc+/eHTx9+gRDhoxEYmKiSpkGDRrj8uVLxcZNwmPCR0SkY1z1ElG1VoXiCxKRWrx+/RpZWVlwdKya71xBx4DcRO5tiYmvkZmZqZIE5nF2dsHRoxFISkqEjY1tqeMraEqhubkFUlJSlG3nxu+Ur5yTk3OZJHz/vd+y9t97NDc3BwCkpCQXeZ23d0csXboYkZGH0a/fQMhkMpw4EQUvr49UErJ//72P4OBluHTpPNLSVKekpqWlqnxtZ2dfosVZ5HI5tm/fjNDQHYiNjVF5H9DS0rIE92jx//eY2495CXy1ajUKbfPJk9x3/n788dsCz3MV0Q8DEz4iIh3zQmEE/ZdpqGxrKnQoRO/MwrM1LDxbCx0G9PVFyMmRq72dt0eMSuvtkay3yeUFxy0SFTwqmDfNTxMKut/C70NW4PGiiMUFJyrF3aO9fUW4u9dHVFRuwnfp0gUkJiaqTOdMSUnBuHEjYGJihqFDR8HBwRGGhoa4c+cWli37Pd9zL2nfrl+/GitXLke3bh9j2LBRsLCwhJ6eHn77bX6BfVkW/ahQ5NY7fvzkIhNDKt+Y8BER6Zi9imqQHLzNffiIBGJtbQ1DQ0mBKyYWdKwgVlbWMDIywuPHj/Kde/z4EYyNjZULrJibmyM1NSVfubi42NIF/lbbEokET58+LrDtkik4eSvK2yNUeSNyABAXF5e/9tJXX2Lt23fEb7/NR0zMM0RGHoG5uQWaNWupPH/58kUkJSVh1qx5aNTozb+zeauIvqtjxyLh4dEUM2b8T+V4amqqsq9Lw8HBEQDw4MF9eHjkX/307TJmZuZo1qxFqdug8oHjsEREREQaJBaL0bRpcxw/HoWXL18ojz99+qTA7Q8Kq6NZs5Y4fjxKJeGJi4vDiRPH0Lx5S+X7e1WqOCI1NRX37t1Vlnvx4gVOnDj2zvE3b94Kx49HISEhXnn84cMHOHfuTInqMDY2RmpqavEF35KXfFy9+ua9sYyMDBw4sC9fWSOj0tdfUt7ePhCJRDh8+CD++isKH33UXmVKZt40x7dH0qRSKUJDt79XuyKRKN/oXGTkEZU+KI1atVxRqVJlbN26Kd+zymundm03VKnigM2bQ5T7A76tuK0sqHzgCB8RERGRhg0ZMgKjRw/B6NFD4e/fA3K5HDt3bkP16jXy7atWmOHDR+PChbP44ouh+OST3NU9Q0N3QCwWY/jwL5TlfHw6Yfny3/H1118iMLAPsrIyERq6A1WrOuHOnfyLfZTE0KEjcfbsaYwePRQBAZ9CJpNh585tcHGpjvv37xZ7vaurGyIiDuD33xfAza0ujI1N4OXVtshrmjdviYoVK+GXX35C374PIRKJsX//HlhZWeP5c9VRPlfXOti5cyvWrl0JR8eqsLa2QZMmzd7pXv+rQgU71K/fEBs3rkd6eprKdE4AqF+/AczNLf5//8Pe0NPTw6FD4XjfGbGtW7fBmjXBmD37B7i7N8C//95DRMRBlQV2SkMsFmPKlGmYPn0KPv+8P7p27Q5raxs8fPgv4uPjMXv2PIjFYkyd+g2mTp2IQYN6o0uX7rC1rYCEhHicO3cGjo6O+O67n97vxkjtmPARERERaZibWx0EBf2GpUsXYeXK5bC3r4hhw0bi4cOHePSoZNMiq1evgSVLgrF8+e9Yv341AKB+/YYYPXqcyh58lpZWmD17Hn7/fSGWLfsNlStXwahRY/HkyeN3Tvhq1qyFBQt+x++/L8SqVStgZ2ePIUNG4OXLFyVK+Pz9P8WdO7cQHr4PW7duQqVKlYtN+PT19TF7dhDmz/8FK1cuh42NLXr16gtzcwvMnv2DStnPPhuC2Nhn2LQpBOnpaWjUyKPMEj4gd/GWq1cvw8rKCh4eqvVaWlph7tyFWLJkEYKDl8Hc3AKdOnVB06bNMXny2Hduc+DAz5GRkYHDhw/i6NEI1K7thrlzF2H58t/fuc5WrbywePEyrFkTjE2bQgAAVatWRY8evZRlmjZtjmXLVmPt2mBs374FmZkZsLW1Q/36DeDv/+k7t02ao6fQ5Bu4avLyZSrk8vJ1G3Z25khIyD9fnoSlLf2SLpUj4sxDtbfTqaULTAzUO/NbW/oE0Ey/vG+fPJk7B2vkdSCp6lTkO3za0i/8XvnwxcU9QqVK+VeiLC/KetGWGTOm4MGDf7FlS2iZ1alrNLWQDpUc++T9FfVvoUikB1vbgvePBPgOHxGRzmmr9wzdW7sIHQaRzsvKylL5+smTxzh9+m80btxEoIiISBtxSicRkY6poZeMqi42QodBpNNycnLQq5c/unTpjipVHBAbG4OwsF0wMDBAv36DhA6PiLQIEz4iIh0TqzCB4nkKnCqaF1+YiNQid6XLljhy5BBevXoJAwMDuLs3xIgRX6Bq1fwbmhMRvSsmfEREOuagwhmSI3e5Dx+RgPT09PDNNzOFDoOIdADf4SMiIiIiItJSTPiIiIiIiIi0FBM+IiIiIiIiLcWEj4iIiIiISEtx0RYiIh3TQe8JKn7UXugwiIiISAOY8BER6RgnvVRUdbQUOgwiIiLSAE7pJCLSMY8VZrj3NEnoMIiIiEgDmPAREemYo4qq2Hn8vtBhEBERkQYw4SMiIiKici82NgZeXk0RHr5X6FDK1NixIzB27Aihw1C6dOkCvLya4tKlC0KHUiYmTRqDBQt+VX4dHr4XXl5Ncffu7WKvDQz0w6xZM9US14oVSzFixGC11P1fTPiIiIiIqNw4c+YUVq1aIXQYKh49eohVq1YgNjZG6FCoFK5cuYRLly6gb9+BQoeST8+efXDv3l2cOnVS7W1x0RYiIiIiKjfOnj2N7ds3Y+jQkSrHK1WqjKNH/4a+vuZ/fX306CHWrAlG48ZNULlyFY23r0mNGnng6NG/YWBgIHQo723z5hA0b97ynfts06adEInUMz5mY2OLNm3aYsuWDfD09FJLG3k4wkdERERE5Z6enh4kEgnEYrHQoWglqVSKnJwciEQiSCQStSU6mvLq1UucOXMK3t4d37kOQ0NDtf6BoX17H1y+fBFxcbFqawPgCB8Rkc7x1XuESj4dhA6DSKetWrUCa9YEY8uWUAQH/4HTp09BIpGgd+9+GDBgMB4/foSFC+fi+vWrsLCwxMiRY9C5c1eVOpKTk7Bq1QqcOHEMiYmvUalSZfTo0Qs9e/ZRlpFKpVi3bhVOnTqJZ8+eQCaToXZtNwwbNgoeHk2V5WJjY9Cz58cYP34yDA0NsWlTCBIS4lGjRi1MmTINderUK/aeShIPAOzYsQW7d+9CbOwzGBgYwsHBEb1790enTr6YNWsmDhzYBwDw8noT38mTF5Qxfv319+ja1Q8AMGvWTPz11zGsXr0R8+f/iqtXL8HKygpDhoxE165++OefG/j99wW4c+c27O0rYsqUqWjWrKWy3ri4WGzYsA4XL57D8+fPYWRkBA+PphgzZoJyVCg8fC9mz/4BADB+/Cjltb/9tlz5DP/++y+EhKzBvXt3IBaL4eHRDGPGTICjY1WVew8L24WNG9fhxYsXqFGjBsaOnVTscwWAAQN6oUKFCli06A+V41KpFH5+HdGuXQdMn/5dqft73LhJkMsV2LVrG54/j8PWrbsRGxuD8eNHqdzf1auXsX37Fty8eQOvX7+CtbUN2rXzxsiRYyCRGCnrzeuPdeu2YMGCX3Hp0gVIJBJ07eqHkSPHqiTrcrkc27Ztwv79e/Ds2VOYmJiibl13jBo1FtWr11CW2bJlI/bvD0NMzDNYWFigXbsOGDlyLExMTIp8ZqdP/w2ZTIamTZsXeD4jIwO//vozoqKOAlCgXbsOGD9+MkxMTJVlAgP90LhxE3zzzUwAuZ/x9evX4Ny504iNjYGengj16zfEqFFjUatWbZX6i/qc52natAUA4OTJ4wgMVP0+KUtM+IiIdExlvXRUrWgudBhEBODbb6ehevUaGD16HI4fj8Ty5UtgZmaODRvW4qOP2sPL6yOEhe3E7Nk/oEGDRsokJCMjA2PHjsCrV68QEPApKlSww+XLF7B4cRBSUpIxZEjuIiBpaWnYu3c3fHw64+OPA5Ceno59+8IwefJYBAevQ61arirxHDy4HxkZGfD37wE9PT1s2rQe33wzFdu2hRU50lHSePbsCcWiRUHo3t0fvXr1RWZmJu7du4ObN2+gUydf+Pv3wMuXL3Hu3Gl8992PJXqGOTk5+PLL8WjSpDlat26D8PA9mDPnR+jrG2Dp0kXw8wuAt3dHbNmyAd99Nx27doUrk4Xo6H9w48Y1+Ph0hp2dPWJjY7B7906MGzcSGzZsh5GRERo2bIxevfpi27bNGDjwc7i4VAMA5f+Gh+/FnDk/wtPTC6NHj0d6ejp27tyKL74YhrVrN8HGxhYAsG/fbsybNxsNGjRCr159ERPzDNOmTYaFhQXs7SsWeY/e3j5Yt24VEhMTYWVlpTx+9uxppKamokOHju/U33v3hiEnJweffBIIPT0RjI0LTqKioo4gKysTn3wSCAsLS0RH/4OdO7chPj4eP//8q0rZnJwcTJ48FvXrN8SYMRNw/vxZbNoUgsqVHfDJJ4HKcrNmzcShQ+Hw9GwDf/8eyM7OxqVLF3D7drQy4fvll59w+PBBdOv2MXr27IunT59g165tePjwARYt+gN6enqFPrMbN67B2toGdnb2BZ6fP/8XWFhYYvjwUXjw4F+Ehe3Cy5cvMG/e4kLrjIl5hr/+Oob27X1QuXIVvH79CmFhuzBu3Ahs2LAdFSrYASj+c57HzMwMDg6OuH79KhM+IiIqO/cVFkh++Ar1XGyEDoXovfy68VK+Y83q2MPbwxFZUhkWbbua73zr+pXh1aAyUtKz8UfojXzn23s4oHmdiniVnIngvTfzne/c3AmNalVA7Ms0rD94G9981jRfmdKoX78hpkyZBgDo2tUPAQFdMH/+L5g27Rt07x6Qe0/NmqNfv0AcOhSOwYOHAQC2bt2IuLg4rF27CVWqOAAAAgI+hbm5BTZsWIfAwD6wsLCAubk5duzYq/I+lp/fJ+jfPxA7dmzFjBn/U4knPj4eW7aEwszMDADg5OSM6dOn4OzZ02jduk2h91HSeE6fPolWrVpj+vTvCqzH3b0BnJ1dcO7c6XwjmoXJzMxE9+7+6NdvEACgfXtvBAR0xU8/fYeFC5cqR3hcXKph8uSxOHEiCr6+3QAAnp5eaN/eR6W+1q3bYtSoz3Hs2FH4+naDg4MjGjVqgm3bNqNZsxYqI2Xp6en47bf56NGjJyZNmqo83qFDRwwc2Atbt27C6NHjkJOTgxUr/oCrax389ttyZfLs4lIdc+fOKjbh69ChE1av/hPHjh1FQMCnyuNRUYdhZWUND49mAFDq/n7xIgFbt4bC0tKqyPZHjx6nMpLn798DDg5V8eefSxEXF4dKlSopz2VmZsLXtxsGDvwcABAQEIghQwZg374wZcJ34cI5HDoUjt69+2PcuDejnP36DYJCoQCQO6oYHr4XP//8K9q1ezMrpU6duvj++69x9uxptGzpWWjMjx49LPLdPUNDCRYuXKrsC1vbCli1agUuX76Ixo2bFHhN9eo1sXnzLpXprp07d0X//oHYty9M+f1Z3Of8bVWqOODhwwfFlnsfH/bkXCIiKrUTCgfs+/uh0GEQEQA/P3/l/5dIJKhZsxZEIhF8fbsrjzs5ucDMzBwxMc+Ux6KijqJRIw+YmJgiMTFR+V+zZi2RnZ2Fmzdzk1mxWKz85V8ulyM5OQkymQxubnVw586tfPH4+HRWJnsA0KBBYwBQabsgJY3HzMwcDx78i0ePHpbySRUtLzkGADs7O9jbV0TFipVUpvPVresOACorbb6dxOTk5CApKRGOjlVhZmZe4PP5r/PnzyI1NRXe3h1V7tvExBQ1a9bG5csXAQC3bt3E69ev8Mknn6qMlHbt6qfyvAvj7OyCGjVqIirqiPJYdnY2Tp48gXbtvJVTJUvb3+3b+xSb7AGqzykjIwOJiYmoX78BFAoF7t7NX6+/fw+Vrxs1aqzyGTp+PApisRiffz4837V5o3ZRUUdhYWGJRo2aqDzbhg09IBaLlc+2MMnJSTA3tyj0vL9/D5W+6NGjJ4DcVWILY2hoqEz2ZDIZkpISYWxsgqpVnXH79pvnUJrPubm5ORITE4st9z44wkdEREQfpGn9PQo9JzEQF3ne3MSwyPM2FkZFnq9sa1rk+ZKqWLGSytdmZmawta2Qb/qkmZkZUlJSlF8/ffoY9+/fRffuqqNTeRITXyv//4ED+7BlywY8evQQOTk5b+6hskOx8VhY5P7CnJKSXOR9lDSe/v0/w4UL59C/fyBcXKqhefOW8PHprEzG3oWxsYkyzjympmb53vHKS6zevpesrEyEhKxFePheJCTEK0eXACA1NbXYtp8+fQwAGDMmf+ICQDnambcoh6Ojk8p5fX19ZZnieHt3xOrVf+LVq5ewsbHF2bOnkJaWlm9RktL0d5UqJVu9Mi4uDqtWLcfJkyfyfRb++5xy+8NS5Zi5ubnKdTExz2BnZ19ksvv06WMkJyeV6DNemLf787+qVlXtC0tLK5ibWyAurvCtN+RyObZv34zQ0B2IjY2BTCZ76/o391yaz7lCgSKnppYFJnxEREREAhGJ8q84WfjqiG9+eVUoFGjRwhN9+vQvsGS1arnvQB06FI5Zs2aiTZt26Nt3IKytbSASibBhw1o8e/Y033ViccFtF/WLc2nicXGphk2bduLUqZM4e/YUjh6NwLZtmzFs2CjldLjSKizmwp7j27eycOE8hIfvRc+efeHuXv//ExA9zJz5dbH3DAByeW6Z77//GVZW1vnOSySS4m+ghLy9OyI4eBmOHYtEjx49ERV1FLa2tmjU6M0fHkrb3yWJTyaTYdKkL5CSkoz+/QfB2dkFRkbGePEiAbNmzcz3nArrj9KSyxWwta2Ab7/9ocDzee/LFcbCwrLYP1SU1vr1q7Fy5XJ06/Yxhg0bBQsLS+jp6eG33+ZDLpcry5Xmc56SklyiUdb3wYSPiHSCDECWVF5sufdR/K8GRERlo0oVR2RlZaJZsxZFljt27CiqVHHA7NnzVEYRVq8u243NSxoPABgbG6NDh47o0KEjcnJy8O23U7FmTTD69RsEQ0NDqHmwQ0Xee3pvv0eWlZWVb9SqsJgcHBwB5L7/9fa7ff9VqVJlALmjVm8naDk5OYiJiUHNmrWKjbVqVSfUqlUbkZGH0b27P06ePIEuXbqpJLbq6O9//72HJ08e45tvZqJLlzdTjc+fP/POdTo4OODChbNISUmBuXnBi4g5ODji8uULaNiwMQwNDUvdhrOzC/7663ih5588eYyGDRsrv05KSkRKSjIqVqxc6DXHjkXCw6NpvnchU1NT8yVtxX3O88TGxsDV1a2Ud1c6fIePiHRCllSOiDMP1fqfTM6Uj4g0o337Drh69TIuXbqQ71xiYqJy1CUvGXh7FOaff27gxo3rgsSTlJSock5fXx/VqtWAQqFATo4UAGBkZAwAKlNY1aWgEdadO7eqTNV7O6b/JoLNm7eEqakpQkLWqEyfzJP3bpabW11YWVlj9+5dKuXCw/ciNbXk9+nt3QnXrl3B/v17kJ6eBm/vTv+5n7Lv74KekUKhwPbtW965zrZt20Mmk2Ht2pUF1g3kfqakUik2bFibr0x2djbS0oqecluvXn0kJr5GXFxcgefDwlT7Yteu7QBQ5EIwIpEo34hmZOQRJCTEqxwryeccyP08PXv2FO7uDYu8l/fFET4iIh3jp/cAlXzffSNaIhJev36D8NdfxzFlyjh06/YxatasjfT0NNy7dxfHjkUiIuI49PX14enZBsePR+Hrr79Eq1Zeym0HXFyqISMjQ+PxTJo0Fra2tnB3bwBbW1s8evQIO3duQ6tWrZX7n7m61gEALFo0Dy1atIJIJIKPT+cyi/Vtnp5eOHQoHKamZnBxqYZ//rmOCxfOqbyPBQA1a9aGWCzGxo3rkJaWCgMDAzRp0gzW1jaYNGkqZs2aiWHDBqFDh46wsLBEXFwsTp48jjZt2mHEiC+gr6+P4cNHY9682Rg/fhS8vX0QGxuD/fv3lvgdPiB3e4YVK5Zg2bLfYW9fEQ0aqCYK6uhvZ2cXODg4YunSRUhIiIepqSmOHYt8r4S8adPm8PHpjK1bN+LJk8do3rwFcnJycOnSRXh7+8DXtxs8PJrCzy8Aq1f/idu3o9GkSXOIRHp48uQxIiOP4H//+6nIEWVPTy+IxWJcuHAW3bv75zufnZ2FSZPGoF07b+W2DM2btypypLZ16zZYsyYYs2f/AHf3Bvj333uIiDiYrw9L8jkHgAsXzgIAvLzalvYRlgoTPiIiHVNBLxOVbU2LL0hE5ZaxsTGWLv0T69atxrFjR7FvXxjMzS3g7OyC0aPHKVdt7NrVD69evURY2C6cO3cGLi7V8L///YSoqCPFrnKojnj8/Xvg8OGD2LZtEzIyMmBnZ4/AwF747LOhyrratPkIPXv2xeHDBxERcQAKhUJtCd+ECV9CJBLh8OEDyMrKRv36DbFo0VJMnjxOpZy1tTWmTv0G69atwi+//ASZTIbfflsOa2sb+Pp2Q4UKdtiwYS02bFiLnJwc2NlVROPGTeDj82YEzt+/B+RyOTZtCsEff/yG6tVr4tdfFyA4eFmJ43VwcISbW13cunUT3bt/nG+xD3X0t76+Pn79dSEWLZqHkJC1kEgM0aZNe3z6aS8MHtz3neoEgO+++xE1a9bC/v17cP78GZiZmaNu3XrKhB8Apk79Bq6ubtizJxTLly+BoaEBKleuAj+/gHx7Cv6XjY0tWrb0RFTU0QITvsmTp+PAgb0IDl4OhUKOLl26Y8KEKUXWOXDg58jIyMDhwwdx9GgEatd2w9y5i7B8+e8q5UryOQdy9zds2LBxkdtHlAU9RUneSAVw7do1LFmyBJcvX0ZOTg6qVq2KwYMHo0ePN8uuHj16FEuWLMG9e/dga2uLwMBAjBo1Kt9KU8nJyZg3bx4OHz6MzMxMNGjQADNmzECdOnX+22yJvHyZqnxptrywszNHQoL6pyJQ6WhLv6T///REdevU0gUmBuqd+a2pPtHEM/Nu7ozIc4/U2sb79smTuXNwW2GFCj16olGtCoWW4/dK6WjT90p5Exf3CJUqOQsdRqH09UXIyVHv+8FUOuyT8keoPrl8+SImTvwCmzfvKtVIqia8fv0Kn37qhx9/nFOiEb6i/i0UifRga1v4iqcl+ul0/Phx9OvXDzk5OZgwYQKmTZsGT09PxMbGqpQZM2YMLC0t8d1338HHxwdLly7FnDlzVOqSy+UYMWIE9u/fjwEDBuCrr77Cy5cvMXDgQDx+/Lgk4RAR0Xs4paiMQ+f47y0REWm3xo2bwMOjKbZs2SB0KPls27YZNWrUVPt0TqAEUzpTUlIwY8YM9OnTB99++22h5ebOnYu6deti1apVymF7U1NT/Pnnnxg4cCBcXFwAAAcPHsTly5exdOlS+Pjk7qvRpUsXdO7cGUuWLMHcuXPL4LaIiIiIiEjXLVy4VOgQCjRy5BiMHDlGI20VO8K3d+9eJCcnY8KECQByV5P57yzQe/fu4d69e+jdu7cy2QOAfv36QS6XIyIiQnns0KFDsLe3R4cOHZTHbGxs0KVLFxw5cgRSqRRERERERET0/ood4Tt9+jSqV6+O48ePY968eYiLi4OFhQV69+6NSZMmQSwW4+bNmwAAd3fV3eMrVqyISpUqKc8DQHR0NOrVq5fvJdP69etj69atePz4MWrUqFEW90ZERO8oJS0b6Wret1BiIEL+xb5JF5RmX0y5Au+05YlIpAcNbudGRFRuFZvwPXr0CHFxcZg+fTqGDRuGunXrIioqCsHBwcjKysI333yDhIQEAICdXf4d7+3s7BAf/2ZvioSEBLRs2TJfOXt7ewBAfHx8qRO+ol5SFJKdXcEbSZKwtKFfnr9Kh6mpRO3tGBkZwM7GRO3taKJPNPHMxGKR2tt43z6JM9SHXpYeDAzFRT7356/S8ff12ELPl4UOzZzU/vni90r59PxVOv6+9KxEZWvaypCWUfrZP+amhtDX18x2w5pqh0qOfVL+sE/ej0gkeuefAcUmfOnp6UhKSsKUKVMwYsQIAECnTp2Qnp6OzZs3Y/To0cjMzAQAlV3j80gkEpV9PzIzMwssl3csr67S4CqdVFLa0i+ZUjnS0rLU306mVO3PS1N9oolnJpOpv4337RNpdg4+wT1U6dy56HrE4nJ/LyVqg98r5VJp+kVuo0COTJ5vZlCx18kVGlkVkCtClj/sk/KHffJ+FAoF5HJ5oT8D3nuVTiMjIwBA9+7dVY77+flBKpXi+vXryjLZ2dn5rs/KylKez6uvoHJ5x94uS0REZc9SLxs2Fvy3lj4MOTI9yGV8v5+IdJdUmg2x+N23Ty824cubplmhgup+TXlfJyUlKcvkTe18W0JCgnK6Zl59b0/xzJN37O2yRERU9m4obHAu+rnQYRCVyNNXciQnv4AsJzvfonFERNpMoVAgOzsLiYkJMDOzeud6ik0V69Wrh1OnTuH58+eoWrWq8nhcXByA3BU2K1asCAC4ceMG6tWrpyzz/PlzxMXFqWyo7ubmhsuXL0OhUKhMz7h27RpMTEzg5OT0zjdDRETFO6+oCMmlZ2hep6LQoRAVKyVThAfxMjjmJEBfrCjxQiwZqfoQaWDVFpFIBLmcU9XKE/ZJ+cM+eXdisT7Mza1hbGz6znUUm/D5+voiODgYO3bswKRJkwDkZpvbt2+HiYkJGjVqBDMzM1SvXh1bt25FYGCgcmuGzZs3QyQSoVOnTir1HTp0CEePHlXuw/fq1SscPHgQHTp0gIGBwTvfDBEREWmflEwRomMAlGLdzU4tq8LEQP2LRGjTu5Xagn1S/rBPhFVswufu7o6AgACsWLECL1++RN26dXH8+HGcPHkSX331FczMcl8QnDp1KkaPHo2hQ4eia9euuHPnDjZu3IjevXujWrVqyvo6d+6MRo0aYerUqRgyZAisra2xefNmyOVyjBs3Tn13SkREREREpGNK9PbfTz/9hMqVK2P37t3YvXs3HB0d8cMPP6BPnz7KMu3bt8eSJUuwZMkS/PTTT7CxscHo0aPxxRdfqNQlFovx559/Yu7cuQgJCUFWVhbq16+PX3/9Fc7OzmV7d0RERERERDqsRAmfoaEhJk6ciIkTJxZZzsfHRzlNsyiWlpaYNWsWZs2aVaIgiYiIiIiIqPTefX1PIiL6IPXSuwuHT7oKHQYRERFpALe8JyLSMaZ6OTA3MRQ6DCIiItIAJnxERDrmsqICTl6LFToMIiIi0gAmfEREOuaKwg5/X2fCR0REpAuY8BEREREREWkpJnxERERERERaiqt0EhEREQBABiBLKld7Owq1t0BERHmY8BERERGA3GQv4sxDtbfj3dxZ7W0QEVEuJnxERDqmv95tOPbyEzoMIiIi0gC+w0dEpGMM9eSQGIiFDoOIiIg0gAkfEZGOOaewR+Slp0KHQURERBrAhI+ISMf8o7DF+eh4ocMgIiIiDWDCR0REREREpKWY8BEREREREWkpJnxERERERERaigkfERERERGRluI+fEREOuZzUTSq9g8QOgwiIiLSAI7wERERERERaSkmfEREOuZvRSUcPPtY6DCIiIhIA5jwERHpmDsKa1y990LoMIiIiEgDmPARERERERFpKS7aQkSCSknLRrpUrvZ2FGpvgYiIiKj8YcJHRIJKz8pBxJmHam/Hu7mz2tsgIiIiKm84pZOISMfoQw4DA/7zT0REpAs4wkdEpGMGim6jaq8eQodBREREGsA/8RIREREREWkpJnxERDrmmKIK9vz9QOgwiIiISAOY8BER6ZgHCktEP3wtdBhERESkAXyHj6iMyQBkqXmbAW4xQEREREQlwYSPqIxlSeVq32aAWwwQERERUUlwSicREREREZGW4ggfEZGOMUYOTIwNhA6DiIiINIAJHxGRjukjuouqPQKFDoOIiIg0gAkfEREJQiTSQzoXOCIiIlIrJnxERDrmsLwqLI7dR2C7GoLGkZ0jR+S5R2ptgwscERGRrmPCR0SkY57CDJJnSUKHQURERBrAVTqJiIiIiIi0FBM+IiIiIiIiLcWEj4iIiIiISEvxHT4iIh1jgWyYWEiEDoOIiIg0gAkfEZGO+VR0H1X9egkdBhEREWkAEz4qlgxAlpr3ypIYiCBWawtERB+2lLRs7ltIRESlxoSPipUllSPizEO1ttGppQtMDPhKKZEmHJA7wezIHfTzqS10KFQK6Vk5av+3mPsWEhFpHyZ8REQ6Jg6mkDxPFToMIiIi0gAOqRAREREREWkpJnxERERERERaigkfERERERGRlnqnhC84OBiurq7w9/fPd+7SpUvo27cvGjZsiNatW+Pnn39GRkZGvnLZ2dmYN28evLy80KBBA/Tq1QunT59+l3CIiKgUbJGJijYmQodBREREGlDqRVsSEhKwbNkymJjk/2UhOjoagwcPRs2aNTF9+nTExcVh9erVePr0KZYvX65Sdvr06YiIiMCgQYPg7OyM0NBQDB8+HCEhIWjcuPG73xERkRYTifTea2l+mQLopvcA9j69iqxHzPX5iYiItEKpE7758+fD3d0dCoUCycnJKucWLFgAKysrhISEwNTUFADg6OiIb7/9FqdPn0arVq0AANeuXcP+/fsxY8YMDB48GAAQEBCA7t27IygoCBs3bnzP2yIi0k7ZOXJEnnv0ztc7J2cCAK4Us7x/x1bV3rkNIiIiKj9KNaXz2rVr2LNnD2bMmJHvXGpqKk6dOoWAgABlsgcA/v7+MDExwYEDB5THDh48CAMDA/Ts2VN5TCKRIDAwEBcvXkR8fPy73AsREZVApKQWrt57IXQYREREpAElTvgUCgV++uknBAQEoE6dOvnO3759Gzk5OXB3d1c5bmhoiDp16iA6Olp5LDo6GtWqVVNJDAGgQYMGUCgUKmWJiKhsJYpMkJohFToMIiIi0oAST+ncvXs37t27h6VLlxZ4PiEhAQBgZ2eX75ydnR2uXLmiUrZixYoFlgNQ6hE+W1uzUpXXFDs7c6FDKBPPX6XD1FSi1jaMjAxgp6FFJNTdL5p4XmKxSO1tAJrpF008L0Azz+xDaEMsFgF6JaunvN9LeWkD0J7vFU09L020o00/V6j02CflD/tEOCVK+FJTUzF//nyMGDEC9vb2BZbJzMx9L8TQ0DDfOYlEojyfV9bAwKDAcgCQlZVVkrCUXr5MhVxevlYYsLMzR0JCitBhlIlMqRxpaaXrk1K3kSnVyPPSRL9o4nnJZOpvA9BQv4jFGrkXTTyzD6ENmUwOKEpWT3m/l/LSBqA93yuael6aaEebfq5Q6bBPyh/2iXqJRHpFDoCVaErnsmXLYGBggM8//7zQMkZGRgByt1v4r6ysLOX5vLJSaf7pRHmJXl7iR0RERERERO+u2BG++Ph4rFu3DhMmTMCLF29e8s/KyoJUKsXTp09hbm6unI6ZN7XzbQkJCSojg3Z2dgVO28y7trBRRCIien8V5KlIMa0qdBhERESkAcUmfC9fvoRUKkVQUBCCgoLyne/QoQOGDx+OkSNHQl9fHzdu3ECnTp2U57OzsxEdHQ0/Pz/lMTc3N4SEhCAtLU1l4ZarV68qzxOpQ0pa9nvtYVYS5WtyMVF+bbPu41F1L6HDICIiIg0oNuFzdHQscKGWRYsWIT09HV9//TVcXFxgbm6OVq1aISwsDCNHjlQmcmFhYUhPT4evr6/yWl9fX6xevRrbt29X7sOXnZ2NXbt2wcPDo8AFXYjKQnpWDiKK2X/sfXk3d1Zr/UREREREJVVswmdubg4fH598x9etWwexWKxybtKkSejTpw8GDhyInj17Ii4uDmvWrEHbtm3h6empLNewYUP4+voiKCgICQkJcHJyQmhoKGJiYjBnzpwyujUiIipIhJEbUu8koEnt/KsqExERkXYp8bYMJVGvXj2sWbMGQUFBmDNnDszMzNCrVy9Mnjw5X9m5c+di0aJFCAsLQ1JSElxdXfHnn3+iSZMmZRkSERH9R6qeBJlZOUKHQURERBrwzglfSEhIgcebNm2KLVu2FHu9RCLBtGnTMG3atHcNgYiIiIiIiIpQom0ZiIiIiIiI6MPDhI+IiIiIiEhLMeEjItIxlWTJsDY3EjoMIiIi0gAmfEREOsYz+wHqulgLHQYRERFpABM+IiIiIiIiLcWEj4hIx4Qb1cX5W/FCh0FEREQawISPiEjHZOoZIFsqEzoMIiIi0gAmfERERERERFqKCR8REREREZGWYsJHRERERESkpZjwERHpGEfZa1SwMhY6DCIiItIAJnxERDqmefZjuFa1EjoMIiIi0gAmfERERERERFqKCR8RkY4JM66PM/88FzoMIiIi0gAmfEREOkYGEWRyudBhEBERkQYw4SMiIiIiItJS+kIHQERE9KETifSQLlXvqKlYodbqiYhISzHhIyIiek/ZOXJEnnuk1jY6tqqm1vqJiEg7cUonEZGOccl5iYo2JkKHQURERBrAhI+ISMd4SJ+ipoOl0GEQERGRBjDhIyIiIiIi0lJM+IiIdMwu44b4+3qs0GEQERGRBjDhIyIiIiIi0lJM+IiIiIiIiLQUEz4iIiIiIiItxYSPiIiIiIhISzHhIyLSMTVzElClgqnQYRAREZEG6AsdABGVXyKRHtKlcrW2IVaotXoqQANpDB5VthA6DCIiItIAJnxEVKjsHDkizz1SaxsdW1VTa/2UnxQi5Mjk0BdzkgcREZG24097IiIds9e4Ps7efC50GERERKQBTPiIiIiIiIi0FBM+IiIiIiIiLcWEj4iIiIiISEsx4SMiIiIiItJSTPiIiHRMHWkcqtqbCR0GERERaQATPiIiHVMn5zmcKpoLHQYRERFpABM+IiIdk6GnjyypTOgwiIiISAOY8BER6ZgDRvVw4Va80GEQERGRBugLHQARAIhEekiXytXejlih9iaIiIiIiMoNJnxULmTnyBF57pHa2+nYqpra2yAiIiIiKi84pZOIiIiIiEhLMeEjIiIiIiLSUkz4iIh0jLs0Bi6VLIQOg4iIiDSA7/AREemY2jkJkNiZCh0GERERaQBH+IiIdEyKngQZWTlCh0FEREQawISPiEjHHDZyw6U7CUKHQURERBrAhI+IiIiIiEhLMeEjIiIiIiLSUsUu2nLt2jWEhobi7NmziImJgZWVFRo3boyJEyfC2dlZpeylS5cwb9483Lx5E2ZmZujSpQumTJkCY2NjlXLZ2dlYvHgxwsLCkJycDDc3N0yaNAmtWrUq27sjIiIi+sDJAGRJ5WptQ2IgglitLRCRUIpN+FauXIlLly7B19cXrq6uSEhIwMaNGxEQEIAdO3agRo0aAIDo6GgMHjwYNWvWxPTp0xEXF4fVq1fj6dOnWL58uUqd06dPR0REBAYNGgRnZ2eEhoZi+PDhCAkJQePGjdVzp0REREQfoCypHBFnHqq1jU4tXWBiwIlfRNqo2IRv8ODBCAoKgqGhofJY165d4efnh+DgYPzyyy8AgAULFsDKygohISEwNc1d7tvR0RHffvstTp8+rRy9u3btGvbv348ZM2Zg8ODBAICAgAB0794dQUFB2LhxY1nfIxERvaVx9lPEO9QWOgwiIiLSgGL/lOPh4aGS7AGAi4sLatWqhfv37wMAUlNTcerUKQQEBCiTPQDw9/eHiYkJDhw4oDx28OBBGBgYoGfPnspjEokEgYGBuHjxIuLj49/7poiIqHDVZC9RycZE6DCIiIhIA95p7F6hUODFixewtrYGANy+fRs5OTlwd3dXKWdoaIg6deogOjpaeSw6OhrVqlVTSQwBoEGDBlAoFCpliYio7L3WM0ZqulToMIiIiEgDip3SWZA9e/bg+fPnmDRpEgAgISF3Pyc7O7t8Ze3s7HDlyhXl1wkJCahYsWKB5QC80wifra1Zqa/RBDs7c6FDKBPPX6XD1FSi1jbEYpHa28ijDfeiqeelqXa05V4+hDbEYhGiJLUhffgKHZo5FVm2vN9LeWlDk+3weZWckZEB7DQ0kq3un/ea+DmsyeelCdryO5g2YZ8Ip9QJ3/379/Hjjz+iSZMm8Pf3BwBkZmYCQL6pn0DudM2883llDQwMCiwHAFlZWaUNCS9fpkIuV5T6OnWyszNHQkKK0GGUiUypHGlppe+X0pDJ1N9GHm24F009L021oy338iG0IZPJAUXJ6inv91Je2tBkO3xeJZeZKdXIz2FN/LzXxM9hTT0vTdCm38G0BftEvUQivSIHwEo1pTMhIQEjR46EpaUlFi9eDJEo93IjIyMAudst/FdWVpbyfF5ZqTT/VKK8RC8v8SMiIiIiIqL3U+IRvpSUFAwfPhwpKSnYvHmzyvTNvP+fN7XzbQkJCbC3t1cpW9C0zbxr3y5LRERERERE765EI3xZWVkYNWoUHj58iBUrVqB69eoq52vXrg19fX3cuHFD5Xh2djaio6NRp04d5TE3Nzc8ePAAaWlpKmWvXr2qPE9ERERERETvr9iETyaTYeLEibhy5QoWL16MRo0a5Stjbm6OVq1aISwsTCWRCwsLQ3p6Onx9fZXHfH19IZVKsX37duWx7Oxs7Nq1Cx4eHgUu6EJERGWnafYj1K5qJXQYREREpAHFTun85ZdfEBkZifbt2yMxMRFhYWHKc6ampvDx8QEATJo0CX369MHAgQPRs2dPxMXFYc2aNWjbti08PT2V1zRs2BC+vr4ICgpCQkICnJycEBoaipiYGMyZM0cNt0hERG9zkiVCYWUsdBhERESkAcUmfLdu3QIAREVFISoqSuWcg4ODMuGrV68e1qxZg6CgIMyZMwdmZmbo1asXJk+enK/OuXPnYtGiRQgLC0NSUhJcXV3x559/okmTJmVxT0REVIQEkSmSUrNgacZFsojeV0paNtKlcrW2Ub7WISeiD02xCV9ISEiJK2vatCm2bNlSbDmJRIJp06Zh2rRpJa6biIjKxl+Smsh88Aqt61cWOhSiD156Vg4izjxUaxvezZ3VWj8RabdSbctAREREREREHw4mfERERERERFqKCR8REREREZGWKvHG60RERESknUQiPbUvPgMAEgMRxGpvhYjexoSPiEjHtMp6gDhnN6HDIKJyJDtHjshzj9TeTqeWLjAx4AQzIk1iwkdEpGMqy5ORbWEkdBhERESkAfwTCxGRjokVWeBVcqbQYRAREZEGMOEjItIxpyXVEP3otdBhEBERkQYw4SMiIiIiItJSTPiIiIiIiIi0FBM+IiIiIiIiLcWEj4iIiIiISEsx4SMi0jFtsu7BvZqN0GEQERGRBnAfvg+cDECWVK7WNhRqrZ2INM1OnoZ0M4nQYRAREZEGMOH7wGVJ5Yg481CtbXg3d1Zr/USkWY/FVkhIzICdlbHQoRAREZGacUonEZGOuWDojDtPEoUOg4iIiDSACR8REREREZGW4pROIiIi0joikR7S1fyOOwCI+aI7EZVzTPiIiIhI62TnyBF57pHa2+nYqpra2yAieh+c0klERERERKSlmPAREemY9pl30LBGBaHDICIiIg1gwkdEpGOsFRkwMzEQOgwiIiLSACZ8REQ65oHYFnGv0oUOg4iIiDSACR8RkY65bOiI+8+ShA6DiIiINIAJHxERERERkZZiwkdERERERKSluA8fEREREVE5JAOQJZWrtQ2JgQhitbZAQmPCR0REREQaIRLpIV3NCUxKWrZa69ekLKkcEWceqrWNTi1dYGLASX/ajAkfEZGO6Zh5C89quwsdBhHpoOwcOSLPPVJrGx+3q8V3lojewu8HIiIdY67IgrGEf+8jIiLSBfyJryYpadlqn7IAAAq1t0BE2uaOvh1eJKTBwc5U6FCIiIhIzZjwqUl6Vo7a51wDgHdzZ7W3QUTa5YZBFWTGJTPhIyJ6R5pYTAXgH/apbDDhIyIiIiIqBU0spgLwD/tUNvgOHxERERERkZZiwkdERERERKSlmPARERERERFpKSZ8REQ6pkvmP2jqZi90GERERKQBXLSFiEjHGCtyIDEQCx0GEZFayOUKZKp5BU2unkkfEiZ8REQ6Jlq/Il48T4FTRXOhQyEiKnNZUhkOq3kFTa6eSR8STukkItIx0QaV8CQ+VegwiIiISAOY8BEREREREWkpJnxERERERERaigkfERERERGRlmLCR0REREREpKW4SicRkY7xy7iOx3UbCB0GERERaQBH+IiIdIwB5NAX859/IiIiXcCf+EREOuaaQRU8iE0WOgwiIiLSAMESvuzsbMybNw9eXl5o0KABevXqhdOnTwsVDhGRzrinb4eYF2lCh0FEREQaINg7fNOnT0dERAQGDRoEZ2dnhIaGYvjw4QgJCUHjxo2FCouIiIiIiD5AMgBZUrla25AYiCBWawtlT5CE79q1a9i/fz9mzJiBwYMHAwACAgLQvXt3BAUFYePGjUKERUREREREH6gsqRwRZx6qtY1OLV1gYvBhvRUnSLQHDx6EgYEBevbsqTwmkUgQGBiIixcvIj4+XoiwiIiIiIiItIogI3zR0dGoVq0aTE1NVY43aNAACoUC0dHRsLe3L3F9IpFeWYf4/kR6MDMxUHsz+mKR2tvRljZy21F/v2jX89LEvfB7RZNtGNlaw9bIFFmWxkXWw++V8teOtvSJptrhz5Xy14am2tGWPtFUO2KxnkZ+l9ZEG2IN9L2mnldpFBePnkKhUGgoFqXu3bujYsWKWLVqlcrxe/fuoVu3bvj5559VRv+IiIiIiIio9ASZ0pmZmQkDg/zZt0QiAQBkZWVpOiQiIiIiIiKtI0jCZ2RkBKlUmu94XqKXl/gRERERERHRuxMk4bOzsytwYZaEhAQAKNX7e0RERERERFQwQRI+Nzc3PHjwAGlpqhv/Xr16VXmeiIiIiIiI3o8gCZ+vry+kUim2b9+uPJadnY1du3bBw8MDFStWFCIsIiIiIiIirSLItgwNGzaEr68vgoKCkJCQACcnJ4SGhiImJgZz5swRIiQiIiIiIiKtI8i2DEDuAi2LFi3C3r17kZSUBFdXV0yePBmenp5ChENERERERKR1BEv4iIiIiIiISL0EeYePiIiIiIiI1I8JHxERERERkZYSZNEWbZWdnY3FixcjLCwMycnJcHNzw6RJk9CqVSuhQ9NZ8fHxWL9+Pa5evYobN24gPT0d69evR4sWLYQOTWddu3YNoaGhOHv2LGJiYmBlZYXGjRtj4sSJcHZ2Fjo8nXX9+nUsX74cN2/exMuXL2Fubg43NzeMGTMGHh4eQodH/y84OBhBQUFwc3NDWFiY0OHopLNnz2LQoEEFngsPD0eNGjU0HBHluXbtGpYsWYLLly8jJycHVatWxeDBg9GjRw+hQ9NJ06dPR2hoaKHnT5w4wZX5NYQJXxmaPn06IiIiMGjQIDg7OyM0NBTDhw9HSEgIGjduLHR4OunBgwcIDg6Gs7MzXF1dcfnyZaFD0nkrV67EpUuX4OvrC1dXVyQkJGDjxo0ICAjAjh07+MuSQJ48eQKZTIaePXvCzs4OKSkp2Lt3LwYMGIDg4GC0bt1a6BB1XkJCApYtWwYTExOhQyEAn332GerVq6dyjL+8Cuf48eMYM2YMmjdvjgkTJkBfXx8PHz5EbGys0KHprN69e+cb9FAoFJg5cyYcHBz4/aJBXLSljFy7dg09e/bEjBkzMHjwYAC5K5F2794d9vb22Lhxo7AB6qjU1FRIpVJYW1vjyJEjGDNmDEf4BHbp0iW4u7vD0NBQeezhw4fw8/NDt27d8MsvvwgYHb0tIyMDPj4+cHd3x4oVK4QOR+dNnz4dMTExUCgUSE5O5gifQPJG+JYuXQofHx+hwyEAKSkp6Ny5M7p27Ypvv/1W6HCoCBcuXED//v0xadIkjBo1SuhwdAbf4SsjBw8ehIGBAXr27Kk8JpFIEBgYiIsXLyI+Pl7A6HSXmZkZrK2thQ6D3uLh4aGS7AGAi4sLatWqhfv37wsUFRXE2NgYNjY2SE5OFjoUnXft2jXs2bMHM2bMEDoUektqaipycnKEDkPn7d27F8nJyZgwYQKA3H7heEb5tG/fPujp6aF79+5Ch6JTmPCVkejoaFSrVg2mpqYqxxs0aACFQoHo6GiBIiMq/xQKBV68eMHkvBxITU3Fq1ev8O+//2LBggW4c+cO30MWmEKhwE8//YSAgADUqVNH6HDo/3311Vdo0qQJGjZsiCFDhuD27dtCh6SzTp8+jerVq+P48eP46KOP0KRJEzRv3hxBQUGQyWRCh0f/TyqV4sCBA2jcuDEcHR2FDken8B2+MpKQkFDgXGQ7OzsA4AgfURH27NmD58+fY9KkSUKHovO+/vprHDp0CABgYGCAPn36cNqNwHbv3o179+5h6dKlQodCyP2+6Ny5M9q2bQtra2vcvn0bq1evRr9+/bBjxw5Uq1ZN6BB1zqNHjxAXF4fp06dj2LBhqFu3LqKiohAcHIysrCx88803QodIAE6ePInExET4+fkJHYrOYcJXRjIzM2FgYJDvuEQiAZD7Ph8R5Xf//n38+OOPaNKkCfz9/YUOR+eNGTMGvXv3RlxcHMLCwpCdnQ2pVJpvGi5pRmpqKubPn48RI0bA3t5e6HAIudPS3165tkOHDvD29sann36KJUuWYP78+QJGp5vS09ORlJSEKVOmYMSIEQCATp06IT09HZs3b8bo0aNhY2MjcJS0b98+GBgYoEuXLkKHonM4pbOMGBkZQSqV5juel+jlJX5E9EZCQgJGjhwJS0tLLF68GCIR/0kSmqurK1q3bo1PP/0Uq1atwj///MP3xgS0bNkyGBgY4PPPPxc6FCqCm5sbWrVqhTNnzggdik4yMjICgHzvhfn5+UEqleL69etChEVvSUtLw9GjR+Hl5cXXNwTA367KiJ2dXYHTNhMSEgCAf5kl+o+UlBQMHz4cKSkpWLlypXL6M5UfBgYG6NChAyIiIpCZmSl0ODonPj4e69atQ79+/fDixQs8ffoUT58+RVZWFqRSKZ4+fYqkpCShw6T/V7lyZfaHQPJ+flSoUEHleN7X7BfhHTlyBBkZGZzOKRAmfGXEzc0NDx48QFpamsrxq1evKs8TUa6srCyMGjUKDx8+xIoVK1C9enWhQ6JCZGZmQqFQ5Pu3jdTv5cuXkEqlCAoKQocOHZT/Xb16Fffv30eHDh0QHBwsdJj0/548ecKRC4Hk7Yf4/PlzleNxcXEAwOmc5cDevXthYmICb29voUPRSUz4yoivry+kUim2b9+uPJadnY1du3bBw8ODm0sS/T+ZTIaJEyfiypUrWLx4MRo1aiR0SATg1atX+Y6lpqbi0KFDqFy5MmxtbQWISrc5Ojpi6dKl+f6rVasWHBwcsHTpUgQEBAgdps4p6HvlwoULOHv2LLy8vASIiHx9fQEAO3bsUB5TKBTYvn07TExM+HNGYK9evcLp06fRsWNHGBsbCx2OTuKiLWWkYcOG8PX1RVBQEBISEuDk5ITQ0FDExMRgzpw5Qoen0/744w8AUO7xFhYWhosXL8LCwgIDBgwQMjSd9MsvvyAyMhLt27dHYmKiyubRpqam3MhYIBMnToREIkHjxo1hZ2eH2NhY7Nq1C3FxcViwYIHQ4ekkc3PzAr8f1q1bB7FYzO8VgUycOBHGxsZo3LgxrK2tcffuXWzduhXW1tYYN26c0OHpJHd3dwQEBGDFihV4+fIl6tati+PHj+PkyZP46quvYGZmJnSIOi08PBw5OTmczikgPQV3piwzWVlZWLRoEfbu3YukpCS4urpi8uTJ8PT0FDo0nebq6lrgcQcHB0RGRmo4Gho4cCDOnTtX4Dn2iXB27NiBsLAw3Lt3D8nJyTA3N0ejRo0wZMgQNG/eXOjw6C0DBw5EcnKyyh9LSHPWr1+PvXv34vHjx0hNTYWNjQ28vLwwbtw4VKlSRejwdFZ2djb++OMP7N69Gy9evICjoyMGDx6MPn36CB2azuvduzeePHmCv/76C2KxWOhwdBITPiIiIiIiIi3Fd/iIiIiIiIi0FBM+IiIiIiIiLcWEj4iIiIiISEsx4SMiIiIiItJSTPiIiIiIiIi0FBM+IiIiIiIiLcWEj4iIiIiISEsx4SMiIiqCq6srpk+fLnQYRERE74QJHxERERERkZbSUygUCqGDICIiKq+ysrIgEolgYGAgdChERESlxoSPiIjoPzIzM6Gvrw99fX2hQyEiInovnNJJREQfhOPHj8PV1RXr168v8Hzv3r3RsmVLSKVSXLt2DdOnT0fnzp3RsGFDNG7cGH369MHhw4fzXTd9+nS4urri1atXmDFjBjw9PdGoUSPExcUBKPgdvvDwcIwaNQrt2rWDu7s7WrRogS+++AK3bt3KV7+3tzcGDhyI+/fvY8SIEWjcuDGaNGmC8ePHIyEhIV/51NRULFy4EF26dEH9+vXRokUL9O3bF/v371cpFx8fj++//14Zg5eXF7777ju8fPmyxM+UiIi0H/90SUREHwQvLy/Y2dlh9+7dGDRokMq5hw8f4sqVKxg4cCAMDAxw+PBh/Pvvv/D19YWDgwMSExMRGhqKsWPHIigoCH5+fvnq//zzz1GhQgV88cUXSE9Ph4mJSaGxbNiwAVZWVujVqxfs7Ozw+PFjbNu2DX379kVoaChcXFxUyj9//hyDBg2Cj48Ppk6dilu3bmHr1q1ITU3F6tWrleWSk5PRr18/3L17F507d0bfvn0hl8tx8+ZNREVFoVu3bgCAmJgY9O7dG1KpFIGBgXBycsKjR4+wefNmnD17Fjt37oS5ufl7PG0iItIWTPiIiOiDIBaL4efnh9WrV+PevXuoWbOm8tzu3bsBAJ988gkAYPTo0ZgyZYrK9QMHDkRAQACWLVtWYMJXq1YtBAUFlSiWlStX5ksIAwIC4O/vj7Vr12LmzJkq5x49eoSFCxeia9euymMikQibNm3Cv//+i+rVqwMAFixYgLt37+LHH39E7969VeqQy+XK///TTz8hJycHu3fvRqVKlZTHfX190bt3b6xduxbjxo0r0b0QEZF245ROIiL6YOQldHkJHgAoFArs2bMHtWvXRr169QBAJRnLyMjA69evkZGRgZYtW+L+/ftITU3NV/fQoUNLHEde/QqFAqmpqXj16hWsra1RrVo1XLt2LV95e3t7lWQPAFq2bAkgNxkEchO68PBw1KhRI1+yB+QmiACQkpKCY8eOwdvbG4aGhnj16pXyPwcHBzg5OeHvv/8u8b0QEZF24wgfERF9MPKSur1792Ly5MkQiUQ4f/48nj17hq+++kpZ7uXLl1i0aBGOHj1a4DttycnJMDMzUzn232mYRbl58yYWL16Mc+fOIT09XeWco6NjvvJVq1bNd8zKygoAkJiYCAB4/fo1kpKS0KZNmyLbfvDgAeRyOXbs2IEdO3YUWKag9oiISDcx4SMiog+Kv78/Zs+ejTNnzsDT0xO7d++GWCzGxx9/DCB31G3IkCG4f/8+Bg0aBHd3d5ibm0MsFmPnzp3Yt2+fyvTIPMbGxiVqPyYmBv3794eZmRlGjx6N6tWrw9jYGHp6epg9e3a+BBDInY5amNIulp1X/uOPP1aOeP6XRCIpVZ1ERKS9mPAREdEHxc/PD/PmzcPu3bvh4eGBQ4cOwdPTE/b29gCA27dv49atWxgzZgzGjx+vcu327dvfu/3Dhw8jPT0dy5YtU07LzJOYmAhDQ8N3qtfa2hqWlpYFrvT5NicnJ+jp6UEqlcLT0/Od2iIiIt3Bd/iIiOiDYmNjgzZt2uDw4cPYu3cvUlNTVUa68t51++/I2Z07dwrclqG08kbr/lv/tm3bCtxmoaREIhG6deuGe/fuFZiY5rVnbW2Njz76CIcPH8aVK1cKLPfq1at3joOIiLQLR/iIiOiD88knnyAyMhK//PILzM3N4ePjozxXo0YN1KpVCytXrkRmZiaqVauGBw8eYOvWrahduzb++eef92q7bdu2MDY2xtSpUzFgwABYWFjg0qVLOHHiBJycnCCTyd657okTJ+LMmTP49ttv8ffff6NJkyZQKBSIjo5GTk4O5s2bBwCYOXMm+vXrhwEDBsDf3x9169aFXC7HkydPcPToUQQEBHCVTiIiAsCEj4iIPkDt2rWDlZUVEhMT0bNnT5V31sRiMVasWIFff/0VoaGhyMjIQK1atfDrr7/i1q1b753wOTk5ITg4GAsWLMDy5cshFovh4eGBkJAQ/PTTT3j27Nk7121paYmtW7di+fLlOHz4MI4cOQJTU1PUqFEDAwYMUJarXLkydu7cieDgYERGRmLPnj2QSCSoXLky2rdvjy5durzXPRIRkfbQU5T2bXEiIiIiIiL6IPAdPiIiIiIiIi3FhI+IiIiIiEhLMeEjIiIiIiLSUkz4iIiIiIiItBQTPiIiIiIiIi3FhI+IiIiIiEhLMeEjIiIiIiLSUkz4iIiIiIiItBQTPiIiIiIiIi3FhI+IiIiIiEhL/R85pn0SBm/mGQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "outcomes = np.array([1, 2, 3, 4, 5, 6])\n",
    "prob = np.array([1/6, 1/6, 1/6, 1/6, 1/6, 1/6])\n",
    "\n",
    "n_repeats = 10000\n",
    "N = 5\n",
    "\n",
    "# compute exp_value & variance\n",
    "truth_exp_val = np.dot(outcomes, prob)\n",
    "truth_var = np.dot(outcomes ** 2, prob) - truth_exp_val**2\n",
    "truth_var\n",
    "\n",
    "list_est_var = list()\n",
    "for _ in range(n_repeats):\n",
    "    # draw observations\n",
    "    observations = choices(outcomes, prob, k=N)\n",
    "    \n",
    "    # compute sampel mean (needed to compute est_var_bias)\n",
    "    est_exp_val = sum(observations) / N\n",
    "    \n",
    "    # compute biased variance estimate\n",
    "    est_var = sum((np.array(observations) - est_exp_val)**2) / N\n",
    "    \n",
    "    # store\n",
    "    list_est_var.append(est_var)\n",
    "    \n",
    "\n",
    "sns.set(font_scale=1.5)\n",
    "plt.hist(list_est_var, bins=30, label='estimated variance (bias)', alpha=.5)\n",
    "plt.axvline(truth_var, label='ground truth variance', color='r')\n",
    "plt.axvline(np.mean(list_est_var), label='mean estimated variance (bias)', color='b', linestyle='--')\n",
    "plt.xlabel('variance')\n",
    "plt.legend()\n",
    "plt.gcf().set_size_inches(15, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca1cf80e",
   "metadata": {},
   "source": [
    "# Punchline:\n",
    "\n",
    "$$ \\hat{\\sigma^2}_{bias} = \\frac{1}{N}\\sum_i (x_i - \\bar{x})^2 $$\n",
    "\n",
    "consistently under estimates our variance.  We'll soon see that:\n",
    "\n",
    "$$ \\hat{\\sigma^2}_{bessel} = \\frac{1}{N-1}\\sum_i (x_i - \\bar{x})^2 $$\n",
    "\n",
    "on average, doesn't over or underestimate variance."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}