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": "code",
   "execution_count": 3,
   "id": "bc7938ee-c442-428a-984b-47975605b850",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.concatenate([data,\n",
    "                       np.array([[10.0, -4.0],\n",
    "                                 [5.0, -5.0]])\n",
    "                      ], axis=0)"
   ]
  },
  {
   "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": 4,
   "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": 5,
   "id": "cabdfc12-3033-4a2f-96ed-4ab1ea8312ab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f8d3c21ed20>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAGwCAYAAABGogSnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABa4UlEQVR4nO3dd1hT1+MG8DesAAJRRAgoCqJVEbdVcY+6xd06Wkdrrbuu/upqBW3ddXyrdbRaO9Rqq7iqpWrde1cRtQ5cGKSABkRZyf39kRKMrCQkZL2f5+EpOTn33kNMycu5Z4gEQRBAREREZCPsTN0AIiIiopLE8ENEREQ2heGHiIiIbArDDxEREdkUhh8iIiKyKQw/REREZFMYfoiIiMimOJi6AeZIqVTi8ePHcHd3h0gkMnVziIiISAuCICA1NRV+fn6wsyu4f4fhJx+PHz+Gv7+/qZtBREREenj48CEqVKhQ4PMMP/lwd3cHoHrxPDw8TNwaIiIiHWVnA6tWAcOHA87Opm5NiUlJSYG/v7/6c7wgDD/5yLnV5eHhwfBDRESWRS4HBgwA9u0Drl0DNm4EbGwIR1FDVjjgmYiIyFrcvQs0baoKPgCwdSvw99+mbZMZYvghIiKyBsePA40bAzExqseensCBA0DduiZtljli+CEiIrJ0P/8MtGsHJCaqHlerBpw5A7Rsadp2mSmGHyIiIkulVAKffQYMHgxkZqrK3noLOHUKqFLFtG0zYxzwTEREZIlevlSFnq1bc8tGjACWLwccHU3XLgvAnh8iIiJL5OgIpKSovrezA5YtU01vZ/ApEsMPERGRJXJwALZsARo1AnbtAsaPt7kp7fribS8iIiJLkZYGlCqV+7h0aeD0aYYeHbHnh4iIyNwJAvDVV0DNmoBMpvkcg4/O2PNDRERkRAqlgLOxyUhITYe3uzMaBXrC3k6HwJKZCYweDaxbp3rcowdw9KhNbVthaAw/RERERhIVLcOs3TGQydPVZb4SZ4SHBaNTiG/RJ0hOBvr0AQ4fzi3r1g0Qiw3fWBvC215ERERGEBUtw6gNFzWCDwDEy9MxasNFREXLCjjyP//8AzRpkht8xGLgl1+AmTN5q6uYGH5MRKEUcOpOEnZejsOpO0lQKAWjXm/o0KEQiUQQiURwdHSEj48P2rdvj++//x5KpVLr8/zwww8oXbq08RpKRGQFFEoBs3bHIL/f7Dlls3bHFPy7/+BBVfC5dUv12NtbFYL69zdCa20Pb3uZQLG7QfXUqVMnrF+/HgqFAk+ePEFUVBTGjx+PrVu3YteuXXBw4NuBiMgQzsYm5+nxeZUAQCZPx9nYZIQGldV8cu1aYNQoIDtb9bhWLWD3bqBSJeM12Maw56eEFbsbtBjEYjGkUinKly+P+vXrY/r06di5cyf++OMP/PDDDwCAJUuWoFatWihVqhT8/f0xevRoPH/+HABw+PBhvP/++5DL5epepIiICADAhg0b0LBhQ7i7u0MqlWLgwIFISEgw2s9CRGTOElILDj6F1jt1Chg+PDf4dO0KnDjB4GNgDD8lqNjdoEbQtm1b1KlTB5GRkQAAOzs7fP3114iOjsaPP/6IgwcP4tNPPwUANG3aFMuWLYOHhwdkMhlkMhk++eQTAEBmZia++OIL/P3339ixYwdiY2MxdOjQEvs5iIjMibe7djOx8tQLDVUtVggAEycCO3cC7u4Gbh3xPkcJKlY3qBFVr14dV65cAQBMmDBBXR4YGIgvvvgCo0aNwsqVK+Hk5ASJRAKRSASpVKpxjg8++ED9feXKlfH111+jUaNGeP78Odzc3Erk5yAiMheNAj3hK3FGvDw93z94RQCkEtW09zwWL1ZtTtqtm7GbabPY81OC9O4GNTJBECD6b+bAoUOH0L59e5QvXx7u7u4YPHgwkpKSkJaWVug5Ll26hB49eqBSpUpwd3dH69atAQAPHjwwdvOJiMyOvZ0I4WHBAFRB51U5j8PDgmF//pxqPI/GwfYMPkbG8FOC9O4GNbLr168jMDAQ9+/fR5cuXRASEoJt27bhwoUL+OabbwAAWVlZBR6flpaGDh06wM3NDRs2bMC5c+ewfft2AKrbYUREtqhTiC9WvVcfUonm73SpxBmr3quPTjHHgFatgH79gAsXTNRK28TbXiWoWN2gRnLw4EFcvXoVEydOxPnz55GdnY3FixfDzk6Vi3/99VeN+k5OTlAoFBplN27cQGJiIubPnw9/f38AwPnz50vmByAiMmOdQnzRPliqucJzQBnYz52jWq8nx6JFwObNpmuojWHPTwnSuhtUl2XPdZCRkYH4+HjExcXh4sWLmDt3Lnr06IFu3bph8ODBCAoKQnZ2NpYvX467d+/i559/xurVqzXOERAQgOfPn+Ovv/5CYmIiXrx4gYoVK8LJyUl93K5du/DFF18Y5WcgIrI09nYihAaVRY+65RFavhTshwzWDD7vvw/89JPpGmiDGH5KWJHdoEZc5ycqKgq+vr4ICAhAp06dcOjQIXz99dfYuXMn7O3tUbduXSxZsgQLFixASEgINm7ciHnz5mmco2nTphg5ciT69euHcuXKYeHChShXrhx++OEH/PbbbwgODsb8+fPx1VdfGe3nICKySAkJQNu2wMaNqsciEbBggWrPLicn07bNxogEQSi5edUWIiUlBRKJBHK5HB4eHka5RrE3uiMiIssRHQ2EhQH37qkeu7oCGzYAvXqZtFnWRtvPb475MZGcblAiIrJyf/4JvP02kJqqeuznp5rhVb++adtlw3jbi4iIyJicnICXL1XfN2gAnD3L4GNiDD9ERETG1KYNsGoV0KcPcPQoUL68qVtk8zjmJx8lMeaHiIis1PPnQKlSqgHNrxKEvGVkUNp+frPnh4iIyFDu3gUaNQLmzMn7HIOP2WD4ISIiMoRjx1TB5/p14PPPga1bTd0iKgDDDxERUXH99BPQrh2QlKR6XL06UK+eadtEBWL4ISIi0pdSCUyfDgwZAuTsgdi+PXDqFBAUZNq2UYEYfggAIBKJsGPHDqNf59tvv4W/vz/s7OywbNkyREREoG7duka/LhGRwb14oVq/59WV8EeOBPbsAUqXNlmzqGic7ZUPa5ztNXToUDx79qzAgBMfH48yZcpALBYbrQ0pKSnw8vLCkiVL0KdPH0gkEiiVSmRkZKBs2bJatZOIyCw8fgx07567G7udHbB0KTBuHAc2mxBXeCadSKVSo1/jwYMHyMrKQteuXeHrm7uHmZubm9GvTURkUEOG5AYfd3dgyxagc2fTtom0ZlG3vebNm4c333wT7u7u8Pb2Rs+ePXHz5k2NOoIgICIiAn5+fnBxcUHr1q1x7do1E7XYcrx62+vevXsQiUSIjIxEmzZt4Orqijp16uDUqVMax5w8eRItW7aEi4sL/P398fHHHyMtLS3f8//www+oVasWAKBy5coQiUS4d++exm2viIgI/Pjjj9i5cydEIhFEIhEOHz5srB+ZiEh/q1cDZcsClSoBJ08y+FgYiwo/R44cwZgxY3D69Gns378f2dnZ6NChg8YH7sKFC7FkyRKsWLEC586dg1QqRfv27ZGas6cKaW3GjBn45JNPcPnyZbzxxhsYMGAAsrOzAQBXr15Fx44d0bt3b1y5cgVbtmzB8ePHMXbs2HzP1a9fPxw4cAAAcPbsWchkMvj7+2vU+eSTT/DOO++gU6dOkMlkkMlkaNq0qXF/SCIifQQFAX/8AZw5A4SEmLo1pCOLuu0VFRWl8Xj9+vXw9vbGhQsX0LJlSwiCgGXLlmHGjBno3bs3AODHH3+Ej48PNm3ahBEjRuR73oyMDGRkZKgfp6Sk6N64JUtUX0WpXx/YtUuzrHt34OLFoo+dNEn1VUI++eQTdO3aFQAwa9Ys1KxZE7dv30b16tWxaNEiDBw4EBMmTAAAVK1aFV9//TVatWqFVatWwdnZWeNcLi4u6nE95cqVy/c2m5ubG1xcXJCRkVEit+GIiLSSmQl89RUwcSLg4pJb/uabpmsTFYtFhZ/XyeVyAICnpycAIDY2FvHx8ejQoYO6jlgsRqtWrXDy5MkCw8+8efMwa9as4jUmJQWIiyu63mu9HQCAf//V7lh9Qlkx1K5dW/19zhidhIQEVK9eHRcuXMDt27exceNGdR1BEKBUKhEbG4saNWqUaFuJiIwiKQno3Vu1J9fffwO//KIa3EwWzWLDjyAImDRpEpo3b46Q/7oc4+PjAQA+Pj4adX18fHD//v0CzzVt2jRMeqVHJSUlJc8tmSJ5eGi3WV25cvmXaXNsCc88c3R0VH8v+m/2glKpVP93xIgR+Pjjj/McV7FixZJpIBGRFhRKAWdjk5GQmg5vd2c0CvSEvZ0WM7Ju3gS6dgXu3FE93rULuHYN+G/8Ilkuiw0/Y8eOxZUrV3D8+PE8z4lem2YoCEKesleJxeLiT/Euzi2p12+DWYD69evj2rVrqFKlikHP6+TkBIVCYdBzEpHtioqWYdbuGMjk6eoyX4kzwsOC0SnEt+AD//oL6NsXePZM9djHB9i5k8HHSlhk3924ceOwa9cuHDp0CBUqVFCX54wTyekBypGQkJCnN8gWyeVyXL58WePrwYMHep1rypQpOHXqFMaMGYPLly/j1q1b2LVrF8aNG1esNgYEBODKlSu4efMmEhMTkZWzYioRkY6iomUYteGiRvABgHh5OkZtuIioaFn+B377LdCxY27wqV0bOHsWaNzYuA2mEmNR4UcQBIwdOxaRkZE4ePAgAgMDNZ4PDAyEVCrF/v371WWZmZk4cuQIZw0BOHz4MOrVq6fxNXPmTL3OVbt2bRw5cgS3bt1CixYtUK9ePXz++eca6/foY/jw4ahWrRoaNmyIcuXK4cSJE8U6HxHZJoVSwKzdMchvFd+cslm7Y6BQvlJDoVD14I8YofoeALp1A44fB3g736pY1ArPo0ePxqZNm7Bz505Uq1ZNXS6RSODy3wj8BQsWYN68eVi/fj2qVq2KuXPn4vDhw7h58ybc3d21uo41rvBMRGRLTt1JwoDvThdZ75fhTRAaVBZISwP69VNtTZFj0iRg4ULA3t6ILSVDssoVnletWgUAaN26tUb5+vXrMXToUADAp59+ipcvX2L06NF4+vQpGjdujH379mkdfIiIyPIlpKYXXenVemJxbm+PgwOwciUwfLiRWkemZlHhR5tOKpFIhIiICERERBi/QUREZJa83Z2LrvRqPQcHYPNmICwMiIgA2rY1XuPI5Cwq/BAREWmjUaAnfCXOiJen5zvuRwSgsliBRoGeuYUSCXDkCDcmtQEWNeCZiIhIG/Z2IoSHBQNQBZ1XiQQB4078gt/XjIS97PFrTzL42AKGHyIiskqdQnyx6r36kEpyb4GJszOxJmoJJh3fCJcnMtX2QunajQ8i68HbXkREZLU6hfiifbAUZ2OTIb/3EM3+7yO4XzmnelIkAgYMUA12JpvC8ENERFbN3k6E0JcyYFhPIGerI1dXYNMmoEcPk7aNTIPhh4iIrNvevUD//kBqqupx+fLA7t1AvXqmbReZDMf8EBGRdRIE4OuvVdPXc4JPw4aqrSoYfGwaww8REVmnEyeA8eMBpVL1uG9f1VR2Pz/TtotMjuGHiIisU/Pmqi0qAGDGDGDLFtVYH7J5HPNDRETWa+FCoHNn4K23TN0SMiPs+SEiIutw7Biwfbtmmb09gw/lwfBDRESW74cfgHbtgIEDgXPnTN0aMnMMP0REZLmUSmDaNOD994GsLNVqzStWmLpVZOY45oeIiCxTWhoweDAQGZlbNmYMsGyZyZpEloHhh4iILE9cnGpfrosXVY/t7ID//Q8YO9a07SKLwPBDRESW5eJF1cKFj//bkd3dHfj1V6BTJ9O2iywGx/wQEZHl2L0baNFCHXzS/PxxecsfUHToaOKGkSVh+CEiIsvh4QFlZhYA4Hz5GmjZex56HnmG5gsOIipaZuLGkaVg+CEiIqNQKAWcupOEnZfjcOpOEhRKodjnjCr7Bqa0H43twa3xbv85SCpVGgAQL0/HqA0XGYBIKyJBEIr/brQyKSkpkEgkkMvl8PDwMHVziIgsTlS0DLN2x0AmT1eXeZZywpc9QtCltq/2J5LLAQ8PQCSCQimg+YKDqnMKAiASaVQVAZBKnHF8SlvY24nyPx9ZNW0/v9nzQ0REBhUVLcOoDRc1gg8AJKdlYvSmi5i3NwaAFj1DN24ADRoAs2cDAM7GJueeU5Q33AgAZPJ0nI1NNvjPRNaFs72IiMhgFEoBs3bHoLBbCmuOxkIpAL9fkWkEJF+JM8LDgtEpxBc4cEC1C7tcDkREADVqIOGNZlq1ISE1vehKZNMYfoiIyGA0emcK8d2x2DxlOeN2fneKRs25MwCFQvVEnTpAkybwznLWqg3e7trVI9vF8ENERAZTnF4XkVKBGYe+R83zO3MLw8KATZsANzc0UgrwlTgjXp6eb89SzpifRoGeereBbAPH/BARkcHo2+vilvEC30V+iWGvBp/Jk1W7tLu5AQDs7UQIDwsGoAo6r8p5HB4WzMHOVCSGHyIiMphGgZ7wLOWk0zHl5QnYuuH/0O6Oajf2LDt7XP58AfDVV4C9vUbdTiG+WPVefUglmiFLKnHGqvfqq8YLERWBt72IiMhg7O1EmNW9Jsb9cknrYxb+sQzVE+8DAOTiUhjZazo+HjKswPqdQnzRPliKs7HJSEhNh7e76lYXe3xIWww/RERkUF5uYp3qT+n0MXb+NAkpzqUwrG8EXgYGFTlux95OhNCgssVpJtkwhh8iIjKo/THxOtV/VFqKwf2+wGOPcnjm4oFVHLdDRsYxP0REpJf8FimMipbh+xP3CjxGnJWBicc2wjkrXWNs0DWfIDhLvTluh0oEe36IiEhnUdEyROy6hviUDHVZGVdHZGQrCzym3POn+DbyS9ST3UStlDg0v/AXLjyUc9wOlTiGHyIi0klUtAwjN1zMU/70RVaBx1RPiMXabbNRIeVfAEDLO+fhcPsfhAYHG62dRAXhbS8iItKaQilgauRVnY5pc+cctm78VB185F4+cDh5AmDwIRNhzw8REWnt9J0kPCukh0eDIOCD87sw49A62Auq22GXfasie9t2NKxb04itJCocww8RERVKoRRwNjYZ8Snp2Hr+oVbHOCiyMevAarx7OUpdtrdaMywcMA1/NWaPD5kWww8RERUoKlqGWbtjtNqsNEepjBdYs30Omt//W122PLQflrZ4Fyv71OegZjI5hh8iIspXVLQMozZczHcT0cK8dBQjw0E1jT3D3gFTOo/HmaadsTIsmNPYySww/BARUR4KpYBZu2N0Dj4AoLSzx8dh/4c12+dgVZtBGD19CBYHlWWPD5kNzvYiIqI8zsYm63SrS/IyVeNxmtgV7/Wfg0bvhqFZVS8GHzIrDD9ERDYsv1WaAeCAlltUiAQlphz+Afu+HwNpSqLGc2VcHTG2bVWDt5mouHjbi4jIRuU3mNlX4ozPuwZj++W4Io93yUzH0j2L0emfUwCAddtmo9egxch0cIQIwLzetdjjQ2aJ4YeIyAYVNJg5Xp6O0Zvyrt78Op/URKzb9gVCntwBAGSL7LC5TgdkOjjCV+KMcA5uJjPG8ENEZGMKG8yszQDnkPjbWLttNqTPkwEAKU6uGNtjCo5WboDPu9bA0GaB7PEhs8bwQ0RkhXIWJsxv01BdBzO/quPNk1j2+2K4ZKs2NH0g8cEHfcNx26siAMDLXczgQ2aP4YeIyMoUNJYn51ZUQqoewUcQMOrMVkw58qO66Fz5YIzoPQPJrhJ1mbe7c7HaTlQSONuLiMiK5Izleb1nRyZPx8gNFxEVLdMroDR6dE0j+Gyr2Qbv9p+jDj4iqAJWo0DPYrWfqCQw/BARWQltFiac9OvfyM5WwsddrNO5z/qHYHWj3gCAhS0HY3LXSch0cNSoEx4WzFteZBFEgiDos4CnVUtJSYFEIoFcLoeHh4epm0NEpJVTd5Iw4LvTWtUViQBdf/vbKRVoEHcd5/xDNMo9Szlibq9anN1FJqft5zd7foiIrIQuY3mKCj7NYy+hy43jGmVKO/s8wQcAPu9Wk8GHLAoHPBMRWQlDDTZ+79JeROxfDYWdPR57lMNlv2qF1pd6cJAzWRb2/BARWYlGgZ4o7eJYdMUC2CsVCD+wBl/uWwkHQQmxIgv9/v6zwPoc5EyWiuGHiMhK2NuJ0KJqWb2Odct4gbXbZuP9C7vVZasb9caMjmPyrZ8zrJmDnMkS8bYXEZGVUCgFnLv3VOfjKsifYO3W2aieeB8AkGVnjxkdxuDXOh3UdexEgPKVcUJSbmFBFozhh4jISpyNTUZ8SoZOx9SPu441kXNQ7sUzAMBTZ3eM6jUNpyvWBpDbw7NiQD2UKSXOd8VoIkvD8ENEZCV0Xbm5843jWPb7YogVWQCAO57lMazPTNzzLK+uwx4eskYMP0REVkLX2V7/upVBzlamJyrVxqie05Hi7IYyrg6YGRYCqQd7eMg6MfwQEVmJpFTdbnmdr1AT0zqNQ8NH1zGz/Uhk26s+Ep6+yIbUwxmhQfoNniYydww/REQW6tWd273cxPh8V3Sh9Uu/TIHc2Q2CKHeib2RIO0SGtMtTV6/NT4ksBMMPEZEFym/n9sJUT4jFuq2zsS2kLZa0HFRkfe7OTtaM6/wQEVmYgnZuL0jb22exdeOnKJ/6Lz4+tQVdrx8rsC4XLiRbwJ4fIiILos3O7WqCgGHndmDGoe9h998Rl32r4qx/zXyrc+FCshUMP0REFuT0nSStenwcFNmYvX81Bv4dpS77vXoLTO4yARmO4nyP4bR2shUMP0REFiIqWoap264WWc8j/TlW7ZiLZvevqMv+13QAljUfoDHYOUdpV0d8M6A+mgSVZY8P2QSGHyIiC5Azzqeo210ByXFYt202gpLjAAAZ9g74tPN47KzZpsBjnr3Igp2diMGHbAbDDxGRmVMoBUyNvKrVOJ/5UcvVwSfRVYKPen2GixVqFHkcp7aTLWH4ISIyM6+u3+Pt7oxTtxPx7EWWVsd+0nUitv80GcmuHhjWNxyPJD5aHcep7WRLGH6IiEzk9ZDTKNAT+2PidVq/53WPJD4Y1O8LPJL44LnYtcj6IqgGOnNqO9kSi1vn5+jRowgLC4Ofnx9EIhF27Nih8bwgCIiIiICfnx9cXFzQunVrXLt2zTSNJSIqQFS0DM0XHMSA705j/ObLGPDdadScGYWROqzf45r5ElMPr4dLpmb9G96B+Qaf10f0cGo72SqLCz9paWmoU6cOVqxYke/zCxcuxJIlS7BixQqcO3cOUqkU7du3R2pqagm3lIgof1HRsnxDTnq2UutzSFMS8dvGKRh5ZhuW7lkMkVD4sX3rl4dUonlrSypxxqr36nNqO9kci7vt1blzZ3Tu3Dnf5wRBwLJlyzBjxgz07t0bAPDjjz/Cx8cHmzZtwogRI0qyqUREeeQMXi6OWrJbWBv5BXyeJwMAmt77G5WT43CnrH++9e1EwNzetWFvJ8pzm409PmSLLC78FCY2Nhbx8fHo0KGDukwsFqNVq1Y4efJkgeEnIyMDGRm5uyGnpKQYva1EZJtWHLyt9eDl/HS+cRxL9iyFS7bqd9YDiQ8+6BteYPABgOEtAuHkoOro507tRBZ426sw8fHxAAAfH83ZDT4+Purn8jNv3jxIJBL1l79/wb9EiIj0pVAKWH8iVr+DBQGjT/2KVTvnq4PP2QrB6Dl4CW57Vcz3EDsRMKJlIKZ1Cda3yURWyap6fnKIRJrduIIg5Cl71bRp0zBp0iT145SUFAYgIjK4s7HJePZS914fp+wszIv6Gn2uHVKXbQtpi2kdxyHTwTFP/VZvlEPLql4YFBqg7vEholxWFX6kUikAVQ+Qr2/uAL6EhIQ8vUGvEovFEIvz3+uGiMhQ9FlI0D0jDeu2zkKjRzHqsoUtB2Nlk7eB1/6oE4mAbwbUR5faHMBMVBir+pMgMDAQUqkU+/fvV5dlZmbiyJEjaNq0qQlbRkQEHIgp+PZ7QdIcnSF3dgcAvHQQY2TPaVgZ+k6e4AMAQ0MrMfgQacHien6eP3+O27dvqx/Hxsbi8uXL8PT0RMWKFTFhwgTMnTsXVatWRdWqVTF37ly4urpi4MCBJmw1EdmazGwlfjx5D+fuJcHVyQE9avvh9yu6hx+lnT3Gh32CFTsXYGnzd3HVt2qBdTvUZPAh0oZIEARttosxG4cPH0abNnk36BsyZAh++OEHCIKAWbNmYc2aNXj69CkaN26Mb775BiEhIVpfIyUlBRKJBHK5HB4eHoZsPhHZgHl7Y/DtsVjo+9u1bNozJJUqrXX9nFWaj09py6nrZNO0/fy2uPBTEhh+iEhf8/bGYM1R/WZ02SsV+OzgWnS9cRw9By/GYw9vrY4TAVyskAjaf35b1ZgfIiJTepmpwLfH9As+bhkvsHbbbLx/YTe8055i7bYv4KgoemaY1EPM4EOkI4sb80NEZI6iomX4v61X9LrVVeFZPNZtm41qiQ8AAJl2DljfoDuy7PNOY3/VxLfewNi2VXiri0hHDD9ERMUUFS3DqA0Xoc8YgvqPruPb7V/C64UcAPDU2R0je03HmYq1CjzGV+KM8LBg9vYQ6Ynhh4ioGBRKARG7YvQKPj2uHcLCP/4HsSIbAHDHszyG9ZmJe57lNepJPcQY0KgiArxKcU8uIgNg+CEiKoYVB28hPkW3xQtFghITj23Ex6e2qMtOVKqNUT2nI8XZTV02tk0VNKvixbBDZGAMP0REeoqKlmHpgVs6H9fkQbRG8NlUpxNmth+JbHvNX8lVfdy4ESmREXC2FxGRHhRKAbN2xxRdMR+nKtXGyiZ9oYQIX7T9ENM7jskTfADA2925uM0konyw54eISEsKpYCzsclISE1HYmoGZHLd9+rKsajlYOyrGorLftXyfd5XohrbQ0SGx/BDRFSInMDz5zUZtl2MQ2p6ts7naHf7DFwz07E7uJW6TBDZFRh8ACA8LJjjfIiMhOGHiKgAUdEyROy6hviUDP1OIAgYdm4HZhz6Hln29ojz8MbFCjUKPaS0qyPm967FaexERsTwQ0SUj6hoGUZuuKj38Q6KbMzevwoD//4TACBWZKPH9cMFhp9STvb4qGUQFy0kKgEMP0REr1EoBUyNvKr38ZKXqVi5cx6a3b+iLlvWbACWNRuYb/2edfywuF9dhh6iEsLwQ0T0mtN3kvDsRdH7auUnIDkO67bNRlByHAAgw94R/9dlAna9Mt7ndRU8XRl8iEoQww8R0WtO3U3U67gmD65g9fa5KJ3+HADwr2tpjOg9AxfLFz7Oh2v5EJUshh8iotfoszlpz2uHsGjvMjgqFQCAG16V8GHfmXgk8Sn0uDKujmhSmeGHqCQx/BARvSIqWoYfTt3T+bgHEimUItWtq0OVG2Bc9yl4LnYt8rh5vWvxlhdRCWP4ISL6T3FmeF2sUAOfdh6POrJbmNN2GBR29oXWL+PqiHmc0k5kEgw/RGSTchYvjE9JR/LzDJR2dcKXe7TfrsIr7SmSXCUQRLm7BO2s2QY7a7Yp9DgRgPHtqmBcuzfY40NkIgw/RGRzoqJlmLU7Ru/tKWrL/sHabV/gt1pvYVGrITod+83A+uhSm709RKbEjU2JyKbsvaK6taVv8Oly4zh+3TQV3mlPMeb0b+h2/ajWx058qyqDD5EZYM8PEdkEhVLA13/9g6//uq3fCQQBY079iv879rO66EyFmjhRqY5Wh/tKnDG2bVX9rk1EBsXwQ0RWLypahqmRV/VeuNApOwvzo75G72uH1GVbQ9phesexyHRwLPTYnFE93KiUyHww/BCRVYuKlmHUhovQY+keAEDZtGdYs30OGsZdV5fNbzUUqxv3AURFhxmpxBnhYcGc1UVkRhh+iMhqKZQCZu2O0Tv4VP33Pr7fNhv+8icAgJcOYkzoNhl/Vmuq1fGfd62Boc0C2eNDZGYYfojIKimUAn44Eav3wGYIAr7ct1IdfOLdPPFhn5mIllYp8lARVD0+DD5E5onhh4iszt4rj/HZzmgkp+k3xgcAIBJhYthk7PxpEuLdyuLDPp/jibuX1odzjA+R+WL4ISKrMm9vDNYcjTXIuR57eGNA/7mI8/DGSydnrY6ReogR0b0mx/gQmTGu80NEVmPvFZnewcc9Iw3hB9bANfOlRvltr4paB5+Jb72BE1PbMfgQmTn2/BCRVVAoBXy2M1qvY/2fxWPd1tl4I+kBKsgTMKLXdCiL2JvrVb6c0UVkURh+iMgqnI1NRnJaps7HNXx0DWsi56Dsy5T/Hseg4rN43PMsX+hxn3etAS93MbzdndEo0JPje4gsCMMPEVmFhFTdZ3X1ij6I+VFfQ6zIBgDc8ayAD/rOxP0yfoUeJ/UQcyYXkQVj+CEiq7D/2hOt64oEJSYd24hxp7aoy45XqoPRPachxdmtyOMjutdk8CGyYAw/RGTxhv90DvtjErSq65yVjsV7lqLrzRPqso11OyH8rZHIti/6V+IHzQI4tofIwjH8EJFF+/1ynNbBxyP9OX769XPUld0CAChEdviy7TCsb9Bdq60qAKB9sFTvthKReWD4ISKLpVAKmBJ5Rev6qWJXPHErC+AWnju5YFz3T3Eo6E2tj5d6iNEo0FOPlhKROWH4ISKLolAKOBubjHj5S2y/9AhpmUqtjxVEdpjQ7RMs+/0rLGnxHm6WC9Dp2hzrQ2QdGH6IyGJERcswa3eM9vt1CQK8nycjwb2suuilkzNG9P5Mp+u6OtljyTt1ONaHyEow/BCRRYiKlmHkhota13dUZOGLfavQ7vZZ9By8BHESb52vKXaww8hWQfi4XVX2+BBZEYYfIjJ7CqWAqZFXta4veZmK1TvmIvSB6pi122YjbMgyrWZzAUBpV0e83zQQY9tWYeghskIMP0Rk9pb/dQvPXmi3Q3tgchzWbZ2Fyk8fAwAy7B2xssnbWgefz7vW4AKGRFaO4YeIzNrvl+Ow7K9bWtUNvf83Vm+fC0lGGgDgX9fS+Kj3Z7hUvnqRx4oASCXODD5ENoDhh4hMKmf2VkJqusY+WQqlgAmbL2L3lXitztPv7z/x5b6VcFQqAADXywXgwz4ztRrrkxN1wsOCGXyIbADDDxGZTH6zt3wlzuhexxdbzj3Cs5dF3+qyUyow9fAP+OjcdnXZX0Fv4uOw/0Oa2FWrdki5KzuRTWH4ISKTiIqWYdSGixBeK5fJ07HmaKzW52ny4KpG8FnbsAfmtvkASjv7Io/tHOKDwaGB3JWdyMbYmboBRGR7FEoBs3bH5Ak++jgZUBfLQ/shW2SH6R3H4Mt2w7UKPgDQKcQXoUFlGXyIbAx7foioxOSM7zlx+1/tFyrUwpIW7+KPas0Q41NZp+O83Z0N1gYishwMP0RUInRenbkAXa8fg6MyGztqtlGXCSI7nYOPr8SZ+3QR2SiGHyIyuoLG9+hEEDDu5GZMPr4RGfYOeCTxxvkKNfU+HWd2EdkujvkhIqMyxPgecXYmlv6+GJOPb1Q9VmSj4z+n9DqXnQhYObAeZ3YR2TD2/BCRUZ2NTS7Wra6yac+wZvscNIy7DgBQQoQFrYdgTaM+ep1vxYD66FKbwYfIljH8EJFRJaTqH3ze+Pce1m37Av7yJwCAF45iTOw2GX++0bTIY10c7fAyS6l+7Mu1fIjoPww/RGRU+s6oanX3AlbsnA/3zJcAgHg3TwzrMxPXpFWKPNbVyR6XZ3bAhftP86wcTUTE8ENERtWgUhnYiQClDoN+3r6yD/OjVsBeUPXcXJFWwfDen+GJu5dWx49oGQQnBzuEBpXVp8lEZOUYfojIqC7cf6pT8AGAW16VkG1nD3uFEn+80RQTu01CuqP2PUgBXtpta0FEtonhh4iM6kCMdhuTvuqyXzVM7joRNRJi8VXLQRBEuk1M5eKFRFQYhh8iMpqoaBnWnbhXZD2f1EQkuHlqhJzfa7TE7zVa6nxNqYeYixcSUaF0Cj8PHz6Ev7+/sdpCRBZKoRRw+k4STt1NhACgtIsTPEs5Yfbuq0Ue2/DRNayJnINfa3fAgtZDi92WiO41ObCZiAqlU/ipXr06Jk2ahKlTp6JUqVLGahMRWZCoaBmmRl7FsxdZOh/bO/ovzItaDrEiG6PObEW0TxD21GihVztKOdlj8Tt1OJWdiIqk0430/fv3Y9++fahatSrWr19vrDYRkYWIipZh5IaLOgcfkaDEJ0d/wpI9SyFWZAMAjgbUw7HAenq1o3c9X1yJ6MjgQ0Ra0Sn8NG3aFGfOnMH8+fMxc+ZM1KtXD4cPHzZS04jInCmUAqZsvaLzcc5Z6fhmx3yMPfWruuznel3w/tsRSHF20/l8ZVwdsejterzVRURa02tvr8GDB+Off/5BWFgYunbtil69euH27duGbhsRmbHxmy9Bnp6t0zHeqUnYsmkauvxzEgCgENkh/K0R+Lz9KCjs7PVqx7zetRh8iEgnem9sKggCOnTogI8++gi7du1CSEgIJk+ejNTUVEO2j4jMiEIp4NSdJMzaFY3fr8h0OrbmkzvY+dMk1Im/BQBIdXLBsD4z8WODMECke3jxlThj9Xv1eauLiHSm04Dn1atX49y5czh37hyuX78Oe3t71K5dG2PGjEHdunWxceNGBAcHY/v27WjYsKGx2kxEJhAVLcOs3TH6bVIqCAg/sAa+z5MAAI88vPFB35n4p1yAXm2Z+FZVjG1blT0+RKQXkSAIWq+96u/vjyZNmqi/GjZsCLFYrFFn7ty52LRpE6Kjow3e2JKSkpICiUQCuVwODw8PUzeHyOSiomUYteEidFyoWYM0JRE7f56ERx7e+Kj3Z0gqVVqv84gASCXOOD6lLcMPEWnQ9vNbp/CjjSdPnsDPzw8KhcKQpy1RDD9EuRRKAc0XHNSvx+c1lZMeIU7ijQwHp2Kf65fhTbh3FxFp0PbzW+8xPwXx9vbGwYMHDX1aIjKRs7HJOgef0i9TMOfPFSiV8UKj/G7ZCgYJPgCQkFr8MEZEtsng21uIRCK0atXK0KclIhPRdW+uykmPsG7bLAQ+lUGamoThvT+DUs+ZXIXh/l1EpC+D9/yYi5UrVyIwMBDOzs5o0KABjh07ZuomEZmlnBlcOy/H4dSdJChe2YI9M1uJLecfan2upvcuY/vPkxH4VDUTrHb8LZRP+deg7RVBNdOL+3cRkb6scmPTLVu2YMKECVi5ciWaNWuGNWvWoHPnzoiJiUHFihVN3Twis5HfDC5fiTPCw4IBAJ/8dgXPM7QbvzfgchRm718FR6Wq/vVyARjWdyYee3gbrL05w5vDw4I52JmI9GbwAc/moHHjxqhfvz5WrVqlLqtRowZ69uyJefPmFXk8BzyTLTDEDC4AsFMqMP3Q9/jw/E512YGgNzE+7P+QJnbV+XzvNwtAh2ApnqZl4Is91/MNZlzbh4jyo+3nt9X1/GRmZuLChQuYOnWqRnmHDh1w8uTJfI/JyMhARkaG+nFKSopR20hkagqlgFm7Y4odfEplvMD/di/CW3fOqcvWNuyBuW0+0Gucz4iWgZjWJVj9uGOIL87GJiMhNR3e7qpbXezxIaLisrrwk5iYCIVCAR8fH41yHx8fxMfnP3Bz3rx5mDVrVkk0j8gs6DOD63WlX6bgl1+mo8a/9wAAWXb2mNl+FH6p20nnc3mWcsSXPULQpbafRrm9nYjT2YnI4Kwu/OQQvbZcviAIecpyTJs2DZMmTVI/TklJgb+/v1HbR2RKhpgmLnd2Q2wZP9T49x7k4lIY1XMaTgbU1erYNm+Uw0etgtijQ0QmYXXhx8vLC/b29nl6eRISEvL0BuUQi8V5VqomsmaGmCYuiOwwqdskKPbaY1nzgbhTVvs/GD5qFcQeHSIyGaub6u7k5IQGDRpg//79GuX79+9H06ZNTdQqIvPSKNATvhJn6NTXIggoL0/QKEp3dMa4HlN0Cj6cpk5EpmZ14QcAJk2ahLVr1+L777/H9evXMXHiRDx48AAjR440ddOIzIK9nQjhYcFaD3gWZ2fif7u/ws6fJqLCM90WPXwdp6kTkalZ3W0vAOjXrx+SkpIwe/ZsyGQyhISEYO/evahUqZKpm0ZkNpRK7ep5pT3Fmsg5aPD4BgBg7bYv0PX9r6HQcTZXGVdHzOtdi9PUicjkrHKdn+LiOj9kjRRKQT1t/E7Cc3x98HaRx7zx7z18v3U2KqSobne9cBRjQrdPsO+NUK2v26mmDwaFBqBJ5bLs8SEio7LZdX6IKK/8VnIuSus757B810K4Z74EAMjcyuLDvjNxzSdIp2v/ee0JetYrz+BDRGbDKsf8EFGunJWctQ4+goD3z+/Eum1fqIPP39Kq6DF4ic7BJ8es3TEae4YREZkSe36IrJiuKzk7KLIRcWAN3rv8h7psT7VmmNx1ItId9ZseLwCQydNxNjaZ09uJyCww/BBZMV1Xcm7y4KpG8Fke2g9LWrwLQVT8TmJDLKxIRGQIvO1FZMUOxOg2Lf14YD0sazYAmXYOmNR1Iha3HGSQ4AMYZmFFIiJDYM8PkRVSKAWcvpuEDWfu63zssmYD8Xv1FrjtVdEgbREBkHJhQyIyIww/RFYmKlqGiF0xiE8p+jZTn6t/AQC21WqXWygSGTT4AFzYkIjMC8MPkRWJipZh5IaLRdYTCUp8cvRnjDn9GzLtHPCwtA/O+ofofd3SLo4Y0jQAW8491AhdUokzwsOCubAhEZkVhh8iK6FQCpgaebXIei6Z6ViyZwk6/3MSAOCkzEaruxeKFX7m91Gt3Pxxu6rqhRS5WzsRmSuGHyIrcfpuEp69yCq0jk9qItZu+wK1ntwBAChEdpjVbjh+ahCm93UnvvWGumfH3k7E6exEZPYYfoisxKI/bxT6fM3421i3bTakz5MBAKlOLhjbYyqOVG5QrOsGeLkW63giopLG8ENkwXJmdS2KuoHLj+QF1uv4z0ks/X0xXLMyAAAPJT4Y1udz/FMuoNht4BR2IrI0DD9EZu7VDUlfHUcTFS3D1MirRd7qeu/iHny5f5X68fnyNTCi1wwklSpdrHZxCjsRWSqGHyIzlt+GpL4SZ3Sv44s1R2O1OsdVaRWkOzjBOTsT24NbY2rnj5Hh4FSsdnEKOxFZMpEgCNxt8DUpKSmQSCSQy+Xw8PAwdXPIRuVsSGqI/0G73DiOoKSHWN60PyDSLax4uTkhI1uJ1PRsdZkvp7ATkRnS9vObPT9EZkjXDUlfVUH+BI/dvaC0s1eX7a3eXO+2LB9QH40CPTmFnYisBsMPkRnSdUPSHM3uXcaqHfOwuXYHzG07rFhteHVMD6ewE5E14camRGZInx3QB17+Az/+OhMeGWn46Nx2dL1+rNjt4JgeIrJG7PkhMkO6TB+3Uyow49D3GHZ+p7psf5VGOBTUUO/rc0wPEVkzhh8iM/HqlPaoaJlWx5TKeIGvdy9Cuzvn1GXfvtkL81sP1Rjzow03sT3eaeiP9sFSjukhIqMoaOmOksbwQ2QG8pvSXpTy8gSs3TYbNf69BwDIsrPHZx1GY0udjnq1Yc17DdGsqpdexxIRFaWgpTtM0cvMMT9EJpYzpV2X4FMv7gZ2/DRJHXyeObth8Dtf6B18Srs4ogkHNBORkRT0ey5eno5RGy5q3dttKAw/RCak15R2QcD0Q9+j3ItnAIC7ZfzQa9BinKpUW+92vN8sgLe5iMgoCvs9l1M2a3cMFMqSW3aQ4YfIhPSa0i4SYWyPTxHv5omTFWuj16DFiPUsr3cbyrg6YmzbqnofT0RUmKJ+zwkAZPJ0nI1NLrE2ccwPkQnpM6UdAJ64e+GdgQsg8/BClr2j3tcXAZjXuxZ7fYjIaLT9Pafv70N9sOeHyIS0mdJe7vlTfLVnKdwyXmiUPyjjW6zg4ytxxqr36nM6OxEZlbZLd+iyxEdxseeHyISepmUW+nz1hFis3TYbFVL+RdkXz/Bhn5lQ6DiF/VV96pVHy2rluEUFEZWYRoGe8JU4I16enu+4n1dXky8pDD9EJqJQCvhk698FPt/mzjks37UQbpkvAQDV/r0PaWoS4iTeel3PV+KMhW/XYeChEmEu67mQ6dnbiRAeFoxRGy5CBGgEoJx3REmvJs/wQ2QiJ28l4kWmIu8TgoAPzu/CjEPrYC8oAQCXfatieO/P8a+bfn8ZicCtKqjkmNN6LmQeOoX4YtV79fO8L6Qmel+IBEEoubllFiIlJQUSiQRyuRweHh6mbg5ZqIL+8s0p/2z7FdxJ1BzH46DIxqwDq/Hu5Sh12Z5qzTC560SkO+p3P7y0qyPm967FDx0qETnrubz+wZITuznOzLYZu0dQ289v9vwQGUFBf/l2r+OLHZce40lqRp5jPNKf45sd89Hi/mV12deh/bC0xbsQRPrPTfhmQH2u3Ewloqj1XERQrefSPljKXkgbZW8nQqgZLKjK8ENkYAX95SuTp2PN0dh8j/F8IcevG6egSvIjAECGvQOmdB6PHTXb6N2OnEGEXLmZSoou67mYwwcg2S5OdScyIL1WbAbw1MUdN7wDAQCJrhIM7D+32MEH4DgfKlnmuJ4LUX7Y80NkQHqt2AxAENlhcpcJyHBwxNJmA/GotLRY7TDVIEKybea4ngtRfhh+iAwoXv5Sq3oiQQn/Z0/woExuOMlwFGNy10nFuv7g0EroHOLLacVkEua4ngtRfnjbi8hAoqJlmP17TJH1XDLTsXr7XOz4eTIqPjXsTsadQ3wRGlSWwYdMImc9FyD31msO3oolc8LwQ2QAUdEyjNxwEU9fZBVaT5qSiN82TUHHW6fh+TIF30Z+CXtlPmv96EgE1Wwy/kVNppaznotUonlrS8rtVMiM8LYXUTEplAKmRl4tsl5I/G2s3TYb0ueqnYtTnFwxr80HxdquAuBf1GR+OoX4on2wlCs8k9li+CEqptN3kvCsiB6fjjdPYtnvi+GSrVrf56HEBx/0mYlb5SoV+/oc3EzmyFzWcyHKD8MPUTGduptY8JOCgNGnf8OnR39SF50vXwMf9f4Mya4Sva858a03EODlyr+oiYj0wPBDVEy3nqTmW+6UnYW5f65A3+i/1GWRNdtgWqdxyHBw0utaZVwdMY9bVRARFQvDD5GW8tuTZmHUdfwZk5Bv/cYPr2oEn0UtBuGb0HcAkfa9NAMb+aO0qyPsRHYIDSqLJpU5k4uIqLgYfoi0kN9eXU72ImQqCl7L+VhgfSxu/i7GnP4Nk7pOwt7qzXW6pq/EGV/0rMWwQ0RkYAw/REUoaK+uwoJPjuVN+2N3jZa451lep2uKwNlbRETGwnV+iAqhy15d717ai7ev7NMsFIl0Dj6+XA+FiMio2PNDVAht9uqyUyrw2cF1+ODCLmTZ2eNhaSlOV6yt87VavVEOI1sFcfYWEZGRseeHqBBF7T7tlvECa7d9gQ8u7AIAOCoVaPTwml7XGtkqiFtTEBGVAPb8EBWisN2nK8ifYO3W2aieeB8AkGVnjxkdxuDXOh10ugY3eyQiKlkMP0QFUCgFZCuUEAF5xvzUj7uObyO/hNcLOQDgmbMbRvaartftLoCDm4mIShLDD1E+8pvanqN7zGEs2vs/iBWqLS3ueJbHsD4zdR7YDABSDzEiutfk4GYiohLE8EP0H4VSwOk7Sdhw5h7+iH6Sb51h53bg84Nr1Y9PVqyNUT2nQe7irvP1Jr71Bsa2rcIeHyKiEsbwQzZPoRSw4uBtrDl6By8yFYXWPVuhJl46iOGSnYFfanfA5x1GI9tet/+NfLkRKRGRSTH8kE2LipZhauTVIndlz3HVtyomdpukGuz8Zi+dtqoYGloJHUN8OZWdiMjEGH7IZu298hijN10qtE5AchwelJZCaWevLouq1kzna7UP9kZEjxCdjyMiIsPjOj9kk/ZekWHsL4UHnzZ3zuH3HydgxqHvi3Wt4S0C8N3gN4t1DiIiMhz2/JDNiYqWYfSmiwVXEAQMO78T0w99D3tBiWHnd+KiX3XsqdFC52v9OPRNtKruXYzWEhGRoTH8kE3J2aurIA6KbMzevxoD/45Sl/1evQUOVGmk1/WepWs3loiIiEoOww/ZlML26vJIf45VO+ai2f0r6rL/NR2AZc0HQBDpd4e4sBWiiYjINBh+yKYUtFdXQHIc1m2bjaDkOABAhr0DPu08HjtrttHrOtyygojIfDH8kFVSKAWcjU1GQmo6vN2d1dPL8+uJafzgKlZvn4sy6akAgERXCT7q9RkuVqih17VzJrFzywoiIvPE8ENWJ7+tKXIWFgyt7KVZWRDwydGf1cHnpldFDOsbjkcSH72vL+UihkREZk0kCMLrezbavJSUFEgkEsjlcnh4eJi6OaSDqGgZRm24mGcj0hz2IkDx2pPlnidj50+T8I9XJYztMQXPxa56XfuDZgFoHyzlIoZERCai7ec3e37IauTM5Coszb8efADgXzdPvP3uQsS7l4XilcUMtcXtKoiILAvDD1mNwmZy5ZCmJGLa4fX4rONopIpLqcvjJLqtxdOksif6vVkRUg9n9vQQEVkYhh+yGgXN5MpRS3YLayO/gM/zZJROT8UHfcN17unxcHbA/N610aU2e3mIiCwVww9ZjcLW1Ol08wSW/r4ELtkZAIDA5DiUe/4U8R5eBR7zOs9Sjjg97S04OXBXGCIiS8bwQ1bjaVpm3kJBwOjTv+HToz+pi85WCMbIXjOQ7CrR6rw5N7Tm9qrF4ENEZAUYfsgi5azjE5+SjuTnGSjt4og5f9zQqOOUnYV5fy5Hn+iD6rJtIW0xreM4ZDo4an0tTl0nIrIuFhV+5syZgz179uDy5ctwcnLCs2fP8tR58OABxowZg4MHD8LFxQUDBw7EV199BScnp5JvMBlFfuv4vK7MCznWbJ+DRo9y9/Fa2HIwVjZ5GxBpNzi5tKsjvhlQH02CynJAMxGRFbGo8JOZmYm3334boaGhWLduXZ7nFQoFunbtinLlyuH48eNISkrCkCFDIAgCli9fboIWk6EVtY4PAHilPcW2Df+HSs/iAQAvHcSY2G0Soqo10/o6IgDze9dCs6rajwkiIiLLYFHhZ9asWQCAH374Id/n9+3bh5iYGDx8+BB+fn4AgMWLF2Po0KGYM2dOgQseZWRkICMjQ/04JSXFsA0ng1AoBUyNvFpo8AGARNfSuOxbDZWexeOJmyc+7P05rvpW1fo6XLeHiMi6WVT4KcqpU6cQEhKiDj4A0LFjR2RkZODChQto0yb/TSrnzZunDlZkvlYcvIVnL7KKrigS4dMu4/HCUYz/NRuo04yuz7vWwNBmgbzNRURkxaxq6kp8fDx8fDT3ZCpTpgycnJwQHx9f4HHTpk2DXC5Xfz18+NDYTSUtKZQCTt1JwvaLj/Dtsbv51rFXKlA56ZFGWYaDE6Z1/lin4AMAXu5iBh8iIitn8vATEREBkUhU6Nf58+e1Pp8on8GsgiDkW55DLBbDw8ND44tMLypahuYLDmLAd6cx8de/kZahyFPHPSMN67bOVo3xefq42NcsbK0gIiKyDia/7TV27Fj079+/0DoBAQFanUsqleLMmTMaZU+fPkVWVlaeHiEyb9oMbK7wLB7rts1GtcQHAIA1kXPQ5f2vodRjfy4RVFPaGwV66tdgIiKyGCYPP15eXvDyMsyMmtDQUMyZMwcymQy+vqrBqvv27YNYLEaDBg0Mcg0yPm02KK3/6Dq+3f4lvF7IAQBPnd0R0X6EXsEnR3hYMG95ERHZAJOHH108ePAAycnJePDgARQKBS5fvgwAqFKlCtzc3NChQwcEBwdj0KBBWLRoEZKTk/HJJ59g+PDhvJVlQU7fTSp0DZ8e1w5h4R//g1iRDQC441kew/rMxD3P8npdT+ohRkT3mpzdRURkIywq/MycORM//vij+nG9evUAAIcOHULr1q1hb2+PPXv2YPTo0WjWrJnGIodkGaKiZZi67Wq+z4kEJSYc34TxJzery05Uqo1RPacjxdlNr+tNfOsNjG1bhT0+REQ2RCQIQlHLpticlJQUSCQSyOVy9hiVgJytKg7ExGPdiXv51hFnZWDx3mXoduOYumxTnU6Y2X4ksu11z/Ds7SEisj7afn5bVM8PWR9ttqoAgMYPo9XBRwkR5rQdhnUNe2i9VcWrZnSpgQ+acy0fIiJbZfKp7mS7cmZ0FRV8AOBo5QZY2HIwnju5YHifz7DuzZ56BZ/Sro4MPkRENo49P2QS2szoet3KJm9jZ3BrxEm89b7u/N61GHyIiGwce37IJM7GJhfc4yMI+PBsJPpfjtIsF4n0Dj52ImDlwHoc40NEROz5IdNISM0/+DgqsjB73yoMuLIPWXb2eFBaipMBdYt9vRUD6qNLbQYfIiJizw+ZyL3EF3nKJC9T8eOv4RhwZR8AwFGpQK0nt4t1HV+JM1a/x+BDRES52PNDJS4qWoZlB/7RKAtMjsO6rbNQ+b/9uTLsHfF/XcZjV3Brnc/v4eyAvg0qoH2wFI0CPTnGh4iINDD8UInKb6Bz6P0rWLVjLkqnPwcA/OtaGiN6z8DF8jV0Pr+HswPOf9YeTg7s1CQiovwx/FCJen2g8zt/78Ocfd/AUanasf2GVyV82HcmHkn024i2b4MKDD5ERFQohh8qMQqlgBO3/1U/HnX6N0w5krtdyaHKDTCu+xQ8F7vqfY32wdJitZGIiKwf/0SmEhEVLUPzBQex4tAdddnRgHp44SgGAHzfoDs+7DOzwODj4lj0W9VX4oxGgZ6GaTAREVkt9vyQ0eWs5Pz6gobXpFUwPuz/4P08GRvrdSn0HGuHvInz955i6WsDpQEgZzhzeFgwBzcTEVGRGH7IqF4d4Fwl8QFiPctDYWevfn5/1SZFnkPqIUaTymXRrIoXqknd8uwFJpU4IzwsmAsYEhGRVhh+yKhyBjh3uXEcS/YswebaHRDRfqRO54joXlPdo9MpxBftg6U4G5uMhNR0eLs7czo7ERHphOGHjEKhFHA2Nhl/XH2MMSe34P+O/QwAGHrxd5zxD8Ef1ZtrdZ6VA+vn6dGxtxMhNKiswdtMRES2geGHDC4qWoZZu2OQlJSK+VFfo/e1Q+rntoa0w19VGmt1ngntqnBlZiIiMjiGHzKInJ6e/THx+P7EPXi+kGNj5By8GRejrrOg1RCsatwXEGl3iyqwnJuxmktERDaM4YeKLaenJ2cQcpXEB/h+6yxUlD8BALx0EGNCt8n4s1pTnc7r7e5s8LYSEREx/FCxvD6NveXdC1ixcwE8MlUbl8a7eeLDPjMRLa2i9TlFUM3g4po9RERkDAw/pLc8+3QJAsad3KIOPtE+QRjW53M8cffS+pxcs4eIiIyNKzyT3l7fpwsiEUb3nIY493L4s2oTvD1wgU7BB1D1+Kx6L+8MLyIiIkNhzw/pLV7+Mk/Zv25l0Oe9RXji7glBpF22/rxrDXi5i7lmDxERlQj2/JDeMm/dwartc+GR/lyjPN7DS+vg4ytxxtBmgehRtzxCg8oy+BARkdGx54f0c+IEeo7sDfGzZLhlvMD7b0cg2173txPH9hARUUljzw/pbsMGoG1biJ8lAwD8UhNR5mWKTqewEwErB9bj2B4iIipxDD+kPaUS+PxzYNAgIDMTAHC2Sn30GvQV/nXTbVr6igH10aW2nzFaSUREVCiGH9LOixdAv37Al1/mlo0Ygae/bkeqs24rMU9oV5XbVhARkckw/FDRZDKgVStg61bVYzs7YNkyYNUqdKxXEaveqw9fiXarMZdxdcS4dlWN11YiIqIicMAzacjZoyshNV019dw5A3ZNGkP06JHq+VKlcHPJt7jVqBW87yajUaAnOoX4on2wVGNvr4LM612LA5yJiMikGH5I7fU9ugCgtIsD5ni+ga6PHuGRRzkM6xuOm3clwN3LAFRT1cPDgtEpxBehQWURGlQWjQI9EbErBvEpued5tR4REZEpiQRBEIquZltSUlIgkUggl8vh4eFh6uaUiNf36HqVODsT0w+tw/Km/ZFYqozGczl9OK+vypynB4mLFxIRkZFp+/nN8JMPWws/CqWA5gsOQiZPh6MiCwHJj3GrXCWtj8/ZiPT4lLYMOEREZDLafn5zwDOp9+iSvEzFT7/OxG+bpiAwOU7r4wUAMnk6zsYmG6+RREREBsLwQ0hITUdgchy2/zwZoQ+uonT6c6zcMQ8iQanzeYiIiMwdBzwTqkSfxfafJ6P0f3t0/etaGtM7jtV6f64c3u7aTXcnIiIyJYYfW7d2LYJHjYIoOxsAcMOrEob1DUecxFvrU+SM+WkUqNsqz0RERKbA2162SqEAPvkEGD5cHXwOVm6It99bpHPwAbhBKRERWQ6GH1v0/DnQqxeweHFu2YQJyIzcDjdvzd6b0q6OKO3qWOCppBLnPNPciYiIzBlve9miEyeA3btV39vbAytWACNHohOA9rXK51mfB4C6zMtNDAhAYloG1+8hIiKLxHV+8mET6/wsWADMm6far+utt0zdGiIiomLjOj9UuE8/Ba5dY/AhIiKbw/Bj7QQB+PJLYPVqzXKRCChf3jRtIiIiMiGO+bFm6enAhx8CGzeqxvZUrgx06GDqVhEREZkUe36sVUIC0LatKvgAgFIJ3Lhh2jYRERGZAfb8WKPoaCAsDLh3T/XY1RXYsEE1vZ2IiMjGMfxYmz/+APr1A1JTVY/9/FTT2uvXN227iIiIzARve1kLQQCWLwe6dcsNPvXrA2fPMvgQERG9guHHWsyaBXz8sWpsDwD07g0cPcoZXURERK9h+LEWvXoBpUqpvp82Dfjtt9zHREREpMYxP9aiTh3gl1+A5GRgyBBTt4aIiMhsMfxYqsuXgZo1AcdXNh0NCzNZc4iIiCwFb3tZop9+Aho3BsaPVw10JiIiIq0x/FgSpRKYMUN1WyszE1i1CoiMNHWriIiILApve1mKFy+AwYOBbdtyy0aOBLp3N12biIiILBDDjyV4/Bjo0QM4f1712M4OWLoUGDdOtUEpERERaY3hx9xduqQayBwXp3rs7g5s3gx06WLadhEREVkojvkxZzt2AM2b5wafSpWAEycYfIiIiIqB4cdcKZXAV1+pxvoAQJMmwJkzQK1apm0XERGRhWP4MVd2dqrBzZUqAQMGAIcOAT4+pm4VERGRxeOYH3Pm4wOcPq36Lwc2ExERGQR7fszFzZtA166q7SleJZUy+BARERkQw485+Osv1ZievXuBt98GsrJM3SIiIiKrxfBjat99B3TqBDx7pnqcmJi394eIiIgMhuHHVBQKYPJk4KOPgOxsVVm3bsDx4xzYTEREZEQMP6aQmgr07AksWZJbNmmSal0fd3dTtYqIiMgmcLZXCVEoBZyNTUbqP3fQbNIHKHXjmuoJBwdg5Upg+HDTNpCIiMhGMPyUgKhoGWbtjoHy0SP8/uMElEp7BgDIcpfAcUck0LataRtIRERkQyzmtte9e/cwbNgwBAYGwsXFBUFBQQgPD0dmZqZGvQcPHiAsLAylSpWCl5cXPv744zx1SlJUtAyjNlyETJ6OJ25lcTSwPgAgtowvOvVbgCjvGiZrGxERkS2ymJ6fGzduQKlUYs2aNahSpQqio6MxfPhwpKWl4auvvgIAKBQKdO3aFeXKlcPx48eRlJSEIUOGQBAELF++vMTbrFAKmLU7BkJOgUiEaR3HIclFgpWhb0Pu4oFZu2PQPlgKezuu5UNERFQSRIIgCEVXM0+LFi3CqlWrcPfuXQDAH3/8gW7duuHhw4fw8/MDAGzevBlDhw5FQkICPDw8tDpvSkoKJBIJ5HK51sfk53RMHGYvikSMT+VC6/0yvAlCg8rqfR0iIiLS/vPbYm575Ucul8PT01P9+NSpUwgJCVEHHwDo2LEjMjIycOHChQLPk5GRgZSUFI2vYnvyBG8M6IHNv0xDUNLDQqsmpKYX/3pERESkFYsNP3fu3MHy5csxcuRIdVl8fDx8Xlsjp0yZMnByckJ8fHyB55o3bx4kEon6y9/fv3iNi44GGjeG55UL8MhIw4qdCyASlAVW93Z3Lt71iIiISGsmDz8REREQiUSFfp0/f17jmMePH6NTp054++238eGHH2o8J8pnHyxBEPItzzFt2jTI5XL118OHhffUFGrvXqBpU+D+fQBAgocXPuk6EYIo70stAuArcUajQM88zxEREZFxmHzA89ixY9G/f/9C6wQEBKi/f/z4Mdq0aYPQ0FB8++23GvWkUinOnDmjUfb06VNkZWXl6RF6lVgshlgs1r3xrxIEYPlyYOJEQPlfL0+DBri28DvE7HsMEYBXB1flRLHwsGAOdiYiIipBJg8/Xl5e8PLy0qpuXFwc2rRpgwYNGmD9+vWws9PsTQkNDcWcOXMgk8ng6+sLANi3bx/EYjEaNGhg8LarZWUB48cDq1bllvXtC/z4I9q4umKVtxSzdsdAJs8d2yOVOCM8LBidQnyN1y4iIiLKw2Jmez1+/BitWrVCxYoV8dNPP8He3l79nFQqBaCa6l63bl34+Phg0aJFSE5OxtChQ9GzZ0+dprrrNNvr2TPgnXeA/ftzy2bMAGbPBl4JZzkrPCekpsPbXXWriz0+REREhqPt57fJe360tW/fPty+fRu3b99GhQoVNJ7LyW/29vbYs2cPRo8ejWbNmsHFxQUDBw5UrwNkFGfOAAcOqL53cgLWrgUGDcpTzd5OxOnsREREZsBien5Kks7r/Hz1FbBgAbB9O9C8ufEbSERERHnYxDo/ZmPyZODaNQYfIiIiC8DwowulEpg+HVixQrNcJAK8vU3TJiIiItKJxYz5Mbm0NGDwYCAyUjWQuUoVoFMnU7eKiIiIdMTwo424OKB7d+Dixdyy/xYxJCIiIsvC8FOUixeBsDDg8WPVY3d3YMsWoHNn07aLiIiI9MIxP4XZvRto0SI3+AQEACdPMvgQERFZMPb8FOa993K/b9pUNZWdA5uJiIgsGnt+tPHuu8BffzH4EBERWQH2/OQjZ93HFEA1tf3TT4HMTNUXERERmaWUlBQAuZ/jBeEKz/l49OgR/P39Td0MIiIi0sPDhw/zbIX1KoaffCiVSjx+/Bju7u5ITU2Fv78/Hj58qN1WF1YsJSWFr8V/+Fpo4uuRi69FLr4Wufha5DLmayEIAlJTU+Hn5wc7u4JH9vC2Vz7s7OzUiVEkUu287uHhYfNv2Bx8LXLxtdDE1yMXX4tcfC1y8bXIZazXQiKRFFmHA56JiIjIpjD8EBERkU1h+CmCWCxGeHg4xGKxqZticnwtcvG10MTXIxdfi1x8LXLxtchlDq8FBzwTERGRTWHPDxEREdkUhh8iIiKyKQw/REREZFMYfoiIiMimMPy85vDhwxCJRPl+nTt3rsDjhg4dmqd+kyZNSrDlxhEQEJDn55o6dWqhxwiCgIiICPj5+cHFxQWtW7fGtWvXSqjFxnHv3j0MGzYMgYGBcHFxQVBQEMLDw5FZxH5v1vK+WLlyJQIDA+Hs7IwGDRrg2LFjhdY/cuQIGjRoAGdnZ1SuXBmrV68uoZYa17x58/Dmm2/C3d0d3t7e6NmzJ27evFnoMQX9Trlx40YJtdo4IiIi8vxMUqm00GOs9X2R3+9JkUiEMWPG5Fvfmt4TR48eRVhYGPz8/CASibBjxw6N5/X9PNi2bRuCg4MhFosRHByM7du3G7TdDD+vadq0KWQymcbXhx9+iICAADRs2LDQYzt16qRx3N69e0uo1cY1e/ZsjZ/rs88+K7T+woULsWTJEqxYsQLnzp2DVCpF+/btkZqaWkItNrwbN25AqVRizZo1uHbtGpYuXYrVq1dj+vTpRR5r6e+LLVu2YMKECZgxYwYuXbqEFi1aoHPnznjw4EG+9WNjY9GlSxe0aNECly5dwvTp0/Hxxx9j27ZtJdxywzty5AjGjBmD06dPY//+/cjOzkaHDh2QlpZW5LE3b97UeB9UrVq1BFpsXDVr1tT4ma5evVpgXWt+X5w7d07jddi/fz8A4O233y70OGt4T6SlpaFOnTpYsWJFvs/r83lw6tQp9OvXD4MGDcLff/+NQYMG4Z133sGZM2cM13CBCpWZmSl4e3sLs2fPLrTekCFDhB49epRMo0pQpUqVhKVLl2pdX6lUClKpVJg/f766LD09XZBIJMLq1auN0ELTWbhwoRAYGFhoHWt4XzRq1EgYOXKkRln16tWFqVOn5lv/008/FapXr65RNmLECKFJkyZGa6OpJCQkCACEI0eOFFjn0KFDAgDh6dOnJdewEhAeHi7UqVNH6/q29L4YP368EBQUJCiVynyft9b3BABh+/bt6sf6fh688847QqdOnTTKOnbsKPTv399gbWXPTxF27dqFxMREDB06tMi6hw8fhre3N9544w0MHz4cCQkJxm9gCViwYAHKli2LunXrYs6cOYXe6omNjUV8fDw6dOigLhOLxWjVqhVOnjxZEs0tMXK5HJ6enkXWs+T3RWZmJi5cuKDx7wkAHTp0KPDf89SpU3nqd+zYEefPn0dWVpbR2moKcrkcALR6H9SrVw++vr5o164dDh06ZOymlYhbt27Bz88PgYGB6N+/P+7evVtgXVt5X2RmZmLDhg344IMP1HtDFsQa3xOv0vfzoKD3iiE/Qxh+irBu3Tp07NgR/v7+hdbr3LkzNm7ciIMHD2Lx4sU4d+4c2rZti4yMjBJqqXGMHz8emzdvxqFDhzB27FgsW7YMo0ePLrB+fHw8AMDHx0ej3MfHR/2cNbhz5w6WL1+OkSNHFlrP0t8XiYmJUCgUOv17xsfH51s/OzsbiYmJRmtrSRMEAZMmTULz5s0REhJSYD1fX198++232LZtGyIjI1GtWjW0a9cOR48eLcHWGl7jxo3x008/4c8//8R3332H+Ph4NG3aFElJSfnWt5X3xY4dO/Ds2bNC/2C21vfE6/T9PCjovWLQzxCD9SGZufDwcAFAoV/nzp3TOObhw4eCnZ2dsHXrVp2v9/jxY8HR0VHYtm2boX4Eg9HntcixdetWAYCQmJiY7/MnTpwQAAiPHz/WKP/www+Fjh07GvxnKS59Xou4uDihSpUqwrBhw3S+njm/L/ITFxcnABBOnjypUf7ll18K1apVy/eYqlWrCnPnztUoO378uABAkMlkRmtrSRs9erRQqVIl4eHDhzof261bNyEsLMwIrTKd58+fCz4+PsLixYvzfd5W3hcdOnQQunXrpvNx1vCewGu3vfT9PHB0dBQ2bdqkUbZhwwZBLBYbrK0OhotR5m3s2LHo379/oXUCAgI0Hq9fvx5ly5ZF9+7ddb6er68vKlWqhFu3bul8rLHp81rkyJmpdPv2bZQtWzbP8zmzPeLj4+Hr66suT0hIyJPkzYGur8Xjx4/Rpk0bhIaG4ttvv9X5eub8vsiPl5cX7O3t8/zFVdi/p1Qqzbe+g4NDvu8ZSzRu3Djs2rULR48eRYUKFXQ+vkmTJtiwYYMRWmY6pUqVQq1atQp8b9vC++L+/fs4cOAAIiMjdT7WGt8T+n4eFPReMeRniM2EHy8vL3h5eWldXxAErF+/HoMHD4ajo6PO10tKSsLDhw81/sHNha6vxasuXboEAAX+XIGBgZBKpdi/fz/q1asHQHUP/MiRI1iwYIF+DTYiXV6LuLg4tGnTBg0aNMD69ethZ6f7XWNzfl/kx8nJCQ0aNMD+/fvRq1cvdfn+/fvRo0ePfI8JDQ3F7t27Ncr27duHhg0b6vX/kjkRBAHjxo3D9u3bcfjwYQQGBup1nkuXLlnMe0BbGRkZuH79Olq0aJHv89b8vsixfv16eHt7o2vXrjofa43vCX0/D0JDQ7F//35MnDhRXbZv3z40bdrUcI0zWB+SlTlw4IAAQIiJicn3+WrVqgmRkZGCIAhCamqqMHnyZOHkyZNCbGyscOjQISE0NFQoX768kJKSUpLNNqiTJ08KS5YsES5duiTcvXtX2LJli+Dn5yd0795do96rr4UgCML8+fMFiUQiREZGClevXhUGDBgg+Pr6WvRrkXOrq23btsKjR48EmUym/nqVNb4vNm/eLDg6Ogrr1q0TYmJihAkTJgilSpUS7t27JwiCIEydOlUYNGiQuv7du3cFV1dXYeLEiUJMTIywbt06wdHRUa/bx+Zm1KhRgkQiEQ4fPqzxHnjx4oW6zuuvx9KlS4Xt27cL//zzjxAdHS1MnTpVAGAxtz4LMnnyZOHw4cPC3bt3hdOnTwvdunUT3N3dbfJ9IQiCoFAohIoVKwpTpkzJ85w1vydSU1OFS5cuCZcuXRIAqD8z7t+/LwiCdp8HgwYN0pg9euLECcHe3l6YP3++cP36dWH+/PmCg4ODcPr0aYO1m+GnAAMGDBCaNm1a4PMAhPXr1wuCIAgvXrwQOnToIJQrV05wdHQUKlasKAwZMkR48OBBCbXWOC5cuCA0btxYkEgkgrOzs1CtWjUhPDxcSEtL06j36mshCKrpjeHh4YJUKhXEYrHQsmVL4erVqyXcesNav359gWOCXmWt74tvvvlGqFSpkuDk5CTUr19fY2r3kCFDhFatWmnUP3z4sFCvXj3ByclJCAgIEFatWlXCLTaOgt4Dr77/X389FixYIAQFBQnOzs5CmTJlhObNmwt79uwp+cYbWL9+/QRfX1/B0dFR8PPzE3r37i1cu3ZN/bwtvS8EQRD+/PNPAYBw8+bNPM9Z83siZ9r+619DhgwRBEG7z4NWrVqp6+f47bffhGrVqgmOjo5C9erVDR4MRYIgCIbrRyIiIiIyb5zqTkRERDaF4YeIiIhsCsMPERER2RSGHyIiIrIpDD9ERERkUxh+iIiIyKYw/BAREZFNYfghIiIim8LwQ0RERDaF4YeIiIhsCsMPERER2RSGHyKyer/88gucnZ0RFxenLvvwww9Ru3ZtyOVyE7aMiEyBG5sSkdUTBAF169ZFixYtsGLFCsyaNQtr167F6dOnUb58eVM3j4hKmIOpG0BEZGwikQhz5sxB37594efnh//97384duwYgw+RjWLPDxHZjPr16+PatWvYt28fWrVqZermEJGJcMwPEdmEP//8Ezdu3IBCoYCPj4+pm0NEJsSeHyKyehcvXkTr1q3xzTffYPPmzXB1dcVvv/1m6mYRkYlwzA8RWbV79+6ha9eumDp1KgYNGoTg4GC8+eabuHDhAho0aGDq5hGRCbDnh4isVnJyMpo1a4aWLVtizZo16vIePXogIyMDUVFRJmwdEZkKww8RERHZFA54JiIiIpvC8ENEREQ2heGHiIiIbArDDxEREdkUhh8iIiKyKQw/REREZFMYfoiIiMimMPwQERGRTWH4ISIiIpvC8ENEREQ2heGHiIiIbMr/A2pG2R5hrbkFAAAAAElFTkSuQmCC",
      "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": 6,
   "id": "c6a1998f-81eb-4a32-8d00-04639ca777c2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.027163347597745746"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ec33c149-7d6f-433b-872a-ca97c0b7d7d6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.8940177705732415"
      ]
     },
     "execution_count": 7,
     "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.12.3"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}