1 Line fits.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1bba0128",
   "metadata": {},
   "source": [
    "# Simple fits\n",
    "\n",
    "Here we are fitting a line from scratch.\n",
    "In the next notebook, we will do fancier fits with neural networks, but let's start with a basic problem and complicate it as we go along.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "23feddde",
   "metadata": {},
   "outputs": [],
   "source": [
    "from typing import Tuple\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb1286f0",
   "metadata": {},
   "source": [
    "We start by generating some fake dataset, which is simple enough that we can visualize the results easily. For this reason, the dataset will contain only two  variables.\n",
    "\n",
    "The simulated example data will be $f(x) = 3 x + \\epsilon$, where $\\epsilon \\sim \\mathcal{N}(\\mu=0, \\sigma=0.5)$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5d457cd8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_data(N: int) -> np.ndarray:\n",
    "    x = 2*np.random.randn(N, 1)\n",
    "    epsilon = 0.5*np.random.randn(N, 1)\n",
    "    z = 3*x + epsilon\n",
    "    return np.concatenate((x, z), axis=1).astype(np.float32)\n",
    "\n",
    "data = generate_data(N=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48433f6f",
   "metadata": {},
   "source": [
    "We can fit this line from scratch, assuming $y = f(x) = \\beta x + \\alpha + \\epsilon$, where $\\epsilon$ is a zero-mean Gaussian noise."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a170344-34ff-4af6-9a6c-daa44b95408e",
   "metadata": {},
   "source": [
    "We can do this in two ways. In both cases, we start assuming the knowledge of $f(x)$ and use the Bayes Theorem:\n",
    "\n",
    "$ p(\\alpha, \\beta|\\mathcal{D}) \\, \\propto p(\\mathcal{D}|\\alpha, \\beta) \\, p(\\alpha) \\, p(\\beta)$\n",
    "\n",
    "We assume $p(\\alpha)$ and $p(\\beta)$ are constants, so we just need to maximize the likelihood $p(\\mathcal{D}|\\alpha, \\beta)$.\n",
    "That likelihood is a Gaussian distribution, because we assumed $y$ is exactly equals to a line plus some Gaussian noise coming from $\\epsilon$.\n",
    "\n",
    "$p(\\alpha,\\beta|\\mathcal{D}) \\propto \\prod \\mathcal{N}(y_i|\\text{mean}=f(x_i), \\text{std. dev.}=\\sigma_\\epsilon)$\n",
    "\n",
    "This leads to minimizing $L = \\sum_i (y_i - \\beta x_i - \\alpha)^2$, as seen in the lecture (simply take $- \\log p$ and to maximize $p$, one has to minimize $-\\log p$).\n",
    "\n",
    "We can try finding out $\\beta$ and $\\alpha$ by taking the derivative of $L$ and setting it to zero:\n",
    "\n",
    "$\\frac{d L}{d\\beta} = 2 \\sum_i (y_i - \\beta x_i - \\alpha) x_i = 0$\n",
    "\n",
    "$\\frac{d L}{d\\alpha} = -2 \\sum_i (y_i - \\beta x_i - \\alpha) = 0$\n",
    "\n",
    "Re-arranging:\n",
    "\n",
    "$2 \\sum_i (y_i x_i - \\beta x_i^2 - \\alpha x_i) = 0$\n",
    "\n",
    "$2 \\sum_i (y_i - \\beta x_i - \\alpha) = 0$\n",
    "\n",
    "This leads to:\n",
    "\n",
    "$ \\beta \\sum_i x_i^2 = \\sum_i (y_i - \\alpha) x_i$\n",
    "\n",
    "$ \\alpha = 1/N \\sum_i y_i - \\beta x_i = \\mathbb{E}[y - \\beta x]$\n",
    "\n",
    "We can then substitute $\\alpha$ in $\\beta$:\n",
    "\n",
    "$ \\beta \\sum_i x_i^2 = \\sum_i y_i x_i - 1/N \\sum_i \\sum_j (y_j - \\beta x_j) x_i$\n",
    "\n",
    "And expand it:\n",
    "\n",
    "$ \\beta (\\sum_i x_i^2 - 1/N \\sum_i \\sum_j x_i x_j) = \\sum_i y_i x_i - 1/N \\sum_i \\sum_j y_j x_i$\n",
    "\n",
    "Taking $x_j$ from a sum that does not depend on $j$, we see that the second term is just the square of the sum over $x_i$:\n",
    "\n",
    "$ \\beta (\\sum_i x_i^2 - 1/N(\\sum_{i} x_i)^2 ) = \\sum_i y_i x_i - 1/N \\sum_i \\sum_j y_j x_i$\n",
    "\n",
    "Dividing by the number of points $N$, several sums reduce to averages, which make this much easier to interpret if we use some standard identities from statistics.\n",
    "For example $1/N \\sum_i x_i^2 - 1/N^2 (\\sum_i x_i)^2$ is a known expression for the variance of $X$. \n",
    "\n",
    "$ \\beta \\text{var}[X] = \\mathbb{E}[XY] - \\mathbb{E}[X]\\mathbb{E}[Y]$\n",
    "\n",
    "And the right-hand-side is a known expression for the covariance between $X$ and $Y$:\n",
    "\n",
    "$ \\beta \\text{var}[X] = \\text{cov}[X,Y]$\n",
    "\n",
    "This leads to:\n",
    "\n",
    "$\\beta = \\frac{\\text{cov}[X,Y]}{\\text{var}[X]}$\n",
    "\n",
    "This, together with the expression from above for $\\alpha$, allows us to easily fit a line in a few lines of code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "30205402",
   "metadata": {},
   "outputs": [],
   "source": [
    "cov = np.cov(data[:,0], data[:,1])[0,1]\n",
    "varX = np.std(data[:,0])**2\n",
    "beta = cov/varX\n",
    "alpha = np.mean(data[:,1] - beta*data[:,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "cabdfc12-3033-4a2f-96ed-4ab1ea8312ab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f62192d0e20>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAGwCAYAAABGogSnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABe4klEQVR4nO3deVSU1R8G8GeGZdiH3UFFATUVd9y3XNLEBSW1zDKX+lmupWYqlbuJWZq526aWlba5F+WaqSglkiJqLuA6CIoygDLAzPv7AxkZGGAWYGDm+ZzDOc2d971z5ZTzdN97v1ckCIIAIiIiIishNvcAiIiIiCoTww8RERFZFYYfIiIisioMP0RERGRVGH6IiIjIqjD8EBERkVVh+CEiIiKrYmvuAVRFarUat2/fhqurK0QikbmHQ0RERHoQBAEZGRmoWbMmxOKS53cYfnS4ffs2/P39zT0MIiIiMsKNGzdQu3btEt9n+NHB1dUVQP4vz83NzcyjISIisnLbtwNjxwK5ufmvPTyAbduA9u21LlMoFPD399d8j5eE4UeHgkddbm5uDD9ERETm1qUL4OQEpKcDHToAP/wAlPKEpqwlKww/REREVLXVqwds2gQcPgwsXQrY25vUHXd7ERERUdWyaxfw8KF2W3g4sGKFycEHYPghIiKiqkKpBCZOBAYNAsaPBwShQj6G4YeIiIjM79o1oGtXYO3a/Ndffw38+WeFfBTDDxEREZnXr78CrVoBf/+d/9rBAfjyS6B79wr5OIYfIiIiMg+VCnj/faB/f+D+/fy2evWA6Gjg1Vcr7GO524uIiIgq3507wEsvAQcPPmkLDwc2bgTc3Sv0oxl+iIiIqHIlJQGdOgFyef5rGxvgww+BadOASjhWiuGHiIiIKledOkDz5vnhx88vv1pz166V9vEMP0RERFTuVGoBMYlpSMnIhq+rA9oFegKApq3mwpVo7b0Q4mUfAzVqVOrYGH6IiIioXEXFyzF/dwLk6dkAgCbJl+FpB5z1b4wHD3M11/k1GY25qWqEVm724W4vIiIiKj9R8XKM3xKbH3wEAcPjovDLlnfw0dYFsElN1bo2OT0b47fEIipeXqljZPgxE5VaQPSVe9gZdwvRV+5Bpa6YKpYFRo8eDZFIBJFIBDs7O9SoUQO9e/fGV199BbVarXc/mzZtgnsFr8InIqLqSaUWMH93AgQAjjnZWLZ3OSJ/Xw2JKheyzDRMOPGj1vUF33zzdydU+PdgYXzsZQZFpwMBwE/qgLlhwQht6ldhnxsaGoqNGzdCpVLhzp07iIqKwltvvYWffvoJu3btgq0t/3UgIiLjxSSmQZ6ejaB7N7F2RyQa3b2meW9j6zAs6T662D0CAHl6NmIS09CxnleljJMzP5VMazqwkMqY+pNIJJDJZKhVqxZCQkLw7rvvYufOnfjtt9+wadMmAMDy5cvRrFkzODs7w9/fHxMmTEBmZiYA4PDhwxgzZgzS09M1s0jz5s0DAHzzzTdo06YNXF1dIZPJ8NJLLyElJaXC/ixERFT1pGRko//5v7Dr66ma4JNp74hJA2dgfq83kGtjV+q9lYXhpxIVng4sylxTfz179kSLFi3wyy+/AADEYjFWrlyJc+fOYfPmzTh48CBmzJgBAOjUqRNWrFgBNzc3yOVyyOVyTJ8+HQCQm5uLhQsX4t9//8WOHTuQlJSE0aNHV9qfg4iIzCwnB+1WLMCaXR/CJecRAOCidx0MGrkcexo/Xebtvq4OFT1CDT7nqEQF04ElMcfUHwA0atQIZ86cAQBMmTJF0x4QEIBFixZh3LhxWLt2Lezt7SGVSiESiSCTybT6eLVQGfKgoCCsXLkSbdu2RWZmJlxcXCrlz0FERGYiCEBoKPwOHdI0/dKkB957diIe2ZceakQAZNInW+ErA2d+KpG+U3qVOfUHAIIgQPS4oub+/fvxzDPPoFatWnB1dcUrr7yCe/fu4eHDh6X2cerUKYSFhaFOnTpwdXVFt27dAADXr1+v8PETEZGZiUQ423MgAEBpY4uIPpMwrf80vYIPAMwNC4aNuOIrOxdg+KlE+k7pVebUHwCcP38egYGBSEpKwoABA9C8eXP8/PPPOHXqFNasWQMAyMnJKfH+rKws9OnTB25ubvj222/x999/Y/v27WXeR0REliEqXo6BmQ3waacXMWTEx/i+ZWixYyrcnezg7qS95kcmdcC6ESEVutlHFz72qkTtAj3hJ3VAcnq2znU/5pj6O3jwIM6ePYupU6fi1KlTUKvVWLZsGcTi/Fz8ww8/aF1vb28PlUql1XbhwgXcu3cPS5Ysgb+/PwDgn3/+qZw/ABERVb6UFGDHDuD117XWs37SdYTOy72c7REd8QxsxKJiVZ8rc8anAMNPJbIRizA3LBjjt8RCBGgFoMqY+lMqlUhOTtba6h4ZGYkBAwZg5MiRiI+PR25uLlatWoWwsDAcO3YM69ev1+ojICAAmZmZOHDgAFq0aAEnJyfUqVMH9vb2WLVqFcaNG4f4+HgsXLiwQv4MRERkZseOAS+8ANy+DXh6IqZVj1LXswLAvawcnLp2Hx3reVXqmtaS8LFXJQtt6od1I0Igk2o/2qqMqb+oqCj4+fkhICAAoaGhOHToEFauXImdO3fCxsYGLVq0wPLly/Hhhx+iadOm+PbbbxEZGanVR6dOnTBu3DgMGzYMPj4+WLp0KXx8fLBp0yb8+OOPCA4OxpIlS/Dxxx9X2J+DiIjMQBCA5cuBbt3ygw8AvPceUu9n6nV7Za9nLY1IEITK21ddTSgUCkilUqSnp8PNza1CPkPXgW/mmPojIiIqU3o6MGYM8Hg9JwCge3fg++8RnWWH4Z+fKLOL78d2qPBZH32/v/nYy0xsxKIqMfVHRERUqrg44PnngcuXn7RFRAALFgC2tminFqrcetay8LEXERER6fbVV0DHjk+Cj7s7sHs3sHgx8PhIpIL1rMCT9asFzLWVvSwMP0RERFTc0qXAa68B2Y/X6rRuDcTGAgMGFLvUnOtZjcE1PzpUxpofIiKiqqbwetTaGfcQMvgZiO7eBSZMyF/sLJHofb851rNyzQ8RERHpLSpejvm7E7S2rQ8MewevNXRGi5kT9eqjuqxn5WMvIiIiC6dSC4i+cg87424h+sq9Ygdo/3bqOs5OmInMlHta7bt9myD8fgCi4uWVOdwKx5kfIiIiC6ZrRsdP6oC5YcEIbeqHg/v+gc+YkXjn1nk0uXMFE8IjNEdTFESkebvOoXewrEotWjYFZ36IiIgsVFS8HOO3xBarwJycno3xW2Lx/QdfoMWgZ9Dm1nkAQK/LMWicmlisn2SFEqsPXi7WXl0x/BAAQCQSYceOHRX+OZ999hn8/f0hFouxYsUKzJs3Dy1btqzwzyUisjaFz9wqSqRW4c2j32HY+6/D65ECAHDTzQdDRyzFed8gnf19sv8/i3n8xd1eOljibq/Ro0fjwYMHJQac5ORkeHh4QFLGSn5TKBQKeHt7Y/ny5RgyZAikUinUajWUSiW8vLz0GicREekn+so9nZWXPR6m49PdH+PppNOatoNBbTBtwDQ8cCz9O89P6oCjM3tW2cdf3O1FBpHJZBX+GdevX0dubi769+8PP78nNR9cXFwq/LOJiKyJSi3g6+ikYu0ht85j9c4PUTPjbv51IjGWdR2BdR2GQhCV/TBInp6NmMS0arGjqzR87EUAtB97JSUlQSQS4ZdffkGPHj3g5OSEFi1aIDo6Wuueo0ePomvXrnB0dIS/vz/efPNNZGVl6ex/06ZNaNasGQAgKCgIIpEISUlJWo+95s2bh82bN2Pnzp0QiUQQiUQ4fPhwRf2RiYgsUlS8HJ2XHMRv8cla7U3uXMG272Zpgk+qkztGDFuItR1f0Cv4FKhKB5Qaq9qFnyNHjiAsLAw1a9bUuU5l9OjRmi/Ogp/Q0FDzDLaae++99zB9+nTExcXhqaeewvDhw5GXlwcAuHLlCkJDQzFkyBCcOXMG27Ztw9GjRzFp0iSdfQ0bNgz79+8HAMTExEAul8Pf31/rmunTp+OFF15AaGgo5HI55HI5OnXqVLF/SCIiC1KwwDlZUTygnPMNwv767QEAJ2s3Qf/Rn+JEQAuDP8PX1aHsi6q4avfYKysrCy1atMCrr76KwYMH67wmNDQUGzdu1LyuyHUsGsuX5/+UJSQE2LVLu23gwPyS4WWZNi3/p5JMnz4d/fv3BwDMnz8fTZo0weXLl9GoUSNERkbi5ZdfxpQpUwAADRo0wMqVK9GtWzesW7cODg7a/3E4Ojpq1vX4+PjofMzm4uICR0dHKJXKSnkMR0RkSUpb4AwAEIkwo99bOOPXAJ+1GwyV2AY9G/rg0IVUACj5voLbUfUOKDVWtQs/ffv2Rd++fUu9RiKRVP6Xp0IB3LpV9nVFZjsAAKmp+t2rUBg+LhM0b95c888Fa3RSUlLQqFEj/Pvvvzhz5gy+/fZbzTWCIECtViMxMRGNGzeu1LESEVm7mMQ0rS3tQ84eQLKrF44FtNS0ZUicsa7D85rXY7vWwwtt/IvVASqqqh5QaqxqF370cfjwYfj6+sLDwwM9e/bEokWLNLMOuiiVSiiVSs1rhTEhw80NqFWr7Ot8fHS36XNvJe88s7Oz0/yz6HHBK7VaDQDIzMzEG2+8gTfffLPYfXXq1KmcARIRkUbBWhxJXg7m7VuP4Wf+wD1HN/QfvRLJbt7FrveTPjl7q3ewTHMmV9Ldh/g+5rrWozNZoaKIlsDiwk9oaCgGDx6MwMBAXLlyBe+++y769u2L6Oho2NjY6LwnMjIS8+fPN+2DTXkkVfQxWDUQEhKChIQE1K9fv1z7tbe3h0qlKtc+iYisga+rA+rcl2Pdjkg0SbkKAPB6pED/C3/hy3bPaV0rgvYsTtEzuSb1rG/WA0ormsWFnxdffFHzz82aNUPz5s1Rr149HD58GM8884zOeyIiIjCtUHBRKBTFFuNagvT0dMTFxWm1eXl5GfVnnTlzJjp06IBJkybhf//7H5ydnZGQkIB9+/Zh9erVRo8xICAAv//+Oy5evAgvLy9IpVKtGSgiItKt3b9HsPfrKXDNzt91+8hWgvf6TMAvTbW/+2RuEswb2KTUWZzqckCpsSwu/BQVFBQEb29vXL58ucTwI5FIKmdRtJkdPnwYrVq10mp77bXX8MUXXxjcV/PmzfHnn3/ivffeQ9euXSEIAurVq4dhw4aZNMaxY8fi8OHDaNOmDTIzM3Ho0CF0797dpD6JiCxabi7w7ruw+fhjuD5uuupZC+PDI3DRJ0Dr0qm9nsKknvUtahbHGNW6wrNIJML27dsRHh5e4jU3b95EnTp1sGPHDgwcOFCvfi2xwjMREVVvKrVQ/FFUshwYNgw4elRznbxPGEa0/x+uKJ8s9fCzsDU7JbHYCs+ZmZm4fPnJ4WqJiYmIi4uDp6cnPD09MX/+fAwZMgQymQxXrlzBjBkzUL9+ffTp08eMoyYiIjKerpPZa7vY4ffPx8H5+uODSO3sgGXL4DdpEv4QYNFrdkxV7WZ+Dh8+jB49ehRrHzVqFNatW4fw8HCcPn0aDx48QM2aNfHss89i4cKFqFGjht6fwZkfIiKqKgoKFxb9shYB6HfhKNbsXJJfRuWHH4AOHcwxxCpD3+/vahd+KgPDDxERVTZdj7UAoMuHB0utwTMmYT9C3x+HNm2esvrZHYYfEzD8EBFRZdL1WMvd0Q69Gvvip9j8Irgtbl/EM5djsPzpV3T2YS3rekqj7/d3tTvbi4iIyJIUPNYqOrvz4FFufvARBLwSuwc/fjsTb0Zvw5CzB3T2k5yejfFbYhEVL6+MYVdrDD9ERERmUtZ5XM7Kh1i5+yMs3Lce9ur8g6UHJRwGdDy0KWiZvzsBKjUf6pSG4YeIiMhMip7HVViD1GvY+fU0DDx/RNP2edtwvDp0LiDSvbZHACBPz0ZMYlpFDNdiVLut7kRERJZif0Kyzvbwc4ew+PfVcMrNP3dSYe+Ed/pNwe8NO+nVb8E5X6Qbww8REZEZ5OSp8cOpG1ptkrwczDnwGV6Oi9K0nfcJwPjwCCR56nEA9mO+rg7lNk5LxPBDRERUyaLi5Xh3ezwysrUPcp51eKNW8NnWrDfm9B4HpZ1+RzCJkH8Ce8E2edKN4YeIiKgSlVS0EABWdxyGvhePwT07E7N7j8ePzXsb3H/h09pJN4YfIiKiSlLW7q57zu4YH/4usu3scd43yKC+WedHf9ztRUREVEkK7+7yyUzDit0fwf2RQuua07UaGRx8pvZqgKMzezL46IkzP0RERCbSeeK6WFSs/fb9hwCADtfPYNWupfDJegD3R5kY8/xcCKKy5yPEIqBwCR/O9hiH4YeIiMgEuo6m8JM6YGALP+z6V67VLoYa40/8jOlHvoGNoAYANExNQk3FXdyS+pb6OSIAq4e3goezhKe1m4jhh4iIyEi/nrmNCd+dLtYuT8/GhiOJWm1u2ZlYvmcZel35W9N2JKAVpoRNR5qTtNTPkTra4sMhzTnDU04YfoiIiIzw6xk5Jn1fPPjo0kx+CWt3LoF/+h0AgBoirOz8IlZ2ehFqsU2Z9zvZ26J3sMyk8dITDD9EREQGioqXY8J3sWVfKAh4Oe43zDnwGSSq/LO50hzdMGXA2zgS1Frvzys4sqJjPS9jh0yFMPwQEREZoGC7uj46XfsXH/yxVvM6tmZDTBw0C3I3H4M/l0dWlB9udSciIjJAaYeRFnW8bgv83KQHAOCr1gMx7KUlRgUfgEdWlCfO/BARERnAoBkYkQjvPzsRext1xcH67Yz6PB5ZUf4480NERKQHlVpA9JV7uHQnQ+f79nm5mLdvPZ79L1qr/ZG9g0nBB+CRFeWNMz9ERERl0FXLp7Ba6SlYs3MJWsr/w+BzhzDAJwDXPcrell60aGHR1zIWMawQDD9EREQ6FFRn3peQjK+OJZV4Xfcr/+CTPcvgkZ0/IyTJy0Hj1MRSw8+kHvXQub4PWtf1wKlr9zVFC4u+ZhHDisHwQ0REVERZMz0AIFarMPXod5gcvU3Tds1dhgnhEThXo16p/Teo4arZtl50+zq3s1c8hh8iIqJCouLlGL8ltsST1wHAO+s+Pt39ETpfO6Np+6NBB0zvNwUKB5cyP4M7t8yL4YeIiKxW0YNHW9f1wPzdCaUGn7Y34rF611LUyEwDAOSJxPiw22h83u45QFT6Iyru3KoaGH6IiMgqRcXLMXdnPO5k5Gja3B1t8eBRXon3OORmY+2OJfB5+AAAcMfFE5MGzsDf/k3L/Dzu3Ko6uNWdiIisTlS8HOO2xGoFHwClBh8AyLZzwIx+bwEAjtdpjv6jP9Ur+AD5Mz7rRoRw51YVwJkfIiKyKiq1gGk//Gv0/YfqtcWIFxbieN3mpR5KKnOTYHi7OgjwdubOrSqG4YeIiKzK8Ut38TBHVfaFgoBhZ/5Am5vn8U6/t7TW8xwNbFXibQXb2Bl2qi6GHyIisio/n75Z5jUOudlY9Mc6DI0/AAA441cf34QM0Kv/wtvYqWpi+CEiIqtS1qxPYNotrN0RicapSZq2uvflevfPbexVH8MPERFVS0W3qZf1mKngeif7ktfp9LtwFB/+9ilccx4BALLsHDArdDJ2B3crczzcxl59MPwQEVG1o6sCs18p52CVVbHZTpWLiEMb8eqpXZq2/7zqYHx4BK54+5c5Hm5jr14YfoiIqFopqQKzPD0b47bE4rXOAegVLNPMBJVVsdlPkYo1O5cg5PZFTdv24O54t88kPLJ/8ghLBMDVwRYt/d1x+sYDZGQ/2RbPA0irF4YfIiKqNlRqocwKzF8eS8KXx5Lg6WyPBWFN8MFv50u9fvpf32iCj9LGFvN7vYHvWoTqrNa8dGhzhDb1M/iRG1UtDD9ERFRtxCSmlXrYaGFpWTmYtPV0mdfNf+Z1tLtxDgAwPjwC8bL6Oq97/elAzcyOjVjEHV3VGMMPERFVGykZ+gWfUgmC1qyOwsEFo4fOQ6qLR6mHku76V44ZoY05w2MBeLwFERFVG6ZuI299MwG7vp4Kn8eHkha44u1f5mns8vRsxCSmlXoNVQ/VLvwcOXIEYWFhqFmzJkQiEXbs2KH1viAImDNnDvz8/ODo6IhevXrh0qVL5hksERGVq3aBnvCTOsDguRdBwGt/78DW7yPQPPkyVu9aChu1HlWeiyiXmScyu2oXfrKystCiRQusWbNG5/tLly7FypUrsX79epw8eRLOzs7o06cPsrP5LywRUXVnIxZhdv/GpS5gLspVmYV1OyIx++AXsHsceAQAzo9r+RiCBQwtQ7Vb89O3b1/07dtX53uCIGDFihV4//33MWjQIADA119/jRo1amDHjh148cUXdd6nVCqhVCo1rxUKRfkPnIiITBYVL8fCvef1vr5xylWs3RGJwEIVmtd2GIplXV+BqpRDSYtiAUPLUu1mfkqTmJiI5ORk9OrVS9MmlUrRvn17REdHl3hfZGQkpFKp5sffv+yCVkREVLkK6vXou9vr+TP7sP2b6Zrgky5xxqtD5mBpt9EGBx+ABQwtSbWb+SlNcnIyAKBGjRpa7TVq1NC8p0tERASmTZumea1QKBiAiIgqmCG1cvSp71NAkqvEgn3rMezsPk3bGVl9TBg0CzfdZQaPkwUMLY9FhR9jSSQSSCQScw+DiMhqGHo8hSH1fbomxWkFn29a9cOinv+D0ta+2LXtAjxwOTULaVk5WuN4sW0dBHg7sYChhbKo8COT5Sf6O3fuwM/vyX88d+7cQcuWLc00KiIiKqyk4yaS07Mxfkss1o0IKRaADNlltb9Be3zTqh+GxB9ARJ9J2Nmkh87rxCJgy/86wEYsYrVmK2NRa34CAwMhk8lw4MABTZtCocDJkyfRsWNHM46MiIiA0h9fFbTN350AlVr7isTUrBL7FOvYsr6w51gMGPVpicEHANQCcOrafU215kEta6FjPS8GHytQ7cJPZmYm4uLiEBcXByB/kXNcXByuX78OkUiEKVOmYNGiRdi1axfOnj2LkSNHombNmggPDzfruImIqOzHVwKKFxP89YwcKw7ortcmU9zFD9/NwsCEw1rtObZ2uOpVu8zxsG6Pdap2j73++ecf9OjxJMkXLFQeNWoUNm3ahBkzZiArKwuvv/46Hjx4gC5duiAqKgoODqzNQERkbvqGjWOXU5GSkY3E1Ex8euCyzmu6JJ7Gp7s/gtcjBYJTriLBNwiXvesYNB7W7bFO1S78dO/eHYJQ8np/kUiEBQsWYMGCBZU4KiIi0oe+YWP1oSslvicS1Jh8fBumHP0O4scPy9Ic3eCQl1PiPbr4sW6P1ap24YeIiKqvguMpktOzDarSXMDjYTpW7FmGbomxmraDQW0wbcA0PHB0M6gv1u2xXtVuzQ8REVVfNmIR5oYFA4DB53O1unUBezZN0QQflUiMpU+PxGtD5xgcfNa+1Ip1e6wYww8REVWq0KZ+WDciBDKpnuttBAGjTu3Gtu9moVZGKgAg1ckdI4YtxNqOL0AQGfZVNuWZBujXvKahwyYLwsdeRERUrkqr3FzwnjJPjaVDmiNBrsCn+y/hYW7JJ6xLszMxKXob7NV5AICY2sGYNHAmUly9DB6bh5MdJj/TwLg/GFkMhh8iIio3JVVunt2/MS6lZGHjsUQ8eJRrUJ/pjq54M2wGvtn2Pr5oG46Pnx6JPBvDv75EACIHN+M6H4JIKG3rlJVSKBSQSqVIT0+Hm5thz5GJiKxVSZWbjWGnykWujZ1WW937t3HNo/THVSJA5+eXdnQGWQ59v78580NERCYz5ODR0kjycjB3/wb4ZqZh7JDZWut5ygo+QH7wmd2/MTyd7ZGWlQNPFwlkbjyygrQx/BARkckMOXi0JHXuy7F25xI0vZNf4+eNk79gfYehBvfj7SrBoJa1TBoLWTaGHyIiMpmpx0T0vnQCy/Z+Ajdl/hlej2wlSHHxMKovVm2msjD8EBGRyYwNHLaqPEw/8jXGxfyiabvqURPjn3sXF30Cil0vFuUfSKqLCICMVZtJDww/RERksvtZSoPv8c24h9W7PkS7mwmatj0Nu2BW3zeRKXEqdv2kHvUR7OeKid+dBqC9sLlgNQ+rNpM+GH6IiMhoKrWAE1fv4d3t8Qbd1/Hav1i56yP4PHwAAMgV2+CDHq9hU+swQKQ7vHSu742O9bywTiwqtp1ext1cZACGHyIiMoqumj76euHMPk3wue3qjUmDZiK2VuMSr/dyttc8zgpt6ofewbISCykSlYXhh4iIDGZqTZ/3n52AZsmXcdvNB2+FTcd9J2mp1y8c1FQr3NiIRehYz/AKz0QAww8RERlApRZw4so9zPr5rEHBxyE3G9l2TxZFZ0mcMHx4JO45uUEttin13jeeDkS/5nycReWHB5sSEZFeouLl6PLhQbz85Un9j6gQBIyI3YsjG8aiVnqK1lupLh5lBp/+zfwQ0S/Y2CET6cTwQ0REZSp4zGXI+h6nnEf4dPfHWLRvHXyz7mPNzkjY5xl2rlfs9ftQlbS3nchIfOxFRERaip7K3rquh8FHV9S/ex3rdkSiwb0bmrZ/agVDXcJOrpLI07MRk5jG9T1Urhh+iIhIQ9cOroJzsvQ1MOEwlkStglNufu2fDHtHvNNvCqIadjZqTKZWjyYqiuGHiIgAlLyDS9/gY5+Xi9kHP8crp3/VtJ33CcD48AgkeRp/1haPq6DyxvBDRERQqQXM22X8qey10+9g7Y5INE++rGn7oVkvzO49Hko7iVF98rgKqigMP0REVqLoWp7ChQFXHfgPyQrjHy81SknSBJ9sW3vM7j0OPzZ/1uj+eFwFVSSGHyIiK6BrLY/f4yMhTl+/jw1HEk3qf3+D9ljfbjD6XIrGhPAInPcNMqk/HldBFUkkCAL3EBahUCgglUqRnp4ONzc3cw+HiMgkJa3lEQFGP+Zyy86EQuKsdQ6XjVoFp9xsZEic9erDwVaM7Dy15rWnsx2ea1kLvYJlPK6CjKLv9zdnfoiILJhKLZS4Td3Y4NP++lms2rUUy7u8jK0tQ598lthGr+Dj7mSHJYOb8XwuMhuGHyIiCxaTmGbUwaO6iAQ13jj5C9458jVsBDXm79+As7L6OCerb1A/a4aHoHMDbwBg/R4yC4YfIiILVl41ctyyM7Fs73L0vhyjaYup3QTJrt4G9eMqsUEHBh4yM4YfIiILVh41cpomX8a6HZHwT78DAFBDhFWdXsSnnV8s82yuooa0rs1HW2R2DD9ERBbsflYOxCLAqOOxBAEv/RuFufs3QKLKAwCkObph6oC38WdQa6PG06cJd2+R+TH8EBFZkMK1fJLuPsSK/f8ZtbDZMScbH/yxBoPPHdK0nfZriInhM3HbzdeosfmxYCFVEQw/REQWQlctH2NJVDlofz1e83pj6zAs7vEqcm3sjOpPBBYspKpDbO4BEBGRaVRqAZ/u/w/jtsSW286uB45umDRoJu47uGLiwJmY3+sNvYOPg632V4uf1AHrRoSwYCFVGSxyqAOLHBJRdREVL8e8XQkmHU0BAHaqXDjmKqFwcNFqd8p5hIf2jnr1IRYBY7sGYkZoY9bvIbNgkUMiIgtXUuVmQ9VKT8GanUuQZe+AV15YqLWDS9/gAwAfDm6G59vWAcD6PVS18bEXEVE1VFrlZkN0v/IP9mx6Cy3l/6HztTN469j3RvelyM4zcTRElYMzP0RE1ZCplZvFahWmHP0Ob0Zv07Rdc5fhj6c6Gt2np4vE6HuJKpPFzfzMmzcPIpFI66dRo0bmHhYRUbkypXKzV9YDfP3DHK3g80eDDggbtQLnatQzul+Zm+kFFYkqg0XO/DRp0gT79+/XvLa1tcg/JhFZMXeJcVvO29w8h9U7P4QsMw0AkCcSY2m3Ufis3WCtE9oNxRo+VJ1YZCqwtbWFTCYz9zCIiCpE5K8J+OxIomE3CQLGxmzHzD83wVZQAwDuuHhi0sAZ+Nu/qdFjKYhLrOFD1YnFPfYCgEuXLqFmzZoICgrCyy+/jOvXr5d6vVKphEKh0PohIqqKIn9NwIYjiYYvdBaJUC/tpib4HK/THP1Hf2pS8AEAGWv4UDVkcXV+fvvtN2RmZqJhw4aQy+WYP38+bt26hfj4eLi6uuq8Z968eZg/f36xdtb5ISJzKnxUha+rA1r6u6PJ3CjjzukCIMlV4udvZ+BQUBt80uUlvQ8lHdyqJgQByMpRoXVddwT7SZH2MIc1fKjK0bfOj8WFn6IePHiAunXrYvny5Xjttdd0XqNUKqFUKjWvFQoF/P39GX6IyGx0HVXh6mCLDAO2k9dOv4Ob0hpabfZ5ucixNWy90CcvtMBzIbUNuofIHFjk8DF3d3c89dRTuHz5conXSCQSSCTcoklEVUNJxQv1DT4OudlY+Md69P3vGAaO/ARXvZ4EF0ODDwCkZeUYfA9RVWaRa34Ky8zMxJUrV+Dnx+fRRFT1qdQC5u0yvnhhYNotbP9mOp6P3w+XnEdYuyMStirTig+yfg9ZGosLP9OnT8eff/6JpKQkHD9+HM899xxsbGwwfPhwcw+NiKhMqw9eMvqcrr4XjmLX5ilonJoEAMiyc8Daji8gz8a0SX7W7yFLY3GPvW7evInhw4fj3r178PHxQZcuXXDixAn4+PiYe2hERCVSqQWsPngJn+y/ZPC9dqpczDq8Ca/9s1PTdsnLH+PC38UVb3+TxuXpbMf6PWRxLC78bN261dxDICIyiCkns/spUrF654doffuCpm1HcDe822eSQYeSlmTRoKbczUUWx+LCDxFRdfLrGTkmfBdr1L2dkuKwatdSeD3Kr02mtLHFwmdex5aWfU2q1lzgjacD0a95TZP7IapqGH6IiMzk1zO3MfH700bfLwLg8SgDAHDTzRfjwyNw1q9BuYxtbNdARPQLLpe+iKoai1vwTERUHUTFyzHhu9MwpdLasYCWWNHlJRyo1xb9R39absEHAL74KxFR8fJy64+oKuHMDxFRJcvfzn7O4Pvq372Oy17+Wo+0VnUaBgAQRIb/v6yDrQjZeSWnr/m7E9A7WMY1P2RxOPNDRFTJYhLTkKxQln1hAUHAa3/vwG8bJ2vt6ALyQ48xwWdAc79Sg48AQJ6ejZjENIP7JqrqGH6IiCqJSi0g+so9rDtccsX5olyUD7F2RyRmH/wCdmoVZh3eiCZ3rhg9BjcHW6x9KQS9g2uUfTGAlAzjag4RVWV87EVEVAl0ndVVlkYpiVi7IxJB929r2j5v9xwu+AQYNQZPZzuciOgFe1sxoq/c0+seX1cWOCTLw/BDRFTOCp/G7u0iQUxiGj49YFjxwqFn92PRH2vhkJd/rla6xBnTBkzDgfrtDR5PwYqdxc81g71t/oR/u0BP+EkdkJyerfMoDREAmdSBBQ7JIjH8EBGVI2NmeAqT5Coxf/8GvHjmD03b2Rr1MD48AjfdZUb1KZM6YG5YMEKbPjnj0EYswtywYIzfEgsRoBWACsLS3LBgLnYmi8TwQ0RUTko6jV1fNRUp+OLnhQhOSdS0bWnZFwufGQulrb1Rfc7u3xijOwfqDDGhTf2wbkRIsbCmKywRWRKGHyKicqBSC5i/2/jT2AFAIXGBQ27+LrCHdhK822cSdjTpYVRfBY+tSgo+BUKb+qF3sEzzmM7XNf9RF2d8yJIx/BARlYOYxDSjH3UVyJQ4YUJ4BD7+dQWm9p+GSz51TepP38dWNmIROtbzMumziKoThh8ionJgzJbwGhl3IRKAZDdvTdsF30AMGLXCpLO5/PjYiqhUDD9ERCZSqQXczTCgaCGAzklx+HT3R7jm7ocXX4pEro3dkzdNCD5TezXApJ4N+NiKqBQMP0REeiq8hb1gbcy+hGSDdneJBDUmHd+GqUe/gxgCvB+mY2L0D1jR5WWTxiYCsOalEPRrztkeorIw/BAR6UHXFnZ3R1s8eJSndx/ujxRYsXsZuiee0rQdCmqNTa3DTB7fmpdaMfgQ6Ynhh4ioDCVtYTck+LS8fRFrdixBrYxUAIBKJMbyLi9jbcfnjTqbq4DMTYJ5A5twfQ+RARh+iIig+5GWjVhk+hZ2QcCo2D147+CXsFfnh6W7TlK8GfYOjge0NGnMpdXwIaKSMfwQkdXT9UirYMeU1NHe6C3sYrUKn+7+GGEX/tK0xdQOxuSBM3DH1buUO/Xj7Sph8CEyAsMPEVm1kh5pJadnY/yWWIzpHGB032qxDVJcnpyNtaHdYHz09Ejk2ZTPX708dJTIOAw/RGS1SnukVdD23cnrJn3Gku6jEZR2E9+3CMUfT3U0qa8CPHSUyDQMP0RktfSpypydp9a7P0leDpomX8ap2sGatlwbO4wZOs+k2j2F8dBRItMZv8WAiKia25+QXG59+T9Ixs9b3sG3295H45Sr2m+WU/AB8md81o0I4e4uIhNw5oeIrFJUvBxfHksql756XTqJ5XuXw02ZBQBYvmc5+o1ZadIW9gIDmvlheLs6uJul5KGjROWE4YeIrE7BWh9T2ahVeOfI1xh38mdN21WPmpg64G2Tg49YBIztGoiIfsFlX0xEBmH4ISKrolIL2HQs0eQT2H0y07B654dof/Ocpm1vw86Y2fctZEqcjO63SU1XDG5VG690DIC9LVcmEFUEhh8ishq66vkYo+O1M1i5eyl8sh4AAHLFNljc41VsbD3Q6PU9nOkhqjwMP0RkFUqq52OoEbF7MX//BtgI+bvAbrt6Y9KgmYit1djoPsNb1sTSoS0400NUSRh+iMjimXxERSHnfQOhFolgIwBHAlphSth0pDlJTepzWNs6DD5ElYjhh4gsnj71fPR1qnYwIru/CldlFlZ1Gga12MbovliskMg8GH6IyGIVHFa6/s/LxnUgCOh9+ST212+ntXvrq7aDTB4bixUSmQ/DDxFZJFMXNzvlPMLi31cjPOFPLOvyMlZ1Hm70WEKb1kBM4n2kZeVo2mSPD05lsUKiysfwQ0QWx9TFzfXvXse6HZFocO8GAGDKse+xp/HTSPSsZXBf7o52WPNSawD5j99SMrJZrJDIzBh+iMhi5OSpsfl4Epbvu2h08BmYcBhLolbBKVcJAMiwd8SMvm8ZFXwAYEznAE3I6VjPy8hREVF5YvghIosQ+WsCPv8rEWojU499Xi5mH/wcr5z+VdN23icAE8IjjA4+Hk52mNSzgXEDIqIKw/BDRNVe5K8J2HAk0ej7a6ffwZodS9Ai+ZKm7cemvTD72XHItnMwqk8RgMjBzfhoi6gKstjCEmvWrEFAQAAcHBzQvn17xMTEmHtIRFQBcvLU+Owv44NPk+TL2LPpLU3wyba1xzt938Q7/acYHXyc7Gx48jpRFWaR4Wfbtm2YNm0a5s6di9jYWLRo0QJ9+vRBSkqKuYdGROVs07FECCZUL7ziVRtyV28AQJK7HwaP+Bg/Nn/WpDFN7f0Ugw9RFSYSBP3/2rhx4wb8/f0rcjzlon379mjbti1Wr14NAFCr1fD398fkyZMxa9asMu9XKBSQSqVIT0+Hm5tbRQ+XiAxQULsnJSMbSXcfYu3hS1DmmVa7OTDtFiYf34q5vcchQ+JsUl8iABcX9WXFZiIz0Pf726A1P40aNcLbb7+NWbNmwcnJ+FOLK1JOTg5OnTqFiIgITZtYLEavXr0QHR2t8x6lUgmlUql5rVAoKnycRGS48jiYtN2NeCS7eOG6x5OZmUTPWpg24O3yGCIGNJcx+BBVcQb9F7pv3z78/vvvaNCgATZt2lRBQzLN3bt3oVKpUKNGDa32GjVqIDk5Wec9kZGRkEqlmp/qMLtFZG0KavcYHXwEAW+c/Anfff8u1u2IhCRXWfY9BnKW2GDFiyHl3i8RlS+Dwk+nTp1w8uRJREZGYvbs2WjdujX++uuvihpbpYmIiEB6errm58aNG+YeEpHVU6kFRF+5h51xt3Ds0l3M22X8waRu2Zn4/JdFiDi8CbaCGk1SrmJE3G/lOl4RgGXPt+DuLqJqwKit7iNHjsTQoUOxZMkS9O3bF6Ghofjoo48QGBhY3uMzmLe3N2xsbHDnzh2t9jt37kAmk+m8RyKRQCKRVMbwiEgP5fF4q0CT5MtYtyMSddKf/J3waacXsbF1mFH9OdnbwN5WjAcPczVtfjyqgqhaManOz7PPPguFQoFVq1Zh7969mDx5MubMmQMXF5fyGp/B7O3t0bp1axw4cADh4eEA8hc8HzhwAJMmTTLbuIhIP6YeTaEhCBj+7++Yt38DJKr8oHLfwRVTB7yNw/XaGNXle/0a4dUuQQB4VAVRdWZQ+Fm/fj3+/vtv/P333zh//jzEYjGaNm2KcePGoUWLFti6dSuCg4Pxyy+/oE0b4/5yKQ/Tpk3DqFGj0KZNG7Rr1w4rVqxAVlYWxowZY7YxEVHZVGrBpMdbBRxzsrHojzUYcu6Qpi3O7ylMCJ+F226+RvXpJ3XAq12CeFQFkQUwKPx88MEHaN++PUaOHIkOHTqgdevWcHR01Lz/+uuvY/HixRg9ejTi4+PLfbD6GjZsGFJTUzFnzhwkJyejZcuWiIqKKrYImoiqltUHLyFZYdqjLkleDn7ZMh2NU5M0bRtbh2Fxj1eRa2NnVJ8iAHPDgjm7Q2QhDKrzo487d+6gZs2aUKlU5dltpWKdH6LKFxUvx7gtseXS17Qj3+DN6G3ItHfErNDJ2NP4aaP7krlJMG9gE67nIaoGKqTOjz58fX1x8ODB8u6WiCyUSi3gxNV7mPXz2XLrc0WXl+CRnYFNrcNwxcv40hVTez2FST3rc8aHyMKU+8yPJeDMD1HliIqXY96uc0hWGF9zp6YiBS1u/4ffGnUpt3F5Odvjg+eacraHqJox28wPEZE+yuMxV7erp7Bi98dwznkEuZsP4mo2NHlcns52iI54hlWaiSwY/+smokqnUguY9Yvxj7nEahWmHfkGG3+cB4/sDNir8zDr8MZyGdvi55ox+BBZOM78EFGlKTiU9K9LKVpFAg3hlfUAn+7+CF2u/atp21e/Pd7uP9WksYlFwOrhIXzURWQFGH6IqFKUR9Xm1jcTsGbnEsgy0wAAeSIxlnYbhc/aDQZEpi1KXj28Ffo1Z/AhsgYMP0RU4Uyu2iwI+N/f2zHr8dlcAJDi7IFJg2Yixr+pSWPj0RRE1ofhh4iMVvAYq6RjHgpvYzdlW+ncA59hzKndmtfRdZrhzbAZSHXxMLrPST3qoXN9Hx5NQWSFGH6ISC9Fg879LCUW7j2v9Rir8CxKeR5O+lPTZ/BSXBQkqlys7vgCPunyMlRiG6P6EgGQSR0wtXdDhh4iK8U6Pzqwzg+RNn2DTEGUeP3pQHx2JNH0w0kLGXL2ANKc3HCoXluj+ygY37oRXNhMZIn0/f5m+NGB4YfoCWPW64hFgNrIv1kccrMx+tQefN7uOaNnd0rC9T1Elo1FDonIZCq1gPm7DT9l3djgE5B2C+t2RKJxahLcszOwpPsY4zoqYmTHuujb1I/re4gIAIscElEpYhLTymXNjj5CLx7D7s1TNKexjzj9K3web2k3Vd+mfuhYz4vBh4gAcOaHiEqRklHxwcdOlYtZhzfhtX92atouefljfHgEUl08Te5f5iZBu0DT+yEiy8HwQ0Ql8nV1qND+ZYq7WLNzCVrfvqBp29m4GyJCJ+GhvWO5fMa8gU0440NEWhh+iKhE7QI94Sd1QHJ6drnu3AKAromxWLH7Y3g9UgAAlDa2WPjM69jSsq/B1ZpFQLHxuTvZYcngZlzcTETFMPwQUYlsxCLMDQvG+C2xOgOGsbpf+Rtf/bQA4sc93nTzxYTwWTjj95TBfYkArHmpFaRO9oi+cg+AgI5B3ujANT5EVAKGHyIqVWhTP7z+dCA+/ysRhQtjiEVAj4Y+iEm8jwxlnkF9Hq/bEmdl9dEi+RIO1GuLaf2nId3R1eCxFd263rm+t8F9EJH1YfgholJFxct1FixUC8CBC6lG9Zlja4eJ4bPQ78JRfN7uOQgiwzaejuxQF32bces6ERmH4YeISmRsnR8tgoBRsXtwvE5zXPKpq2m+Ka2Bz9oPMbi7tS+F8PR1IjIJww8RlejElXsm1flxUT7Eh799iv4Xj+GyZ20MHPWJSbu4PJzs0KepzOj7iYgAFjkkohJExcsx8btYo+9vmJqEXZunoP/FYwCA+mk38czlGJPGdP9hLmISy6fwIRFZL878EFExxpznVdiQswew6I+1cMxTAgAUEmdM6z8N+xu0N3lslVF4kYgsG8MPEQHIX98Tk5iGZEU23t9x1qjgI8nLwbx96zH8zB+atrM16mFCeARuuJfP46qKLrxIRJaP4YfIihUEnv0JydgedwtpWblG91XnvhzrdkSiScpVTdu3LUOx4JnXobS1N3msIgAyqQOPqiAikzH8EFmpqHg55u9OKJeDS12VWdjxzdvwfFyt+aGdBO89OxHbm/Y0uW8gP/gAwNywYG5tJyKTccEzkRUqWNNTXie2Z0ic8Vm7wQCAK561Ef7KsnILPkD+jM+6ESE8qoKIygVnfoishGZNT/ojzN55rtzP6trQfjDyxGJ83yIUWRKncunT3ckOa4aH8KgKIipXDD9EVqA8H3EBQKekONS/dwNftw7TtAkiMb54PPujD6mDLZ5+yhu7zyQXe68g5iwZ3AydG/DICiIqXww/RBbO1G3rhYkENSZG/4Bpf30LQSTCfz51caJOc6P6WjuiNTrX90b/5sWDmazImV1EROWJ4YfIgpXL8RSPuT9S4JM9y9Dj6qn8BkHA0LMHjAo/flIHdAjyApB/cGrvYBliEtOQkpENX1cHntlFRBWK4YfIgsUkppXLo64Wty9izc4lqK3IP8hUDRGWd30Zazq+YFR/RXdt2YhF6FjPy+RxEhHpg+GHyIKZXA1ZEDAydg/eP/gl7NV5AIC7TlK8FfYOjgW0NLg7Z3sbLHuhBR9nEZFZMfwQWbCkuw+NvtdZ+RCRv6/GwPNHNG1/1wrGpEEzcMfV8EXIDrZinJ7zLOxtWWGDiMyL4YfIQuXkqbE5Osno+5fvXY4+l05oXn/W9jks7TYKeTbG/bWx4sWWDD5EVCUw/BBZoKh4Od7dftak4yo+enokuiadRp7IBu/0m4LfG3Yyqh+ZmwTzBjbhoy4iqjIsLvwEBATg2rVrWm2RkZGYNWuWmUZEVLnKa2v7Ze86mDhoFq561sI1j5oG3duithSvdgnkzi0iqpIsLvwAwIIFCzB27FjNa1dXVzOOhqjiFa7evHDveYODj/+DZEw6vg2zn52AHFs7Tfuhem2NGs/AFjUxqGUto+4lIqpoFhl+XF1dIZPJzD0MokphavXmXpdOYtne5ZAqs6C0tcecZ8ebNB6xCHilY4BJfRARVSSLDD9LlizBwoULUadOHbz00kuYOnUqbG1L/qMqlUoolUrNa4VCURnDJDJKwSxPSkY2ku5m4ZP9l4zqx0atwvQj32D8yZ80bV2STsNVmYUMibPR4xvbNZALm4moSrO48PPmm28iJCQEnp6eOH78OCIiIiCXy7F8+fIS74mMjMT8+fMrcZREximvM7p8MtOwetdStL8Rr2n79alOmNnvLRODTwAi+gWbNDYiooomEgShvA93LnezZs3Chx9+WOo158+fR6NGjYq1f/XVV3jjjTeQmZkJiUSi815dMz/+/v5IT0+Hm5ubaYMnKifltZC5w/UzWLVrKXyyHgAAcsU2iOz+Kr5qMxAQGb8wedXwVghrYdjCaCKi8qRQKCCVSsv8/q4W4Sc1NRX37t0r9ZqgoCDY29sXaz937hyaNm2KCxcuoGHDhnp9nr6/PKKKVnQhc1pWjtF9iQQ1xp38GdOPfAMbQQ0AuO3qjUkDZyK2dmOj+/XjIaREVEXo+/1dLR57+fj4wMfHx6h74+LiIBaL4evrW86jIqpY5fWIq8DQswcw88/NmtdHAlphSth0pDlJDeqne0NvvPF0fR5CSkTVVrUIP/qKjo7GyZMn0aNHD7i6uiI6OhpTp07FiBEj4OHhYe7hEemtvB5xFfZL054YGn8AbW+cw6edh2NVp2FQi20M7ueNp+vzEFIiqtYsKvxIJBJs3boV8+bNg1KpRGBgIKZOnYpp06aZe2hEelOpBczfnVCuwQcAVGIbTBo4Aw1Tr+FoYCuj+vCT5s/0EBFVZxYVfkJCQnDixImyLySqwmIS00x+1OWU8wjz923A5tYDEC+rr2lPdfFEqotx4UUEYG5YMB9xEVG1x2IcRFVMSoZpwafe3RvY+fU0PB+/H+t2RMItO9PkMflJHbBuRAgXNRORRbComR8iS+Dr6mD0vQMT/kRk1Co45+YHKPdHCjRKTUKMf1Oj+nOR2GDDiDboUM+LMz5EZDEYfoiqEJVagFotwN3RDg8e6X8iu31eLt4/+AVGnt6raTvvE4AJ4RFI9DT+jK2lQ1qgcwNvo+8nIqqKGH6Iqghjt7bXTr+DNTuWoEXyk2Mufmr6DN5/djyy7YyfRXrj6UD0a87HXERkeRh+iMyooIjh/oRkfHksyeD7e1z5G5/sWQb3x+t6lDZ2mN17PH5o3tvoas2eznZYNKgp+jVntWYiskwMP0RmEhUvx7xdCUhWGLfA2TfjHtZvXwyJKv/x2DV3GSaER+BcjXoG9+XmYIuhrWujd7CMRQuJyOIx/BCZQVS8HOO2xJrUR4qrFxb1fA0L963H7w064J1+U6BwcDG4H09nO5yI6MWT2InIajD8EFUylVrArF/Olktf37Tqj9tuPjhQr51Rj7lEABY/14zBh4isCv/GI6pEKrWAGT/9iwcP9d/JBQAQBLx+8mdM/etb7XaRCAfqtzcq+LB2DxFZK878EFWSqHg5Zv1y1uDg45adiY9/XYFnL+VXLz9d8ykcrtfW4M/v0dAHXep7w9PZHjKpI9f2EJHVYvghqgTGHlTaJPky1u5cgroPkjVtjVKvGRx+RAAuJGfgi1FtGXiIyOox/BBVMKMOKhUEvPjv75i/f4NmN9d9B1dMHfA2DtdrY/AYBADy9GzEJKbxRHYisnoMP0QVpKCGz7HLqQYVLnTMycaifWsxJP6gpi3O7ylMHDQLt6S+Jo3J1HPDiIgsAcMPUQUwtlpz0L2bWLsjEo3uXtO0bQoZgMU9XkOOrZ3J4zLl3DAiIkvB8ENUzoxd3wMAi39frQk+WXYOmBU6GbuDu5k8JhEAmdQB7QI9Te6LiKi641Z3onJk1PqeQmb0fQsKiTP+86qDgSM/KbfgAwBzw4K52JmICJz5ISpXMYlphj3qEgStGj3XPfzwygsL8J93XTyyL59HVDKpA+aGBbOeDxHRYww/ROXIkAXF3a6ewriTP+HVIXO1gs6/NRvq3YcI0Jpl8nS2w4KwJvBydUBKRjZ8XR1Yz4eIqAiGHyITFOzoKgga3i6SMu8Rq1V469hWTD6+FWIIWPTHGrzdf5pBVZqd7W3w2Stt0DbQE6eu3WfQISIyAMMPkZF07eiSuUng7mSH9Ie5Otf9eGU9wIrdH6PrtThNmzQ7ExJVLpS29np/9ott/SEWi2AjFrFuDxGRgUSCIBi7NtNiKRQKSKVSpKenw83NzdzDoSqopB1dRR9DFdb6ZgJW7/wQfpn3AAAqkRgfPT0SG9oPhiDSb++BWASoC32AH9fzEBFp6Pv9zd1eRAYqbUeXzuAjCHgtZju2fh+hCT6pzu54+cVFWN9hqN7BB9AOPgCQnJ6N8VtiERUv17sPIiJrx8deRAYyZEeXqzILH/26AqH/RWvaTvg3xeSBM5DqYnrNHQH5s03zdyegd7CM632IiPTA8ENUgqKLmQsWExuyo2vI2QNawWdth6FY1vUVqMQ25TZOnttFRGQYhh8iHXQtZi5YX2PIERGbWw9Aj6v/oOXti5g64G0crN+uIoYLgOd2ERHpi+GHqIiSFjMXrK9Z9WIruEhskanMK3avSFBrreERRGJMGfA2nHMe4aa7zKjxOEtskKVUlXkdz+0iItIPFzwTFVLWYmYBwOStp3UGn4C0W9j59TS0uXlOq/2+k9To4AMAHw1pDj+pA0pazSNC/qwUz+0iItIPww9RIfosZtYVjPpcPI5dm6eiefJlrNn5Ibyz7ps8Fj+pA9aPCEG/5jUxNywYAIoFIJ7bRURkOD72IirE0HUztqo8zDq8Ef/7Z6emLcPeCS7Kh7jr7GFQXx5OdhjdKQAB3s7FqjWHNvXDuhEhxYsqss4PEZHBGH6ICjFk3YxMcRerd32INrfOa9p2NX4as0In46G9o8GfPSesCZ5rVavE90Ob+qF3sEznDjQiItIfww9RIa3resDDyQ73H+aWel2XxNP4dPdH8HqkAADkiG2x4Jmx2NKqn0FndBUmcys7ePE4CyIi0zH8ED1WsL29tOAjEtSYfHwbphz9DuLHq39uuvlg4qBZBp3GXhQXLBMRVR6GH7J6KrWA1Qcv4ZP9l8q8tu59Ocaf+EkTfA4GtcG0AdPwwNH4M+BE4IJlIqLKxN1eZNWi4uXovOSgXsEHAJI8a+H9ZydAJRJj6dMj8drQOSYFHz+pA9aNCOGCZSKiSsSZH7JaJRUz1CIIsBHUWsdR/NzsGcT5PYUr3v5Gf7a7ox3WvByCDkFenPEhIqpknPkhq1RaMcMCLsqHWL3zQ7x38Mti7xkbfESPf5YMaYbO9b0ZfIiIzIAzP2SVyipm2DA1CWt3RKJe2i0AwD+1g/Froy4mfy7r8hARmV+1Cj8ffPAB9u7di7i4ONjb2+PBgwfFrrl+/TrGjx+PQ4cOwcXFBaNGjUJkZCRsbavVH5UqQOFT2i/dySjxusHxB/DB72vhmKcEACgkzlDa2hn9uWHNZegVLGNdHiKiKqJaJYKcnBw8//zz6NixI778svijCJVKhf79+0Mmk+H48eOQy+UYOXIk7OzssHjxYjOMmKoKXae0FyXJy8Hc/Rvw0r+/a9ria9TD+PAI3DDibC6xCBjbNRAR/YKNGjMREVUMkSAIpa73rIo2bdqEKVOmFJv5+e233zBgwADcvn0bNWrUAACsX78eM2fORGpqKuzt7XX2p1QqoVQqNa8VCgX8/f2Rnp4ONzfjd/JQ1fDrGTkmfBdb6jV17suxducSNL1zRdP2XYs+mN/rDShtdf97U5KX2/kjyMcFr3QMgL0tl9UREVUWhUIBqVRa5ve3Rf3NHB0djWbNmmmCDwD06dMHCoUC586dK/G+yMhISKVSzY+/v/G7eMh8VGoB0VfuYWfcLURfuQeVWsCvZ25j0velB5/el05gz+YpmuDzyFaCt/tNxbuhkw0OPm88HYgPBjfHa12DGHyIiKqoavXYqyzJyclawQeA5nVycnKJ90VERGDatGma1wUzP1R96Hqs5e5khwdlHFMhEtR4/eQvcFNmAQCueNbChPAIXPQJMOjzRQBef5qPuIiIqgOz/6/prFmzIBKJSv25cOFChY5BIpHAzc1N64eqj4J6PUXX85QVfABAEIkxeeAM3HN0w56GXTBo5CcGBx8AWDm8FYMPEVE1YfaZn7fffhujR48u9ZqgoCC9+pLJZIiJidFqu3PnjuY9sjwqtYB5u0qv11OUJC9H63FWsps3Bo5agVtuPkYdSioCsPjX8+jXzI87uYiIqgGzhx8fHx/4+PiUS18dO3bEBx98gJSUFPj6+gIA9u3bBzc3NwQH8//KLdHqg5eQrCh5B1dhIkGNCdE/YvC5QwgfuQwZEmfNe7ekvkaPQQAgT89GTGIaT1wnIqoGzP7YyxDXr19HXFwcrl+/DpVKhbi4OMTFxSEzMxMA8OyzzyI4OBivvPIK/v33X/z+++94//33MXHiREgkEjOPnspbVLxc7zO53B8p8NVP8/HOX9+gXtpNLP31U6CcNzqmZOgXwoiIyLzMPvNjiDlz5mDz5s2a161atQIAHDp0CN27d4eNjQ327NmD8ePHo2PHjnB2dsaoUaOwYMECcw2ZKkjB8RT6aHH7ItbsXILailQAgBoiJPgGQgQBAsrvMZWvq0O59UVERBWnWtb5qWj61gkg84m+cg/DPz9R+kWCgFdO78XsA1/AXp0HALjn6IY3B87AsYCW5TYWEfKPrTg6syfX/BARmZG+39/VauaHqEBZj5icch5hSdQqDDx/RNP2T63GmDRwJpLdvMttHAVRZ25YMIMPEVE1wfBD1VJpj5jq372O9dsXo37aTU3b523D8WG30cizKd9/5XlQKRFR9cPwQ9VSu0BPyNwcdO706pp0WhN8FPZOeKffFPzesJPefbs72uLBo7xi7SLk7+ya2qsBArydeVApEVE1xfBDVVrhk9gLh419CcnIzlPpvGdj64Foe+McAh7IMT48Atc8aur1WQURZsmQ5gBQrGI0Z3mIiCwDFzzrwAXPVYOuIyv8pA4Y2MIPG44katqclQ+RJXHSutdZ+RB5Yhso7XSXOBgaUgtHL9/TmjnyKxJuSgpeRERUNen7/c3wowPDj/kVHFlR9F/OgkdPBZ65fBIf712BNwe+g78CQ8rst/DOLAAMN0REFoS7vajaKqjhoyuVF7TZqFWYfuQbjD/5EwDg090fo//oTyF3K7lauK6dWazITERkfRh+qMqJSUwrdkhpYT6Z97Fq14focCNe03bSvykyizz6KoprdoiICGD4oSqg6Nqa5PRHJV7b/vpZrNq1FL5Z9wEAuWIbLOk+Bl+2GaTzUFJPZzvMHtAEMjc+1iIionwMP2RWuhY1ezrbF7tOJKjxxslf8M6Rr2EjqAEAchcvTBw0C7G1G5fY/6JBTdGvuX67vYiIyDow/JDZlLSoOS0rR+u1W3Ymlu1djt6XYzRtRwJaYUrYdKQ5SUv9DA9nHmhLRETaGH7ILEpb1FyU18N0dLh+FkD+oaQrO7+IlZ1ehFpsU+a9PGmdiIiKEpt7AGSdylrUXFiiZy3M6PsW7jm6YfTz87Ciy8t6BR+AJ60TEVFxnPkhs9iXkFzie4452VCLRFoFCn9r1AVHA1shQ+KsV/8F9XzaBXqaOlQiIrIwnPmhShcVL8dXx5J0vlfv3g3s/Hoa5u/fUOw9Q4IPwJPWiYhIN878UIUrvJXd20WCebvO6bwuLOFPLIlaBefcbDx17zr+qR2Mn5r1MvjzWM+HiIhKw/BDFSoqXo55uxJ0nr5ewD4vF+8d+gKjYvdq2i5410VszUYGf97UXg0wqWcDzvgQEVGJ+NiLKkxUvBzjtsSWGnxqpafgh+9maAWfn5v2RPjIZbjqVdugzxMB2Pr3DWOHS0REVoIzP1QhVGoBs345W+o13a/8jRV7lsE9OxMAoLSxw9xeb2Briz46qzWXRQAgT89GTGIaz+wiIqISMfxQhTh++S4ePMzV+Z5YrcLUo99hcvQ2Tds1dxkmhEfgXI16Jn82a/sQEVFpGH6o3EXFy/H2j/+Wek2r2xc0//xHgw6Y3m8KMhxcIALgLLFFpjLP6M9nbR8iIioN1/xQuSo4siJLqSrxGrXYBm+FvYObbr74oPureP2596BwcEENNwnWjQjB0iHNjfpsEQA/1vYhIqIycOaHyk2JR1YIAmpk3sMdV29N0z1ndzzzv3VahQwLKvT0a+6HN24GYsORRL0/m7V9iIhIX5z5oXKj68gKt+xMbNj+AX755h14PEzXek87+AB3FNkYvyUWUfFyRPQLxtqXWsHT2U6vz5ZJHbBuRAhr+xARUZk480MGK1y00Nc1/zGTjVhUbKFxkztXsHZHJOo+yD/KYvne5RgzdF6JO7kE5M/gzN+dgN7BMvRrXhN9mvo9KZDoLAFEwN1MpdY/Fx4DERFRWRh+yCBR8XLM352gNcPj97iismahsSBg2Jk/sGDfekhU+Tu+Hji44OuQAWVuYS+6Xd1GLOK2dSIiKlcMP6S3gsXMRdf0yNOzMW5LLFYNb4UAR2DSTyswNP6A5v04vwaYOCgCt6S+en8Wt6sTEVFFYfghvZS4mLmQFWv3YNsfH6NG0iVN2+aQ/vigx/+QY6vf2p0C3K5OREQVheGH9KJrMXNhoReP4aNfV8A15xEA4KG9A2b1mYxdwd0M+hwR8hcvc7s6ERFVFIYfKlXB4ubf4uWlXifNztQEn6u+dVF7/14Md5KhzuVUrD50Ra/P4nZ1IiKqDAw/VCJdi5tLsq35s2h7MwE2ahXe7TMJXznJ0LGeF9oFeuLn2FtITs8u9ZEZkD/jMzcsmNvViYioQjH8kE4lLW4uEJB2C0metZ40iESYFToZeWIbQCTCvoRkzW6tuWHBGL8lFiJAq7+C1692DkDvYBm3qxMRUaVgkUMqprTFzWK1ClP/+hYHvhiPZ/+L1novz8ZWs5V9Z9xtqNT5PYQ29cO6ESGQSbUXMcukDlg/IgRzwppoghIREVFF48wPFVPS4mbPh+n4dNdH6HotDgDw8a8r8KysAZLdvItdey8rR1OrB8gPQL2DZTqLIxIREVUmhh/SolILOHoptVh7yM3zWLNzCfwy7+VfJxJjXYehuONa8q6sorV6WLCQiIiqAoYf0oiKl2PWL2fx4GHuk0ZBwGv/7MSswxthp84/qT3V2R2TB87AiTqln77OWj1ERFQVVas1Px988AE6deoEJycnuLu767xGJBIV+9m6dWvlDrQaioqXY9yWWK3g46rMwtodkZh98AtN8Dnh3xT9Rq9ETN3mKOmBlQj5R16wVg8REVVF1Sr85OTk4Pnnn8f48eNLvW7jxo2Qy+Wan/Dw8MoZYDWVk6fGu9vjtdrq372OXZunoN9/xzVtazsMxcsvfoC7Lp4Y2zUQAIoFINbqISKiqq5aPfaaP38+AGDTpk2lXufu7g6ZTFYJI6r+ouLleHf7WaRl5Wq1P7JzgPujTABAusQZ0wZMw4H67eHpbIfFzzVDaFM/tKrjUawOEGv1EBFRVVetwo++Jk6ciP/9738ICgrCuHHjMGbMGIhKOU1cqVRCqVRqXisUisoYptmVVsvnltQXUwe8jalHv8XEQTNx0z0/TM4e0EQTbLiDi4iIqiOLCz8LFixAz5494eTkhD/++AMTJkxAZmYm3nzzzRLviYyM1MwqWYuitXzq3r+Nu07uyJI4aa45XK8N/gwKgSB68nTU10Wi1Q93cBERUXVj9jU/s2bN0rlIufDPhQsX9O5v9uzZ6Ny5M1q1aoWZM2dixowZ+Oijj0q9JyIiAunp6ZqfGzdumPrHqvIK1/Lpc/E4dm+agiVRqwBBex6ocPABgLd//BdRZZzzRUREVJWZfebn7bffxujRo0u9JigoyOj+27dvj4ULF0KpVEIikei8RiKRlPieJVKpBRy7fBe2qjzM/HMTxv69AwAQduEv/BnUGj8161XivXcU2Ri/JRbrRoRwXQ8REVVLZg8/Pj4+8PHxqbD+4+Li4OHhYVXhpjQFh5Wqb97E9zuXou2tBM17uxo/jV8bdi71fgH5O7rm705A72AZ1/cQEVG1Y/bwY4jr168jLS0N169fh0qlQlxcHACgfv36cHFxwe7du3Hnzh106NABDg4O2LdvHxYvXozp06ebd+BVRMEC505Jcfh090fwfpgOAMgR22LhM//DN636a87mKo0AQJ6erXV8BRERUXVRrcLPnDlzsHnzZs3rVq1aAQAOHTqE7t27w87ODmvWrMHUqVMhCALq16+P5cuXY+zYseYacpWhUgtYsDMek459j6lHv4P48VLnm24+mDhoFv6t2dDgPoseX0FERFQdiARB0LXT2aopFApIpVKkp6fDzc3N3MMpFyfPXkf2oOfQLTFW03YoqDWmDngbDxy1/4wiQOf296K+H9uBMz9ERFRl6Pv9Xa1mfsh4ybkiONnYAcg/lHR5l5extuPzxXZzAWUHHxHyixny+AoiIqqOzL7VnSqHr9QJb/efitiaDTFi2EKs6TRMZ/ApisdXEBGRpeHMj6XKyACuXgVatAAAtAv0hHMNbwwZ8TEEPRY1F/BwtkdaVo7mNY+vICKi6o7hxxLFxwNDhwL37wOnTwM1a8JGLMLcsGCM3xKr95oeAJjdvzFkUkceX0FERBaDj70szTffAO3aARcvAikpwBtvaN4KbeqHdSNCIJM66N2dTOqIjvW8MKhlLXSs58XgQ0RE1R7Dj6XIzs4POiNHAo8e5be1bAmsWKF1WWhTPxyd2RPf/q893B3tSuxOBMCPi5qJiMgCMfxYgqtXgU6dgM8+e9I2dixw/DhQr16xy23EInSu740lQ5pBBC5qJiIi68LwU93t3AkhJCR/bQ8AlYMj1F9tzA9Cjo6l3lrSYzCZ1IFndxERkcXigufqbN48YP58zUzNVY+aGP/cu1Ak+2NuvFyv8BLa1A+9g2WISUzjomYiIrIKDD/V2L+OPmjx+J/3NuyMmX3fQqbECaJ0w05etxGLWKmZiIisBsNPNaVSCxinboSxrQfihnsNbGw9UHMoKU9eJyIiKhnX/FQXajUQFaV5GZOYBnl6Nhb0eh0b2wwqdhp74ZPXiYiI6AmGn+ogLQ0YOBDo2xf49lsA+p+ozpPXiYiItDH8VHX//AOEhAB79+a/Hj8eSEuDr6t+hQr1vY6IiMhaMPxUVYIArFsHdO4MXLuW3+btDfz8M+DpiXaBnvCTOhSr0VOARQqJiIh044LnSqJSCzhx9R6ir9wDIKBjkDc6lHRcRGZmfrXm77570tapE7BtG1C7NgCUelYXixQSERGVTCQIgr5nXFoNhUIBqVSK9PR0uLm5mdxfVLwcs345iwcPc7Xa3Z3ssGRwM+3t6AkJ+YeSnj//pG3aNGDJEsCu+HEUUfFyzN+dAHn6k7U9fjx5nYiIrJC+398MPzqUZ/iJipdj3JbYUq9ZX1CPZ98+4LnngKys/DdcXYGNG4EhQ0q9X6UWWKSQiIisnr7f33zsVYFUagHzdiWUed28Xefy6/E0a5YfeLKygObNgZ9+Aho0KPN+FikkIiLSHxc8V6CYxDQkK8reap6sUObX45HJgK1bgddeA06c0Cv4EBERkWEYfipQWTV2ul09BemjDO1ru3UDvviizENJiYiIyDgMPxWopBo7NmoVZvy5CZt/nItle5dDJKhZj4eIiKiScM1PBWoX6AmZm4PWoy+fzPtYuXspOl4/CwDodeVvvHjjH7QLHGCuYRIREVkVzvxUIBuxCPMGBmtet7sRj72b3tQEn1yxDRb2/B+6RbzB3VlERESVhDM/FSy0qR/Wv9QSF6bPxaT9G2ErqAEAchcvRAx7Dy9OeZH1eIiIiCoRw09Fu38fobMnIHTfLk1TUquOSFn7Jb5s14gzPkRERJWM4acipaYC7dsDiYn5r0UiYPZsBMyZgwAbG/OOjYiIyEox/FQkb2+gQ4f88OPlBWzZAoSGmntUREREVo3hpyKJRMBnnwG2tsCiRUCdOuYeERERkdVj+ClPFy4At28DPXs+aXNxAb7+2nxjIiIiIi3c6l5etm4F2rTJP5E9KcncoyEiIqISMPyYSqkEJk8Ghg/PP5D0/n1g9mxzj4qIiIhKwMdeprh2DXjhBSAm5knbK68A69aZb0xERERUKs78GOu334CQkCfBRyLJX9y8eTPg7GzesREREVGJGH4MpVLlP9bq1w9IS8tvCwwEjh8Hxo7N3+FFREREVRYfexnqxReBn3568nrgQGDTJsDDw2xDIiIiIv0x/OggCAIAQKFQFH9zyJD88CMWA/PmAW++mT/bo+taIiIiqjQF39sF3+MlEQllXWGFbt68CX9/f3MPg4iIiIxw48YN1K5du8T3GX50UKvVuH37NlxdXSGq4DU8CoUC/v7+uHHjBtzc3Cr0s6oz/p70w9+T/vi70g9/T/rh70l/Ffm7EgQBGRkZqFmzJsTikpc187GXDmKxuNTEWBHc3Nz4H4we+HvSD39P+uPvSj/8PemHvyf9VdTvSiqVlnkNd3sRERGRVWH4ISIiIqvC8GNmEokEc+fOhUQiMfdQqjT+nvTD35P++LvSD39P+uHvSX9V4XfFBc9ERERkVTjzQ0RERFaF4YeIiIisCsMPERERWRWGHyIiIrIqDD9VzN69e9G+fXs4OjrCw8MD4eHh5h5SlaVUKtGyZUuIRCLExcWZezhVTlJSEl577TUEBgbC0dER9erVw9y5c5GTk2PuoZndmjVrEBAQAAcHB7Rv3x4xMTHmHlKVExkZibZt28LV1RW+vr4IDw/HxYsXzT2sKm/JkiUQiUSYMmWKuYdS5dy6dQsjRoyAl5cXHB0d0axZM/zzzz9mGQvDTxXy888/45VXXsGYMWPw77//4tixY3jppZfMPawqa8aMGahZs6a5h1FlXbhwAWq1Ghs2bMC5c+fwySefYP369Xj33XfNPTSz2rZtG6ZNm4a5c+ciNjYWLVq0QJ8+fZCSkmLuoVUpf/75JyZOnIgTJ05g3759yM3NxbPPPousrCxzD63K+vvvv7FhwwY0b97c3EOpcu7fv4/OnTvDzs4Ov/32GxISErBs2TJ4eHiYZ0ACVQm5ublCrVq1hC+++MLcQ6kWfv31V6FRo0bCuXPnBADC6dOnzT2kamHp0qVCYGCguYdhVu3atRMmTpyoea1SqYSaNWsKkZGRZhxV1ZeSkiIAEP78809zD6VKysjIEBo0aCDs27dP6Natm/DWW2+Ze0hVysyZM4UuXbqYexganPmpImJjY3Hr1i2IxWK0atUKfn5+6Nu3L+Lj4809tCrnzp07GDt2LL755hs4OTmZezjVSnp6Ojw9Pc09DLPJycnBqVOn0KtXL02bWCxGr169EB0dbcaRVX3p6ekAYNX//pRm4sSJ6N+/v9a/W/TErl270KZNGzz//PPw9fVFq1at8Pnnn5ttPAw/VcTVq1cBAPPmzcP777+PPXv2wMPDA927d0daWpqZR1d1CIKA0aNHY9y4cWjTpo25h1OtXL58GatWrcIbb7xh7qGYzd27d6FSqVCjRg2t9ho1aiA5OdlMo6r61Go1pkyZgs6dO6Np06bmHk6Vs3XrVsTGxiIyMtLcQ6myrl69inXr1qFBgwb4/fffMX78eLz55pvYvHmzWcbD8FPBZs2aBZFIVOpPwdoMAHjvvfcwZMgQtG7dGhs3boRIJMKPP/5o5j9FxdP397Rq1SpkZGQgIiLC3EM2G31/V4XdunULoaGheP755zF27FgzjZyqq4kTJyI+Ph5bt24191CqnBs3buCtt97Ct99+CwcHB3MPp8pSq9UICQnB4sWL0apVK7z++usYO3Ys1q9fb5bx2JrlU63I22+/jdGjR5d6TVBQEORyOQAgODhY0y6RSBAUFITr169X5BCrBH1/TwcPHkR0dHSxM2HatGmDl19+2Wz/F1GZ9P1dFbh9+zZ69OiBTp064bPPPqvg0VVt3t7esLGxwZ07d7Ta79y5A5lMZqZRVW2TJk3Cnj17cOTIEdSuXdvcw6lyTp06hZSUFISEhGjaVCoVjhw5gtWrV0OpVMLGxsaMI6wa/Pz8tL7fAKBx48b4+eefzTIehp8K5uPjAx8fnzKva926NSQSCS5evIguXboAAHJzc5GUlIS6detW9DDNTt/f08qVK7Fo0SLN69u3b6NPnz7Ytm0b2rdvX5FDrDL0/V0B+TM+PXr00MwkisXWPdlrb2+P1q1b48CBA5oyEmq1GgcOHMCkSZPMO7gqRhAETJ48Gdu3b8fhw4cRGBho7iFVSc888wzOnj2r1TZmzBg0atQIM2fOZPB5rHPnzsVKJfz3339m+35j+Kki3NzcMG7cOMydOxf+/v6oW7cuPvroIwDA888/b+bRVR116tTReu3i4gIAqFevHv+vtIhbt26he/fuqFu3Lj7++GOkpqZq3rPmWY5p06Zh1KhRaNOmDdq1a4cVK1YgKysLY8aMMffQqpSJEyfiu+++w86dO+Hq6qpZEyWVSuHo6Gjm0VUdrq6uxdZBOTs7w8vLi+ujCpk6dSo6deqExYsX44UXXkBMTAw+++wzs81GM/xUIR999BFsbW3xyiuv4NGjR2jfvj0OHjxovjoIVK3t27cPly9fxuXLl4sFQ0EQzDQq8xs2bBhSU1MxZ84cJCcno2XLloiKiiq2CNrarVu3DgDQvXt3rfaNGzeW+diVqKi2bdti+/btiIiIwIIFCxAYGIgVK1bg5ZdfNst4RII1/y1IREREVse6FwAQERGR1WH4ISIiIqvC8ENERERWheGHiIiIrArDDxEREVkVhh8iIiKyKgw/REREZFUYfoiIiMiqMPwQERGRVWH4ISIiIqvC8ENERERWheGHiCze999/D0dHR8jlck3bmDFj0Lx5c6Snp5txZERkDjzYlIgsniAIaNmyJZ5++mmsWrUKc+fOxVdffYUTJ06gVq1a5h4eEVUyW3MPgIiooolEInzwwQcYOnQoZDIZVq1ahb/++ovBh8hKceaHiKxGSEgIzp07hz/++APdunUz93CIyEy45oeIrEJUVBQuXLgAlUqFGjVqmHs4RGRGnPkhIosXGxuL7t27Y8OGDdi0aRPc3Nzw448/mntYRGQmXPNDRBYtKSkJ/fv3x7vvvovhw4cjKCgIHTt2RGxsLEJCQsw9PCIyA878EJHFSktLQ6dOndC9e3esX79e096/f3+oVCpERUWZcXREZC4MP0RERGRVuOCZiIiIrArDDxEREVkVhh8iIiKyKgw/REREZFUYfoiIiMiqMPwQERGRVWH4ISIiIqvC8ENERERWheGHiIiIrArDDxEREVkVhh8iIiKyKv8HoFaJxhSZwGYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(data[:,0], data[:,1], label=\"Data\")\n",
    "plt.axline((0, alpha), slope=beta, c='r', lw=2, ls='--', label=\"Line fit\")\n",
    "plt.xlabel(r\"$x$\")\n",
    "plt.ylabel(r\"$y$\")\n",
    "plt.legend(frameon=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd74e2d4-b75b-4235-8541-d5565ec4d7a8",
   "metadata": {},
   "source": [
    "Try comparing the fit result (values of $\\alpha$ and $\\beta$) with the simulated values in the beginning of the notebook. Are they close? Which criteria should you use to establish that they are close?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c6a1998f-81eb-4a32-8d00-04639ca777c2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.007394708"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ec33c149-7d6f-433b-872a-ca97c0b7d7d6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.0088975646258267"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91f0020f-6be7-4d3f-898c-2b062f825959",
   "metadata": {},
   "source": [
    "# Why is this a terrible idea?\n",
    "\n",
    "The previous calculation assumes that $p(\\alpha)$ and $p(\\beta)$, that is, the priors, are uniform. That means that any value of $\\alpha$ and $\\beta$ is acceptable. As we have discussed in the previous lecture, this leads to a fit that is very vulnerable to outliers. In the next notebook we will try to fix this by changing the priors to give a higher weight to values of $\\alpha$ and $\\beta$ close to zero, by setting $p(\\alpha)$ and $p(\\beta)$ to a Gaussian distribution, making an assumption on $\\alpha$ and $\\beta$ before looking at data. This method is called *regularization* and it starts from the assumption that, if the data has been previously *standardized*, we tend not to get wildly large values of a single parameter.\n",
    "\n",
    "Nevertheless, the calculation above can be changed to consider a Gaussian $p(\\alpha)$ and $p(\\beta)$ without using the methods described in the next notebook. Feel free to try redoing the calculation at home considering that.\n",
    "\n",
    "There are many other regularization techniques that make different assumptions on the parameters before looking at data. The choice of this method is strongly dependent on the problem at hand and it applies not only to line fitting, but fits in general.\n",
    "\n",
    "Try for example to add a data point in the variable `data` that is very far off to simulate an outlier. What happens with the fit? Are the fit parameters still close?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cc1b58d-08fc-4f83-8aa4-e2979aa92c89",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "2bdcc1c0-427c-4f6c-83d3-2f39a498cfa1",
   "metadata": {},
   "source": [
    "Another problem with the procedure above is that it does not calculate the errors on $\\alpha$ and $\\beta$. We will also cover this in the next notebook.\n",
    "\n",
    "One way of estimating uncertainties from the previous calculations is by making the approximation that the posterior probability $p(\\alpha, \\beta|\\mathcal{D})$ is approximately Gaussian around the best values of $\\alpha$ and $\\beta$ (this is called the *Laplace approximation*). In the example above, this is exactly true, but depending on the likelihood and priors this is not always so, although the approximation tends to be relatively good in several common cases. If we assume the posterior is Gaussian, the variance of the posterior can be taken as a measure of the uncertainty squared of $\\alpha$ and $\\beta$. How could we calculate that variance?\n",
    "\n",
    "One theorem called the Cramer Rao Lower Bound shows that the covariance matrix of the \n",
    "It can be shown, furthermore, that the covariance matrix of the posterior, $C$ satisfies:\n",
    "\n",
    "$C - I^{-1}$ is semi-positive definite,\n",
    "\n",
    "where $I_{ij} = -\\mathbb{E}[\\frac{\\partial^2 L}{\\partial \\theta_i \\partial \\theta_j}]$ (this is called the *Fisher information matrix*) and $L$ is the negative log-posterior, as above.\n",
    "\n",
    "This means, within the context of the calculation above, that the variance of the posterior could be said to have a *lower bound* of at least $-(\\frac{\\partial^2 L}{\\partial^2 \\theta})^{-1}$, calculated at the optimal values of the parameters.\n",
    "\n",
    "By taking this second derivative we could get estimates to the variances of the parameters. This is often how curve-fitting algorithms establish approximate uncertainties for the parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0812975c-a8fc-4d9e-98bc-d3851e75be70",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.12"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}