{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <span style=\"color:#F72585\"><center>Modelos Lineales Generalizados</center></span>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<figure>\n",
    "<center>\n",
    "<img src=\"../Imagenes/moon.jpg\" width=\"600\" height=\"600\" align=\"center\" /> \n",
    "</center>   \n",
    "</figure>\n",
    "\n",
    "<a href=\"https://commons.wikimedia.org/wiki/File:Moon_in_Blue_Sky.jpg\">Keshani Kaveesha</a>, <a href=\"https://creativecommons.org/licenses/by-sa/4.0\">CC BY-SA 4.0</a>, via Wikimedia Commons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color:#4361EE\">Referencias</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. [Dobson y Baernett, An introduction to Generalizer Linear Modesl](http://library.lol/main/472B57FA461867F6CFB4334BFED60010)\n",
    "1. [Nelder y Wedderburn, Generalized  Linear Models](http://www.medicine.mcgill.ca/epidemiology/hanley/bios601/Likelihood/NelderWedderburn1972.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color:#4361EE\">Introducción</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El proceso de modelación estadística clásicamente incluye los siguientes pasos:\n",
    "\n",
    "1. Especificar los modelos en dos partes: ecuaciones que relacionan la respuesta con variables exploratorias y la distribución de probabilidad de la   variable respuesta.\n",
    "1. Estimar los parámetros usados en los modelos.\n",
    "1. Chequear que tan bien los modelos ajustan a los datos reales.\n",
    "1. Hacer inferencias; Por ejemplo calculando intervalos de confianza o de credibilidad, pruebas de hipótesis acerca de los parámetros y de los modelos.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color:#4361EE\">Familia exponencial de distribuciones</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La familia exponencial de distribuciones es un conjunto de distribuciones muy especial. Muchas propiedades que poseen la hacen muy atrayente para  la modelación  estadística. Algunas de las distribuciones más notables de la probabilidad pertenecen a ella, lo que facilita el desarrollo de modelos de regresión relativamente sencillos.\n",
    "\n",
    "Consideremos una variable aleatoria $Y$ cuya función de probabilidad depende de un parámetro $\\theta$. La distribución pertenece a la familia exponencial si y solo sí puede escribirse en la forma\n",
    "\n",
    "$$\n",
    "f(y|\\theta) =  \\exp\\left[a(y)b(\\theta) + c(\\theta) + d(y) \\right],\n",
    "$$\n",
    "\n",
    "en  donde $a, b, c, d$ son funciones conocidas.\n",
    "\n",
    "* Si $a(y) = y$, se dice que la expresión de la densidad está en `forma canónica`.\n",
    "* $b(\\theta)$ se llama parámetro canónico.\n",
    "* Si la distribución tiene otros parámetros, estos se denominan parámetros de ruido.\n",
    "\n",
    "\n",
    "Los siguientes son los tres ejemplos más notables de distribuciones de la familia exponencial.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:#4CC9F0\">Distribución normal</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Arriba hemos escrito la forma como comúnmente se conoce esta distribución. Se deja al lector verificar que esta densidad se puede reescribir en la forma. \n",
    "\n",
    "$$\n",
    "f(y|\\theta, \\sigma^2) = \\exp \\left[-\\frac{y^2}{2\\sigma^2}  + \\frac{y\\theta}{\\sigma^2} - \\frac{\\theta^2}{2\\sigma^2} -  \\frac{1}{2}\\log (2\\pi \\sigma^2) \\right].\n",
    "$$\n",
    "\n",
    "En este caso, el parámetro de interés de $\\theta$ y $\\sigma^2$ es un parámetro de ruido. Se tiene que:\n",
    "\n",
    "* La densidad está en su forma canónica, pues si  $a(y) = y$,\n",
    "* El parámetro natural es $b(\\theta) = \\frac{\\theta}{\\sigma^2}$,\n",
    "* $c(\\theta) = - \\frac{\\theta^2}{2\\sigma^2} - \\frac{1}{2}\\log (2\\pi \\sigma^2) $, \n",
    "* $d(y) =  -\\frac{y^2}{2\\sigma^2}$.\n",
    "\n",
    "Se puede verificar que si una variable aleatoria $Y$ tiene distribución $\\text{N}(\\theta, \\sigma^2)$, entonces\n",
    "\n",
    "* $E[Y] = \\theta$,\n",
    "* $Var[Y] = \\sigma^2$.\n",
    "\n",
    "El siguiente fragmento de \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGyCAYAAADH859HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABftElEQVR4nO3de1xUdf4/8NdwmwGVEUS5BCKSCt4RFMHwLt5y1fol1YrZauVmrWhtZZql237Jti231utWslYqlZnuLqaYlhqUSuD9lje8gAgGgzeu5/fHp7kxgAwMnLm8no/HefCZM2fOvA8DM+/5XBWSJEkgIiIismJOcgdAREREdC9MWIiIiMjqMWEhIiIiq8eEhYiIiKweExYiIiKyekxYiIiIyOoxYSEiIiKrx4SFiIiIrB4TFiIiIrJ6TFhINg8//DB69eqFkpISuUMhohZ28eJFeHt7Y9GiRXKHQjaCCYsNS0lJgUKhqHV78cUXZY1t+vTp6NSpU533L1u2DPv378e2bdugVqubPZ4LFy5AoVAgJSXFLp6nIRQKBd544417Hqf9O7pw4YLFnvter/+9dOrUCdOnT7dYPEOHDsXQoUMtdr661HbdDX0dDKWlpZn9mNqeS/vaHjx40OxzNVV9111eXo4pU6Zg4sSJWLJkiUWfV3vNKpUKFy9eNLl/6NCh6Nmzp8n+4uJi+Pj4YOPGjUb7CwoKMH36dPj4+MDDwwMxMTH49ttvGxTL5cuXkZSUhCFDhqBt27Z1vjdUVFQgNDQUy5Yta9B5HZWL3AFQ061duxZhYWFG+wICAmSKRnjttdcwZ86cWu/78ccf8de//hW7d+9GYGBgC0fmODIzM/n7tQKNeR3S0tKwfPlys5MWa3rN64vlhRdegJeXF/71r3812/OXlZVh4cKF+OSTTxp0/OLFixEQEICEhASjc4wYMQLFxcX4xz/+gQ4dOmD58uUYM2YMdu7ciSFDhtR7zl9++QWfffYZ+vbti3HjxmHDhg21Hufq6opFixZh7ty5SExMRLt27Rp+oQ6ECYsd6NmzJ6KiouQOw0hoaGid9w0cOBDXr19vwWgc08CBA+UOgdD8r4MkSbh79y7c3d2t6jWvL5YPPvig2Z9/zJgxWL9+PV588UX06dOn3mNv3LiB1atX47333oNCodDt/+ijj3D06FFkZGQgJiYGADBs2DD06dMHL730En766ad6zzt48GDde93BgwfrTFgA4LHHHsO8efOwevVqvPrqqw29TIfCJiE7V1e1bM0qd2016u7du/HHP/4RPj4+aNeuHR566CFcvXrV5PHr169HTEwMWrdujdatW6Nv37746KOPdPfXVjV+9+5dzJ8/HyEhIXBzc8N9992H2bNno7i42CS2Bx98EN988w369esHd3d3hIWF4eOPP27QNV+9ehVTpkxBmzZtoFarkZCQgPz8/FqPPXjwIH73u9/B29sbKpUKERER+Pzzz2V5HnNeg127dmHo0KFo164d3N3d0bFjRzz88MO4ffu27pjaXvsff/wRgwYNgkqlQkBAAObPn4+KigqTeFNTUxEfHw9/f3+4u7sjPDwcr7zyCm7dumVybEpKCrp16walUonw8HCsW7euQb8/QFSFv/TSS/Dz84OHhwceeOAB7N+/v9Zj8/Pz8cwzzyAwMBBubm4ICQnB4sWLUVlZ2eDnM7R48WJER0fD29sbnp6e6NevHz766CM0dAH7hl53zdfh9u3bePHFFxESEgKVSgVvb29ERUXpPsymT5+O5cuX6x6r3bRNdgqFAs899xxWrVqF8PBwKJVK/Pvf/671ubR+/fVXPPnkk/D29karVq0wYcIEnDt3zuiYuprhamtOKy4uxgsvvIDOnTtDqVSiQ4cOGDduHE6ePFnndQPA0aNHMXHiRHh5eUGlUqFv37662LW+++47KBQKbNiwAQsWLEBAQAA8PT0xcuRInDp1qtbfcW1eeukltGvXDi+//PI9j01JSUFlZaVR7QoAbN68Gd26ddMlKwDg4uKCqVOnYv/+/bhy5Uq953VyavhHrJubGxISErBmzZoG/w06Gtaw2IGqqiqTN20Xl8a9tDNnzsT48eOxfv16XLp0CX/+858xdepU7Nq1S3fMokWL8Je//AUPPfQQXnjhBajVahw9erTW9mItSZIwadIkfPvtt5g/fz7i4uJw+PBhvP7668jMzERmZiaUSqXu+EOHDuGFF17AK6+8Al9fX3z44YeYMWMG7r//fgwePLjO57lz5w5GjhyJq1evIjk5GV27dsX//vc/kzciANi9ezfGjBmD6OhorFq1Cmq1Ghs3bkRCQgJu375dbx+K5nyee70GFy5cwPjx4xEXF4ePP/4Ybdu2xZUrV/DNN9+gvLwcHh4etcZ8/PhxjBgxAp06dUJKSgo8PDywYsUKrF+/3uTYM2fOYNy4cUhKSkKrVq1w8uRJLF26FPv37zf6W0hJScGTTz6JiRMn4u9//ztKSkrwxhtvoKysrEFv1k899RTWrVuHF198EaNGjcLRo0fx0EMPobS01Oi4/Px8DBgwAE5OTli0aBFCQ0ORmZmJN998ExcuXMDatWvv+Vw1XbhwAc888ww6duwIQCRzzz//PK5cuXLPjqBNue558+bhk08+wZtvvomIiAjcunULR48eRVFREQDRnHrr1i18+eWXyMzM1D3O399fV/7666+xd+9eLFq0CH5+fujQoUO9zzljxgyMGjVK9ze1cOFCDB06FIcPH0bbtm3rfWxNpaWleOCBB3DhwgW8/PLLiI6Oxs2bN7Fnzx7k5eWZNE9rnTp1CrGxsejQoQPef/99tGvXDp9++immT5+Oa9eu4aWXXjI6/tVXX8WgQYPw4YcfQqPR4OWXX8aECRNw4sQJODs73zPONm3aYOHChZgzZw527dqF4cOH13ns//73P0RERJj8Lo4ePYq4uDiT43v37g0AOHbsGO677757xtJQQ4cOxcqVK3H06FH06tXLYue1GxLZrLVr10oAat0qKiokSZIkANLrr79u8tjg4GDpiSeeMDnXs88+a3Tc22+/LQGQ8vLyJEmSpHPnzknOzs7S73//+3pje+KJJ6Tg4GDd7W+++UYCIL399ttGx6WmpkoApDVr1hjFplKppIsXL+r23blzR/L29paeeeaZep935cqVEgBpy5YtRvufeuopCYC0du1a3b6wsDApIiJC97vSevDBByV/f3+pqqqqRZ+noa/Bl19+KQGQcnJy6v1d1HztExISJHd3dyk/P1+3r7KyUgoLC5MASOfPn6/1PNXV1VJFRYX0/fffSwCkQ4cOSZIkSVVVVVJAQIDUr18/qbq6Wnf8hQsXJFdXV6PXvzYnTpyQAEhz58412v/ZZ59JAIz+Pp955hmpdevWRn8TkiRJ77zzjgRAOnbsWL3PNWTIEGnIkCF13l9VVSVVVFRIS5Yskdq1a2d0PbUda85113wdevbsKU2aNKneeGfPni3V9fYMQFKr1dKNGzdqvc/wubR/U5MnTzY67ocffpAASG+++aZuX833BK2av7slS5ZIAKT09PR6r6FmLI8++qikVCql3Nxco+PGjh0reXh4SMXFxZIkSdLu3bslANK4ceOMjvv8888lAFJmZma9z6u95gMHDkhlZWVS586dpaioKN1rNWTIEKlHjx5Gj/Hw8JBmzZplci5XV9da33MyMjIkANL69evrjcXQgQMHTN4bajpz5owEQFq5cmWDz+tI2CRkB9atW4cDBw4YbY2tYfnd735ndFv7TUJbe5Keno6qqirMnj3brPNqv5XXrE145JFH0KpVK5Ne93379tV98wUAlUqFrl271luLA4jajDZt2phcx+OPP250+5dffsHJkyfx+9//HgBQWVmp28aNG4e8vLx6q5+b83nu9Rr07dsXbm5uePrpp/Hvf//bpGq/vphHjBgBX19f3T5nZ+daa4XOnTuHxx9/HH5+fnB2doarq6uug+GJEycAiG/MV69exeOPP27U7h8cHIzY2NgGxQNA97vRmjJlisnf73//+18MGzYMAQEBRr/DsWPHAgC+//77hvwKjOzatQsjR46EWq3WXeOiRYtQVFSEgoKCOh/X1OseMGAAtm3bhldeeQXfffcd7ty5Y3bsw4cPh5eXV4OPr/k7jo2NRXBwsO41MMe2bdvQtWtXjBw50qzH7dq1CyNGjEBQUJDR/unTp+P27dtGtUnAvf8PGsLNzQ1vvvkmDh48WGdTb3FxMW7fvl1nLZXha2zOfY2hjeFeTU2OigmLHQgPD0dUVJTR1lg1e6drm2m0b6raDmTmjkQoKiqCi4sL2rdvb7RfoVDAz89PVx1eVxzaWO715l5UVGT0gazl5+dndPvatWsAgBdffBGurq5G27PPPgsAKCwslOV57vUahIaGYufOnejQoQNmz56N0NBQhIaG4h//+Eed8WpjrhlfbTHfvHkTcXFx+Omnn/Dmm2/iu+++w4EDB/DVV18ZxaF9zRpyzrriqe1YFxcXk9/BtWvX8J///Mfkd9ijRw8A9b9Wtdm/fz/i4+MBAP/617/www8/4MCBA1iwYAEA1Pt31tTrfv/99/Hyyy/j66+/xrBhw+Dt7Y1JkybhzJkzDY7fsHmoIeqKteb/XUNcv369USORioqKao1bO6LxXu8BNf8PGurRRx9Fv379sGDBglr7a2nPp1KpTO5r165drb+jGzduAAC8vb3NiuVetDE0Jol1BOzDYueUSiXKyspM9jfmjQqALuG4fPmyyTel+rRr1w6VlZW4fv26UdIiSRLy8/PRv3//RsVT2/PU1mmzZmdYHx8fAMD8+fPx0EMP1Xqubt26yf48dYmLi0NcXByqqqpw8OBBfPDBB0hKSoKvry8effTROmOurVNwzX27du3C1atX8d133xkN26zZOVr7gdKQc9YVj/ZYw34AlZWVJn+fPj4+6N27N/7617/Wei5zh/Fv3LgRrq6u+O9//2v0QfX111+bFXdNDbnuVq1aYfHixVi8eDGuXbumq22ZMGGCUafV+pj7zb6uWO+//37dbZVKVet7RWFhoe7vGBDvAZcvXzbr+QHxe8vLyzPZr+1QbvgclqRQKLB06VKMGjUKa9asqTUuQJ+EGOrVqxeOHDlisl+7r7b5XJpCG0Nz/S5sHWtY7FynTp1w+PBho327du3CzZs3G3W++Ph4ODs7Y+XKlWY9bsSIEQCATz/91Gj/pk2bcOvWLd39TTVs2DCUlpZi69atRvtrdizt1q0bunTpgkOHDpnUTmm3Nm3ayP489+Ls7Izo6GjdqJKff/653pi//fZbXa0PIDpsp6amGh2n/TA07AQNAKtXrza5Nn9/f2zYsMFoVMPFixeRkZFxz9i1I08+++wzo/2ff/65SSfyBx98EEePHkVoaGitv0NzExaFQgEXFxejzpt37txp0JwdTb1uQ76+vpg+fToee+wxnDp1SjfKq7G1CXWp+TvOyMjAxYsXjUb/1PZecfr0aZMmy7Fjx+L06dNGna8bYsSIEbpk2NC6devg4eHRrEOyR44ciVGjRmHJkiUm731ubm7o3Lkzzp49a/K4yZMn4+TJk0bDlysrK/Hpp58iOjra4vNdaZt3u3fvbtHz2gvWsNi5xMREvPbaa1i0aBGGDBmC48eP45///GejZ5ft1KkTXn31VfzlL3/BnTt38Nhjj0GtVuP48eMoLCzE4sWLa33cqFGjMHr0aLz88svQaDQYNGiQbpRQREQEEhMTm3KZOtOmTcN7772HadOm4a9//Su6dOmCtLQ0bN++3eTY1atXY+zYsRg9ejSmT5+O++67Dzdu3MCJEyfw888/44svvpD9eWqzatUq7Nq1C+PHj0fHjh1x9+5d3ZDv+voVLFy4EFu3bsXw4cOxaNEieHh4YPny5SZDlWNjY+Hl5YVZs2bh9ddfh6urKz777DMcOnTI6DgnJyf85S9/wcyZMzF58mQ89dRTKC4uxhtvvNGgppHw8HBMnToVy5Ytg6urK0aOHImjR4/inXfegaenp9GxS5YsQXp6OmJjY/GnP/0J3bp1w927d3HhwgWkpaVh1apVZjVTjB8/Hu+++y4ef/xxPP300ygqKsI777xjkqTVpqnXHR0djQcffBC9e/eGl5cXTpw4gU8++QQxMTG6EV7aESJLly7F2LFj4ezsjN69e8PNza3B12jo4MGDmDlzJh555BFcunQJCxYswH333adrlgTEe8XUqVPx7LPP4uGHH8bFixfx9ttvmzTjJiUlITU1FRMnTsQrr7yCAQMG4M6dO/j+++/x4IMPYtiwYbXG8Prrr+v6Ii1atAje3t747LPP8L///Q9vv/12s894vXTpUkRGRqKgoEDXlKg1dOhQbNu2zeQxf/jDH7B8+XI88sgjeOutt9ChQwesWLECp06dws6dO42OfeONN7B48WLs3r3bKBH88ssvAeiTkYMHD6J169YAgP/3//6f0Tl+/PFHODs71zsS0qHJ3OmXmsCwN3xdysrKpJdeekkKCgqS3N3dpSFDhkg5OTl1jhKqeS5tj/3du3cb7V+3bp3Uv39/SaVSSa1bt5YiIiKMer/XHCUkSWKkz8svvywFBwdLrq6ukr+/v/THP/5R+vXXX42OCw4OlsaPH29yLfca6aF1+fJl6eGHH5Zat24ttWnTRnr44Yd1vfpr9tA/dOiQNGXKFKlDhw6Sq6ur5OfnJw0fPlxatWpViz9PQ1+DzMxMafLkyVJwcLCkVCqldu3aSUOGDJG2bt1q9DjUMkLshx9+kAYOHCgplUrJz89P+vOf/yytWbPGZJRQRkaGFBMTI3l4eEjt27eXZs6cKf3888+1XtuHH34odenSRXJzc5O6du0qffzxx7W+/rUpKyuTXnjhBalDhw6SSqWSBg4cKGVmZtY6YuX69evSn/70JykkJERydXWVvL29pcjISGnBggXSzZs3632e2v52Pv74Y6lbt26SUqmUOnfuLCUnJ0sfffRRvSOmGnPdNV+HV155RYqKipK8vLx0zz137lypsLDQ6Pcyc+ZMqX379pJCoTCKCYA0e/bsWmOq+Vzav6kdO3ZIiYmJUtu2bSV3d3dp3Lhx0pkzZ4weW11dLb399ttS586dJZVKJUVFRUm7du2q9Xf366+/SnPmzJE6duwoubq6Sh06dJDGjx8vnTx5ss5YJEmSjhw5Ik2YMEFSq9WSm5ub1KdPH5O/J+3f+xdffGG0//z58/ccZWN4zbW9Lz7++OMSAJNRQt9++60EQNq/f7/JY/Lz86Vp06ZJ3t7eur/R2kZIvfDCC5JCoZBOnDhhtB91jOSs7eM3Li5OmjBhQr3X58gUksQZaoiIyLH17t0bgwYNMru5W2vAgAEIDg42u8ZU6+zZs+jSpQu2b9+OUaNGNeoc9o4JCxERObxvvvkGkydPxpkzZ8weBaXRaNC+fXvk5OQgPDy8Uc//5JNP4vLly0hPT2/U4x0BO90SEZHDGzNmDP72t7/h/PnzZj/W09MTZWVljU5WKisrERoaqus8T7VjDQsRERFZPdawEBERkdVjwkJERERWjwkLERERWT27mTiuuroaV69eRZs2bSy+IBURERE1D0mSUFpaioCAADg51V2PYjcJy9WrV81a24aIiIisx6VLl+odUm43CYt2PZZLly6ZTOtNRERE1kmj0SAoKOie66rZTcKibQby9PRkwkJERGRj7tWdg51uiYiIyOoxYSEiIiKrx4SFiIiIrB4TFiIiIrJ6TFiIiIjI6jFhISIiIqvHhIWIiIisHhMWIiIisnqNSlhWrFiBkJAQqFQqREZGYu/evQ163A8//AAXFxf07dvX5L5Nmzahe/fuUCqV6N69OzZv3tyY0IiIiMgOmZ2wpKamIikpCQsWLEB2djbi4uIwduxY5Obm1vu4kpISTJs2DSNGjDC5LzMzEwkJCUhMTMShQ4eQmJiIKVOm4KeffjI3PCIiIrJDCkmSJHMeEB0djX79+mHlypW6feHh4Zg0aRKSk5PrfNyjjz6KLl26wNnZGV9//TVycnJ09yUkJECj0WDbtm26fWPGjIGXlxc2bNjQoLg0Gg3UajVKSko4NT8REZGNaOjnt1k1LOXl5cjKykJ8fLzR/vj4eGRkZNT5uLVr1+Ls2bN4/fXXa70/MzPT5JyjR4+u95xlZWXQaDRGGxEREdknsxY/LCwsRFVVFXx9fY32+/r6Ij8/v9bHnDlzBq+88gr27t0LF5fany4/P9+scwJAcnIyFi9ebE74RGQjCgqAb74Bjh0DqqqM72vfHhg+HIiMBJw4bIDIYTRqteaaKypKklTrKotVVVV4/PHHsXjxYnTt2tUi59SaP38+5s2bp7utXZ6aiGzXlSvAtm3AwYNAXY3VpaXAuXPA1q3AmDFAdDRQx3chIrIjZv2b+/j4wNnZ2aTmo6CgwKSGBABKS0tx8OBBZGdn47nnngMAVFdXQ5IkuLi4YMeOHRg+fDj8/PwafE4tpVIJpVJpTvhEZKVKS4HPPgOys433K5VAmzb625WVQHGxKBcUAOvWAf/5DzBlCtCvX4uFS0QyMCthcXNzQ2RkJNLT0zF58mTd/vT0dEycONHkeE9PTxw5csRo34oVK7Br1y58+eWXCAkJAQDExMQgPT0dc+fO1R23Y8cOxMbGmnUxRGR7fv0VeO894No1/b42bYARI4ChQwF3d/1+SQJOnQLS0sRP7ePXrAEefxwYPLhFQyeiFmR2Req8efOQmJiIqKgoxMTEYM2aNcjNzcWsWbMAiKaaK1euYN26dXByckLPnj2NHt+hQweoVCqj/XPmzMHgwYOxdOlSTJw4EVu2bMHOnTuxb9++Jl4eEVmz69dFslJUJG57eopmnrg4wM3N9HiFAggLE9u5c8B//yv6uUiSqKEpKwNGjWrZayCilmF2wpKQkICioiIsWbIEeXl56NmzJ9LS0hAcHAwAyMvLu+ecLDXFxsZi48aNWLhwIV577TWEhoYiNTUV0dHR5oZHRDYiLw9YtkzfxNOhA5CUBLRr17DHd+4MPP88sHkzsH272PfllyJpGT9eJDdEZD/MnofFWnEeFiLbcemSSFZu3hS3AwJEsqJWm38uSRIddbds0e8bNQp4+GEmLUS2oFnmYSEiaqriYuNkpWNH4IUXGpesACIpGTcOeOQR/b70dLERkf1gwkJELUaSgJQUfbISGgrMmwe0bt30c48cCUydqq9V+fprwMzWaSKyYkxYiKjFfPstcOKEKKvVwLPPGo8Caqq4OEA7aXZVFfDhh6JPCxHZPiYsRNQiLl0SHWS1nnzSMjUrNf3ud8BvYwBw7ZroiEtEto8JCxE1u/Jy4KOPxMRvgOgUGx7ePM/l4gLMmKEfFr1nD3DoUPM8FxG1HCYsRNTsNm0Sw5gBICgImDSpeZ/P1xdISNDf/ve/9cOnicg2MWEhomZ19Cjw3Xei7Ooqaj9aYu2fQYOAiAhRvnVLJC32MYkDkWNiwkJEzaaqCvj8c/3tKVMAf/+WeW6FAkhMBNq2FbePHwcOH26Z5yYiy2PCQkTNZt8+/RpBXbqIUTwtqVUr4NFH9be/+gqorm7ZGIjIMpiwEFGzuHtXrKSsJdfMs337ivleACA/H/jhh5aPgYiajgkLETWL9HSgtFSUIyOB3xZnb3EKhUiWtLZu5dwsRLaICQsRWVxJCbBjhyg7OTX/qKB7CQ3Vd8DVaICdO+WNh4jMx4SFiCzuP/8Rc68AwJAhYiVmuU2eLJInQKzurNHIGw8RmYcJCxFZVF6e6GwLACoVMH68vPFo+foCgweLclkZ8L//yRsPEZmHCQsRWdTmzfr5TsaMAdq0kTceQ+PHA0qlKO/Zox/BRETWjwkLEVnMxYv6afDbtgVGjJA1HBOenvrFEaurWctCZEuYsBCRxaSn68vjxunX87Emo0aJ+VkA4MAB4Ndf5Y2HiBqGCQsRWcSNG0BWlii3bg3ExsobT12USmDoUFGurgZ27ZI1HCJqICYsRGQR336rn0V26FCxbpC1GjpUv57Rnj1ikjsism5MWIioye7c0Y8McnHR12BYK09PIDpalO/e1cdORNaLCQsRNdnevfpaipgY6xoZVJdRo/TlXbu4xhCRtWPCQkRNUlVl3A9k5Ej5YjGHvz/Qs6coFxUBP/8sbzxEVD8mLETUJFlZ+pE2vXsDfn7yxmMOw1qWHTv088cQkfVhwkJEjSZJxkOZDRMAW9CtGxAUJMoXLwK//CJvPERUNyYsRNRop08DubmiHBwMdOkibzzmUiiMkyzD5IuIrAsTFiJqtG+/1ZdHjRIJgK2JigK8vET50CHg+nV54yGi2jFhIaJGKSkBjhwR5bZtgX79ZA2n0ZydgWHD9Ld/+EG+WIiobkxYiKhRMjP1Q4FjY8UHv62KiQGcfns3zMjgEGcia8SEhYjMJknGNRGDBskXiyV4eooRToCoOTp6VN54iMgUExYiMtuZM0BBgSiHhQE+PvLGYwmGSRdnviWyPkxYiMhshh/oDzwgXxyW1LOn6IsDiL45JSWyhkNENTQqYVmxYgVCQkKgUqkQGRmJvXv31nnsvn37MGjQILRr1w7u7u4ICwvDe++9Z3RMSkoKFAqFyXaXK5IRWZ3bt/Wzwnp4AH37yhqOxTg5ib4sgOjD8uOP8sZDRMZczH1AamoqkpKSsGLFCgwaNAirV6/G2LFjcfz4cXTs2NHk+FatWuG5555D79690apVK+zbtw/PPPMMWrVqhaefflp3nKenJ06dOmX0WJVK1YhLIqLmtH8/UFEhytHR1r0qs7kGDQK2bRPlffuA+HjbHKpNZI/MrmF59913MWPGDMycORPh4eFYtmwZgoKCsHLlylqPj4iIwGOPPYYePXqgU6dOmDp1KkaPHm1SK6NQKODn52e0EZH1Mexsay/NQVrt24s+OYDoo8OZb4msh1kJS3l5ObKyshAfH2+0Pz4+HhkZGQ06R3Z2NjIyMjBkyBCj/Tdv3kRwcDACAwPx4IMPIjs7u97zlJWVQaPRGG1E1Lxyc41ntg0MlDee5sDOt0TWyayEpbCwEFVVVfD19TXa7+vri/z8/HofGxgYCKVSiaioKMyePRszZ87U3RcWFoaUlBRs3boVGzZsgEqlwqBBg3DmzJk6z5ecnAy1Wq3bgrQLghBRs7Hn2hWtiAjRNwcQCzvevi1vPEQkNKrTraJGo64kSSb7atq7dy8OHjyIVatWYdmyZdiwYYPuvoEDB2Lq1Kno06cP4uLi8Pnnn6Nr16744IMP6jzf/PnzUVJSotsuXbrUmEshogaqqAB++kmUXV2B/v3ljae5uLqKvjmAuOYDB+SNh4gEszrd+vj4wNnZ2aQ2paCgwKTWpaaQkBAAQK9evXDt2jW88cYbeOyxx2o91snJCf3796+3hkWpVEKpVJoTPhE1QU4OcOeOKEdFAe7usobTrB54ANi9W5QzMoAaLdhEJAOzaljc3NwQGRmJ9BpLmqanpyM2NrbB55EkCWVlZfXen5OTA39/f3PCI6JmZFjTMHCgfHG0hMBAff+cCxe4ICKRNTB7WPO8efOQmJiIqKgoxMTEYM2aNcjNzcWsWbMAiKaaK1euYN26dQCA5cuXo2PHjgj7rev9vn378M477+D555/XnXPx4sUYOHAgunTpAo1Gg/fffx85OTlYvny5Ja6RiJro9m39dPVqNdC1q7zxtIT+/YHLl0X5wAFg3Dh54yFydGYnLAkJCSgqKsKSJUuQl5eHnj17Ii0tDcHBwQCAvLw85GqHEQCorq7G/Pnzcf78ebi4uCA0NBRvvfUWnnnmGd0xxcXFePrpp5Gfnw+1Wo2IiAjs2bMHAwYMsMAlElFTZWcDVVWiHBWlXyjQnvXvD2zeLMpMWIjkp5AkSZI7CEvQaDRQq9UoKSmBp6en3OEQ2ZX33gNOnhTlV14BfuuSZveWLgXOnRPlRYuA++6TNx4ie9TQz28H+J5ERE2h0QDaSah9fIBOnWQNp0UZVvJytBCRvJiwEFG9Dh4EtPWwAwY41lT1kZH66z1wQP97IKKWx4SFiOplWLNgr3Ov1MXTUz9Vf2GhGDFERPJgwkJEdSoq0vfhuO8+ICBA3njkYJiksVmISD5MWIioTo5cu6IVEQG4/Dae8sABoLpa3niIHBUTFiKqExMWsa5Qjx6irNEAp0/LGw+Ro2LCQkS1ysvTT5zWubMYIeSoOFqISH5MWIioVvv368uOWrui1bs3oF267OefgcpKeeMhckRMWIjIhCSJ4cyAGNYbFSVvPHJzcwP69BHl27eBEyfkjYfIETFhISITeXlAQYEod+kihvc6OsOkLTtbvjiIHBUTFiIy8fPP+nK/fvLFYU26d9c3Cx06xNFCRC2NCQsRmTCsQejbV7YwrIqrq3600M2bwC+/yBsPkaNhwkJERgoL9aODOnUCvLxkDceqGNY2GdZCEVHzY8JCREYMa1fYHGSsVy/9JHLZ2VxbiKglMWEhIiNsDqqbSqVfW6i4GLh4UdZwiBwKExYi0ikpAc6eFeWAAMDXV954rBGbhYjkwYSFiHRycvTliAjZwrBqvXuLuWkANgsRtSQmLESkY9gcxISldm3aiLlpADFXTV6evPEQOQomLEQEALh1Czh1SpR9fIDAQHnjsWaGyRwnkSNqGUxYiAgAcPiwfjK0iAh9sweZYsJC1PKYsBARADYHmcPLS8xRAwCXLom5a4ioeTFhISKUlQHHj4uypyfQubO88dgCw6TOsLMyETUPJixEhGPHgIoKUe7bl81BDWGYsHB4M1HzY8JCRDh0SF/mZHEN4+sL+PuL8rlzYn0hImo+TFiIHFx1NXD0qCgrlUC3bvLGY0t69xY/JUn/OySi5sGEhcjBnT+vrx3o3l2/Vg7dmzZhAYxrqYjI8piwEDk4ww9aww9gurfOnYFWrUT5+HGgslLeeIjsGRMWIgd3+LD4qVCI1Yip4Zyc9L+zu3eBM2fkjYfInjFhIXJghYX6qeVDQsS082QeNgsRtQwmLEQOTFu7ArA5qLF69ACcnUX58GEuhkjUXJiwEDkww4SlTx/54rBlKhXQtasoFxVxMUSi5tKohGXFihUICQmBSqVCZGQk9u7dW+ex+/btw6BBg9CuXTu4u7sjLCwM7733nslxmzZtQvfu3aFUKtG9e3ds3ry5MaERUQPdvQucPi3K7drp5xQh87FZiKj5mZ2wpKamIikpCQsWLEB2djbi4uIwduxY5Obm1np8q1at8Nxzz2HPnj04ceIEFi5ciIULF2LNmjW6YzIzM5GQkIDExEQcOnQIiYmJmDJlCn766afGXxkR1evYMaCqSpR79+bstk1hmLAY1loRkeUoJMm8Ftfo6Gj069cPK1eu1O0LDw/HpEmTkJyc3KBzPPTQQ2jVqhU++eQTAEBCQgI0Gg22bdumO2bMmDHw8vLChg0bGnROjUYDtVqNkpISeHp6mnFFRI5p7Vrgxx9FOSkJCA+XNRybt3gxcPWqSPz+9jd2YCZqqIZ+fptVw1JeXo6srCzEx8cb7Y+Pj0dGRkaDzpGdnY2MjAwMGTJEty8zM9PknKNHj27wOYnIPNXVwJEjoqxUAl26yBuPPTCc9Vb7uyUiyzErYSksLERVVRV8fX2N9vv6+iI/P7/exwYGBkKpVCIqKgqzZ8/GzJkzdffl5+ebfc6ysjJoNBqjjYga5tw54NYtUe7Rg7PbWoJhsxATFiLLa1SnW0WNxm5Jkkz21bR3714cPHgQq1atwrJly0yaesw9Z3JyMtRqtW4LCgoy8yqIHBeHM1teSAjQurUoHzvGWW+JLM2shMXHxwfOzs4mNR8FBQUmNSQ1hYSEoFevXnjqqacwd+5cvPHGG7r7/Pz8zD7n/PnzUVJSotsuXbpkzqUQOTTD2W179pQ3Fnvh5KT/XZaV6UdgEZFlmJWwuLm5ITIyEunp6Ub709PTERsb2+DzSJKEsrIy3e2YmBiTc+7YsaPecyqVSnh6ehptRHRvhnOFcHZbyzKcy4arNxNZltkt1/PmzUNiYiKioqIQExODNWvWIDc3F7NmzQIgaj6uXLmCdevWAQCWL1+Ojh07IiwsDICYl+Wdd97B888/rzvnnDlzMHjwYCxduhQTJ07Eli1bsHPnTuzbt88S10hEBgw/SLl2kGWFh4uaFm2n5ilT5I6IyH6YnbAkJCSgqKgIS5YsQV5eHnr27Im0tDQEBwcDAPLy8ozmZKmursb8+fNx/vx5uLi4IDQ0FG+99RaeeeYZ3TGxsbHYuHEjFi5ciNdeew2hoaFITU1FdHS0BS6RiAwxYWk+7u7A/feL5qCCArF16CB3VET2wex5WKwV52EhureKCmDePKC8HPD0BN5+mxPGWdr27cBXX4nyo48Cw4bJGw+RtWuWeViIyLadOSOSFUB0EGWyYnmGnZg5vJnIcpiwEDkQw+Ygjg5qHgEBQNu2onz6tD5BJKKmYcJC5EC0CYuTE6fiby4Khb5vUEUFhzcTWQoTFiIHcf06cO2aKIeGAh4e8sZjz9gsRGR5TFiIHASbg1pOWBjg7CzKR4+K9YWIqGmYsBA5CCYsLUel0i8oWVgohjcTUdMwYSFyABUVwKlToty2LXDffbKG4xB69NCX2SxE1HRMWIgcwKlTImkBOJy5pRhOysdp+omajgkLkQNgc1DL8/MD2rUT5TNnxIKIRNR4TFiI7Jwk6ZskOJy55SgU+mahykrg5El54yGydUxYiOxcQYHo+AmIjqAqlbzxOBI2CxFZDhMWIjvH5iD5dOsGuPy2xCyHNxM1DRMWIjt37Ji+zISlZSmV+uHNN27oJ+4jIvMxYSGyY4ZTw7dtC/j7yxqOQzIc3myYPBKReZiwENmx06f1w5l79OBwZjkwYSGyDCYsRHbs+HF9mc1B8vD3B7y8RNkwgSQi8zBhIbJj2m/0CoVY34ZankIBdO8uyly9majxmLAQ2akbN4C8PFHu3JmrM8vJsHaLzUJEjcOEhchOGX4wGvajoJYXFiYm7QOMm+mIqOGYsBDZKSYs1sPDAwgJEeW8PKCoSN54iGwRExYiO1RVBZw4IcqtWgEdO8obDxknjaxlITIfExYiO3T+PHD3rih3765vjiD5cHgzUdPwbYzIDrE5yPp07ChquwBR+1VVJW88RLaGCQuRHTJMWLRDakleTk761+LuXVELRkQNx4SFyM6UlgIXL4pyYCCgVssbD+mxWYio8ZiwENkZww6dbA6yLoa1XUxYiMzDhIXIznA6fuulVgNBQaJ88aKoDSOihmHCQmRHJEn/zV2pFDPcknUxrGXh8GaihmPCQmRHLl/Wf2vv1g1wcZE3HjJlWOvFhIWo4ZiwENkRDme2fp07i9ovQCQskiRvPES2ggkLkR1hwmL9XFxE7RcAaDSiVoyI7q1RCcuKFSsQEhIClUqFyMhI7N27t85jv/rqK4waNQrt27eHp6cnYmJisH37dqNjUlJSoFAoTLa72qk6ieieysqAs2dF2ccHaN9e3niobuzHQmQ+sxOW1NRUJCUlYcGCBcjOzkZcXBzGjh2L3NzcWo/fs2cPRo0ahbS0NGRlZWHYsGGYMGECsrOzjY7z9PREXl6e0aZSqRp3VUQO6NQp/eyprF2xbpyPhch8ZnfJe/fddzFjxgzMnDkTALBs2TJs374dK1euRHJyssnxy5YtM7r9f//3f9iyZQv+85//ICIiQrdfoVDAz8/P3HCI6DdsDrId7duLWrDCQuCXX0TtmLZfCxHVzqwalvLycmRlZSE+Pt5of3x8PDIyMhp0jurqapSWlsLb29to/82bNxEcHIzAwEA8+OCDJjUwRFQ/bdOCk5O+jwRZJ4VC3yxUVQWcPi1vPES2wKyEpbCwEFVVVfD19TXa7+vri/z8/Aad4+9//ztu3bqFKVOm6PaFhYUhJSUFW7duxYYNG6BSqTBo0CCcOXOmzvOUlZVBo9EYbUSOqrAQKCgQ5dBQgK2p1o/NQkTmadQsDQqFwui2JEkm+2qzYcMGvPHGG9iyZQs6dOig2z9w4EAMHDhQd3vQoEHo168fPvjgA7z//vu1nis5ORmLFy9uTPhEdofNQbYnLEzUhlVXs+MtUUOYVcPi4+MDZ2dnk9qUgoICk1qXmlJTUzFjxgx8/vnnGDlyZP1BOTmhf//+9dawzJ8/HyUlJbrt0qVLDb8QIjtj+IHH1Zltg0qln4n42jWgqEjeeIisnVkJi5ubGyIjI5Genm60Pz09HbGxsXU+bsOGDZg+fTrWr1+P8ePH3/N5JElCTk4O/P396zxGqVTC09PTaCNyRFVVwMmToty6NdCxo7zxUMOxWYio4cwe1jxv3jx8+OGH+Pjjj3HixAnMnTsXubm5mDVrFgBR8zFt2jTd8Rs2bMC0adPw97//HQMHDkR+fj7y8/NRUlKiO2bx4sXYvn07zp07h5ycHMyYMQM5OTm6cxJR3c6dA7RTFoWHiw6dZBs4HwtRw5ndhyUhIQFFRUVYsmQJ8vLy0LNnT6SlpSE4OBgAkJeXZzQny+rVq1FZWYnZs2dj9uzZuv1PPPEEUlJSAADFxcV4+umnkZ+fD7VajYiICOzZswcDBgxo4uUR2T/DDzr2X7EtHTsCrVoBt26JWrLqatGvhYhMKSTJPlay0Gg0UKvVKCkpYfMQOZT/+z/g4kVRfvttQK2WNx4yz4cfAgcOiPJLL4lRXkSOpKGf38zliWzYzZuAtkIzMJDJii0ybBZiPxaiujFhIbJhJ07oV/vl6CDbxISFqGGYsBDZMMMPOCYstqltWyAgQJQvXhT9WYjIFBMWIhslSfoOt66uwP33yxsPNZ62s7QkiVozIjLFhIXIRl29CmhnB+jWTSQtZJs4HwvRvTFhIbJRbA6yH/ffr084jx/X90siIj0mLEQ2ivOv2A9XV6BrV1EuLgby8mQNh8gqMWEhskHl5YB2qS1vb+AeS3mRDTBMOjnrLZEpJixENuj0aaCyUpR79OB0/PaAw5uJ6seEhcgGcXVm++PnB3h5ifKZM0BFhbzxEFkbJixENkj7DVyhAMLC5I2FLEOh0CefFRWiFo2I9JiwENmYGzeA/HxRDgkBPDzkjYcsh/1YiOrGhIXIxnB0kP0KD9f3R2I/FiJjTFiIbAznX7FfHh6i1gwQQ5t//VXeeIisCRMWIhtSXQ2cPCnKHh5Ap06yhkPNwDAJZbMQkR4TFiIbcuECcPu2KIeHA078D7Y7nKafqHZ8uyOyIYYfYOy/Yp86ddJ3pD5xQtSqERETFiKbwvlX7J+Tk36o+u3bolaNiJiwENmM27eB8+dF2d9fP8kY2R8ObyYyxYSFyEacOKFfxZfNQfaN/ViITDFhIbIR7L/iOLy8RC0aIGrVbt2SNx4ia8CEhcgGSJI+YXF1Bbp0kTcean7apFSSRO0akaNjwkJkA65eBYqLRblbN5G0kH3r2VNfZrMQERMWIpvA5iDHc//9+sT02DF9/yUiR8WEhcgGMGFxPK6uojYNAEpKgCtX5I2HSG5MWIisXFkZcOaMKPv4AB06yBsPtRw2CxHpMWEhsnKnTgFVVaLco4d+NV+yfxzeTKTHhIXIyrE5yHG1by9q1QDgl1+Au3fljYdITkxYiKyYJAFHj4qys7O+TwM5BoVC3yxUVSVq24gcFRMWIit2/TpQWCjK998PqFTyxkMtj81CRAITFiIrxuYg6tZN1K4BoraNw5vJUTFhIbJi2uYggAmLo1Iq9TMbFxUBBQXyxkMkl0YlLCtWrEBISAhUKhUiIyOxd+/eOo/96quvMGrUKLRv3x6enp6IiYnB9u3bTY7btGkTunfvDqVSie7du2Pz5s2NCY3IblRU6PssqNXAfffJGw/Jh81CRI1IWFJTU5GUlIQFCxYgOzsbcXFxGDt2LHJzc2s9fs+ePRg1ahTS0tKQlZWFYcOGYcKECcjOztYdk5mZiYSEBCQmJuLQoUNITEzElClT8NNPPzX+yohs3C+/iKQF4HBmR8eEhQhQSJJ5LaLR0dHo168fVq5cqdsXHh6OSZMmITk5uUHn6NGjBxISErBo0SIAQEJCAjQaDbZt26Y7ZsyYMfDy8sKGDRsadE6NRgO1Wo2SkhJ4enqacUVE1umLL4CdO0X5qaeAqCh54yH5SBLwyitiPSlXV+C997ieFNmPhn5+m1XDUl5ejqysLMTHxxvtj4+PR0ZGRoPOUV1djdLSUnh7e+v2ZWZmmpxz9OjR9Z6zrKwMGo3GaCOyJ9pv0goFEB4ubywkL4VCX8tSUQGcPi1vPERyMCthKSwsRFVVFXx9fY32+/r6Ij8/v0Hn+Pvf/45bt25hypQpun35+flmnzM5ORlqtVq3BQUFmXElRNatqAjIyxPlzp2BVq3kjYfkZzhNv2FnbCJH0ahOt4oajemSJJnsq82GDRvwxhtvIDU1FR1qLIhi7jnnz5+PkpIS3Xbp0iUzroDIuhl+IBl+UJHjCg8HnH57x2bCQo7IxZyDfXx84OzsbFLzUVBQYFJDUlNqaipmzJiBL774AiNHjjS6z8/Pz+xzKpVKKJVKc8InshmGH0i9eskXB1kPd3cxeeDp02Joc0EBF8Ikx2JWDYubmxsiIyORnp5utD89PR2xsbF1Pm7Dhg2YPn061q9fj/Hjx5vcHxMTY3LOHTt21HtOIntVUQGcOCHKajUQGChvPGQ92CxEjsysGhYAmDdvHhITExEVFYWYmBisWbMGubm5mDVrFgDRVHPlyhWsW7cOgEhWpk2bhn/84x8YOHCgribF3d0darUaADBnzhwMHjwYS5cuxcSJE7Flyxbs3LkT+/bts9R1EtmM06f1w5l79uRwZtLr1Qv46itRPnoUGD5c3niIWpLZfVgSEhKwbNkyLFmyBH379sWePXuQlpaG4OBgAEBeXp7RnCyrV69GZWUlZs+eDX9/f902Z84c3TGxsbHYuHEj1q5di969eyMlJQWpqamIjo62wCUS2Rb2X6G6+PsDXl6ifOoUUF4ubzxELcnseVisFedhIXvx2muif4KTE/Duu6LvApHWp58C2snFn3uOfZzI9jXLPCxE1Ly0nSkB0cGSyQrVZJigsB8LORImLERWhM1BdC9hYfrVm48c4erN5DiYsBBZEQ5npntRKoGuXUW5qAho4JydRDaPCQuRlSgv16/O7OUlOlgS1caw9o2LIZKjYMJCZCVOnQIqK0WZw5mpPoYJy5Ej8sVB1JKYsBBZCTYHUUP5+gI+PqJ85gxw96688RC1BCYsRFZAkvTflJ2dRcdKorooFPpalqoq4ORJeeMhaglMWIisQH6+6EAJiA6VXCaL7oXDm8nRMGEhsgKG/RA4nJkaomtXwNVVlDm8mRwBExYiK3D4sL7cu7d8cZDtcHPTNx0WFwOXLskaDlGzY8JCJLNbt4CzZ0XZ1xfo0EHeeMh2GCa3hkkvkT1iwkIks6NHgepqUWbtCpnDsB8LExayd0xYiGRm2H+FCQuZw8sLCAoS5YsXRdMQkb1iwkIko6oq/QgPDw8gNFTeeMj2GCa5nESO7BkTFiIZnT0L3Lkjyj166Be1I2ooJizkKJiwEMno0CF9mc1B1BjBwYCnpygfPw5UVMgbD1FzYcJCJCPtN2InJ1HDQmQuhULf+baiQr+AJpG9YcJCJJNr18QGiL4rrVrJGw/ZLsPaOcNaOyJ7woSFSCacLI4sJTwccHERZc56S/aKCQuRTAwTlj595IuDbJ9SqZ/19tdfgcuX5Y2HqDkwYSGSwe3bwC+/iHKHDpzdlpqOk8iRvWPCQiSDY8eMZ7dVKOSNh2wfp+kne8eEhUgGhh8oht+MiRrL2xsIDBTlCxcAjUbWcIgsjgkLUQsznN1WpQK6dJE3HrIfrGUhe8aEhaiFnT4t+rAAonaFs9uSpRh23s7Oli8OoubAhIWoheXk6MsREbKFQXYoOBho21aUT54E7t6VNRwii2LCQtSCJEmfsLi4cHZbsiyFAujbV5QrK0XnbiJ7wYSFqAXl5gLFxaIcFib6sBBZkjZhAYxr84hsHRMWohZk2K/A8IOFyFK6dgXc3UX5yBFR00JkD5iwELUg7TdehYKz21LzcHbWjxa6c0d08iayB0xYiFrItWtAXp4od+4MeHrKGw/ZLzYLkT1qVMKyYsUKhISEQKVSITIyEnv37q3z2Ly8PDz++OPo1q0bnJyckJSUZHJMSkoKFAqFyXaXXdzJjhiuosvmIGpOPXroF0M8dIiLIZJ9MDthSU1NRVJSEhYsWIDs7GzExcVh7NixyM3NrfX4srIytG/fHgsWLECfeurAPT09kZeXZ7Sp2COR7IjhN10mLNSclEqge3dRLi4GLl6UNRwiizA7YXn33XcxY8YMzJw5E+Hh4Vi2bBmCgoKwcuXKWo/v1KkT/vGPf2DatGlQq9V1nlehUMDPz89oI7IXGg1w7pwoBwRwsUNqfoZJMSeRI3tgVsJSXl6OrKwsxMfHG+2Pj49HRkZGkwK5efMmgoODERgYiAcffBDZ9/gPKysrg0ajMdqIrJVhtTxrV6glGC6qyX4sZA/MSlgKCwtRVVUFX19fo/2+vr7Iz89vdBBhYWFISUnB1q1bsWHDBqhUKgwaNAhnzpyp8zHJyclQq9W6LSgoqNHPT9Tc2BxELa1NG+D++0U5P19sRLasUZ1uFdq0/TeSJJnsM8fAgQMxdepU9OnTB3Fxcfj888/RtWtXfPDBB3U+Zv78+SgpKdFtly5davTzEzWnu3fFNOkA4OUFdOwobzzkOAyTY8NO30S2yKyExcfHB87Ozia1KQUFBSa1Lk0KyskJ/fv3r7eGRalUwtPT02gjskZHj+on7+rbV19NT9TcDMc5/PyzfHEQWYJZCYubmxsiIyORnp5utD89PR2xsbEWC0qSJOTk5MDf399i5ySSS1aWvszmIGpJ7dsDgYGifOECcOOGrOEQNYmLuQ+YN28eEhMTERUVhZiYGKxZswa5ubmYNWsWANFUc+XKFaxbt073mJzfGvBv3ryJ69evIycnB25ubuj+27i7xYsXY+DAgejSpQs0Gg3ef/995OTkYPny5Ra4RCL5lJWJ6dEB0aega1d54yHHExkJXL4syj//DIwcKW88RI1ldsKSkJCAoqIiLFmyBHl5eejZsyfS0tIQHBwMQEwUV3NOloiICF05KysL69evR3BwMC5cuAAAKC4uxtNPP438/Hyo1WpERERgz549GDBgQBMujUh+R44AFRWiHBEBOHFuaWphkZHAli2inJXFhIVsl0KS7GMORI1GA7VajZKSEvZnIauxerW+78DcuWKFZqKWtmQJcOWKKL/1luj8TWQtGvr5ze97RM3EsDmodWs2B5F8IiP1ZcM+VUS2hAkLUTM5epTNQWQdmLCQPeBbKFEzMfxgMPzAIGppfn5iSQhALBHx66/yxkPUGExYiJpBeblxc1C3bvLGQ2SYNHNOFrJFTFiImsHRoyJpAdgcRNaBzUJk6/g2StQM2BxE1sbfX2wAcPYsm4XI9jBhIbKw8nLg8GFRbtWKzUFkPdgsRLaMCQuRhbE5iKwVm4XIlvGtlMjC2BxE1iogwLhZqLhY1nCIzMKEhciCDCeLY3MQWSPWspCtYsJCZEGHDomkBQD69QOcneWNh6imqCh9+aef5IuDyFxMWIgsyPADIDpavjiI6uLvDwQFifLFi8C1a/LGQ9RQTFiILKS0FDh+XJS9vID775c3HqK6GCbT+/fLFweROZiwEFlIVhZQXS3KAwYACoW88RDVpX9//d/n/v2AJMkbD1FDMGEhshDDb6oDBsgXB9G9tG2r7xBeUCCahoisHRMWIgsoLBTDRAExdDQwUN54iO7FMKlm51uyBUxYiCzAsHaFnW3JFkREAC4uonzggL45k8haMWEhaiJJMk5Y+veXLxaihvLwAHr1EuXSUuDkSXnjIboXJixETXT5MpCXJ8r33w+0aydvPEQNZVgbyGYhsnZMWIiaiM1BZKt69gTc3UU5O1u/BhaRNWLCQtQE1dX6hMXJiWsHkW1xdRUzMgNihmbtKuNE1ogJC1ETnDmjX0CuZ0+xfhCRLWGzENkKJixETfDjj/oy514hW9Sli5iXBQCOHhUdcImsERMWoka6exc4eFCU3d2Bvn1lDYeoUZycgIEDRbm6mrUsZL2YsBA1UlaWvpPigAGiPwCRLYqN1Zf37eNU/WSdmLAQNdIPP+jLgwbJFwdRU/n66hfrzMvjVP1knZiwEDVCfr5+Kv777gM6dpQ3HqKmMky6DZNxImvBhIWoETIy9OVBg7gyM9m+yEhAqRTl/fs5JwtZHyYsRGaqqgIyM0XZ2Zmjg8g+KJVAVJQo370rJpIjsiZMWIjMdOwYoNGIcp8+QJs28sZDZClsFiJrxoSFyEzsbEv2qnNn0QEXAE6dAgoL5Y2HyFCjEpYVK1YgJCQEKpUKkZGR2Lt3b53H5uXl4fHHH0e3bt3g5OSEpKSkWo/btGkTunfvDqVSie7du2Pz5s2NCY2oWWk0+unL27YFuneXNRwii1IojIc4G/bVIpKb2QlLamoqkpKSsGDBAmRnZyMuLg5jx45Fbm5urceXlZWhffv2WLBgAfr06VPrMZmZmUhISEBiYiIOHTqExMRETJkyBT9xBiOyMj/9JCbXAoCYGDHpFpE9GThQ34k8I0P/904kN4UkmTdFUHR0NPr164eVK1fq9oWHh2PSpElITk6u97FDhw5F3759sWzZMqP9CQkJ0Gg02LZtm27fmDFj4OXlhQ0bNjQoLo1GA7VajZKSEnh6ejb8gogaSJKAxYvFPBUA8Je/AB06yBsTUXP45z+BI0dE+U9/Anr0kDcesm8N/fw26/theXk5srKyEB8fb7Q/Pj4eGU2oO8zMzDQ55+jRo+s9Z1lZGTQajdFG1Jx++UWfrHTpwmSF7Jdh36w9e+SLg8iQWQlLYWEhqqqq4KvtlfUbX19f5OfnNzqI/Px8s8+ZnJwMtVqt24KCghr9/EQN8d13+vLgwbKFQdTsevfWL4h46BDw66+yhkMEoJGdbhU1ZsmSJMlkX3Ofc/78+SgpKdFtly5datLzE9WnpAT4+WdRbtMG6NdP3niImpOzMxAXJ8qSxFoWsg5mJSw+Pj5wdnY2qfkoKCgwqSExh5+fn9nnVCqV8PT0NNqImsvevfrOhw88ALi4yBsPUXN74AF9p/K9e4HKSnnjITIrYXFzc0NkZCTS09ON9qenpyPWcCycmWJiYkzOuWPHjiadk8hSqqrEGzYgRk+wOYgcQdu2QESEKJeWcuZbkp/Z3xPnzZuHxMREREVFISYmBmvWrEFubi5mzZoFQDTVXLlyBevWrdM9JicnBwBw8+ZNXL9+HTk5OXBzc0P33yaxmDNnDgYPHoylS5di4sSJ2LJlC3bu3Il9+/ZZ4BKJmubQIaC4WJT79AG8vWUNh6jFDB0KZGWJ8u7dQP/+soZDDs7shCUhIQFFRUVYsmQJ8vLy0LNnT6SlpSE4OBiAmCiu5pwsEdo0HUBWVhbWr1+P4OBgXLhwAQAQGxuLjRs3YuHChXjttdcQGhqK1NRUREdHN+HSiCzDsLPt0KFyRUHU8rp0AQICgKtXxerkly8DgYFyR0WOyux5WKwV52Gh5pCXB7zxhij7+op5WLgyMzmS774DtNNhxcUBU6fKGg7ZoWaZh4XI0RjWrgwZwmSFHM/AgYBKJco//QTcvi1vPOS4mLAQ1eHuXeDHH0XZzU1MxU/kaFQqkbQAQHm5/n+CqKUxYSGqw08/iaQFAKKjAQ8PeeMhkoth363vvhNzsxC1NCYsRLWorgZ27tTfHjJEvliI5ObvD3TrJsrXrunXGSJqSUxYiGpx6BBQUCDKYWEAV34gRzdypL68Y4d8cZDjYsJCVIMkGb8h11iXk8gh9eolaloA4MwZ4Px5eeMhx8OEhaiGs2eBc+dE+b77gN/mNyRyaAoFMGqU/jZrWailMWEhqqFm7QqHMhMJ0dGAWi3K2dn6ZlOilsCEhchAfr7ovwKItVSiomQNh8iquLgAw4eLsiQBNZaAI2pWTFiIDBi+AY8YwVWZiWoaPBhQKkU5M1MsjEjUEpiwEP1Go9FPiqVScVVmotp4eIgp+gGgosJ4Nmii5sSEheg3u3YBlZWiPHiwfjpyIjI2YgTg9Nunx+7dQFmZvPGQY2DCQgQxo+3334uys7N4Qyai2nl7A/37i/KtW0BGhrzxkGNgwkIEUbuiXdRtwADR4ZaI6mY4P9E334jmIaLmxISFHN6dO/rOtk5OwLhx8sZDZAsCA4E+fUS5uBjYu1fWcMgBMGEhh/ftt/ralYEDgQ4d5I2HyFZMmKAvb9smVnMmai5MWMih3b5tXLsyfry88RDZkqAgoF8/UdZogD175I2H7BsTFnJoO3eKDrcAEBsL+PjIGw+RrZkwQT8b9DffcMQQNR8mLOSwbt0SzUGAGBnEvitE5gsIACIjRbm0VD/ajsjSmLCQw0pP19euDBoEtGsnbzxEturBB41rWbT/V0SWxISFHFJpqRjKDIjp98eOlTceIlvm7288L8vu3fLGQ/aJCQs5pO3b9W3tDzwgJsIiosYzrGXZsUM/8o7IUpiwkMMpLNR/A2TtCpFl+PoC0dGifPs2kJYmbzxkf5iwkMPZtEm/ZtCoUZzVlshSJk4EXF1FedcuoKBA3njIvjBhIYdy5gzw88+i7OkJjBkjbzxE9sTbW3wJAICqKuCrr+SNh+wLExZyGJIEfPGF/vbEiVyRmcjSxowRXwYAIDsbOH1a3njIfjBhIYfx44/AxYuiHBgoJoojIstSKoFJk/S3P/8cqK6WLRyyI0xYyCGUlQGbN+tvP/KImIqfiCwvJkZM2w8Aly6JLwtETcW3bHII27cDJSWi3KcPEBYmbzxE9szJCZgyRX9782ZO2U9Nx4SF7F5RkZgXAhBvpA8/LG88RI6ga1egb19R1mjEas5ETcGEheyaJAHr1wMVFeL2sGFivggian4PPyzW6QJELeeVK/LGQ7atUQnLihUrEBISApVKhcjISOzdu7fe47///ntERkZCpVKhc+fOWLVqldH9KSkpUCgUJttdLkhBTXTgAHD0qCi3bStWliWiltGhg37qgOpq4JNP2AGXGs/shCU1NRVJSUlYsGABsrOzERcXh7FjxyI3N7fW48+fP49x48YhLi4O2dnZePXVV/GnP/0JmzZtMjrO09MTeXl5RpuKY06pCW7dEiMUtB57DHB3ly8eIkc0dqy+VvP8ea7mTI1ndsLy7rvvYsaMGZg5cybCw8OxbNkyBAUFYeXKlbUev2rVKnTs2BHLli1DeHg4Zs6ciT/84Q945513jI5TKBTw8/Mz2oia4osvxCKHANCvn749nYhajqsrkJiov715M3DjhnzxkO0yK2EpLy9HVlYW4uPjjfbHx8cjIyOj1sdkZmaaHD969GgcPHgQFdqOBQBu3ryJ4OBgBAYG4sEHH0R2dna9sZSVlUGj0RhtRFonTgCZmaLs7g48+qi88RA5si5dgMGDRbmsTPQrkyR5YyLbY1bCUlhYiKqqKvjW6LXo6+uL/Pz8Wh+Tn59f6/GVlZUoLCwEAISFhSElJQVbt27Fhg0boFKpMGjQIJw5c6bOWJKTk6FWq3VbkHbQPzm88nLg00/1tx9+GFCr5YuHiICHHtL/Hx45AmRlyRsP2Z5GdbpVaNcQ/40kSSb77nW84f6BAwdi6tSp6NOnD+Li4vD555+ja9eu+OCDD+o85/z581FSUqLbLl261JhLITu0datYkRkQ3+weeEDeeIhI1HQ+9pj+9saN+iZbooYwK2Hx8fGBs7OzSW1KQUGBSS2Klp+fX63Hu7i4oF27drUH5eSE/v3711vDolQq4enpabQRnTgBpKeLsouLaDuvJ5cmohYUEaHvS1ZaCqxbx6YhajizEhY3NzdERkYiXfuJ8Jv09HTE1rEwS0xMjMnxO3bsQFRUFFy165DXIEkScnJy4O/vb0545OBKS4GPP9bfnjSJc64QWZvHHwdatxblw4eB776TNRyyIWY3Cc2bNw8ffvghPv74Y5w4cQJz585Fbm4uZs2aBUA01UybNk13/KxZs3Dx4kXMmzcPJ06cwMcff4yPPvoIL774ou6YxYsXY/v27Th37hxycnIwY8YM5OTk6M5JdC+SBKSkiBk1AaBHD2DkSFlDIqJaqNXAk0/qb3/5JXD5snzxkO1wMfcBCQkJKCoqwpIlS5CXl4eePXsiLS0NwcHBAIC8vDyjOVlCQkKQlpaGuXPnYvny5QgICMD777+Phw3mRy8uLsbTTz+N/Px8qNVqREREYM+ePRgwYIAFLpEcwa5d+gni2rQBpk9nUxCRterZExgxAvj2W6CyEvjXv4AFCwA3N7kjI2umkCT7aEHUaDRQq9UoKSlhfxYHc+kS8NZb4o0PAP70J1HDQkTWq7JS/N9qx0vExQFTp8obE8mjoZ/fXEuIbFpZmfh2pk1WRo5kskJkC1xcgJkz9bUqe/cCP/8sb0xk3ZiwkM2SJGDtWuDaNXG7Y0dg8mR5YyKihvPzM57UMSWFCyRS3ZiwkM3auhXQToisUolvay5m98oiIjnFxgL9+4tyWRmwfDnnZ6HaMWEhm3TgAJCWJsoKBfD00xzCTGSLFApg2jTgt3EbKCoCVq3SN/MSaTFhIZtz4QLw73/rbz/yCPutENkyNzfg2Wf1U/f/8gvXGyJTTFjIphQXAytWANp1MwcNAoYPlzUkIrKAtm1F0qKdT/SHH8SwZyItJixkM27fBv75T6CkRNzu0kXMmsn5VojsQ6dOwBNP6G9/+SVHDpEeExayCXfvAh98oJ+zoV074Jln2MmWyN707w+MGyfKkgR8+KF+UkhybExYyOpVVIiRA+fOidtt2ojJ4dq0kTcuImoev/udGD0EAFVVwMqVwMmT8sZE8mPCQlatslK8WZ0+LW57eABJSWL+BiKyTwqFWGk9KkrcrqwUfdfOnpU3LpIXExayWlVVojr42DFxW6UC5swBAgPljYuImp+TE/CHPwB9+ojbZWXA++8DFy/KGxfJhwkLWaWyMvGNSjsxnKsr8NxzolMeETkGZ2fgqaeA8HBx++5d4N13gVOn5I2L5MGEhaxOaal4U9J2tHNxAWbPFqOCiMixuLqK4c7a//+7d0VNy8GD8sZFLY8JC1mVwkLg7bfF5HAA4O4umoG037CIyPG4uQHPPw/06iVuV1aK5uLdu+WNi1oWExayGpcuAUuXAgUF4nbbtsCf/wx07SprWERkBZRKUdMyaJC4LUnAxo3AV19xRlxHwVksyCr8+CPw6af6GWz9/cXQZW9veeMiIuvh5CRGD6nV+rXEtm8H8vKAJ58UowjJfrGGhWRVUSESlbVr9clKaKioWWGyQkQ1KRTAxInGs1wfPgy8+SZHENk7Jiwkm8JC0QS0d69+X1wcMHcu0KqVfHERkfUbMkT0a9G+VxQVif5ve/awicheKSTJPl5ajUYDtVqNkpISeHp6yh0O1UOSgP37Rfvz7dtin6sr8PvfAzEx8sZGRLbl11+BNWv0M2EDYsK5Rx/lbNi2oqGf30xYqEX9+qtoAjJcG6RDB7EuECeEI6LGqKwUnW8NV3du3VokLVFRXCDV2jFhIasiSaKq9quvxDwKWv37A1OnillsiYia4uBB4LPP9DW3ANC7t6i9bdtWtrDoHhr6+c1RQtTszp0Ty8QbrgOiVotOc337yhYWEdmZqCgxDcKGDcDPP4t9hw+LtcjGjweGDRPNz2SbWMNCzebaNeDrr/VvHFoPPAA8/DCHIBJR8/n5Z5G4aDT6fV5ewKRJwIABYog0WQc2CZFsfv0V+OYb0QRUXa3f7+sralXCwuSLjYgcx+3bwKZNwA8/GI8cCgwUQ6N79WL/FmvAhIVa3OXLQHq6GAFkmKh4egITJoiaFX6rIaKWduWK6D9n2NkfEBNUxseLGhcXdpCQDRMWahHV1cCJE8DOncDx48b3KZXizWDUKFEmIpLTqVOixqXmBHNqNTB8uJj2n0OhWx4TFmpW164BmZliKy42vq9VKzGp0/Dh/OcnIusiScChQ8COHcYDAQBRA9yrFxAbK346O8sTo6NhwkIWV1go/tGzskz/0QHAxwcYOVL8s7NGhYis3blzInHJyTGdHbdNGzHqqG9foEsXJi/NiQkLNVl1NZCbK4YFHjok+qjU5OQE9OwpkpQ+fdhHhYhsT0GB6Jj744+mNcaAGNHYs6d4jwsP59IhlsaEhcxWXS06p506JbbTp40neTMUECCSlOho0amWiMjWafvkZWSIWpfKStNjFArgvvuAbt3E1qULp2hoKiYsVK+qKvGtIjdXdEC7cEGUtSsm16ZTJ1E92qeP6F3P4YBEZK/u3hWjinJyxM87d+o+tn178f4YHCy2++5jLYw5mjVhWbFiBf72t78hLy8PPXr0wLJlyxAXF1fn8d9//z3mzZuHY8eOISAgAC+99BJmzZpldMymTZvw2muv4ezZswgNDcVf//pXTJ48ucExMWExJUli0qTCQrEVFABXrwJ5eaLTrOHQ49q0aSNmjQwLE9Nbc2prInJElZXAmTMicTl1SjSP3+uT09NTfLELCAD8/EQfPx8foF07zrZbU7NNzZ+amoqkpCSsWLECgwYNwurVqzF27FgcP34cHTt2NDn+/PnzGDduHJ566il8+umn+OGHH/Dss8+iffv2ePjhhwEAmZmZSEhIwF/+8hdMnjwZmzdvxpQpU7Bv3z5ER0ebG6Ldq64Gbt4UW2mpfisuNt6KiuqvManJx0d8S7j/flHVyVoUIiIxR0t4uNgAMSHdmTMieTl/Hrh0yfS9VqMR26lTpudr21bMutu2rX5Tq8WXxNat9T/d3PgebMjsGpbo6Gj069cPK1eu1O0LDw/HpEmTkJycbHL8yy+/jK1bt+LEiRO6fbNmzcKhQ4eQmZkJAEhISIBGo8G2bdt0x4wZMwZeXl7YsGFDg+KSs4ZFksRWVSW26mp9WbtVVoqfFRWirP1ZXi7K5eViKyvTb3fv6rfbt/VbXf1KGsrFRcw66+8vqi61VZmswiQiMl9Vlai5vnhRNK1ra7JLS5t2XhcX0T/G3V28P7u7ixGYKpXxTzc3/ebqqt9cXMRPZ2dRdnERZcPNyUlszs7yJUfNUsNSXl6OrKwsvPLKK0b74+PjkZGRUetjMjMzER8fb7Rv9OjR+Oijj1BRUQFXV1dkZmZi7ty5JscsW7bMnPAsrqgIeOstkYBokxJtubpav1ljLyA3N1H12K6daF/VVkf6+4vbHM1DRGQZzs5iuv/AQDH5nFZpqUhcrl8XzfLXr4vPlcJC4zWO6lJZqa+paSnaBMbJSSQw2p/a8gMPiPWY5GBWwlJYWIiqqir4+voa7ff19UV+fn6tj8nPz6/1+MrKShQWFsLf37/OY+o6JwCUlZWhrKxMd1vTTK9oS/6h1Mcw0/bwMK421JYNqxhVKlYlEhHJqU0bfV/AmrTJiLYJv6TEuJn/5k3L1qw3lPaLeF3Ky1smjto0avUERY1PQkmSTPbd6/ia+809Z3JyMhYvXtzgmBvD2VnUUBhmmYY/a2aihlVshlVuhtVwhtV02p+G1XlubqZVfiqVOI4JCBGRfXBxAby9xdYQ1dVipJK2u4D2p7Y7geFm2O2gosK4W4L2p7b7gmEXBsMWhKoqfcuCYeuCnF0HzEpYfHx84OzsbFLzUVBQYFJDouXn51fr8S4uLmjXrl29x9R1TgCYP38+5s2bp7ut0WgQFBRkzuXcU9u2wP/9n0VPSUREZDYnJ5EsOHJfQ7N6Mri5uSEyMhLp6elG+9PT0xEbG1vrY2JiYkyO37FjB6KiouD629iuuo6p65wAoFQq4enpabQRERGRfTK7SWjevHlITExEVFQUYmJisGbNGuTm5urmVZk/fz6uXLmCdevWARAjgv75z39i3rx5eOqpp5CZmYmPPvrIaPTPnDlzMHjwYCxduhQTJ07Eli1bsHPnTuzbt89Cl0lERES2zOyEJSEhAUVFRViyZAny8vLQs2dPpKWlITg4GACQl5eH3Nxc3fEhISFIS0vD3LlzsXz5cgQEBOD999/XzcECALGxsdi4cSMWLlyI1157DaGhoUhNTeUcLERERASAU/MTERGRjBr6+c3ZOIiIiMjqMWEhIiIiq8eEhYiIiKweExYiIiKyekxYiIiIyOoxYSEiIiKrx4SFiIiIrB4TFiIiIrJ6TFiIiIjI6pk9Nb+10k7Yq9FoZI6EiIiIGkr7uX2vifftJmEpLS0FAAQFBckcCREREZmrtLQUarW6zvvtZi2h6upqXL16FW3atIFCobDYeTUaDYKCgnDp0iW7XaPI3q+R12f77P0aeX22z96vsTmvT5IklJaWIiAgAE5OdfdUsZsaFicnJwQGBjbb+T09Pe3yj9CQvV8jr8/22fs18vpsn71fY3NdX301K1rsdEtERERWjwkLERERWT0mLPegVCrx+uuvQ6lUyh1Ks7H3a+T12T57v0Zen+2z92u0huuzm063REREZL9Yw0JERERWjwkLERERWT0mLERERGT1mLAQERGR1WPC0khlZWXo27cvFAoFcnJy5A7HYn73u9+hY8eOUKlU8Pf3R2JiIq5evSp3WBZx4cIFzJgxAyEhIXB3d0doaChef/11lJeXyx2aRf31r39FbGwsPDw80LZtW7nDabIVK1YgJCQEKpUKkZGR2Lt3r9whWcyePXswYcIEBAQEQKFQ4Ouvv5Y7JItKTk5G//790aZNG3To0AGTJk3CqVOn5A7LYlauXInevXvrJlOLiYnBtm3b5A6r2SQnJ0OhUCApKUmW52fC0kgvvfQSAgIC5A7D4oYNG4bPP/8cp06dwqZNm3D27Fn8v//3/+QOyyJOnjyJ6upqrF69GseOHcN7772HVatW4dVXX5U7NIsqLy/HI488gj/+8Y9yh9JkqampSEpKwoIFC5CdnY24uDiMHTsWubm5codmEbdu3UKfPn3wz3/+U+5QmsX333+P2bNn48cff0R6ejoqKysRHx+PW7duyR2aRQQGBuKtt97CwYMHcfDgQQwfPhwTJ07EsWPH5A7N4g4cOIA1a9agd+/e8gUhkdnS0tKksLAw6dixYxIAKTs7W+6Qms2WLVskhUIhlZeXyx1Ks3j77belkJAQucNoFmvXrpXUarXcYTTJgAEDpFmzZhntCwsLk1555RWZImo+AKTNmzfLHUazKigokABI33//vdyhNBsvLy/pww8/lDsMiyotLZW6dOkipaenS0OGDJHmzJkjSxysYTHTtWvX8NRTT+GTTz6Bh4eH3OE0qxs3buCzzz5DbGwsXF1d5Q6nWZSUlMDb21vuMKgW5eXlyMrKQnx8vNH++Ph4ZGRkyBQVNUVJSQkA2OX/XFVVFTZu3Ihbt24hJiZG7nAsavbs2Rg/fjxGjhwpaxxMWMwgSRKmT5+OWbNmISoqSu5wms3LL7+MVq1aoV27dsjNzcWWLVvkDqlZnD17Fh988AFmzZoldyhUi8LCQlRVVcHX19dov6+vL/Lz82WKihpLkiTMmzcPDzzwAHr27Cl3OBZz5MgRtG7dGkqlErNmzcLmzZvRvXt3ucOymI0bNyIrKwvJyclyh8KEBQDeeOMNKBSKereDBw/igw8+gEajwfz58+UO2SwNvT6tP//5z8jOzsaOHTvg7OyMadOmQbLiCZHNvT4AuHr1KsaMGYNHHnkEM2fOlCnyhmvMNdoLhUJhdFuSJJN9ZP2ee+45HD58GBs2bJA7FIvq1q0bcnJy8OOPP+KPf/wjnnjiCRw/flzusCzi0qVLmDNnDj777DOoVCq5w+HU/ID4JldYWFjvMZ06dcKjjz6K//znP0ZvllVVVXB2dsbvf/97/Pvf/27uUBuloddX2x/k5cuXERQUhIyMDKut5jT3+q5evYphw4YhOjoaKSkpcHKy/ry9Ma9hSkoKkpKSUFxc3MzRNY/y8nJ4eHjgiy++wOTJk3X758yZg5ycHHz//fcyRmd5CoUCmzdvxqRJk+QOxeKef/55fP3119izZw9CQkLkDqdZjRw5EqGhoVi9erXcoTTZ119/jcmTJ8PZ2Vm3r6qqCgqFAk5OTigrKzO6r7m5tNgzWTEfHx/4+Pjc87j3338fb775pu721atXMXr0aKSmpiI6Oro5Q2yShl5fbbT5bFlZmSVDsihzru/KlSsYNmwYIiMjsXbtWptIVoCmvYa2ys3NDZGRkUhPTzdKWNLT0zFx4kQZI6OGkiQJzz//PDZv3ozvvvvO7pMVQFyzNb9fmmPEiBE4cuSI0b4nn3wSYWFhePnll1s0WQGYsJilY8eORrdbt24NAAgNDUVgYKAcIVnU/v37sX//fjzwwAPw8vLCuXPnsGjRIoSGhlpt7Yo5rl69iqFDh6Jjx4545513cP36dd19fn5+MkZmWbm5ubhx4wZyc3NRVVWlmyfo/vvv1/3N2op58+YhMTERUVFRiImJwZo1a5Cbm2s3/Y5u3ryJX375RXf7/PnzyMnJgbe3t8n7jS2aPXs21q9fjy1btqBNmza6vkdqtRru7u4yR9d0r776KsaOHYugoCCUlpZi48aN+O677/DNN9/IHZpFtGnTxqS/kbZ/oyz9kGQZm2Qnzp8/b1fDmg8fPiwNGzZM8vb2lpRKpdSpUydp1qxZ0uXLl+UOzSLWrl0rAah1sydPPPFErde4e/duuUNrlOXLl0vBwcGSm5ub1K9fP7saErt79+5aX6snnnhC7tAsoq7/t7Vr18odmkX84Q9/0P1ttm/fXhoxYoS0Y8cOucNqVnIOa2YfFiIiIrJ6ttGAT0RERA6NCQsRERFZPSYsREREZPWYsBAREZHVY8JCREREVo8JCxEREVk9JixERERk9ZiwEBERkdVjwkJERERWjwkLERERWT0mLERERGT1mLAQERGR1fv/+t7SpngoWAQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import norm\n",
    "\n",
    "mu = 0\n",
    "sigma = 1\n",
    "normal = norm(mu, sigma)\n",
    "\n",
    "x = np.linspace(-4,  4, 100)\n",
    "y = normal.pdf(x)\n",
    "plt.plot(x, y, 'b', lw=2, alpha=0.6)\n",
    "plt.title('Función de densidad de la distribución N(0,1)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Con el siguiente código calculamos la media varianza y algunos cuantiles de la distribución $\\text{N}(0,1)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Distribución Normal(0,1); media = 0.0, varianza = 1.0\n",
      "Distribución Normal(0,1); probabilidades=[0.025 0.5   0.975]; cuantiles=[-1.96  0.    1.96]\n"
     ]
    }
   ],
   "source": [
    "# media y varianza de la distribución Normal(0,1)\n",
    "import numpy as np\n",
    "from scipy.stats import norm\n",
    "\n",
    "mu = 0\n",
    "sigma = 1\n",
    "normal = norm(mu, sigma)\n",
    "media, var = normal.stats(moments=\"mv\")\n",
    "print('Distribución Normal(0,1); media = {}, varianza = {}'.format(media, var, curtosis ))\n",
    "\n",
    "# calcula algunos cuantiles de la distribución\n",
    "prob = np.array([0.025,  0.5, 0.975])\n",
    "cuantiles = normal.ppf(prob)\n",
    "print('Distribución Normal(0,1); probabilidades={}; cuantiles={}'.format(prob, np.round(cuantiles,3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:#4CC9F0\">Distribución de Poisson</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La distribución de Poisson es una distribución discreta usada en problemas de conteos. La función de probabildiad tiene la forma\n",
    "\n",
    "$$\n",
    "p(y|\\theta) = \\frac{e^{-\\theta} \\theta^y }{y!}, \\theta >0, \\hspace{1mm} y =0,1,2, \\ldots.\n",
    "$$\n",
    "\n",
    "Denotamos esta densidad como $\\text{Pois}(\\theta)$ . En foma de la familia exponencial, la distribución puede escribirse como\n",
    "\n",
    "$$\n",
    "p(y|\\theta) = \\exp \\left[y\\log \\theta - \\theta  - log y! \\right]\n",
    "$$\n",
    "\n",
    "Entonces se tiene que:\n",
    "\n",
    "* la distribución está en forma canónica,: $a(y) = y$,\n",
    "* el parámetro natural es $b(\\theta)=log (\\theta)$, \n",
    "* $c(\\theta) = - \\theta$,\n",
    "* $d(y) = - \\log y!$.\n",
    "\n",
    "Se puede verificar que si una variable aleatoria $Y$ tiene distribución $\\text{Pois}(\\theta)$, entonces\n",
    "\n",
    "* $E[Y] = \\theta$,\n",
    "* $Var[Y] = \\theta$.\n",
    "\n",
    "El siguiente fragmento de código dibuja la función de probabilidad $\\text{Pois}(10)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAHACAYAAABDKXcJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA7MElEQVR4nO3deVyVdf7//yeyndwXEMUQ1BYxRB0owzQrFROzxqVMS63UhnBS4NMiqWNZaZpjpAnklmOpWS6TJp8SU3M7mQvaomN9JhVGYRQrUEnZrt8f/jhfj4ByWOTo9bjfbtdtPO/zfr+v13XOaXz6vq5zHRfDMAwBAACYWK2aLgAAAKCmEYgAAIDpEYgAAIDpEYgAAIDpEYgAAIDpEYgAAIDpEYgAAIDpEYgAAIDpEYgAAIDpEYgAXPdiY2PVokULpaen13QpAK5TBCKgBixevFguLi62zWKxqFmzZrr//vs1bdo0nTx5ssSYV199VS4uLg7tJzc3V6+++qq2bNni0LjS9hUQEKCHHnrIoXmqwtWOe82aNVq0aJH+93//V35+flW674CAALv3qW7duurcubOWLFni8FxbtmyRi4uLw+9FTZoyZYratWunoqIiW9uSJUv0+OOP6/bbb1etWrUUEBBQ5vizZ88qOjpavr6+slgs6tixoz7++OMS/e69915FR0dXwxEADjAAXHMffPCBIcn44IMPDKvVamzdutVYuXKlER0dbTRo0MBo3LixkZKSYjcmPT3dsFqtDu3n1KlThiRj8uTJDo0rbV/+/v5G3759HZqnKlzpuP/9738b3t7exsaNG6tl3/7+/sY999xjWK1Ww2q1Gp9++qnRuXNnQ5KRkJDg0FzZ2dmG1Wo1srOzq6XWqnb8+HGjTp06xqeffmrX3rNnTyMoKMh48sknjVtuucXw9/cvc45evXoZDRs2NJKSkoxNmzYZo0aNMiQZS5cuteu3ZcsWw93d3fjXv/5VHYcClAuBCKgBxYFo9+7dJZ47duyY4efnZ9SrV8/IzMys1H4cDUTnzp0r87maCkQ1qbRj/u2334z69esbt9xySw1VdW289NJLRosWLYzCwkK79ksf9+3bt8xAtH79ekOSsWzZMrv2Xr16Gb6+vkZBQYFde1BQkDF69OiqKR6oAE6ZAU6mZcuW+vvf/64zZ87o/ffft7WXdupo06ZNuu+++9SkSRPddNNNatmypQYOHKjc3FwdPXpU3t7ekqTXXnvNdtrnqaeesptv3759GjRokBo1aqQ2bdqUua9ia9asUXBwsCwWi1q3bq3Zs2fbPV98OvDo0aN27WWdMvriiy/Uo0cPNWjQQLVr11ZgYKCmTZt2xeMuKirSjBkz1LZtW3l6eqpp06YaPny4/vOf/9j1u++++xQUFKTdu3erW7duql27tlq3bq233nrL7jSQIxo2bKjbb79dx44ds7Vt375dPXr0UL169VS7dm116dJF69evv+rx//LLL3r88cfl6+srT09P+fj4qEePHtq/f7+tz5Xe42K//vqroqKi1KJFC3l4eKh169aaMGGCLly4YFeDi4uL/vrXv+rDDz9UYGCgateurQ4dOujzzz+365eXl6eFCxdq6NChqlXL/q+Jyx+XZc2aNapbt64effRRu/ann35aJ06c0K5du+zahw0bpmXLlunMmTPlmh+oagQiwAlFRETI1dVVW7duLbPP0aNH1bdvX3l4eGjRokX64osv9NZbb6lOnTrKy8tT8+bN9cUXX0iSRo4cKavVKqvVqkmTJtnNM2DAAN1yyy369NNPlZSUdMW69u/fr+joaMXExGjNmjXq0qWLxo0bp5kzZ1boOBcuXKiIiAgVFRUpKSlJ69at09ixY0sEm8s999xzevnll9WrVy+tXbtWr7/+ur744gt16dJFWVlZdn0zMzP1xBNP6Mknn9TatWvVp08fxcXF6aOPPqpQzfn5+Tp27JgtbH799dd64IEHlJ2drYULF2r58uWqV6+e+vXrpxUrVlxxroiICO3du1czZsxQSkqKEhMT1alTJ/3++++Srv4eS9L58+d1//33a8mSJYqNjdX69ev15JNPasaMGRowYECJfa5fv17vvfeepkyZolWrVqlx48bq37+/fvnlF1ufXbt26fTp07r//vsr9BpJ0g8//KDAwEC5ubnZtQcHB9uev9R9992nc+fOXVfXWOEGU9NLVIAZXemUWTEfHx8jMDDQ9njy5MnGpf/Jrly50pBk7N+/v8w5rnTKrHi+v/3tb2U+dyl/f3/DxcWlxP569epl1K9f33a6rfjYjhw5Ytdv8+bNhiRj8+bNhmEYxpkzZ4z69esbXbt2NYqKiso8hstrOXTokCHJiIqKsuu3a9cuQ5Lxyiuv2Nq6d+9uSDJ27dpl17ddu3ZG7969y9znpcccERFh5OfnG/n5+caRI0eMESNGGJKMF1980TAMw7j77ruNpk2bGmfOnLGNKygoMIKCgoybb77ZdmyXH39WVpYhyYiPjy9z/+V5j5OSkgxJxieffGLXPn36dEOSsWHDBlubJMPHx8fIycmxtWVmZhq1atUypk2bVmLs1U7ZXumU2a233lrqa3zixAlDkjF16lS79ry8PMPFxcV4+eWXr7hPoLqwQgQ4KcMwrvh8x44d5eHhoWeffVb/+Mc/7P6F74iBAweWu+8dd9yhDh062LUNHTpUOTk52rdvn0P73blzp3JychQVFeXQt+c2b94sSbZTf8XuuusuBQYG6quvvrJrb9asme666y67tuDgYLtTXleSnJwsd3d3ubu7q1WrVvrkk0/0/PPP64033tC5c+e0a9cuDRo0SHXr1rWNcXV11bBhw/Sf//xHhw8fLnXexo0bq02bNnr77bc1a9YspaamljiNV573eNOmTapTp44GDRpk1178+lz+etx///2qV6+e7bGPj4+aNm1q93qcOHFCLi4u8vLyKtdrVJYrva+XP+fu7q6GDRvq+PHjldonUFEEIsAJnTt3TqdPn5avr2+Zfdq0aaONGzeqadOmGjNmjNq0aaM2bdro3XffdWhfzZs3L3ffZs2aldl2+vRph/Z76tQpSdLNN9/s0Lji/ZRWt6+vb4k6mjRpUqKfp6en/vjjj3Ltr2vXrtq9e7f27NmjgwcP6vfff9fs2bPl4eGh3377TYZhlFnLpfVezsXFRV999ZV69+6tGTNm6E9/+pO8vb01duxY23U05XmPT58+rWbNmpUIGE2bNpWbm1uFXo8//vhD7u7ucnV1LddrVJomTZqUeuy//vqrpIuB8HIWi6Xc7wtQ1QhEgBNav369CgsLdd99912xX7du3bRu3TplZ2frm2++UVhYmKKjo0u910tZHFmdyczMLLOt+C9ai8UiSSUu6L382p7ia3Cudr3Q5Yr3k5GRUeK5EydOVHpV43INGjRQaGioQkJCFBgYKA8PD9tzjRo1Uq1atcqsRdIV6/H399fChQuVmZmpw4cPKyYmRgkJCXrxxRdtfa72Hjdp0kT//e9/S6wonjx5UgUFBRV6Pby8vJSXl6dz5845PLZY+/btdejQIRUUFNi1f//995KkoKCgEmN+++23Kn//gPIiEAFOJi0tTS+88IIaNGigv/zlL+Ua4+rqqs6dO2vu3LmSZDt95enpKUlV9q/uH3/8UQcOHLBrW7ZsmerVq6c//elPkmS7Ud93331n12/t2rV2j7t06aIGDRooKSnpqqcHL/XAAw9IUomLonfv3q1Dhw6pR48e5Z6rsurUqaPOnTtr9erVdq9xUVGRPvroI91888267bbbyjXXbbfdpokTJ6p9+/alnn4s6z3u0aOHzp49q3/+8592/YtvHlmR16Nt27aSpH//+98Ojy3Wv39/nT17VqtWrbJr/8c//iFfX1917tzZrv3EiRM6f/682rVrV+F9ApXhdvUuAKrLDz/8oIKCAhUUFOjkyZPatm2bPvjgA7m6umrNmjW2VZTSJCUladOmTerbt69atmyp8+fPa9GiRZKknj17SpLq1asnf39/ffbZZ+rRo4caN24sLy+vK95d+Ep8fX318MMP69VXX1Xz5s310UcfKSUlRdOnT1ft2rUlSXfeeaduv/12vfDCCyooKFCjRo20Zs0abd++3W6uunXr6u9//7tGjRqlnj17avTo0fLx8dH//d//6cCBA3rvvfdKreH222/Xs88+qzlz5qhWrVrq06ePjh49qkmTJsnPz08xMTEVOraKmjZtmnr16qX7779fL7zwgjw8PJSQkKAffvhBy5cvL3MF7rvvvtNf//pXPfroo7r11lvl4eGhTZs26bvvvtP48eMlle89Hj58uObOnasRI0bo6NGjat++vbZv366pU6cqIiLC1s8RxSuT33zzje1bYcUOHjyogwcPSrq4Opibm6uVK1dKktq1a2cLNH369FGvXr303HPPKScnR7fccouWL1+uL774Qh999FGJ03HffPONJFXqm21ApdTsNd2AORV/E6t48/DwMJo2bWp0797dmDp1qnHy5MkSYy7/tpXVajX69+9v+Pv7G56enkaTJk2M7t27G2vXrrUbt3HjRqNTp06Gp6enIckYMWKE3XynTp266r4M4//dpHDlypXGHXfcYXh4eBgBAQHGrFmzSoz/6aefjPDwcKN+/fqGt7e38fzzz9tu1Ff8LatiycnJRvfu3Y06deoYtWvXNtq1a2dMnz79irUUFhYa06dPN2677TbD3d3d8PLyMp588kkjPT3drl/37t2NO+64o0R9I0aMuOIdli8/5qvZtm2b8cADDxh16tQxbrrpJuPuu+821q1bZ9fn8m+Z/fe//zWeeuopo23btkadOnWMunXrGsHBwcY777xju2lhed/j06dPG5GRkUbz5s0NNzc3w9/f34iLizPOnz9v10+SMWbMmFKPs/hzUaxbt25GREREib7F70dp2+XfZjxz5owxduxYo1mzZoaHh4cRHBxsLF++vNTXcNiwYUb79u1LfQ64FlwMw4G1agCAKaxatUqDBw/WsWPH1KJFi2rdV05Ojnx9ffXOO+9o9OjR1bovoCwEIgBACYZhqEuXLgoJCSnz9GVVee2117RixQp99913JW7kCFwrXFQNACjBxcVF8+fPl6+vb4V/5qS86tevr8WLFxOGUKNYIQIAAKbHChEAADA9AhEAADA9TtiWQ1FRkU6cOKF69eo5dFdfAABQcwzD0JkzZ+Tr66tata68BkQgKocTJ07Iz8+vpssAAAAVkJ6eftXfTSQQlUPxL0Onp6erfv36NVwNAAAoj5ycHPn5+dn+Hr8SAlE5FJ8mq1+/PoEIAIDrTHkud+GiagAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHpuNV0AYHavftLb8TGPfWn789NrHnR4/Af9v3B4DADcyFghAgAApkcgAgAApkcgAgAApkcgAgAApkcgAgAApkcgAgAApkcgAgAApkcgAgAApkcgAgAApud0gSghIUGtWrWSxWJRSEiItm3bVmbfjIwMDR06VLfffrtq1aql6OjoEn3mz5+vbt26qVGjRmrUqJF69uypb7/9thqPAAAAXG+cKhCtWLFC0dHRmjBhglJTU9WtWzf16dNHaWlppfa/cOGCvL29NWHCBHXo0KHUPlu2bNGQIUO0efNmWa1WtWzZUuHh4Tp+/Hh1HgoAALiOOFUgmjVrlkaOHKlRo0YpMDBQ8fHx8vPzU2JiYqn9AwIC9O6772r48OFq0KBBqX2WLl2qqKgodezYUW3bttX8+fNVVFSkr776qjoPBQAAXEecJhDl5eVp7969Cg8Pt2sPDw/Xzp07q2w/ubm5ys/PV+PGjcvsc+HCBeXk5NhtAADgxuU0gSgrK0uFhYXy8fGxa/fx8VFmZmaV7Wf8+PFq0aKFevbsWWafadOmqUGDBrbNz8+vyvYPAACcj9MEomIuLi52jw3DKNFWUTNmzNDy5cu1evVqWSyWMvvFxcUpOzvbtqWnp1fJ/gEAgHNyq+kCinl5ecnV1bXEatDJkydLrBpVxMyZMzV16lRt3LhRwcHBV+zr6ekpT0/PSu8TAABcH5xmhcjDw0MhISFKSUmxa09JSVGXLl0qNffbb7+t119/XV988YVCQ0MrNRcAALjxOM0KkSTFxsZq2LBhCg0NVVhYmObNm6e0tDRFRkZKungq6/jx41qyZIltzP79+yVJZ8+e1alTp7R//355eHioXbt2ki6eJps0aZKWLVumgIAA2wpU3bp1Vbdu3Wt7gAAAwCk5VSAaPHiwTp8+rSlTpigjI0NBQUFKTk6Wv7+/pIs3Yrz8nkSdOnWy/Xnv3r1atmyZ/P39dfToUUkXb/SYl5enQYMG2Y2bPHmyXn311Wo9HgAAcH1wqkAkSVFRUYqKiir1ucWLF5doMwzjivMVByMAAICyOM01RAAAADWFQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEzP6X7LDMC11+ezMQ6P+d9H5lZDJQBQM1ghAgAApkcgAgAApscpM6ASZi/t7fCYsU98WQ2VAAAqgxUiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgek4XiBISEtSqVStZLBaFhIRo27ZtZfbNyMjQ0KFDdfvtt6tWrVqKjo4utd+qVavUrl07eXp6ql27dlqzZk01VQ8AAK5HThWIVqxYoejoaE2YMEGpqanq1q2b+vTpo7S0tFL7X7hwQd7e3powYYI6dOhQah+r1arBgwdr2LBhOnDggIYNG6bHHntMu3btqs5DAQAA1xGnCkSzZs3SyJEjNWrUKAUGBio+Pl5+fn5KTEwstX9AQIDeffddDR8+XA0aNCi1T3x8vHr16qW4uDi1bdtWcXFx6tGjh+Lj46vxSAAAwPXEaQJRXl6e9u7dq/DwcLv28PBw7dy5s8LzWq3WEnP27t37inNeuHBBOTk5dhsAALhxOU0gysrKUmFhoXx8fOzafXx8lJmZWeF5MzMzHZ5z2rRpatCggW3z8/Or8P4BAIDzc5pAVMzFxcXusWEYJdqqe864uDhlZ2fbtvT09ErtHwAAODe3mi6gmJeXl1xdXUus3Jw8ebLECo8jmjVr5vCcnp6e8vT0rPA+AQDA9cVpVog8PDwUEhKilJQUu/aUlBR16dKlwvOGhYWVmHPDhg2VmhMAANxYnGaFSJJiY2M1bNgwhYaGKiwsTPPmzVNaWpoiIyMlXTyVdfz4cS1ZssQ2Zv/+/ZKks2fP6tSpU9q/f788PDzUrl07SdK4ceN07733avr06XrkkUf02WefaePGjdq+ffs1Pz4AAOCcnCoQDR48WKdPn9aUKVOUkZGhoKAgJScny9/fX9LFGzFefk+iTp062f68d+9eLVu2TP7+/jp69KgkqUuXLvr44481ceJETZo0SW3atNGKFSvUuXPna3ZcAADAuTlVIJKkqKgoRUVFlfrc4sWLS7QZhnHVOQcNGqRBgwZVtjQAAHCDcppriAAAAGoKgQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJie092HCMD1J2LNGxUal9x/YhVXAgAVwwoRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPX7tHqb2wT/CHR7z9IgN1VAJAKAmsUIEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMz+kCUUJCglq1aiWLxaKQkBBt27btiv2//vprhYSEyGKxqHXr1kpKSirRJz4+Xrfffrtuuukm+fn5KSYmRufPn6+uQwAAANcZpwpEK1asUHR0tCZMmKDU1FR169ZNffr0UVpaWqn9jxw5ooiICHXr1k2pqal65ZVXNHbsWK1atcrWZ+nSpRo/frwmT56sQ4cOaeHChVqxYoXi4uKu1WEBAAAn51bTBVxq1qxZGjlypEaNGiXp4srOl19+qcTERE2bNq1E/6SkJLVs2VLx8fGSpMDAQO3Zs0czZ87UwIEDJUlWq1X33HOPhg4dKkkKCAjQkCFD9O23316bgwIAAE7PaVaI8vLytHfvXoWHh9u1h4eHa+fOnaWOsVqtJfr37t1be/bsUX5+viSpa9eu2rt3ry0A/fLLL0pOTlbfvn3LrOXChQvKycmx2wAAwI3LaVaIsrKyVFhYKB8fH7t2Hx8fZWZmljomMzOz1P4FBQXKyspS8+bN9fjjj+vUqVPq2rWrDMNQQUGBnnvuOY0fP77MWqZNm6bXXnut8gcFAACuC06zQlTMxcXF7rFhGCXartb/0vYtW7bozTffVEJCgvbt26fVq1fr888/1+uvv17mnHFxccrOzrZt6enpFT0cAABwHXCaFSIvLy+5urqWWA06efJkiVWgYs2aNSu1v5ubm5o0aSJJmjRpkoYNG2a7Lql9+/Y6d+6cnn32WU2YMEG1apXMhJ6envL09KyKwwIAANcBp1kh8vDwUEhIiFJSUuzaU1JS1KVLl1LHhIWFlei/YcMGhYaGyt3dXZKUm5tbIvS4urrKMAzbahIAADA3pwlEkhQbG6sFCxZo0aJFOnTokGJiYpSWlqbIyEhJF09lDR8+3NY/MjJSx44dU2xsrA4dOqRFixZp4cKFeuGFF2x9+vXrp8TERH388cc6cuSIUlJSNGnSJD388MNydXW95scIAACcj9OcMpOkwYMH6/Tp05oyZYoyMjIUFBSk5ORk+fv7S5IyMjLs7knUqlUrJScnKyYmRnPnzpWvr69mz55t+8q9JE2cOFEuLi6aOHGijh8/Lm9vb/Xr109vvvnmNT8+AADgnJwqEElSVFSUoqKiSn1u8eLFJdq6d++uffv2lTmfm5ubJk+erMmTJ1dViQAA4AbjVKfMAAAAagKBCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmJ7T/do9AHPqu3q2w2PWDxhbDZUAMCNWiAAAgOlVaoUoPz9fmZmZys3Nlbe3txo3blxVdQEAAFwzDq8QnT17Vu+//77uu+8+NWjQQAEBAWrXrp28vb3l7++v0aNHa/fu3dVRKwAAQLVwKBC98847CggI0Pz58/XAAw9o9erV2r9/vw4fPiyr1arJkyeroKBAvXr10oMPPqiff/65uuoGAACoMg6dMtu5c6c2b96s9u3bl/r8XXfdpWeeeUaJiYlatGiRvv76a916661VUigAAEB1cSgQffrpp7Y/p6eny8/Pr9R+FotFUVFRlasMAADgGqnwRdX+/v5q1KiROnTooA4dOqhjx47q0KGDLly4oLlz52rJkiVVWScAAEC1qXAg+uWXX7R//37t379fqampWrlypU6cOCFJql+/fpUVCAAAUN0qHIgCAgIUEBCgP//5z7Y2q9WqESNGaPr06VVRGwAAwDVRpTdmDAsL07vvvqs33nijKqcFAACoVhUORPn5+aW233rrrfrxxx8rXBAAAMC1VuFTZnXq1FG7du3UqVMndezYUZ06dZKvr6/mzJmj8PDwqqwRAACgWlU4EG3atEkHDhzQgQMHtHTpUr3yyiv6448/JEnh4eGaMGGCgoODFRwcrMDAwCorGAAAoKpVOBB17dpVXbt2tT0uKirS4cOHbd8827t3rxYtWqSTJ0+qsLCwSooFAACoDpX6cddL1apVS4GBgQoMDNSQIUNs7f/973+rahcAAADVokq/ZVYaHx+f6t4FAABApTi0QtSqVSu5uLg4vJPo6GiNHTvW4XEAAADXgkOBaPHixRXaSUBAQIXGAQAAXAsOBaLu3btXVx0AAAA1psLXEG3cuLHM595///2KTgsAAHDNVTgQ9e3bV//zP/+jvLw8W9upU6fUr18/xcXFVUlxAAAA10KFA9HWrVu1bt063Xnnnfrxxx+1fv16BQUF6ezZszpw4EBV1ggAAFCtKhyIOnfurNTUVAUHByskJET9+/fX//zP/2jTpk3y8/OryhoBAACqVaXuQ3T48GHt3r1bN998s9zc3PSvf/1Lubm5VVUbAADANVHhQPTWW28pLCxMvXr10g8//KDdu3fbVoysVmtV1ggAAFCtKhyI3n33Xf3zn//UnDlzZLFYdMcdd+jbb7/VgAEDdN9991VhiQAAANWrwr9l9v3338vLy8uuzd3dXW+//bYeeuihShcGAABwrVQ4EF0ehi7FDRxxLaz84EGHxwx6+otqqAQAcL1z6JRZWlqaQ5MfP37cof4AAAA1waFAdOedd2r06NH69ttvy+yTnZ2t+fPnKygoSKtXr650gQAAANXNoVNmhw4d0tSpU/Xggw/K3d1doaGh8vX1lcVi0W+//aaDBw/qxx9/VGhoqN5++2316dOnuuoGAACoMg6tEDVu3FgzZ87UiRMnlJiYqNtuu01ZWVn6+eefJUlPPPGE9u7dqx07dlQ4DCUkJKhVq1ayWCwKCQnRtm3brtj/66+/VkhIiCwWi1q3bq2kpKQSfX7//XeNGTNGzZs3l8ViUWBgoJKTkytUHwAAuPFU6KJqi8WiAQMGaMCAAZIkwzAkSS4uLpUqZsWKFYqOjlZCQoLuuecevf/+++rTp48OHjyoli1bluh/5MgRRUREaPTo0froo4+0Y8cORUVFydvbWwMHDpQk5eXlqVevXmratKlWrlypm2++Wenp6apXr16lagUAADeOSt2peuHChQoKCpLFYpHFYlFQUJAWLFhQ4flmzZqlkSNHatSoUQoMDFR8fLz8/PyUmJhYav+kpCS1bNlS8fHxCgwM1KhRo/TMM89o5syZtj6LFi3Sr7/+qn/+85+655575O/vr65du6pDhw4VrhMAANxYKhyIJk2apHHjxqlfv3769NNP9emnn6pfv36KiYnRxIkTHZ4vLy9Pe/fuVXh4uF17eHi4du7cWeoYq9Vaon/v3r21Z88e5efnS5LWrl2rsLAwjRkzRj4+PgoKCtLUqVNVWFhYZi0XLlxQTk6O3QYAAG5cFb4PUWJioubPn68hQ4bY2h5++GEFBwfr+eef1xtvvOHQfFlZWSosLJSPj49du4+PjzIzM0sdk5mZWWr/goICZWVlqXnz5vrll1+0adMmPfHEE0pOTtbPP/+sMWPGqKCgQH/7299KnXfatGl67bXXHKofAABcvyq8QlRYWKjQ0NAS7SEhISooKKhwQZdfh2QYxhWvTSqt/6XtRUVFatq0qebNm6eQkBA9/vjjmjBhQpmn4SQpLi5O2dnZti09Pb2ihwMAAK4DFQ5ETz75ZKmhYt68eXriiSccns/Ly0uurq4lVoNOnjxZYhWoWLNmzUrt7+bmpiZNmkiSmjdvrttuu02urq62PoGBgcrMzFReXl6p83p6eqp+/fp2GwAAuHFV+JSZdPGi6g0bNujuu++WJH3zzTdKT0/X8OHDFRsba+s3a9asq87l4eGhkJAQpaSkqH///rb2lJQUPfLII6WOCQsL07p16+zaNmzYoNDQULm7u0uS7rnnHi1btkxFRUWqVeti/vvpp5/UvHlzeXh4OHbAAADghlThQPTDDz/oT3/6kyTp3//+tyTJ29tb3t7e+uGHH2z9HPkqfmxsrIYNG6bQ0FCFhYVp3rx5SktLU2RkpKSLp7KOHz+uJUuWSJIiIyP13nvvKTY2VqNHj5bVatXChQu1fPly25zPPfec5syZo3Hjxun555/Xzz//rKlTp2rs2LEVPXQAAHCDqXAg2rx5c1XWIUkaPHiwTp8+rSlTpigjI0NBQUFKTk6Wv7+/JCkjI8Pu99RatWql5ORkxcTEaO7cufL19dXs2bNt9yCSJD8/P23YsEExMTEKDg5WixYtNG7cOL388stVXj8AALg+VeqUWXWIiopSVFRUqc8tXry4RFv37t21b9++K84ZFhamb775pirKAwAAN6BK3ZgRAADgRkAgAgAApkcgAgAApkcgAgAApkcgAgAApkcgAgAApkcgAgAApkcgAgAApkcgAgAApkcgAgAApkcgAgAApkcgAgAApkcgAgAApkcgAgAApkcgAgAApudW0wUAQFXou2pBhcatHziqiisBcD1ihQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJgegQgAAJie0wWihIQEtWrVShaLRSEhIdq2bdsV+3/99dcKCQmRxWJR69atlZSUVGbfjz/+WC4uLvrzn/9cxVUDAIDrmVMFohUrVig6OloTJkxQamqqunXrpj59+igtLa3U/keOHFFERIS6deum1NRUvfLKKxo7dqxWrVpVou+xY8f0wgsvqFu3btV9GAAA4DrjVIFo1qxZGjlypEaNGqXAwEDFx8fLz89PiYmJpfZPSkpSy5YtFR8fr8DAQI0aNUrPPPOMZs6cadevsLBQTzzxhF577TW1bt36WhwKAAC4jjhNIMrLy9PevXsVHh5u1x4eHq6dO3eWOsZqtZbo37t3b+3Zs0f5+fm2tilTpsjb21sjR44sVy0XLlxQTk6O3QYAAG5cThOIsrKyVFhYKB8fH7t2Hx8fZWZmljomMzOz1P4FBQXKysqSJO3YsUMLFy7U/Pnzy13LtGnT1KBBA9vm5+fn4NEAAIDridMEomIuLi52jw3DKNF2tf7F7WfOnNGTTz6p+fPny8vLq9w1xMXFKTs727alp6c7cAQAAOB641bTBRTz8vKSq6tridWgkydPllgFKtasWbNS+7u5ualJkyb68ccfdfToUfXr18/2fFFRkSTJzc1Nhw8fVps2bUrM6+npKU9Pz8oeEgAAuE44zQqRh4eHQkJClJKSYteekpKiLl26lDomLCysRP8NGzYoNDRU7u7uatu2rb7//nvt37/ftj388MO6//77tX//fk6FAQAASU60QiRJsbGxGjZsmEJDQxUWFqZ58+YpLS1NkZGRki6eyjp+/LiWLFkiSYqMjNR7772n2NhYjR49WlarVQsXLtTy5cslSRaLRUFBQXb7aNiwoSSVaAcAAOblVIFo8ODBOn36tKZMmaKMjAwFBQUpOTlZ/v7+kqSMjAy7exK1atVKycnJiomJ0dy5c+Xr66vZs2dr4MCBNXUIAADgOuRUgUiSoqKiFBUVVepzixcvLtHWvXt37du3r9zzlzYHAAAwN6cLRDCPLxdGODym98jkaqgEAGB2TnNRNQAAQE0hEAEAANMjEAEAANMjEAEAANMjEAEAANMjEAEAANMjEAEAANPjPkQA8P97aOVSh8d8PuiJaqgEwLXGChEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9t5ouAABuFP1Wrq7QuHWDBlRxJQAc5XQrRAkJCWrVqpUsFotCQkK0bdu2K/b/+uuvFRISIovFotatWyspKcnu+fnz56tbt25q1KiRGjVqpJ49e+rbb7+tzkMAAADXGacKRCtWrFB0dLQmTJig1NRUdevWTX369FFaWlqp/Y8cOaKIiAh169ZNqampeuWVVzR27FitWrXK1mfLli0aMmSINm/eLKvVqpYtWyo8PFzHjx+/VocFAACcnFMFolmzZmnkyJEaNWqUAgMDFR8fLz8/PyUmJpbaPykpSS1btlR8fLwCAwM1atQoPfPMM5o5c6atz9KlSxUVFaWOHTuqbdu2mj9/voqKivTVV19dq8MCAABOzmkCUV5envbu3avw8HC79vDwcO3cubPUMVartUT/3r17a8+ePcrPzy91TG5urvLz89W4ceMya7lw4YJycnLsNgAAcONymkCUlZWlwsJC+fj42LX7+PgoMzOz1DGZmZml9i8oKFBWVlapY8aPH68WLVqoZ8+eZdYybdo0NWjQwLb5+fk5eDQAAOB64jSBqJiLi4vdY8MwSrRdrX9p7ZI0Y8YMLV++XKtXr5bFYilzzri4OGVnZ9u29PR0Rw4BAABcZ5zma/deXl5ydXUtsRp08uTJEqtAxZo1a1Zqfzc3NzVp0sSufebMmZo6dao2btyo4ODgK9bi6ekpT0/PChwFAAC4HjlNIPLw8FBISIhSUlLUv39/W3tKSooeeeSRUseEhYVp3bp1dm0bNmxQaGio3N3dbW1vv/223njjDX355ZcKDQ2tngMwGeu8hxweE/bs59VQCQAAledUp8xiY2O1YMECLVq0SIcOHVJMTIzS0tIUGRkp6eKprOHDh9v6R0ZG6tixY4qNjdWhQ4e0aNEiLVy4UC+88IKtz4wZMzRx4kQtWrRIAQEByszMVGZmps6ePXvNjw8AADgnp1khkqTBgwfr9OnTmjJlijIyMhQUFKTk5GT5+/tLkjIyMuzuSdSqVSslJycrJiZGc+fOla+vr2bPnq2BAwfa+iQkJCgvL0+DBg2y29fkyZP16quvXpPjAgAAzs2pApEkRUVFKSoqqtTnFi9eXKKte/fu2rdvX5nzHT16tIoqAwAANyqnOmUGAABQEwhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9Jzu1+4BwMweWfmFw2M+G/RgNVQCmAsrRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPTcaroA1Ix/zX3E4TFtx3xWDZUAAFDzCEQAcIPpv2q7w2PWDOxaDZUA1w9OmQEAANMjEAEAANMjEAEAANMjEAEAANMjEAEAANMjEAEAANMjEAEAANMjEAEAANPjxowAADuPrvrO4TGfDgyuhkqAa4cVIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHpcVH0dOjE31uExvmNmVUMlAADcGAhEAIAqN3ZNusNjZvf3q4ZKgPJxulNmCQkJatWqlSwWi0JCQrRt27Yr9v/6668VEhIii8Wi1q1bKykpqUSfVatWqV27dvL09FS7du20Zs2a6iofAABch5xqhWjFihWKjo5WQkKC7rnnHr3//vvq06ePDh48qJYtW5bof+TIEUVERGj06NH66KOPtGPHDkVFRcnb21sDBw6UJFmtVg0ePFivv/66+vfvrzVr1uixxx7T9u3b1blz52t9iACAcpi75r8OjxnT36caKoFZONUK0axZszRy5EiNGjVKgYGBio+Pl5+fnxITE0vtn5SUpJYtWyo+Pl6BgYEaNWqUnnnmGc2cOdPWJz4+Xr169VJcXJzatm2ruLg49ejRQ/Hx8dfoqAAAgLNzmhWivLw87d27V+PHj7drDw8P186dO0sdY7VaFR4ebtfWu3dvLVy4UPn5+XJ3d5fValVMTEyJPlcKRBcuXNCFCxdsj7OzsyVJOTk5jhxSqU4tLHlKrzy8R0ba/nzmjwtX6Fm6y2s/+0d+peY4V8nxVTFH7h8Fla7hj0rOcT638jVcqOQceVVQQ0FuXqXmyM897/D4qpjDfvwfVVBD7jUfX3KOc5UaXxVz5OeerXQNeblnKjXHHxUaf5Pd47X/PO3wHA//uYntz1tXZDk8/t7BXg6PQfUp/kwZhnH1zoaTOH78uCHJ2LFjh137m2++adx2222ljrn11luNN998065tx44dhiTjxIkThmEYhru7u7F06VK7PkuXLjU8PDzKrGXy5MmGJDY2NjY2NrYbYEtPT79qDnGaFaJiLi4udo8NwyjRdrX+l7c7OmdcXJxiY//fV9uLior066+/qkmTJmWOy8nJkZ+fn9LT01W/fv0y5y5LZcdTAzVQg/POQQ3UQA01U4NhGDpz5ox8fX2vOpfTBCIvLy+5uroqMzPTrv3kyZPy8Sn9QrlmzZqV2t/NzU1NmjS5Yp+y5pQkT09PeXp62rU1bNiwXMdRv379Cr+xVTGeGqiBGpx3DmqgBmq49jU0aNCgXHM4zUXVHh4eCgkJUUpKil17SkqKunTpUuqYsLCwEv03bNig0NBQubu7X7FPWXMCAADzcZoVIkmKjY3VsGHDFBoaqrCwMM2bN09paWmKjLx4QXFcXJyOHz+uJUuWSJIiIyP13nvvKTY2VqNHj5bVatXChQu1fPly25zjxo3Tvffeq+nTp+uRRx7RZ599po0bN2r79u01cowAAMD5OFUgGjx4sE6fPq0pU6YoIyNDQUFBSk5Olr+/vyQpIyNDaWlptv6tWrVScnKyYmJiNHfuXPn6+mr27Nm2exBJUpcuXfTxxx9r4sSJmjRpktq0aaMVK1ZU+T2IPD09NXny5BKn2q7VeGqgBmpw3jmogRqowTlruJSLYZTnu2gAAAA3Lqe5hggAAKCmEIgAAIDpEYgAAIDpEYgAAIDpEYgAAIDpEYgq4eOPP9aQIUP0xBNPaOjQoXb3P6qomTNnOtR/06ZNeu6557R//35J0rx58xze5+eff66XX35Z27Zt06OPPqpVq1Y5NH7Dhg3asGGDvvzyS/Xr108bNmxwuIatW7fq6NGjeuqppzR06FBt3brVofH/+Mc/tHHjRg0cOFCPP/645s6d63ANXbt21fz583X2rOM/bClJM2bM0FNPPaXFixfrscce00svveTwHDt27NDgwYPVt29fDRs2TD/99FOFatm5c6dWrFihHTt2VGj8pU6dOlXpOS79sWRHZWRkVHr/u3btqvDYM2fOVPg1yMu7+KO5mzdv1rp162yPy+vnn3+u0H4vlZubq8OHD6uoqEhr166t0Ot54sQJ7d69W1lZjv/YaTE+k/b4TFbuMylVzefyUk51H6LrzZYtW+xC0JgxYzRkyBCH5hg6dKjtz4ZhaO/evXrhhRfKPf69997TokWL9Oabb+q3336zBSNHLFiwQIsXL1Z4eLisVqueffZZu3s5Xc0rr7yixx9/XN7e3jp37lyJn0opj+XLlysvL08zZsxQw4YNNWLECN17773lHr9nzx7t2bPHFubGjRvncA2BgYHy8fHRs88+q/r16+upp57S3XffXe7xP/30kxYvXqyHHnpIn3/+uZ577jmHa/jwww+1fPlyTZo0SePHj9dzzz2njz76yKE5xowZozvuuEMtW7bUgQMHtHTpUiUkJJR7/KX/B2kYhiZOnKj333/foRpeeuklpaenq23btpo8ebJiYmIcqmHMmDHKz8/XTTfdpMLCQrm4uGjOnDkO1TBkyBC5uLjY/rsKDQ3VsmXLyj1++vTp8vPz07p169SwYUM1btxYb775pkM1vPzyy2rYsKG8vLzUuHFjRUZGatGiReUeHxERofbt2+vhhx/W4MGDddNNN1190GWefvppderUSd9++60GDRqk559/XitXriz3+LfeekuZmZlKT09Xo0aNFBoaartZbnnxmbyIz+RFlf1MSlXzubwcK0SVkJ+fry+//FIHDx7Uhg0bdP78eYfnqF27tpYtW6Zly5Zp+fLl6tGjh0Pjvby81LBhQ7399tvauHFjhf7V0bx5czVs2FDPP/+8XF1dVbt2bYfGb926VWfOnFHt2rXVtm1bDR8+3OEaDh48qMzMTDVt2lQeHh4O/65NgwYNdPbsWc2bN08rV67UuXPnHK7Bzc1NDz/8sJYtW6bJkydr8+bNDo3PysrS4sWLlZeXp127dunXX391uIbff/9dJ0+e1K+//qp69eqpXr16Ds9Rq1YtRUVF6aGHHlJUVJTDNywLDg5WZGSkbXN0tU66+HuBxZ/nl156SY7e7uzcuXOaN2+e/vOf/+i9995zeLwkdezYUREREbY6HPmLR5KOHj2qTZs2afny5UpMTNTp06cdrqGgoEDZ2dm2fyw5+rnu2bOnPvnkE9WvX19PP/20xowZ43ANDRs21Pjx4yVd/AfYlX7HsTTp6emKj49XQECAFixYoNTUVIdrqOxnsn379nwmxWfyUlXxubwcK0SVMGfOHK1evVqpqany8/Nz+F8LkjRhwgS7x46m/X79+tmNvfnmmx2uYdiwYXb/62goq127tl577TWlpKSU+0f0LjdlyhS5uLjYHj/44IMOjX/99df12Wef6fDhw6pbt26F3osRI0bY/ty8eXPFxcU5NH7OnDk6fPiwPv74Y33wwQf629/+5nANUVFRmjp1qm2VcOTIkQ7Pccstt+iZZ56Rl5eXsrKyFBwc7ND4QYMG6Y033rA9fueddxyuofi97Nq1qwzD0GOPPabExMRyjy/+C3Pq1KkO77vYyy+/rG+++Ubjxo1Tdna2w+P37dtnF0hzc3MdniM8PFyJiYm655571KJFC4dWPYu5ublpwIABGjBggI4fP+7weIvFoqFDh+qOO+5QZGSkw3+RZ2dn66233lKdOnUkSa6urg7XUNnP5KOPPspnUlJqaqrq1q1re2zWz6RUNZ/LEgwAN5y8vDwjMzPTyMvLM95+++1KzVWR8Tk5OXaPZ8yYUakaXn/99QqPzc3NNbZu3Vrp16Gix1BYWFij78WlsrKyHD6O3NxcY9++fZWq4auvvjKefvppY+fOnUZeXp7x/vvvOzw+MjLS2L9/v2EYhsPjDcMwPvvsMyMyMtJITU01DMMw3njjjQrVUDz+rbfecriGTZs2GZGRkYbVajW2bt1aodfhL3/5i62GadOmOVzDV199ZYwePdqwWq2Vei+Ka6jIe3HpcWRlZRmJiYkOz7F+/XqjV69etjqSkpIcnuNyrBABN5jKXpdW2fGS9Je//KXEHC+++GKlapg4caJDNVTH6+DIMVRXDY6+F5U9jktXKStaw6XXOp4/f97hax0vHf/rr79W6FrJxYsX211v6ejKxuXXax47dszhGubMmWObIyIiolKvw2+//Wb3254VmeOPP/6odA0VeS8unSMiIkLfffedw3MsWLBAn3zyia2OAwcOODzH5QhEwA2mdu3aWrBgge2xoxd3V3Y8NVDD5S691nHChAkOX+tY2fHUcGPVUFVzlFDpNSYATuWXX36xe3z69OlrOp4aqOFya9eutXuckJBwTcdTw41VQ1XNcTl+7R4AAJgeX7sHAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACYEqnTp1Ss2bN7H6bateuXfLw8NCGDRtqsDIANYH7EAEwreTkZP35z3/Wzp071bZtW3Xq1El9+/ZVfHx8TZcG4BojEAEwtTFjxmjjxo268847deDAAe3evVsWi6WmywJwjRGIAJjaH3/8oaCgIKWnp2vPnj0KDg6u6ZIA1ACuIQJgar/88otOnDihoqKiCv2COYAbAytEAEwrLy9Pd911lzp27Ki2bdtq1qxZ+v777+Xj41PTpQG4xghEAEzrxRdf1MqVK3XgwAHVrVtX999/v+rVq6fPP/+8pksDcI1xygyAKW3ZskXx8fH68MMPVb9+fdWqVUsffvihtm/frsTExJouD8A1xgoRAAAwPVaIAACA6RGIAACA6RGIAACA6RGIAACA6RGIAACA6RGIAACA6RGIAACA6RGIAACA6RGIAACA6RGIAACA6RGIAACA6RGIAACA6f1/Oed/mTcGDGIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import poisson \n",
    "\n",
    "mu = 10\n",
    "pois = poisson(mu)\n",
    "import pandas as pd\n",
    "x = np.arange(30)\n",
    "prob_p = pois.pmf(x)\n",
    "data = pd.DataFrame(zip(x,prob_p))\n",
    "data.columns=['x','p']\n",
    "ax = sns.barplot(data=data, x='x', y='p')\n",
    "ax.set(xlabel='x', ylabel='p(x|$\\theta$)', title='Distribución Poisson({})'.format(mu))\n",
    "ax.set_xticklabels(x, rotation=90, size=5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Con el siguiente código calculamos la media varianza y algunos cuantiles de la distribución $\\text{Pois}(10)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Distribución Poisson(10); media = 10.0, varianza 10.0\n",
      "Distribución Poisson(10); probabilidades=[0.25 0.5  0.75]; cuantiles=[ 8. 10. 12.]\n"
     ]
    }
   ],
   "source": [
    "# media y varianza de la distribución Poisson(17.5)\n",
    "import numpy as np\n",
    "from scipy.stats import poisson \n",
    "\n",
    "mu = 10\n",
    "pois = poisson(mu)\n",
    "media, var = pois.stats(moments=\"mv\")\n",
    "print('Distribución Poisson(10); media = {}, varianza {}'.format(media, var))\n",
    "\n",
    "# calcula algunos cuantiles de la distribución\n",
    "prob = np.array([0.25, 0.5, 0.75])\n",
    "cuantiles = pois.ppf(prob)\n",
    "print('Distribución Poisson(10); probabilidades={}; cuantiles={}'.format(prob, cuantiles))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:#4CC9F0\">Distribución Binomial</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La distribución Binomial surge en problemas con ensayos binarios que pueden resultar en éxito (1) o fallo (0). Si se realizan $n$ ensayos (de Bernoulli) de forma independiente y se define la variable aleatoria $y$ como el  número de éxitos, entonces se dice que $Y$ tiene distribución Binomial con parámetro $\\theta$, cuya densidad es dada por:\n",
    "\n",
    "$$\n",
    "p(y|\\theta) = \\binom{n}{y}\\theta^y(1-\\theta)^{n-y}, \\hspace{1mm} y=1,\\ldots,n,\n",
    "$$\n",
    "\n",
    "en donde \n",
    "\n",
    "$$\n",
    "\\binom{n}{y} =\\frac{n!}{y!(n-y)!}\n",
    "$$\n",
    "\n",
    "es el número de grupos diferentes de tamaño $y$, que se pueden organizar si se tienen en total $n$ elementos. Esta distribución se denota $\\text{Binom}(n,\\theta)$. La densidad se puede reescribir como\n",
    "\n",
    "$$\n",
    "p(y|\\theta) = \\exp\\left[ y \\log \\theta - y\\log(1-\\theta) + n\\log  (1-\\theta) + \\log \\binom{n}{y} \\right],\n",
    "$$\n",
    "\n",
    "de donde se tiene que:\n",
    "\n",
    "* la densidad está en forma canónica: $a(y) =y$,\n",
    "* parámetro natural: $b(\\theta) = \\log \\theta - log (1-\\theta) = \\log\\left[\\frac{\\theta}{1-\\theta} \\right]$,\n",
    "* $c(\\theta) = n \\log(1-\\theta)$,\n",
    "* $d(y) = \\log \\binom{n}{y}$.\n",
    "\n",
    "Se puede verificar que si una variable aleatoria $Y$ tiene distribución $\\text{Binom}(n,\\theta)$, entonces se tiene que\n",
    "\n",
    "* $E[Y] = n\\theta$.\n",
    "* $Var[Y] = n\\theta(1-\\theta)$. \n",
    "\n",
    "El siguiente fragmento de código muestra la función de probabilidad de la distribución $\\text{Binom}(4,0.25)$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAHGCAYAAACVcJQUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABIrElEQVR4nO3deVxVdeL/8feV7bpBrqCpQKYC4YobGJqZuKdZSZZojZaOViJTX4csU2vCpVG0FLOxGJtEMjWtNKUyl8RKA7NlHGfKcBRyKUVNQeD8/nC4P6/3goAiV8/r+Xicx8P7uZ/zOZ/PPXd58zmLFsMwDAEAAJhYtaruAAAAQFUjEAEAANMjEAEAANMjEAEAANMjEAEAANMjEAEAANMjEAEAANMjEAEAANMjEAEAANMjEJnQvffeq9atW+vkyZNV3RUAN7jz58+ra9euuuuuu3T+/Pmq7g5QIgJRBSQnJ8tisThdnnrqqSrt28MPP6yAgIASn09MTNSXX36pDRs2yMfHp9L7c+DAAVksFiUnJ1f6tipDQECABg4ceFXbtFgsevzxxy9b77PPPpPFYtFnn31mK5s2bZosFotdvTvuuEN33HHHZdurjH1R/Fk4cOBAhda/3Pu1vJy9PpXB2bjLuh8u9v3332vatGnlfv0u3Vbxvn355ZfL1c7VcLlxP/300zp79qxWr14tDw+Pq7bd4jFfvHh7e6tt27ZKTExUYWFhufrpCq7kM1raujNmzFBISIiKioqcrvvLL7+oXr16slgsevfdd8u8zRUrVqhdu3ayWq1q3LixYmNjdfr06cuu969//UtPPfWUwsLCdNNNN6lu3brq1q2b022X9nubk5Njq3f+/Hk1b95ciYmJZe7/pdwrvCb05ptvKigoyK6scePGVdSbC5577jlNnDjR6XM7d+7UX/7yF23evFlNmjS5xj1DeXXo0EHp6ekKCQkptd6iRYuuUY9Qmorsh++//17Tp0/XHXfcUa5g6Er7vLS+rF69WmvWrFF6erq8vb0rZftPPPGEHnzwQUnSiRMntG7dOk2aNEkHDx7UX//61zL101U0atRI6enpat68+VVr8/Dhw5o9e7aSk5NVrZrzOZAJEybIarWWq923335bI0aM0JgxYzRv3jz961//0uTJk/X9999r06ZNpa67adMmffjhh4qJiVGnTp1UUFCg1NRU3X///Zo+fbqmTp3qsI6z39t69erZ/u3h4aGpU6dq0qRJiomJsXuurAhEVyA0NFQdO3as6m7YKe2D1LVrVx09evQa9sb1/P7776pRo0ZVd6NMvL291bVr18vWu1xgwrVxLfZD8fvXlfZ5aX0ZOnSohg4dWqnbb9asmd3npG/fvvr222+VkpJiF4hc6TUriZeXV5k+8+Uxf/583XTTTSXuh1WrVmnjxo1auHChRo0aVaY2CwsL9fTTTysqKkqvv/66JKlnz56qXbu2HnroIW3YsEH9+vUrcf0HHnhAEyZMsJvN7devn44dO6ZZs2Zp8uTJ8vLyslunLL+3w4cPV1xcnF577TU988wzZRrLxThkVkksFoumTZvmUB4QEKCHH37Y9rh4OnDz5s364x//qPr166tevXoaOnSoDh8+7LD+8uXLFR4erlq1aqlWrVpq166dli5danve2SGIc+fOKT4+XoGBgfL09NTNN9+sCRMm6MSJEw59GzhwoD766CN16NBB1atXV1BQkN54440yjfnw4cMaNmyYateuLR8fH0VHR9tNaV5s165duvvuu1W3bl1ZrVa1b99e77zzzmW3UTwtPHv2bP3lL39Rs2bNZLVa1bFjR33yySd2dYsPn3z99de67777VKdOHVtgLOtrUmzNmjVq06aNrFarbrnlFi1YsMDu+XPnzulPf/qT2rVrJx8fH9WtW1fh4eFau3ZtiWN57bXX1LJlS3l5eSkkJEQrVqywe97ZITNnnB0KKOu+2LVrlx544AEFBASoevXqCggI0PDhw/Xzzz871N25c6e6detmmx6Pj48v1zkhycnJatWqlby8vBQcHKxly5Y5rZefn68XX3xRQUFB8vLyUoMGDfTII49UOMynpqYqKipKjRo1UvXq1RUcHKw///nPOnPmTJnWL+u4ne2HpKQktW3bVrVq1VLt2rUVFBRk+6JOTk7W/fffL+nCj0nxYYDiQx533HGHQkNDtXXrVkVERKhGjRr6wx/+UOK2JKmoqOiyn4uSDlM6O9xYVFSkV155Re3atVP16tV10003qWvXrlq3bl2p4/711181fvx43XzzzfL09NQtt9yiKVOmKC8vz65e8eHjt956S8HBwapRo4batm2rDz74wKF/5eHj4+NweK60w4xz585VYGCgatWqpfDwcO3cudOhzXXr1ik8PFw1atRQ7dq11bt3b6Wnp9vVKX4Nv/nmG91///2274K4uDgVFBRo37596tu3r2rXrq2AgADNnj3bbn1nh73+/e9/65FHHlGLFi1Uo0YN3XzzzRo0aJD27t172dchPz9fS5cu1YMPPuh0dujXX3/VhAkTbO+Zstq5c6eys7P1yCOP2JXff//9qlWrltasWVPq+vXr13d6aLtz5876/fff9euvv5a5Lxfz9PRUdHS0lixZIsMwyr0+gegKFBYWqqCgwG6pqDFjxsjDw0PLly/X7Nmz9dlnn2nEiBF2daZOnaqHHnpIjRs3VnJystasWaNRo0Y5/fEqZhiGhgwZopdfflkxMTH68MMPFRcXp7///e+68847Hb6g9uzZoz/96U+aNGmS1q5dqzZt2mj06NHaunVrqf0/e/as7rrrLm3atEkJCQlauXKl/Pz8FB0d7VB38+bN6tatm06cOKHFixdr7dq1ateunaKjo8t87PzVV1/VRx99pMTERP3jH/9QtWrV1K9fP4cvKOnCX6m33nqrVq5cqcWLF5f7NcnMzFRsbKwmTZqkNWvWKCIiQhMnTrQ7XyMvL0+//vqrnnrqKb333ntKSUnR7bffrqFDhzr94V+3bp0WLFigGTNm6N1335W/v7+GDx9eruP3JSnPvjhw4IBatWqlxMREbdy4UbNmzVJ2drY6deqkY8eO2ep9//336tWrl06cOKHk5GQtXrxYGRkZevHFF8vUp+TkZD3yyCMKDg7WqlWr9Oyzz+qFF17Qp59+alevqKhIgwcP1syZM/Xggw/qww8/1MyZM5WWlqY77rhDZ8+eLffrsX//fvXv319Lly7VRx99pNjYWL3zzjsaNGjQZde9knGvWLFC48ePV48ePbRmzRq99957mjRpki2IDRgwQC+99JIkaeHChUpPT1d6eroGDBhgayM7O1sjRozQgw8+qPXr12v8+PGlbrM8n4uyePjhhzVx4kR16tRJqampWrFihe6+++5Sz3k6d+6cevbsqWXLlikuLk4ffvihRowYodmzZzudpfjwww/16quvasaMGVq1apXq1q2re+65Rz/++GOZ+lhUVGT7/j1+/LjeeOMNffTRR4qJiSnT+gsXLlRaWpoSExP19ttv68yZM+rfv7/dRSfLly/X4MGD5e3trZSUFC1dulS//fab7rjjDm3fvt2hzWHDhqlt27ZatWqVHn30Uc2bN0+TJk3SkCFDNGDAAK1Zs0Z33nmnJk+erNWrV5fav8OHD6tevXqaOXOmPvroIy1cuFDu7u7q0qWL9u3bV+q6X3zxhY4fP66ePXs6ff7JJ59UYGBgmc5pvNi3334rSWrTpo1duYeHh4KCgmzPl9fmzZvVoEEDNWzY0OG5gQMHys3NTXXr1tXQoUNL3MYdd9yhn3/+uWJ9MFBub775piHJ6XL+/HnDMAxDkvH88887rOvv72+MGjXKoa3x48fb1Zs9e7YhycjOzjYMwzB+/PFHw83NzXjooYdK7duoUaMMf39/2+OPPvrIkGTMnj3brl5qaqohyViyZIld36xWq/Hzzz/bys6ePWvUrVvXGDt2bKnbTUpKMiQZa9eutSt/9NFHDUnGm2++aSsLCgoy2rdvb3utig0cONBo1KiRUVhYWOJ2fvrpJ0OS0bhxY+Ps2bO28tzcXKNu3brGXXfdZSt7/vnnDUnG1KlT7doo72tisViMzMxMu7q9e/c2vL29jTNnzjjtZ0FBgXH+/Hlj9OjRRvv27e2ek2RUr17dyMnJsasfFBRk3HrrrbayzZs3G5KMzZs3O4zpYj169DB69Ohhe1yefeGs36dPnzZq1qxpzJ8/31YeHR1dYp8lGT/99FOJbRYWFhqNGzc2OnToYBQVFdnKDxw4YHh4eNi9X1NSUgxJxqpVq+za+OqrrwxJxqJFi0rcjmE4f30uVlRUZJw/f97YsmWLIcnYs2dPqe2VZ9yX7ofHH3/cuOmmm0ptf+XKlQ77+OL2JBmffPKJ0+cu3lZ5PheXfkcUu/S127p1qyHJmDJlSqljuLQvixcvNiQZ77zzjl29WbNmGZKMTZs22cokGb6+vkZubq6tLCcnx6hWrZqRkJBQ6naLx+xsefjhh42CgoJS+1m8fuvWre3qfvnll4YkIyUlxTCM///+bd26td1306lTp4yGDRsaERERtrLi1/Cvf/2r3bbbtWtnSDJWr15tKzt//rzRoEEDY+jQoQ59utxnND8/32jRooUxadKkUtctfs0vfv8W++CDDwwPDw9j7969hmH8/++blStXlrjtYn/5y1/sfp8uFhUVZbRs2fKybVzq9ddfNyTZfe8YhmFs2LDBmDJlivH+++8bW7ZsMV599VWjSZMmRs2aNR2+lw3DMPbv329IMpKSksrdB2aIrsCyZcv01Vdf2S3u7hU7Levuu++2e1ycvItnf9LS0lRYWKgJEyaUq93iv8AvPkwnXZjarFmzpsN0ert27eymTq1Wq1q2bFnqLJR0IdnXrl3bYRzFJzsW+/e//61//vOfeuihhyTJbnatf//+ys7OvuxfPdKFWZ+LTwKsXbu2Bg0apK1btzpcXXLvvffaPS7va3Lbbbepbdu2DuPKzc3V119/bStbuXKlunXrplq1asnd3V0eHh5aunSpfvjhB4f+9+rVS76+vrbHbm5uio6O1r///W/997//vez4S1PWfSFJp0+f1uTJk3XrrbfK3d1d7u7uqlWrls6cOWPX782bN5fY58vZt2+fDh8+rAcffNBumtzf318RERF2dT/44APddNNNGjRokN17o127dvLz87vs4UNnfvzxRz344IPy8/OTm5ubPDw81KNHD0lyum8udiXj7ty5s06cOKHhw4dr7dq1djNuZVWnTh3deeedZa5fns/F5WzYsEGSKvSdU7NmTd1333125cWft0s/X8XnnhTz9fVVw4YNL/udU2zixIm279/NmzfrpZde0jvvvKPhw4eXaf0BAwbIzc3N9vjS797i929MTIzdYadatWrp3nvv1c6dO/X777/btXnplanBwcGyWCx259W4u7vr1ltvvew4CwoK9NJLLykkJESenp5yd3eXp6en9u/ff9n37+HDh2WxWFS/fn278pMnT2rs2LGaPHmyQkNDS22jNCVd0VneKz03bNigCRMm6L777tMTTzxh91zfvn314osvauDAgerevbsmTJigbdu2yWKxOD35unh26dChQ+Xqg8RJ1VckODj4qp1UfekZ8cUnlBUfIig+f6K8V4cdP35c7u7uatCggV25xWKRn5+fjh8/Xmo/ivtyuUMVx48ft/vRKObn52f3+JdffpEkPfXUUyXeoqAsPxyXtltclp+fr9OnT9vdUqBRo0YOfS3Pa1LStorbki5cTTNs2DDdf//9evrpp+Xn5yd3d3clJSU5PQfrcm1eyVWAZd0X0oWQ9Mknn+i5555Tp06d5O3tLYvFov79+9vt8+PHj5fa58v1p6S6fn5+dodffvnlF504cUKenp5O2ypvqDh9+rQiIyNltVr14osvqmXLlqpRo4YOHjyooUOHlul9XdFxx8TEqKCgQK+//rruvfdeFRUVqVOnTnrxxRfVu3fvMvX/0vfu5ZTnc3E5R48elZubW5nGerHi1+zSH8WGDRvK3d39qn3nFGvSpInd9/Add9whi8Wi+Ph4bdy4UX369Cl1/ct99xb319m+aNy4sYqKivTbb7/ZXaxRt25du3qenp6qUaOGw5Vcnp6eys3NLbV/cXFxWrhwoSZPnqwePXqoTp06qlatmsaMGXPZ1+js2bPy8PCwC3ySNGXKFHl4eOjxxx+3nTdZfLn877//rhMnTsjHx6fEYFP8mjn7rvn1118dxl+ajRs3aujQoerdu7fefvvtMoWpgIAA3X777U7P9Sp+jStyeJ1AVEm8vLwczkWR5PBlUFbFP97//e9/1bRp0zKvV69ePRUUFOjo0aN2AcAwDOXk5KhTp04V6o+z7Xz55ZcO5ZeeyFv8l0p8fHyJVz20atXqsttzdoJwTk6OPD09VatWLbvySz9g5X1NStpWcVuS9I9//EOBgYFKTU21256z90BZ26yosu6LkydP6oMPPtDzzz+vP//5z7by4vOhLm2ztD5frj8l1XX2/qhXr54++ugjp21dPJNQFp9++qkOHz6szz77zDYrJKnEk+cvdSXjlqRHHnlEjzzyiM6cOaOtW7fq+eef18CBA/Wvf/1L/v7+l12/vH9pl+VzYbVanb4vLw2bDRo0UGFhoXJycsoVzOrVq6cvvvhChmHY9f/IkSMqKChwmK2oDMWzPHv27LlsILqc4vdvdna2w3OHDx9WtWrVVKdOnSvaRmn+8Y9/aOTIkbbzzYodO3ZMN910U6nr1q9fX/n5+Tpz5oxq1qxpK//222914MABp2G3+Eqz3377rcT2W7duLUnau3ev3dV7BQUF+uc//1nm2bmNGzdqyJAh6tGjh1atWlXiH0LOGIZR4onikir0PuOQWSUJCAjQN998Y1f26aeflummVc5ERUXJzc1NSUlJ5VqvV69eki58qC62atUqnTlzxvb8lerZs6dOnTpld/WJdOFkxIu1atVKLVq00J49e9SxY0enS1l+9FavXq1z587ZHp86dUrvv/++IiMjHf4aulR5X5PvvvtOe/bscRhX7dq11aFDB0kXfrg8PT3tfgBycnJKvMrsk08+sc2WSRdO0E9NTVXz5s2v+B5RZd0XFotFhmE4XN76t7/9zeHwSs+ePUvs8+W0atVKjRo1UkpKit2VHz///LN27NhhV3fgwIE6fvy4CgsLnb43yhKWLx2jJIcxvvbaa2Va/0rGfbGaNWuqX79+mjJlivLz8/Xdd9/Z9asif806U5bPRUBAgI4cOWI3pvz8fG3cuNGureLDOxX5zjl9+rTee+89u/Liiwuu1ndOaTIzMyXJ6cm55dWqVSvdfPPNWr58ud3798yZM1q1apXtyrPKYrFYHN6/H374YZkOCRXft+c///mPXXliYqI2b95st8ybN0/ShSvlNm/e7PCH5cW6dOmiRo0aOVwE8+677+r06dNlutXCpk2bNGTIEN1+++167733HMZYmp9++kmff/6501sUFJ+MX5HbLDBDVEliYmL03HPPaerUqerRo4e+//57vfrqqxW+O3RAQICeeeYZvfDCCzp79qyGDx8uHx8fff/99zp27JimT5/udL3evXurT58+mjx5snJzc9WtWzd98803ev7559W+ffsyX4lxOSNHjtS8efM0cuRI/eUvf1GLFi20fv16hy9Z6cKPUb9+/dSnTx89/PDDuvnmm/Xrr7/qhx9+0Ndff62VK1dedntubm7q3bu34uLiVFRUpFmzZik3N7fE1+Fi5X1NGjdurLvvvlvTpk1To0aN9I9//ENpaWmaNWuW7Ytw4MCBWr16tcaPH6/77rtPBw8e1AsvvKBGjRpp//79Dn2oX7++7rzzTj333HOqWbOmFi1apH/+858Ol95XRFn3hbe3t7p37645c+aofv36CggI0JYtW7R06VKHvwyfffZZrVu3TnfeeaemTp2qGjVqaOHChWW6dL1atWp64YUXNGbMGN1zzz169NFHdeLECU2bNs3hL9QHHnhAb7/9tvr376+JEyeqc+fO8vDw0H//+19t3rxZgwcP1j333FPm1yIiIkJ16tTRuHHj9Pzzz8vDw0Nvv/22Q8AtyZWM+9FHH1X16tXVrVs3NWrUSDk5OUpISJCPj49tFrL4/I0lS5aodu3aslqtCgwMrPAsYVk+F9HR0Zo6daoeeOABPf300zp37pwWLFjgEIIjIyMVExOjF198Ub/88osGDhwoLy8vZWRkqEaNGg7nehQbOXKk7Z42Bw4cUOvWrbV9+3a99NJL6t+/v+66664Kja0kWVlZtkMnZ86cUXp6uhISEuTv739V7oFUrVo1zZ49Ww899JAGDhyosWPHKi8vT3PmzNGJEyc0c+bMK95GaQYOHKjk5GQFBQWpTZs22r17t+bMmVOmP5yKbzOwc+dOuyvC2rVrV+I6t912m8NtFCwWi3r06GE7h8/NzU2zZ89WTEyMxo4dq+HDh2v//v36v//7P/Xu3Vt9+/a1rbtlyxb16tVLU6dOtZ3zs337dg0ZMkR+fn565plnbAG2WEhIiO1GnnfddZe6d++uNm3ayNvbW3v37tXs2bNlsVj0wgsvOPR/586dcnNzU/fu3S/7+jgo92nYsF0Z9tVXX5VYJy8vz/i///s/o2nTpkb16tWNHj16GJmZmSVeZXZpW86uMDIMw1i2bJnRqVMnw2q1GrVq1TLat29vd1WBsytIzp49a0yePNnw9/c3PDw8jEaNGhl//OMfjd9++82unr+/vzFgwACHsVx6dUZJ/vvf/xr33nuvUatWLaN27drGvffea+zYscPpVRN79uwxhg0bZjRs2NDw8PAw/Pz8jDvvvNNYvHhxqdsovpJi1qxZxvTp040mTZoYnp6eRvv27Y2NGzfa1S2+4uPo0aMO7ZT3NXn33XeN2267zfD09DQCAgKMuXPnOrQ5c+ZMIyAgwPDy8jKCg4ON119/3elVT5KMCRMmGIsWLTKaN29ueHh4GEFBQcbbb79tV6+iV5kZRtn3RXG9OnXqGLVr1zb69u1rfPvttw7vU8MwjM8//9zo2rWr4eXlZfj5+RlPP/20sWTJksteZVbsb3/7m9GiRQvD09PTaNmypfHGG284fb+eP3/eePnll422bdva3udBQUHG2LFjjf3795e6DWevz44dO4zw8HCjRo0aRoMGDYwxY8YYX3/99WWv5invuC/dD3//+9+Nnj17Gr6+voanp6fRuHFjY9iwYcY333xj135iYqIRGBhouLm52fWpR48exm233ea0TyVdMVWWz4VhGMb69euNdu3aGdWrVzduueUW49VXX3X62hUWFhrz5s0zQkNDDU9PT8PHx8cIDw833n///RL7YhiGcfz4cWPcuHFGo0aNDHd3d8Pf39+Ij483zp07Z1ev+LNwKWfvv0s5u8rMarUaLVu2NGJjYx2ugCrpNZszZ45D23JylfB7771ndOnSxbBarUbNmjWNXr16GZ9//rldnZK+c0aNGmXUrFnTYTuX7mNnV4r99ttvxujRo42GDRsaNWrUMG6//XZj27ZtJY7n0vd0ZGSk0b9/f4dtX6qkq8xOnTplSDIeeOABh3WWL19utGnTxvD09DT8/PyMJ5980jh16pTTdi9+PYtfp5KWi7/zYmNjjZCQEKN27dqGu7u70bhxY2PEiBHGvn37nI4jMjLSGDRo0GXH64zFMCpw9yKgihw4cECBgYGaM2dOlf+/cQDg6latWqXo6Gj9/PPPuvnmm8u9/vr16zVw4EDt2bPHdu6Qq/rPf/6jFi1aaOPGjWW+cOFinEMEAMANaujQoerUqZMSEhIqtP7mzZv1wAMPuHwYkqQXX3xRvXr1qlAYkjiHCACAG5bFYtHrr7+udevWqaioqMT/4LUkc+bMqaSeXV0FBQVq3ry54uPjK9wGh8wAAIDpccgMAACYHoEIAACYHoEIAACYHidVl0FRUZEOHz6s2rVrl/tW+gAAoGoYhqFTp06pcePGlz2hnEBUBocPHy7X/x8GAABcx8GDBy97d28CURkU/99aBw8etN1OHAAAuLbc3Fw1bdq0TP9HJoGoDIoPk3l7exOIAAC4zpTldBdOqgYAAKZHIAIAAKZHIAIAAKZHIAIAAKZHIAIAAKZHIAIAAKZHIAIAAKZHIAIAAKZHIAIAAKZHIAIAAKZHIAIAAKZHIAIAAKZHIAIAAKZHIAIAAKZHIAIAAKbnXtUdAG4UWTNaV3UX8D/Npu6t6i4AuM4wQwQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEzP5QLRokWLFBgYKKvVqrCwMG3btq1M633++edyd3dXu3btHJ5btWqVQkJC5OXlpZCQEK1Zs+Yq9xoAAFzPXCoQpaamKjY2VlOmTFFGRoYiIyPVr18/ZWVllbreyZMnNXLkSPXq1cvhufT0dEVHRysmJkZ79uxRTEyMhg0bpi+++KKyhgEAAK4zFsMwjKruRLEuXbqoQ4cOSkpKspUFBwdryJAhSkhIKHG9Bx54QC1atJCbm5vee+89ZWZm2p6Ljo5Wbm6uNmzYYCvr27ev6tSpo5SUlDL1Kzc3Vz4+Pjp58qS8vb3LPzCYQtaM1lXdBfxPs6l7q7oLAFxAeX6/XWaGKD8/X7t371ZUVJRdeVRUlHbs2FHiem+++ab+85//6Pnnn3f6fHp6ukObffr0KbXNvLw85ebm2i0AAODG5TKB6NixYyosLJSvr69dua+vr3Jycpyus3//fv35z3/W22+/LXd3d6d1cnJyytWmJCUkJMjHx8e2NG3atJyjAQAA1xOXCUTFLBaL3WPDMBzKJKmwsFAPPvigpk+frpYtW16VNovFx8fr5MmTtuXgwYPlGAEAALjeOJ9WqQL169eXm5ubw8zNkSNHHGZ4JOnUqVPatWuXMjIy9Pjjj0uSioqKZBiG3N3dtWnTJt15553y8/Mrc5vFvLy85OXldRVGBQAArgcuM0Pk6empsLAwpaWl2ZWnpaUpIiLCob63t7f27t2rzMxM2zJu3Di1atVKmZmZ6tKliyQpPDzcoc1NmzY5bRMAAJiTy8wQSVJcXJxiYmLUsWNHhYeHa8mSJcrKytK4ceMkXTiUdejQIS1btkzVqlVTaGio3foNGzaU1Wq1K584caK6d++uWbNmafDgwVq7dq0+/vhjbd++/ZqODQAAuC6XCkTR0dE6fvy4ZsyYoezsbIWGhmr9+vXy9/eXJGVnZ1/2nkSXioiI0IoVK/Tss8/queeeU/PmzZWammqbQQIAAHCp+xC5Ku5DhLLgPkSug/sQAZCu0/sQAQAAVBUCEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD2XC0SLFi1SYGCgrFarwsLCtG3bthLrbt++Xd26dVO9evVUvXp1BQUFad68eXZ1kpOTZbFYHJZz585V9lAAAMB1wr2qO3Cx1NRUxcbGatGiRerWrZtee+019evXT99//72aNWvmUL9mzZp6/PHH1aZNG9WsWVPbt2/X2LFjVbNmTT322GO2et7e3tq3b5/dulartdLHAwAArg8WwzCMqu5EsS5duqhDhw5KSkqylQUHB2vIkCFKSEgoUxtDhw5VzZo19dZbb0m6MEMUGxurEydOVLhfubm58vHx0cmTJ+Xt7V3hdnBjy5rRuqq7gP9pNnVvVXcBgAsoz++3yxwyy8/P1+7duxUVFWVXHhUVpR07dpSpjYyMDO3YsUM9evSwKz99+rT8/f3VpEkTDRw4UBkZGVet3wAA4PrnMofMjh07psLCQvn6+tqV+/r6Kicnp9R1mzRpoqNHj6qgoEDTpk3TmDFjbM8FBQUpOTlZrVu3Vm5urubPn69u3bppz549atGihdP28vLylJeXZ3ucm5t7BSMDAACuzmUCUTGLxWL32DAMh7JLbdu2TadPn9bOnTv15z//WbfeequGDx8uSeratau6du1qq9utWzd16NBBr7zyihYsWOC0vYSEBE2fPv0KRwIAAK4XLhOI6tevLzc3N4fZoCNHjjjMGl0qMDBQktS6dWv98ssvmjZtmi0QXapatWrq1KmT9u/fX2J78fHxiouLsz3Ozc1V06ZNyzoUAABwnXGZc4g8PT0VFhamtLQ0u/K0tDRFRESUuR3DMOwOdzl7PjMzU40aNSqxjpeXl7y9ve0WAABw43KZGSJJiouLU0xMjDp27Kjw8HAtWbJEWVlZGjdunKQLMzeHDh3SsmXLJEkLFy5Us2bNFBQUJOnCfYlefvllPfHEE7Y2p0+frq5du6pFixbKzc3VggULlJmZqYULF177AQIAAJfkUoEoOjpax48f14wZM5Sdna3Q0FCtX79e/v7+kqTs7GxlZWXZ6hcVFSk+Pl4//fST3N3d1bx5c82cOVNjx4611Tlx4oQee+wx5eTkyMfHR+3bt9fWrVvVuXPnaz4+AADgmlzqPkSuivsQoSy4D5Hr4D5EAKTr9D5EAAAAVYVABAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATM+9qjtwIwt7ellVdwH/s3vOyKruAgDAhTFDBAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATI9ABAAATM/lAtGiRYsUGBgoq9WqsLAwbdu2rcS627dvV7du3VSvXj1Vr15dQUFBmjdvnkO9VatWKSQkRF5eXgoJCdGaNWsqcwgAAOA641KBKDU1VbGxsZoyZYoyMjIUGRmpfv36KSsry2n9mjVr6vHHH9fWrVv1ww8/6Nlnn9Wzzz6rJUuW2Oqkp6crOjpaMTEx2rNnj2JiYjRs2DB98cUX12pYAADAxVkMwzCquhPFunTpog4dOigpKclWFhwcrCFDhighIaFMbQwdOlQ1a9bUW2+9JUmKjo5Wbm6uNmzYYKvTt29f1alTRykpKWVqMzc3Vz4+Pjp58qS8vb3LPB7+6w7XcS3+646sGa0rfRsom2ZT91Z1FwC4gPL8frvMDFF+fr52796tqKgou/KoqCjt2LGjTG1kZGRox44d6tGjh60sPT3doc0+ffqU2mZeXp5yc3PtFgAAcONymUB07NgxFRYWytfX167c19dXOTk5pa7bpEkTeXl5qWPHjpowYYLGjBljey4nJ6fcbSYkJMjHx8e2NG3atAIjAgAA1wuXCUTFLBaL3WPDMBzKLrVt2zbt2rVLixcvVmJiosOhsPK2GR8fr5MnT9qWgwcPlnMUAADgeuJe1R0oVr9+fbm5uTnM3Bw5csRhhudSgYGBkqTWrVvrl19+0bRp0zR8+HBJkp+fX7nb9PLykpeXV0WGAQAArkMuM0Pk6empsLAwpaWl2ZWnpaUpIiKizO0YhqG8vDzb4/DwcIc2N23aVK42AQDAjc1lZogkKS4uTjExMerYsaPCw8O1ZMkSZWVlady4cZIuHMo6dOiQli27cPXWwoUL1axZMwUFBUm6cF+il19+WU888YStzYkTJ6p79+6aNWuWBg8erLVr1+rjjz/W9u3br/0AAQCAS3KpQBQdHa3jx49rxowZys7OVmhoqNavXy9/f39JUnZ2tt09iYqKihQfH6+ffvpJ7u7uat68uWbOnKmxY8fa6kRERGjFihV69tln9dxzz6l58+ZKTU1Vly5drvn4AACAa3Kp+xC5Ku5DdP3jPkTmwn2IAEjX6X2IAAAAqgqBCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmB6BCAAAmJ7LBaJFixYpMDBQVqtVYWFh2rZtW4l1V69erd69e6tBgwby9vZWeHi4Nm7caFcnOTlZFovFYTl37lxlDwUAAFwnXCoQpaamKjY2VlOmTFFGRoYiIyPVr18/ZWVlOa2/detW9e7dW+vXr9fu3bvVs2dPDRo0SBkZGXb1vL29lZ2dbbdYrdZrMSQAAHAdcK/qDlxs7ty5Gj16tMaMGSNJSkxM1MaNG5WUlKSEhASH+omJiXaPX3rpJa1du1bvv/++2rdvbyu3WCzy8/Or1L4DAIDrl8vMEOXn52v37t2KioqyK4+KitKOHTvK1EZRUZFOnTqlunXr2pWfPn1a/v7+atKkiQYOHOgwg3SpvLw85ebm2i0AAODG5TKB6NixYyosLJSvr69dua+vr3JycsrUxl//+ledOXNGw4YNs5UFBQUpOTlZ69atU0pKiqxWq7p166b9+/eX2E5CQoJ8fHxsS9OmTSs2KAAAcF1wmUBUzGKx2D02DMOhzJmUlBRNmzZNqampatiwoa28a9euGjFihNq2bavIyEi98847atmypV555ZUS24qPj9fJkydty8GDBys+IAAA4PJc5hyi+vXry83NzWE26MiRIw6zRpdKTU3V6NGjtXLlSt11112l1q1WrZo6depU6gyRl5eXvLy8yt55AABwXXOZGSJPT0+FhYUpLS3NrjwtLU0RERElrpeSkqKHH35Yy5cv14ABAy67HcMwlJmZqUaNGl1xnwEAwI3BZWaIJCkuLk4xMTHq2LGjwsPDtWTJEmVlZWncuHGSLhzKOnTokJYtWybpQhgaOXKk5s+fr65du9pml6pXry4fHx9J0vTp09W1a1e1aNFCubm5WrBggTIzM7Vw4cKqGSQAAHA5VxSIzp8/r5ycHP3+++9q0KCBw9Vd5RUdHa3jx49rxowZys7OVmhoqNavXy9/f39JUnZ2tt09iV577TUVFBRowoQJmjBhgq181KhRSk5OliSdOHFCjz32mHJycuTj46P27dtr69at6ty58xX1FQAA3DgshmEY5Vnh9OnTevvtt5WSkqIvv/xSeXl5tueaNGmiqKgoPfbYY+rUqdNV72xVyc3NlY+Pj06ePClvb+8yrxf29LJK7BXKY/eckZW+jawZrSt9GyibZlP3VnUXALiA8vx+l+sconnz5ikgIECvv/667rzzTq1evVqZmZnat2+f0tPT9fzzz6ugoEC9e/dW3759Sz1xGQAAwFWU65DZjh07tHnzZrVu7fwv4c6dO+sPf/iDkpKS9MYbb2jLli1q0aLFVekoAABAZSlXIFq5cqXt3wcPHizxhoVWq1Xjx4+/sp4BAABcIxU+qdrf31916tRR27Zt1bZtW7Vr105t27ZVXl6eFi5caLsSDAAAwNVVOBD9+OOPyszMVGZmpjIyMvTuu+/q8OHDklSuE48BAACqWoUDUUBAgAICAjRkyBBbWXp6ukaNGqVZs2Zdjb4BAABcE1f1TtXh4eGaP3++XnzxxavZLAAAQKWqcCA6f/680/IWLVrou+++q3CHAAAArrUKHzKrWbOmQkJC1L59e7Vr107t27dX48aN9corrygqKupq9hEAAKBSVTgQffrpp9qzZ4/27Nmjt99+W88884zOnj0rSYqKitKUKVPUpk0btWnTRsHBwVetwwAAAFdbhQPR7bffrttvv932uKioSPv27bNdebZ792698cYbOnLkiAoLC69KZwEAACrDVfvf7qtVq6bg4GAFBwdr+PDhtvJffvnlam0CAACgUlzVq8yc8fX1rexNAAAAXJFyzRAFBgbKYrGUeyOxsbF68skny70eAADAtVCuQJScnFyhjQQEBFRoPQAAgGuhXIGoR48eldUPAACAKlPhc4g+/vjjEp977bXXKtosAADANVfhQDRgwAD96U9/Un5+vq3s6NGjGjRokOLj469K5wAAAK6FCgeirVu36v3331enTp303Xff6cMPP1RoaKhOnz6tPXv2XM0+AgAAVKoKB6IuXbooIyNDbdq0UVhYmO655x796U9/0qeffqqmTZtezT4CAABUqiu6D9G+ffv01VdfqUmTJnJ3d9c///lP/f7771erbwAAANdEhQPRzJkzFR4ert69e+vbb7/VV199ZZsxSk9Pv5p9BAAAqFQVDkTz58/Xe++9p1deeUVWq1W33XabvvzySw0dOlR33HHHVewiAABA5arw/2W2d+9e1a9f367Mw8NDc+bM0cCBA6+4YwAAANdKhWeILg1DF+MGjgAA4HpSrkCUlZVVrsYPHTpUrvoAAABVoVyHzDp16qS7775bjz76qDp37uy0zsmTJ/XOO+9o/vz5Gjt2rJ544omr0lEAcCXdXulW1V3A/3z+xOdV3QXcAMoViH744Qe99NJL6tu3rzw8PNSxY0c1btxYVqtVv/32m77//nt999136tixo+bMmaN+/fpVVr8BAACumnIdMqtbt65efvllHT58WElJSWrZsqWOHTum/fv3S5Ieeugh7d69W59//jlhCAAAXDcqdJWZ1WrV0KFDNXToUEmSYRiSJIvFcvV6BgAAcI1c0Z2qly5dqtDQUFmtVlmtVoWGhupvf/vb1eobAADANVHhQPTcc89p4sSJGjRokFauXKmVK1dq0KBBmjRpkp599tkKd2jRokUKDAyU1WpVWFiYtm3bVmLd1atXq3fv3mrQoIG8vb0VHh6ujRs3OtRbtWqVQkJC5OXlpZCQEK1Zs6bC/QMAADeeCgeipKQkvf7660pISNDdd9+tu+++WwkJCVqyZIkWL15coTZTU1MVGxurKVOmKCMjQ5GRkerXr1+Jl/tv3bpVvXv31vr167V792717NlTgwYNUkZGhq1Oenq6oqOjFRMToz179igmJkbDhg3TF198UaE+AgCAG4/FKD4BqJzq1KmjL7/8Ui1atLAr/9e//qXOnTvrxIkT5W6zS5cu6tChg5KSkmxlwcHBGjJkiBISEsrUxm233abo6GhNnTpVkhQdHa3c3Fxt2LDBVqdv376qU6eOUlJSytRmbm6ufHx8dPLkSXl7e5d5PGFPLytzXVSu3XNGVvo2sma0rvRtoGyaTd1b6dvgsnvXwWX3KEl5fr8rPEM0YsQIu+BSbMmSJXrooYfK3V5+fr52796tqKgou/KoqCjt2LGjTG0UFRXp1KlTqlu3rq0sPT3doc0+ffqU2mZeXp5yc3PtFgAAcOOq8P9lJl04qXrTpk3q2rWrJGnnzp06ePCgRo4cqbi4OFu9uXPnXratY8eOqbCwUL6+vnblvr6+ysnJKVN//vrXv+rMmTMaNmyYrSwnJ6fcbSYkJGj69Oll2iYAALj+VTgQffvtt+rQoYMk6T//+Y8kqUGDBmrQoIG+/fZbW73yXop/aX3DMMrURkpKiqZNm6a1a9eqYcOGV9RmfHy8XaDLzc1V06ZNy9J9AABwHapwINq8efPV7Ifq168vNzc3h5mbI0eOOMzwXCo1NVWjR4/WypUrddddd9k95+fnV+42vby85OXlVc4RAACA69UV3YfoavL09FRYWJjS0tLsytPS0hQREVHieikpKXr44Ye1fPlyDRgwwOH58PBwhzY3bdpUapsAAMBcrugcoqstLi5OMTEx6tixo8LDw7VkyRJlZWVp3Lhxki4cyjp06JCWLbtw9VZKSopGjhyp+fPnq2vXrraZoOrVq8vHx0eSNHHiRHXv3l2zZs3S4MGDtXbtWn388cfavn171QwSAAC4HJeZIZIuXCKfmJioGTNmqF27dtq6davWr18vf39/SVJ2drbdPYlee+01FRQUaMKECWrUqJFtmThxoq1ORESEVqxYoTfffFNt2rRRcnKyUlNT1aVLl2s+PgAA4JpcaoZIksaPH6/x48c7fS45Odnu8WeffVamNu+77z7dd999V9gzAABwo3KpGSIAAICqQCACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACm53KBaNGiRQoMDJTValVYWJi2bdtWYt3s7Gw9+OCDatWqlapVq6bY2FiHOsnJybJYLA7LuXPnKnEUAADgeuJSgSg1NVWxsbGaMmWKMjIyFBkZqX79+ikrK8tp/by8PDVo0EBTpkxR27ZtS2zX29tb2dnZdovVaq2sYQAAgOuMSwWiuXPnavTo0RozZoyCg4OVmJiopk2bKikpyWn9gIAAzZ8/XyNHjpSPj0+J7VosFvn5+dktAAAAxVwmEOXn52v37t2KioqyK4+KitKOHTuuqO3Tp0/L399fTZo00cCBA5WRkXFF7QEAgBuLywSiY8eOqbCwUL6+vnblvr6+ysnJqXC7QUFBSk5O1rp165SSkiKr1apu3bpp//79Ja6Tl5en3NxcuwUAANy4XCYQFbNYLHaPDcNwKCuPrl27asSIEWrbtq0iIyP1zjvvqGXLlnrllVdKXCchIUE+Pj62pWnTphXePgAAcH0uE4jq168vNzc3h9mgI0eOOMwaXYlq1aqpU6dOpc4QxcfH6+TJk7bl4MGDV237AADA9bhMIPL09FRYWJjS0tLsytPS0hQREXHVtmMYhjIzM9WoUaMS63h5ecnb29tuAQAANy73qu7AxeLi4hQTE6OOHTsqPDxcS5YsUVZWlsaNGyfpwszNoUOHtGzZMts6mZmZki6cOH306FFlZmbK09NTISEhkqTp06era9euatGihXJzc7VgwQJlZmZq4cKF13x8AADANblUIIqOjtbx48c1Y8YMZWdnKzQ0VOvXr5e/v7+kCzdivPSeRO3bt7f9e/fu3Vq+fLn8/f114MABSdKJEyf02GOPKScnRz4+Pmrfvr22bt2qzp07X7NxAQAA1+ZSgUiSxo8fr/Hjxzt9Ljk52aHMMIxS25s3b57mzZt3NboGAABuUC5zDhEAAEBVIRABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTIxABAADTc7lAtGjRIgUGBspqtSosLEzbtm0rsW52drYefPBBtWrVStWqVVNsbKzTeqtWrVJISIi8vLwUEhKiNWvWVFLvAQDA9cilAlFqaqpiY2M1ZcoUZWRkKDIyUv369VNWVpbT+nl5eWrQoIGmTJmitm3bOq2Tnp6u6OhoxcTEaM+ePYqJidGwYcP0xRdfVOZQAADAdcSlAtHcuXM1evRojRkzRsHBwUpMTFTTpk2VlJTktH5AQIDmz5+vkSNHysfHx2mdxMRE9e7dW/Hx8QoKClJ8fLx69eqlxMTEShwJAAC4nrhMIMrPz9fu3bsVFRVlVx4VFaUdO3ZUuN309HSHNvv06VNqm3l5ecrNzbVbAADAjctlAtGxY8dUWFgoX19fu3JfX1/l5ORUuN2cnJxyt5mQkCAfHx/b0rRp0wpvHwAAuD6XCUTFLBaL3WPDMBzKKrvN+Ph4nTx50rYcPHjwirYPAABcm3tVd6BY/fr15ebm5jBzc+TIEYcZnvLw8/Mrd5teXl7y8vKq8DYBAMD1xWVmiDw9PRUWFqa0tDS78rS0NEVERFS43fDwcIc2N23adEVtAgCAG4vLzBBJUlxcnGJiYtSxY0eFh4dryZIlysrK0rhx4yRdOJR16NAhLVu2zLZOZmamJOn06dM6evSoMjMz5enpqZCQEEnSxIkT1b17d82aNUuDBw/W2rVr9fHHH2v79u3XfHwAAMA1uVQgio6O1vHjxzVjxgxlZ2crNDRU69evl7+/v6QLN2K89J5E7du3t/179+7dWr58ufz9/XXgwAFJUkREhFasWKFnn31Wzz33nJo3b67U1FR16dLlmo0LAAC4NpcKRJI0fvx4jR8/3ulzycnJDmWGYVy2zfvuu0/33XfflXYNAADcoFzmHCIAAICqQiACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACm53KBaNGiRQoMDJTValVYWJi2bdtWav0tW7YoLCxMVqtVt9xyixYvXmz3fHJysiwWi8Ny7ty5yhwGAAC4jrhUIEpNTVVsbKymTJmijIwMRUZGql+/fsrKynJa/6efflL//v0VGRmpjIwMPfPMM3ryySe1atUqu3re3t7Kzs62W6xW67UYEgAAuA64V3UHLjZ37lyNHj1aY8aMkSQlJiZq48aNSkpKUkJCgkP9xYsXq1mzZkpMTJQkBQcHa9euXXr55Zd177332upZLBb5+fldkzEAAIDrj8vMEOXn52v37t2KioqyK4+KitKOHTucrpOenu5Qv0+fPtq1a5fOnz9vKzt9+rT8/f3VpEkTDRw4UBkZGaX2JS8vT7m5uXYLAAC4cblMIDp27JgKCwvl6+trV+7r66ucnByn6+Tk5DitX1BQoGPHjkmSgoKClJycrHXr1iklJUVWq1XdunXT/v37S+xLQkKCfHx8bEvTpk2vcHQAAMCVuUwgKmaxWOweG4bhUHa5+heXd+3aVSNGjFDbtm0VGRmpd955Ry1bttQrr7xSYpvx8fE6efKkbTl48GBFhwMAAK4DLnMOUf369eXm5uYwG3TkyBGHWaBifn5+Tuu7u7urXr16TtepVq2aOnXqVOoMkZeXl7y8vMo5AgAAcL1ymUDk6empsLAwpaWl6Z577rGVp6WlafDgwU7XCQ8P1/vvv29XtmnTJnXs2FEeHh5O1zEMQ5mZmWrduvXV6zwA4Ia2pXuPqu4C/qfH1i2V0q5LHTKLi4vT3/72N73xxhv64YcfNGnSJGVlZWncuHGSLhzKGjlypK3+uHHj9PPPPysuLk4//PCD3njjDS1dulRPPfWUrc706dO1ceNG/fjjj8rMzNTo0aOVmZlpaxMAAMBlZogkKTo6WsePH9eMGTOUnZ2t0NBQrV+/Xv7+/pKk7Oxsu3sSBQYGav369Zo0aZIWLlyoxo0ba8GCBXaX3J84cUKPPfaYcnJy5OPjo/bt22vr1q3q3LnzNR8fAABwTS4ViCRp/PjxGj9+vNPnkpOTHcp69Oihr7/+usT25s2bp3nz5l2t7gEAgBuQSx0yAwAAqAoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHoEIgAAYHouF4gWLVqkwMBAWa1WhYWFadu2baXW37Jli8LCwmS1WnXLLbdo8eLFDnVWrVqlkJAQeXl5KSQkRGvWrKms7gMAgOuQSwWi1NRUxcbGasqUKcrIyFBkZKT69eunrKwsp/V/+ukn9e/fX5GRkcrIyNAzzzyjJ598UqtWrbLVSU9PV3R0tGJiYrRnzx7FxMRo2LBh+uKLL67VsAAAgItzqUA0d+5cjR49WmPGjFFwcLASExPVtGlTJSUlOa2/ePFiNWvWTImJiQoODtaYMWP0hz/8QS+//LKtTmJionr37q34+HgFBQUpPj5evXr1UmJi4jUaFQAAcHUuE4jy8/O1e/duRUVF2ZVHRUVpx44dTtdJT093qN+nTx/t2rVL58+fL7VOSW0CAADzca/qDhQ7duyYCgsL5evra1fu6+urnJwcp+vk5OQ4rV9QUKBjx46pUaNGJdYpqU1JysvLU15enu3xyZMnJUm5ubnlGlNh3tly1UflKe++q4hT5worfRsom2uxvwvOFlT6NlA212J/nylgf7uK8uzv4rqGYVy2rssEomIWi8XusWEYDmWXq39peXnbTEhI0PTp0x3KmzZtWnLH4dJ8XhlX1V3AtZTgU9U9wDXkM5n9bSo+5d/fp06dks9l1nOZQFS/fn25ubk5zNwcOXLEYYanmJ+fn9P67u7uqlevXql1SmpTkuLj4xUXF2d7XFRUpF9//VX16tUrNUjdaHJzc9W0aVMdPHhQ3t7eVd0dVDL2t7mwv83FrPvbMAydOnVKjRs3vmxdlwlEnp6eCgsLU1pamu655x5beVpamgYPHux0nfDwcL3//vt2ZZs2bVLHjh3l4eFhq5OWlqZJkybZ1YmIiCixL15eXvLy8rIru+mmm8o7pBuGt7e3qT5AZsf+Nhf2t7mYcX9fbmaomMsEIkmKi4tTTEyMOnbsqPDwcC1ZskRZWVkaN+7C4Y74+HgdOnRIy5YtkySNGzdOr776quLi4vToo48qPT1dS5cuVUpKiq3NiRMnqnv37po1a5YGDx6stWvX6uOPP9b27durZIwAAMD1uFQgio6O1vHjxzVjxgxlZ2crNDRU69evl7+/vyQpOzvb7p5EgYGBWr9+vSZNmqSFCxeqcePGWrBgge69915bnYiICK1YsULPPvusnnvuOTVv3lypqanq0qXLNR8fAABwTRajLKdew5Ty8vKUkJCg+Ph4h0OIuPGwv82F/W0u7O/LIxABAADTc5kbMwIAAFQVAhEAADA9AhEAADA9AhEAADA9AhFKtGjRIgUGBspqtSosLEzbtm2r6i6hEmzdulWDBg1S48aNZbFY9N5771V1l1CJEhIS1KlTJ9WuXVsNGzbUkCFDtG/fvqruFipJUlKS2rRpY7shY3h4uDZs2FDV3XJJBCI4lZqaqtjYWE2ZMkUZGRmKjIxUv3797O4DhRvDmTNn1LZtW7366qtV3RVcA1u2bNGECRO0c+dOpaWlqaCgQFFRUTpz5kxVdw2VoEmTJpo5c6Z27dqlXbt26c4779TgwYP13XffVXXXXA6X3cOpLl26qEOHDkpKSrKVBQcHa8iQIUpISKjCnqEyWSwWrVmzRkOGDKnqruAaOXr0qBo2bKgtW7aoe/fuVd0dXAN169bVnDlzNHr06KruikthhggO8vPztXv3bkVFRdmVR0VFaceOHVXUKwCV4eTJk5Iu/EjixlZYWKgVK1bozJkzCg8Pr+ruuByX+q874BqOHTumwsJC+fr62pX7+voqJyeninoF4GozDENxcXG6/fbbFRoaWtXdQSXZu3evwsPDde7cOdWqVUtr1qxRSEhIVXfL5RCIUCKLxWL32DAMhzIA16/HH39c33zzDf/Z9Q2uVatWyszM1IkTJ7Rq1SqNGjVKW7ZsIRRdgkAEB/Xr15ebm5vDbNCRI0ccZo0AXJ+eeOIJrVu3Tlu3blWTJk2qujuoRJ6enrr11lslSR07dtRXX32l+fPn67XXXqvinrkWziGCA09PT4WFhSktLc2uPC0tTREREVXUKwBXg2EYevzxx7V69Wp9+umnCgwMrOou4RozDEN5eXlV3Q2XwwwRnIqLi1NMTIw6duyo8PBwLVmyRFlZWRo3blxVdw1X2enTp/Xvf//b9vinn35SZmam6tatq2bNmlVhz1AZJkyYoOXLl2vt2rWqXbu2bSbYx8dH1atXr+Le4Wp75pln1K9fPzVt2lSnTp3SihUr9Nlnn+mjjz6q6q65HC67R4kWLVqk2bNnKzs7W6GhoZo3bx6X5d6APvvsM/Xs2dOhfNSoUUpOTr72HUKlKuk8wDfffFMPP/zwte0MKt3o0aP1ySefKDs7Wz4+PmrTpo0mT56s3r17V3XXXA6BCAAAmB7nEAEAANMjEAEAANMjEAEAANMjEAEAANMjEAEAANMjEAEAANMjEAEAANMjEAEAANMjEAEAANMjEAEAANMjEAEwpaNHj8rPz08vvfSSreyLL76Qp6enNm3aVIU9A1AV+L/MAJjW+vXrNWTIEO3YsUNBQUFq3769BgwYoMTExKruGoBrjEAEwNQmTJigjz/+WJ06ddKePXv01VdfyWq1VnW3AFxjBCIApnb27FmFhobq4MGD2rVrl9q0aVPVXQJQBTiHCICp/fjjjzp8+LCKior0888/V3V3AFQRZogAmFZ+fr46d+6sdu3aKSgoSHPnztXevXvl6+tb1V0DcI0RiACY1tNPP613331Xe/bsUa1atdSzZ0/Vrl1bH3zwQVV3DcA1xiEzAKb02WefKTExUW+99Za8vb1VrVo1vfXWW9q+fbuSkpKqunsArjFmiAAAgOkxQwQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEzv/wHgEyqtXbT+bgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import binom\n",
    "\n",
    "mu = 0.25\n",
    "n=4\n",
    "binomial = binom(n, mu)\n",
    "x = np.arange(4)\n",
    "prob_p = binomial.pmf(x)\n",
    "data = pd.DataFrame(zip(x,prob_p))\n",
    "data.columns=['x','p']\n",
    "ax = sns.barplot(data=data, x='x', y='p')\n",
    "ax.set(xlabel='x', ylabel='p(x|$\\theta$)', title='Función de probabilidad de la distribución Binomial({},{})'.format(n,mu))\n",
    "ax.set_xticklabels(x, rotation=0, size=10)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Con el siguiente código calculamos la media varianza y algunos cuantiles de la distribución $\\text{Binom(4, 0.25)}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Distribución Binomial(4, 0.25); media = 1.0, varianza 0.75\n",
      "Distribución Binomial(4, 0.25); probabilidades=[0.25 0.5  0.75]; cuantiles=[0. 1. 2.]\n"
     ]
    }
   ],
   "source": [
    "# media y varianza de la distribución Binomial(4, 0.25)\n",
    "import numpy as np\n",
    "from scipy.stats import binom\n",
    "\n",
    "mu = 0.25\n",
    "n=4\n",
    "binomial = binom(n, mu)\n",
    "\n",
    "media, var = binomial.stats(moments=\"mv\")\n",
    "print('Distribución Binomial(4, 0.25); media = {}, varianza {}'.format(media, var))\n",
    "\n",
    "# calcula algunos cuantiles de la distribución\n",
    "prob = np.array([0.25, 0.5, 0.75])\n",
    "cuantiles = binomial.ppf(prob)\n",
    "print('Distribución Binomial(4, 0.25); probabilidades={}; cuantiles={}'.format(prob, cuantiles))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:#4CC9F0\">Propiedades de la familia exponencial</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recordemos que la función de probabilidad (densidad) de una distribución en la familia exponencial se expresa como\n",
    "$$\n",
    "f(y|\\theta) = \\exp\\left[a(y)b(\\theta) + c(\\theta) + d(y) \\right].\n",
    "$$\n",
    "\n",
    "Adicionalmente la verosimilitud, si no se consideran distribuciones a priori para los parámetros, es la función dada por\n",
    "\n",
    "$$\n",
    "f(\\theta|y) = \\exp\\left[a(y)b(\\theta) + c(\\theta) + d(y) \\right],\n",
    "$$\n",
    " en donde la $y$ ha sido observada y en consecuencia se conoce y $\\theta$ es ahora la variable en la función.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <span style=\"color:#4CC9F0\">Media y varianza de la variable respuesta</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si $Y$ es una variable aleatoria con distribución en la familia exponencial, se puede comprobar que:\n",
    "\n",
    "* $E[a(Y)] = - \\frac{c'(\\theta)}{b'(\\theta)}$\n",
    "* $Var[a(Y)] = \\frac{b''(\\theta)c'(\\theta) -c''(\\theta)b'(\\theta) }{[b'(\\theta)]^3}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color:#4361EE\">Modelos Lineales Generalizados (GLM) clásicos</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En los modelos lineales generalizados clásicos o frecuentista no asumimos una distribución a priori constante. Por ejemplo, es usual tomar como distribución a priori dela parámetros $f(\\theta) = 1$, que realmente no es una densidad propia, es decir no suma uno necesariamente. Entonces, la posterior no necesariamente resulta siendo una distribución, por lo que la tratamos simplemente como una función de los parámetros.\n",
    "\n",
    "Un modelo lineal generalizado GLM(por su sigla en inglés), se define de la configura de la siguiente forma. Asumimos que se tienen variables aleatorias independientes $Y_1, \\ldots, Y_n$, cada una con distribución en la familia exponencial tal que\n",
    "\n",
    "* `Componente aleatoria`: $f(y_i|\\theta) = \\exp \\left[y_ib(\\theta_i) + c(\\theta_i) + d(y_i) \\right]$. \n",
    "\n",
    "* Todas las distribuciones tiene la misma forma Normal, Binomial, Poisson,..., pero probablemente distintos parámetros. Denotemos $\\mathbf{y}= (y_1, \\ldots, y_n)$, $\\boldsymbol{\\theta} = (\\theta_1, \\ldots, \\theta_n)$. Entonces la función de densidad conjunta es dada por\n",
    "$$\n",
    "f(\\mathbf{y}|\\boldsymbol{\\theta}) = \\prod_{i=1}^n \\exp \\left[y_ib(\\theta_i) + c(\\theta_i) + d(y_i)  \\right] = \\exp \\left[\\sum_{i=1}^n y_ib(\\theta_i) + \\sum_{i=1}^nc(\\theta_i) + \\sum_{i=1}^nd(y_i)  \\right]\n",
    "$$\n",
    "\n",
    "* `Componente sistemática`: Se supone que las esperanzas de las $E[Y_i]$  son dadas por $E[Y_i]= \\mu_i$. Antes hemos mencionado que  se tiene que $E[Y] = - \\frac{c'(\\theta)}{b'(\\theta)}$. De esta forma las esperanzas $\\mu_i$ se relacionan con los parámetros $\\theta_i$, mediante la función invertible $h(\\theta)=   - \\frac{c'(\\theta)}{b'(\\theta)}$. Es decir, tenemos que \n",
    "$$\n",
    "\\mu_i = h(\\theta_i).\n",
    "$$\n",
    "Los parámetros $\\mu_i$ están relacionados con parámetros estructurales $\\boldsymbol{\\beta}= (\\beta_1,\\ldots, \\beta_p)$, con $p<n$ en la forma\n",
    "$$\n",
    "g(\\mu_i) = \\mathbf{x}_i^t\\boldsymbol{\\beta},\n",
    "$$\n",
    "en donde $g$ es una función invertible, llamada `función de enlace` y en donde $\\mathbf{x}$ representa un vector de variables observadas.\n",
    "\n",
    "La idea central es que las $\\mathbf{x}_i$ son realizaciones de un vector aleatorio $X=(x_1, \\ldots, x_n)$ que se supone permite predecir adecuadamente a los valores observados $y_i$.\n",
    "\n",
    "Es común escribir $\\hat{y_i} = \\mu_i$, en el sentido que los valores $\\mu_i$ son predictores de las $y_i$. Además observe que\n",
    "\n",
    "$$\n",
    "\\hat{y_i} = g^{-1}(\\mathbf{x}_i^t\\boldsymbol{\\beta}) =  g^{-1}(x_1\\beta_1, \\ldots, x_p\\beta_p).\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color:#4361EE\">Estimación puntual </span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Supongamos que se tienen $n$ variables aleatorias independientes $Y=(Y_1, \\ldots, Y_n)$. Además supondremos que cada una de las variables aleatorias están indexadas por parámetros $\\theta_i$. Es decir $Y_i\\sim f(\\theta_i)$. Además vamos a suponer que existe un parámetro estructural de tamaño $p$, denotado $\\boldsymbol{\\beta}$, y variables aleatorias $X_i$ con realizaciones $x_{ij}$ tal que\n",
    "\n",
    "$$\n",
    "g(\\theta_i) = \\beta_1 x_{1i} + \\beta_2 x_{2i} + \\ldots + \\beta_{p} x_{pi}.\n",
    "$$\n",
    "\n",
    "La función $g$ se denomina función de enlace. La ecuación anterior pone en  comunicación el parámetro privilegiado $\\theta_i$ con una expresión lineal (componente lineal). La función $g$ es necesaria para permitir que la componente lineal del modelo puede ser transportado a la escala del parámetro $\\theta_i$.\n",
    "\n",
    "En el modelo lineal clásico la función de enlace es la función identidad. Es decir, en ese caso se tiene\n",
    "\n",
    "$$\n",
    "\\theta_i= \\beta_1 x_{1i} + \\beta_2 x_{2i} + \\ldots + \\beta_{p} x_{pi} + b.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:#4CC9F0\">Función de densidad conjunta</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La función de probabilidad conjunta si la variables son independientes es dada por\n",
    "\n",
    "$$\n",
    "f(y|\\boldsymbol{\\beta}, \\mathbf{X}) = \\prod_{i=1}^n f(y_i|\\theta_i) = \\prod_{i=1}^n f(y_i|\\boldsymbol{\\beta}, \\mathbf{X}) = \\exp \\left[\\sum_{i=1}^n y_ib(\\theta_i) + \\sum_{i=1}^nc(\\theta_i) + \\sum_{i=1}^nd(y_i)  \\right]\n",
    "$$\n",
    "\n",
    "En donde $\\mathbf{X}$ es la matriz completa de observaciones, que se denomina matriz de diseño, que para efectos prácticos se considera constante. Las columnas de la matriz de diseño representan a las variables explicativas. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:#4CC9F0\">Función de log-verosimilitud</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La función de log-verosimilitud  es dada por\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "L &= \\sum_{i=1}^{n} L_i =\\log \\prod_{i=1}^n f(y_i|\\boldsymbol{\\beta},\\mathbf{X}) = \\left[\\sum_{i=1}^n y_ib(\\theta_i) + \\sum_{i=1}^nc(\\theta_i) + \\sum_{i=1}^nd(y_i)  \\right]\\\\\n",
    "\\mu_i &= h (\\theta_i)\\\\\n",
    "g(\\mu_i) &= \\beta_1 x_{1i} + \\beta_2 x_{2i} + \\ldots + \\beta_{p} x_{pi} = \\eta_i.\n",
    "\\end{align}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:#4CC9F0\">Gradiente de la función de log-verosimilitud</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para estimar $\\boldsymbol{\\beta}$ se require calcular el gradiente de la función de log-verosimilitud $\\nabla_{\\boldsymbol{\\beta}}L$. Usando la regla de la cadena para derivadas se tiene que\n",
    "\n",
    "$$\n",
    "\\frac{\\partial L}{\\partial \\beta_j} = U_j= \\sum_{i=1}^n \\frac{\\partial l_i}{\\partial\\theta_i}\\frac{\\partial \\theta_i}{\\partial \\mu_i}\\frac{\\partial \\mu_i}{\\partial \\beta_j}.\n",
    "$$\n",
    "\n",
    "\n",
    "Puede verificarse directamente que\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\frac{\\partial l_i}{\\partial \\theta_i}  &= b'(\\theta_i)(y_i-\\mu_i), \\\\\n",
    "\\frac{\\partial \\theta_i}{\\partial \\mu_i}  &= \\frac{1}{\\frac{\\partial \\mu_i}{\\partial \\theta_i}}\\\\\n",
    "\\frac{\\partial \\mu_i}{\\partial \\theta_i}  &= b'(\\theta_i)var(Y_i)\\\\\n",
    "\\frac{\\partial \\mu_i}{\\partial \\beta_j}  &=  \\frac{\\partial \\mu_i}{\\partial \\eta_i}x_{ij}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Esto es todo lo que se requiere para escribir un algoritmo de optimización basado en el gradiente."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:#4CC9F0\">El método de Fisher-Scoring</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El método de Fisher-Scoring es un algoritmo muy rápido que se implementa usualmente para este problema e estimar los parámetros $\\boldsymbol{\\beta}$, por lo rápido. En problemas con muchos datos y muchas variables puede ser más lento que un método basado en el gradiente e incluso prohibitivo.\n",
    "\n",
    "El método está basado en la `matriz de información`, la cual es la matriz de covarianza de variables score. Concretamente, se tiene que \n",
    "\n",
    "$$\n",
    "U_j = \\sum_{i=1}^n \\left[ \\frac{(y_i-\\mu_i)}{Var(Y_i)}x_{ij}\\frac{\\partial \\mu_i}{\\partial \\eta_i} \\right]\n",
    "$$ \n",
    "\n",
    "La matriz de información tiene como elemento $ij$\n",
    "\n",
    "$$\n",
    "\\mathcal{I}_{jk} = E(U_iU_j) = \\sum_{i=1}^n \\frac{x_{ij}x_{ik}}{Var(Y_i)}\\left( \\frac{\\partial\\mu_i}{\\partial \\eta_j}  \\right)^2.\n",
    "$$\n",
    "\n",
    "El algoritmo es una variación del método de Newton y  tiene como ecuación de actualización\n",
    "\n",
    "$$\n",
    "\\mathbf{b}^{(m+1)} = \\mathbf{b}^{(m)} + \\left[\\mathcal{I}^{(m)} \\right]^{-1}U^{(m)}\n",
    "$$\n",
    "\n",
    "en donde $\\mathbf{b}^{(m)}$ es  la estimación del parámetro $\\boldsymbol{\\beta}$ en el paso $m$.\n"
   ]
  }
 ],
 "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}