06_RegresionBayesiana.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modelos Lineales - Parte 1. Regresión"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Índice<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Distribución-de-los-parámetros\" data-toc-modified-id=\"Distribución-de-los-parámetros-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Distribución de los parámetros</a></span><ul class=\"toc-item\"><li><span><a href=\"#Prior-de-los-parámetros\" data-toc-modified-id=\"Prior-de-los-parámetros-1.1\"><span class=\"toc-item-num\">1.1&nbsp;&nbsp;</span>Prior de los parámetros</a></span></li><li><span><a href=\"#Máximo-a-posteriori-(MAP)\" data-toc-modified-id=\"Máximo-a-posteriori-(MAP)-1.2\"><span class=\"toc-item-num\">1.2&nbsp;&nbsp;</span>Máximo a posteriori (MAP)</a></span></li><li><span><a href=\"#Distribución-posterior\" data-toc-modified-id=\"Distribución-posterior-1.3\"><span class=\"toc-item-num\">1.3&nbsp;&nbsp;</span>Distribución posterior</a></span></li><li><span><a href=\"#¿Por-qué-distribuciones-normales?\" data-toc-modified-id=\"¿Por-qué-distribuciones-normales?-1.4\"><span class=\"toc-item-num\">1.4&nbsp;&nbsp;</span>¿Por qué distribuciones normales?</a></span></li></ul></li><li><span><a href=\"#Distribición-predictiva\" data-toc-modified-id=\"Distribición-predictiva-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Distribición predictiva</a></span><ul class=\"toc-item\"><li><span><a href=\"#Predicciones-semi-bayesianas\" data-toc-modified-id=\"Predicciones-semi-bayesianas-2.1\"><span class=\"toc-item-num\">2.1&nbsp;&nbsp;</span>Predicciones semi-bayesianas</a></span></li><li><span><a href=\"#La-función-predictiva-posterior\" data-toc-modified-id=\"La-función-predictiva-posterior-2.2\"><span class=\"toc-item-num\">2.2&nbsp;&nbsp;</span>La función predictiva posterior</a></span></li><li><span><a href=\"#Ejemplo\" data-toc-modified-id=\"Ejemplo-2.3\"><span class=\"toc-item-num\">2.3&nbsp;&nbsp;</span>Ejemplo</a></span></li></ul></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Arrancamos con una celda preparatoria, como de costumbre."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 319,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To support both python 2 and python 3\n",
    "from __future__ import division, print_function, unicode_literals\n",
    "\n",
    "# Common imports\n",
    "import numpy as np\n",
    "import os\n",
    "import sys\n",
    "\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# to make this notebook's output stable across runs\n",
    "np.random.seed(42)\n",
    "\n",
    "# To plot pretty figures\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "mpl.rc('axes', labelsize=14)\n",
    "mpl.rc('xtick', labelsize=12)\n",
    "mpl.rc('ytick', labelsize=12)\n",
    "\n",
    "# Where to save the figures\n",
    "PROJECT_ROOT_DIR = \".\"\n",
    "CHAPTER_ID = \"05_ModelosLineales\"\n",
    "IMAGES_PATH = os.path.join(PROJECT_ROOT_DIR, \"plots\", CHAPTER_ID)\n",
    "\n",
    "def save_fig(fig_id, tight_layout=True, fig_extension=\"png\", resolution=300):\n",
    "    os.makedirs(IMAGES_PATH, exist_ok=True)\n",
    "    path = os.path.join(IMAGES_PATH, fig_id + \".\" + fig_extension)\n",
    "    print(\"Saving figure\", fig_id)\n",
    "    if tight_layout:\n",
    "        plt.tight_layout()\n",
    "    plt.savefig(path, format=fig_extension, dpi=resolution)\n",
    "\n",
    "# Ignore useless warnings (see SciPy issue #5998)\n",
    "# import warnings\n",
    "# warnings.filterwarnings(action=\"ignore\", message=\"^internal gelsd\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regresión lineal bayesiana."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hasta ahora, nuestro acercamiento a la regresión lineal se centró en encontrar los parámetros que minimizaban una determinada función de error.\n",
    "\n",
    "En el caso de 1) errores normales, 2) independientes, 3) con la misma varianza, vimos que esa función de error correspondía al logaritmo de la verosimilitud, y que los parámetros que encontrábamos eran estimadores de máxima verosimilitud, que tienen una serie de propiedades muy buenas.\n",
    "\n",
    "La maximización de la verosimilitud tiene un límite, sin embargo: la complejidad (y flexibilidad) del modelo no puede ser muy grande, porque si no, el modelo tiende a ajustar. Hay que *elegir*, entonces, la complejidad del modelo, teniendo en cuenta el conjunto de datos disponible. Vimos que dejar de lado un conjunto de validación (o hacer CV) puede ayudarnos a tomar esta decisión de manera criteriosa.\n",
    "\n",
    "Por otro lado, vimos que podemos modificar la función de error agregando un término de regularización, y controlar la complejidad efectiva del modelo a través de un nuevo hiperparámetro. En este caso, sin embargo, no es trivial encontrar una función de verosimilitud asociada con la función a minimizar."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Distribución de los parámetros"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\newcommand{\\vv}[1]{\\boldsymbol{#1}}$\n",
    "$\\newcommand{\\om}[0]{\\boldsymbol{\\omega}}$\n",
    "$\\newcommand{\\norm}[0]{\\mathcal{N}}$\n",
    "$\\newcommand{\\b}[1]{\\mathrm{\\mathbf{#1}}}$\n",
    "$\\newcommand{\\T}{^\\mathrm{T}}$\n",
    "\n",
    "En la visión bayesiana, los parámetros del modelo también tienen distribuciones de probabilidad, que nos indican nuestro grado de conocimiento del valor de ese parámetro.\n",
    "\n",
    "Para empezar con el *approach* bayesiano, lo primero que tenemos que hacer es decidir un prior para los parámetros.\n",
    "\n",
    "Pero primero, recordemos que, en general, el objetivo de la regresión es encontrar los parámetros $\\vv{\\omega}$ de la distribución *condicional* $p(t | x, \\om, \\beta)$. Desde el primer momento, hicimos una suposición <font size=5>enorme</font>:\n",
    "\n",
    "$$\n",
    "p(t | x, \\om, \\beta) = \\norm(t| y(x, \\om), \\beta^{-1})\\;\\;.\n",
    "$$\n",
    "\n",
    "Es decir, que la distribución condicional de los targets es una normal, centrada en la predicción del modelo, $y(x, \\om)$, que puede expresarse, en el caso de modelos lineales, como:\n",
    "\n",
    "$$\n",
    "y(x, \\om) = \\sum_{j=0}^M \\omega_j \\phi_j(x)\\;\\;,\n",
    "$$\n",
    "\n",
    "donde supusimos $\\phi_0(x) = 1$. Esta sumatoria puede expresarse de manera vectorial, definiendo el vector de las funciones de base, $\\boldsymbol{\\phi} = (\\phi_0, \\ldots, \\phi_M)\\T$, y el de los parámetros, $\\om = (\\omega_0, \\ldots, \\omega_M)\\T$:\n",
    "\n",
    "$$\n",
    "y(x, \\om) = \\om\\T \\boldsymbol{\\phi}(x)\\;\\;.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Con todo esto, encontramos los valores de los parámetros que maximizaban la verosimilitud, $\\vv{\\om}_\\mathrm{ML}$ y $\\beta_\\mathrm{ML}$:\n",
    "\n",
    "$$\n",
    "\\om_\\mathrm{ML} = (\\boldsymbol{\\Phi}\\T \\boldsymbol{\\Phi})^{-1}\\,\\boldsymbol{\\Phi}\\T\\,\\b{t}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\frac{1}{\\beta_{ML}} = \\frac{1}{N}\\sum_{i=1}^N\\left\\{y(x_i, \\boldsymbol{\\omega}_{ML}) - t_i\\right\\}^2\\;\\;,\n",
    "$$\n",
    "\n",
    "donde $\\boldsymbol{\\Phi}$ es la *matriz de diseño*, cuya fila $i$-ésima es simplemente $\\boldsymbol{\\phi}(x_i)$.\n",
    "\n",
    "En la filosofía bayesiana, esto es solo parte de la historia. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prior de los parámetros"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fijemos un prior para los parámetros. Como la verosimilitud es normal, una elección natural (y simple) es elegir priors normales para $\\om$:\n",
    "\n",
    "$$\n",
    "p(\\om | \\gamma) = \\norm(\\om | 0, \\gamma^{-1}\\mathrm{\\mathbf{I}})\\;\\;,\n",
    "$$\n",
    "que depende de un único hiperparámetro escalar, $\\gamma$ que define la precisión de los parámetros (la inversa de la varianza.\n",
    "\n",
    "En el ejemplo del ajuste polinomial, con grado $M$ (en realidad, siempre que haya $M$ funciones de base más un parámetro de sesgo $\\omega_0$) esto es:\n",
    "\n",
    "$$\n",
    "p(\\om | \\gamma) = \\left(\\frac{\\gamma}{2\\pi}\\right)^{(M+1)/2} \\exp{\\left\\{-\\frac{\\gamma}{2}\\om^T\\om\\right\\}}\\;\\;.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Máximo a posteriori (MAP)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para encontrar la distribución posterior de los parámetros, usamos el teorema de Bayes.\n",
    "\n",
    "$$\n",
    "p(\\om | \\b{t}, \\beta, \\gamma) = \\frac{p(\\b{t} | \\om, \\beta, \\gamma)\\,p(\\om | \\gamma)}{p(\\b{t} | \\beta, \\gamma)}\\;\\;,\n",
    "$$ \n",
    "\n",
    "donde escribí de forma explícita todos los elementos condicionales (los datos, $\\b{t}$ --recordemos que cada $t_i$ siempre viene asociado con un $x_i$, que ya no escribo para no abundar--, y los hiperparámetros $\\beta$ y $\\gamma$). Recuerde que el numerador es una constante con respecto a $\\om$, de manera que podemos obviarla para buscar la maximización de la posterior.\n",
    "\n",
    "Entonces, tenemos\n",
    "\n",
    "$$\n",
    "p(\\om | \\b{t}, \\beta, \\gamma) \\propto p(\\b{t} | \\om, \\beta, \\gamma)\\,p(\\om | \\gamma)\\;\\;.\n",
    "$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En este punto ponemos nuestras decisiones de prior y de verosimilitud:\n",
    "\n",
    "$$\n",
    "p(\\om | \\b{t}, \\beta, \\gamma) \\propto \\prod_{i=1}^N \\norm(t_i| y(x_i, \\om), \\beta^{-1})\\; \\left(\\frac{\\gamma}{2\\pi}\\right)^{(M+1)/2} \\exp{\\left\\{-\\frac{\\gamma}{2}\\om^T\\om\\right\\}}\\;\\;.\n",
    "$$\n",
    "\n",
    "Como de costumbre, es más fácil maximizar el logaritmo de la posterior:\n",
    "\n",
    "$$\n",
    "\\ln p(\\om | \\b{t}, \\beta, \\gamma) = -\\frac{\\beta}{2}\\sum_{i=1}^N \\left\\{t_i - \\om\\T\\boldsymbol{\\phi}(x_i)\\right\\}^2 - \\frac{\\gamma}{2}\\om\\T\\om + \\text{const}\\;\\;,\n",
    "$$\n",
    "\n",
    "que es equivalente a la función de error de la regresión *ridge*, con $\\lambda = \\gamma/\\beta$. O sea, que (con esta elección de verosimilitud y prior), encontramos naturalmente una regresión regularizada. Entonces, siguiendo la misma línea de razonamiento que en su momento, obtenemos que el vector de parámetros que maximiza el valor de la distribución posterior es:\n",
    "\n",
    "$$\n",
    "\\om_\\mathrm{MAP} = \\left(\\frac{\\gamma}{\\beta}\\mathrm{\\mathbf{I}} + \\boldsymbol{\\Phi}\\T \\boldsymbol{\\Phi}\\right)^{-1}\\,\\boldsymbol{\\Phi}\\T\\,\\b{t}\n",
    "$$\n",
    "\n",
    "***\n",
    "**Una nota**. Otra vez vemos que el problema del sobreajuste está íntimamente vinculado con la búsqueda de maximizar la verosimilitud. Es lo mismo que pasaba con la moneda: si hacemos tres tiradas y sacamos caras en todas, podemos calcular, via ML, que $\\mu=1$, y el conjunto de entrenamiento (tres caras) estará perfectamente ajustado, pero el modelo va a fallar de forma miserable con datos que no haya visto antes. En cambio, si vemos la posterior de $\\mu$, encontramos que existe posibilidad de obtener ceca en tiradas ulteriores.\n",
    "***\n",
    "\n",
    "Con otras elecciones de prior, podemos llegar a Lasso u otro tipo de regularizaciones (ver ec. 3.56 del [Bishop](https://www.microsoft.com/en-us/research/uploads/prod/2006/01/Bishop-Pattern-Recognition-and-Machine-Learning-2006.pdf))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Distribución posterior"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ya tenemos el valor que maximiza la posterior. Pero, ¿cómo es la distribución posterior completa? \n",
    "\n",
    "Para encontrar esto, usamos el hecho de que la distribución normal es el prior conjugado para una verosimilitud normal (si suponemos conocida la matriz de precisión, $\\beta$). Por lo tanto, la posterior también será una distribución multinormal.\n",
    "\n",
    "$$\n",
    "p(\\om | \\b{t}, \\beta, \\gamma) = \\norm(\\om | \\mathbf{m}_N, \\mathbf{S}_N)\\;\\;,\n",
    "$$\n",
    "\n",
    "donde $\\mathbf{m}_N$ y $\\mathbf{S}_N$ son el vector de valores medios y la matriz de covarianza, respectivamente, y el subíndice hace alusión al tamaño de los datos tenidos en cuenta."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para el prior que elegimos arriba, \n",
    "\n",
    "$$\n",
    "p(\\om | \\gamma) = \\norm(\\om | \\mathbf{0}, \\gamma^{-1} \\mathbf{I}) = \\left(\\frac{\\gamma}{2\\pi}\\right)^{(M+1)/2} \\exp{\\left\\{-\\frac{\\gamma}{2}\\om^T\\om\\right\\}}\\;\\;,\n",
    "$$\n",
    "\n",
    "obtenemos, usando propiedades de la distribución normal, que:\n",
    "\n",
    "$$\n",
    "\\boxed{\\begin{array}{lll}\n",
    "\\mathbf{m}_N &=& \\beta\\;\\mathbf{S}_N \\boldsymbol{\\Phi}\\T \\b{t}\\\\\n",
    "\\mathbf{S}^{-1}_N &=& \\gamma\\;\\mathbf{I} + \\beta\\;\\boldsymbol{\\Phi}\\T \\boldsymbol{\\Phi}\\;\\;.\n",
    "\\end{array}}\n",
    "$$\n",
    "\n",
    "Noten que si metemos el resultado de $\\mathbf{S}_N$ en al ecuación de la media, recuperamos la ecuación de arriba para $\\om_\\mathrm{MAP}$. Esto es razonable, como la normal es simétrica, el valor que maximiza la posterior es también el valor medio."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En un ejercicio vamos a ver cómo pueden usarse estas fórmulas para un aprendizaje secuencial de la posterior. Acá, veamos cómo se compara lo que sale de acá con los estimadores de máxima verosimilitud que vimos para la ordenada al origen ($\\omega_0$) y pendiente ($\\omega_1$) de una regresión lineal hace unas clases (la visión frecuentista).\n",
    "\n",
    "Si tenemos un modelo lineal, $y(x, \\om) = \\omega_0 + \\omega_1\\,x$, la matriz de diseño y su transpuesta son\n",
    "\n",
    "$$\n",
    "\\boldsymbol{\\Phi} = \\begin{pmatrix}\n",
    "                    1 & x_1\\\\\n",
    "                    1 & x_2\\\\\n",
    "                    \\vdots& \\vdots\\\\\n",
    "                    1 & x_N\\\\\\end{pmatrix}\n",
    "$$                    \n",
    "\n",
    "$$\n",
    "\\boldsymbol{\\Phi}^T = \\begin{pmatrix}\n",
    "                    1 & 1 & \\cdots & 1\\\\\n",
    "                    x_1 & x_2 & \\cdots & x_N\\\\\n",
    "   \\end{pmatrix}\\;\\;.\n",
    "$$\n",
    "\n",
    "Entonces, la matriz de precisión queda\n",
    "\n",
    "$$\n",
    "\\boldsymbol{S}_N^{-1} = \\beta \\begin{pmatrix}N + \\gamma/\\beta & N\\bar{x}\\\\\n",
    "                                       N\\bar{x} & \\sum_i x_i^2 + \\gamma/\\beta\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "si tenemos un prior que se va ensanchando hasta el infinito: es decir, que no tenemos prior, y en términos matemáticos, $\\gamma \\rightarrow 0$, queda\n",
    "\n",
    "$$\n",
    "\\boldsymbol{S}_N^{-1} = N\\beta \\begin{pmatrix}1 & \\bar{x}\\\\\n",
    "                                       \\bar{x} & \\frac{1}{N}\\sum_i x_i^2\n",
    "\\end{pmatrix}\\;\\;,\n",
    "$$\n",
    "\n",
    "que es la misma expresión que teníamos para la covarianza de los estimadores de ML."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pensemos ahora en prior más general (pero siempre normal):\n",
    "\n",
    "$$\n",
    "p(\\om | \\mathbf{m}_0, \\mathbf{S}_0) = \\norm(\\om | \\mathbf{m}_0, \\mathbf{S}_0)\\;\\;,\n",
    "$$\n",
    "\n",
    "donde ahora la matiz de covarianza es general. Entre otras cosas, esto permite acomodar constantes de regularización diferentes para cada parámetro, que fue una pregunta de hace unas clases; además, podemos pensar en una matriz no diagonal, mucho más flexible. \n",
    "\n",
    "Como sea, el valor medio y matriz de precisión de la posterior en este caso es:\n",
    "\n",
    "$$\n",
    "\\boxed{\\begin{array}{lll}\n",
    "\\mathbf{m}_N &=& \\mathbf{S}_N\\; \\left(\\mathbf{S}^{-1}_0 \\mathbf{m}_0 + \\beta\\;\\boldsymbol{\\Phi}\\T \\b{t}\\right)\\\\\n",
    "\\mathbf{S}^{-1}_N &=& \\mathbf{S}^{-1}_0 + \\beta\\;\\boldsymbol{\\Phi}\\T \\boldsymbol{\\Phi}\\;\\;.\n",
    "\\end{array}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Con esto, tenemos una descripción completa de la distribución del vector de parámetros $\\om\\T = (\\omega_0, \\ldots, \\omega_M)$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ¿Por qué distribuciones normales?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A esta altura del partido, es lícito preguntarse por qué estamos insistiendo tanto con las distribuciones normales, y dándole un rol tan central.\n",
    "\n",
    "Las razones son múltiples, pero una de las más potentes es posiblemente el Teorema Central del Límite (ver  [el notebook correspondiente](06a_TeoremaCentralLimite.ipynb)) y las diapositivas de la clase de hoy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Distribición predictiva"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hace un tiempo hablamos de la prevalencia de la visión predictiva vs. la visión inferencista o explicativa. De forma que todo lo de arriba, referido a los parámetros solo es una parte de la visión bayesiana, relacionada con cómo encontrar los mejores parámetros y su distribución.\n",
    "\n",
    "Si queremos hacer predicciones, tenemos dos caminos: la forma semi-bayesiana, en la que fijamos los valores de los parámetros al valor MAP, o la forma *fully bayesian* en la que usamos la regla del producto y de la suma de forma sistemática."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Predicciones semi-bayesianas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En este approach, la predicción consiste en tomar la verosimilitud y evaluarla en el valor MAP del vector de parámetros $\\om$:\n",
    "\n",
    "$$\n",
    "p(t | x, \\om_\\mathrm{MAP}, \\beta) = \\norm(t| y(x, \\om_\\mathrm{MAP}), \\beta^{-1})\\;\\;,\n",
    "$$\n",
    "\n",
    "donde, recordemos, suponemos que $\\beta$ es conocida (ya vamos a ver cómo se puede *inferir* en esta visión). Con esta distribución, las predicciones sobre $t$ tienen naturaleza probabilística, como queda más o menos claro en esta figura del Bishop:\n",
    "\n",
    "<img width=400px src=\"images/Figure1.16.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### La función predictiva posterior"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sin embargo, haber fijado el parámetro $\\om$ a un valor, cuando en realidad contamos con toda su distribución (la posterior, $p(\\om | \\b{t})$) no es totalmente bayesiano. Mi abuela decía que fijar parámetros hace llorar al niño Bayes.\n",
    "\n",
    "Entonces, ¿cómo hacemos para tener en cuenta *toda* la información disponible a la hora de hacer predicciones? ¿Cómo hacemos para ser completamente bayesianos.\n",
    "\n",
    "Primero, pensemos qué es lo que queremos realmente. \n",
    "\n",
    "***\n",
    "**Pregunta**: ¿Alguna idea? ¿Qué es predecir?\n",
    "\n",
    "¿Les parece razonable pensar que queremos saber cómo será la distribución de una nueva medición $t^\\prime$, pero una vez que hayamos considerado los datos ya obtenidos, $\\b{t}$?\n",
    "\n",
    "Y si es así, ¿cómo escribo eso? ¡¿Y cómo lo calculo?!\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El objetivo, entonces, es encontrar\n",
    "\n",
    "$$\n",
    "p(t^\\prime | \\b{t}, \\beta, \\gamma)\\;\\;,\n",
    "$$\n",
    "\n",
    "donde condicioné a los valores de los hiperparámetros. Mal no puede hacer.\n",
    "\n",
    "Así escrito esto no lo sé calcular, pero podemos escribir la distribución conjunta con una variable extra, y usar la regla de la suma para marginalizar:\n",
    "\n",
    "$$\n",
    "p(t^\\prime | \\b{t}, \\beta, \\gamma) = \\int p(t^\\prime, \\om | \\b{t}, \\beta, \\gamma)\\;\\mathrm{d}\\om\\;\\;.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En otras palabras, marginalizo sobre $\\om$. Ahora, el integrando puedo expresarlo de una manera más simpática, usando la regla del producto:\n",
    "\n",
    "$$\n",
    "p(t^\\prime, \\om | \\b{t}, \\beta, \\gamma) = p(t^\\prime | \\om, \\b{t}, \\beta, \\gamma)\\;p(\\om | \\b{t}, \\beta, \\gamma)\\;\\;.\n",
    "$$\n",
    "\n",
    "Veamos los dos términos. En el primero, como estamos condicionando con respecto a $\\om$, la distribución de $t^\\prime$ se vuelve independiente de $\\b{t}$ (estamos usando fuertemente la hipótesis de independencia). Entonces, el primer término es $p(t^\\prime | \\om, \\beta, \\gamma)$, donde reconocemos a la verosimilitud de un dato (es independiente de $\\gamma$):\n",
    "\n",
    "$$\n",
    "p(t^\\prime | \\om, \\beta) = \\norm(t^\\prime| y(x, \\om), \\beta^{-1})\\;\\;.\n",
    "$$\n",
    "\n",
    "El segundo término, $p(\\om | \\b{t}, \\beta, \\gamma)$ no es ni más ni menos que la posterior del vector de parámetros."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Podemos entonces interpretar a la integral de arriba como el valor medio de la verosimilitud con respecto a la distribución posterior de $\\om$\n",
    "\n",
    "***\n",
    "**Recordatorio**. El valor medio de una función $f(x)$ con respecto a la función de distribución $p(x)$ es\n",
    "\n",
    "$$\n",
    "<f(x)>_{p(x)} = \\int f(x)\\,p(x)\\;\\mathrm{d}x\\;\\;.\n",
    "$$\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gracias a algunas propiedades adicionales de las distribuiones normales, podemos resolver analíticamente la integral que da lugar a la posterior predictiva:\n",
    "\n",
    "$$\n",
    "p(t^\\prime | \\b{t}, \\beta, \\gamma) = \\norm(t^\\prime | \\mathbf{m}_N\\T \\boldsymbol{\\phi}(x), \\sigma_N^2)\\;\\;,\n",
    "$$\n",
    "\n",
    "donde la variaza es\n",
    "$$\n",
    "\\sigma_N^2 = \\frac{1}{\\beta} + \\boldsymbol{\\phi}(x)\\T \\mathbf{S}_N \\boldsymbol{\\phi}(x)\\;\\;,\n",
    "$$\n",
    "que consiste en una contribución debida a la incertidumbre en los parámetros, y otra que viene de la incerteza de los datos."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ejemplo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Veamos esto en acción con nuestro ya conocido dataset sinusoidal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 354,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Array para plotear\n",
    "xx = np.linspace(0, 1, 100).reshape([-1, 1])\n",
    "\n",
    "def plot_data_sine(x, t, ax=None):\n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "    ax.plot(x, t, 'ob', mfc='None', ms=10)\n",
    "    ax.plot(xx, np.sin(2*np.pi * xx), 'g-', lw=2, alpha=0.7, label='Ground Truth')\n",
    "    ax.set_xlabel('x')\n",
    "    ax.set_ylabel('t')\n",
    "    ax.legend(loc=0)\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 355,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAERCAYAAAC+ZEqyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2dd3hU1daH351QQieFcLEENIJSBAuogKKC9TogVhQVUARRgoIVUK54RRALlgt6r6iIJSifimJBBWxUFZQioJAIIiAQeg1JyPr+2JkwGWYmM5mSmWS9z3OeyZyz55y9z0z27+y91l7LiAiKoiiKEgxx5V0BRVEUJfZRMVEURVGCRsVEURRFCRoVE0VRFCVoVEwURVGUoFExURRFUYJGxURRFEUJmqgTE2NMdWPMa8aYP40xe40xS4wxl/soP8QYs9kYs8cY87oxpnok66soiqJEoZgAVYC/gPOBesAjwFRjTBP3gsaYS4GhQBegMXAi8FikKqooiqJYTCysgDfGLAMeE5EP3PZnAutEZHjR+y7AOyLyD1/nS0lJkSZNmoSruoqiKBWSxYsXbxORBp6OVYl0ZQLFGNMQaAas8HC4JfCxy/ulQENjTLKIbHc7T3+gP0BaWhqLFi0KU40VRVEqJsaYP70di8ZprmKMMVWBd4DJIvKbhyK1gd0u751/13EvKCKviEhbEWnboIFHYVUURVHKSNSKiTEmDngLyAMyvBTbB9R1ee/8e28Yq6YoiqK4EZViYowxwGtAQ+AaEcn3UnQF0MblfRtgi/sUl6IoihJeotVm8jLQHLhIRA76KPcm8IYx5h1gE9bz643wV0+JJrKzYfx4yMyEbdsgJQV69oSMDEhPL+/aKf6Qn5/Phg0byM3NLe+qKEBCQgLHHXccVatW9fszUScmxpjGwB3AIWCzHaRA0b45wEqghYisF5EvjDFPAd8ANYAPgEcjX2ulvJgxA3r1gn79YP58aNwY/vwTXnsNzjkH3nwTLve6SkmJFjZs2ECdOnVo0qQJLv/zSjkgImzfvp0NGzZwwgkn+P25mHANDjVt27YV9eaKfbKzrWBMnw7t2x99fMEC6NYNFi7UEUq0s2rVKk455RQVkihBRPjtt99o3rx5if3GmMUi0tbTZ6LSZqIo/jB+vB2ReBISsPtvvx0mTIhsvZSyoUISPZTlu1AxUWKWzEzo29d3mdtvt+UURQkvKiZKzLJtm7WR+CItzZZTlNLYsmULPXv25MQTT+TMM8+kffv2TJs2LaJ1WLduHa1atSqxb/ny5Zx22mmcdtppJCUlccIJJ3Daaadx0UUX+X3OTJcnqjfeeIOMDG+rLcqOiokSs6SkWGO7L9avt+UUxRciQvfu3enUqRN//PEHixcv5t1332XDhg1HlS0oKIho3U499VSWLFnCkiVL6NatG08//TRLlixh1qxZftXJXUzChYqJErP07Gm9tnzx6qu2nKL44uuvv6ZatWoMGDCgeF/jxo0ZNGgQYJ/mu3XrRufOnenSpQs7duyge/futG7dmnPOOYdly5YBMHLkSJ555pnic7Rq1Yp169axbt06mjdvTr9+/WjZsiWXXHIJBw/aVQ+LFy+mTZs2tGnThgkBGPguuOACBg8eTNu2bXnhhRfo06cP77//fvHx2rVrAzB06FDmzJnDaaedxnPPPQfApk2buOyyy2jatCkPPvhgGe9aSaLONVhR/CUjw3pzde3q3Zvr1VetN5cSO3Sd0jUs5/3kxk+8HluxYgVnnHGGz8///PPPLFu2jKSkJAYNGsTpp5/ORx99xNdff02vXr1YsmSJz8+vWbOGKVOmMHHiRK6//no++OADbr75Zm699VbGjx9Pp06deOCBBwJqU15eXnGcwT59+ngs8+STT/LMM8/w6aefAlYYlyxZwi+//EL16tU5+eSTGTRoEMcff3xA13ZHRyZKzJKebteRdOsGw4ZZV+H8fPs6bJjd/+ab6hasBM7AgQNp06YN7dq1K9538cUXk5SUBMDcuXO55ZZbAOjcuTPbt29nz549Ps/ptHUAnHnmmaxbt45du3axa9cuOnXqBFB8Tn/p0aNHQOWddOnShXr16pGQkECLFi34s7T5Yj/QkYkS01x+uR15TJgAHTuWXAGv60tiE18jiHDRsmVLPvjgSIaLCRMmsG3bNtq2PbKkolatWqWep0qVKhQWFha/d13RX736kbx98fHxxdNcweBaJ9drFxYWkpeX5/Vz7nUJhR1IRyZKzJOeDuPGwebNUFBgX8eNUyFR/Kdz587k5uby8ssvF+87cOCA1/LnnXce77zzDgDffvstKSkp1K1blyZNmvDzzz8Ddlps7dq1Pq9bv3596tevz9y5cwGKz1kWmjRpwuLFiwGYPn06+fk2pGGdOnXYuzf8sW9VTBRFqfQYY/joo4/47rvvOOGEEzjrrLPo3bs3Y8eO9Vh+5MiRLF68mNatWzN06FAmT54MwDXXXMOOHTto2bIl48ePp1mzZqVee9KkSQwcOJDTTjuNYCKS9OvXj++++442bdqwYMGC4lFL69atiY+Pp02bNsUG+HCg4VQURSl3Vq1adVToDqV88fSdaDgVRVEUJayomCiKoihBo2KixCTZ2TBkCDRsCPHx9nXIELtfiU0q45R7tFKW70LFRIk5ZsywixVr1LA5TA4dsq81atj9M2aUdw2VQElISGD79u0qKFGAM59JQkJCQJ9TA7wSU2gOk4qJZlqMLrxlWvRlgFcxUWKKIUPsCGT0aO9lhg2zo5Vx4yJXL0WpDKg3l1Jh0BwmihKdqJgoMYXmMFGU6ETFRIkpNIeJokQnUSkmxpgMY8wiY8whY8wbPsr1McYcNsbsc9kuiFxNlUijOUwUJTqJ1qjBm4BRwKVAjVLKLhCRc8NfJSUa0BwmihKdROXIREQ+FJGPgO3lXRclutAcJooSnUSlmATI6caYbcaY1caYEcYYj6MtY0z/oqmzRTk5OZGuoxJCnDlMDh2yOUxq1LCvhw7Z/ZdfXt41VJTKR6yLyfdAKyAVuAa4EfCY91JEXhGRtiLStkGDBhGsYujRUCKaw0RRoo2YFhMR+UNE1opIoYgsB/4NXFve9QonGkpEUZRoJFoN8GVFAFPelQgX2dnQq9fRoUTS0+2K8K5dNZSIoijlQ1SOTIwxVYwxCUA8EG+MSfBkCzHGXG6MaVj09ynACODjyNY2cowfD/36efZiArv/9tttPnRFUZRIEpViAjwCHASGAjcX/f2IMSataC1JWlG5LsAyY8x+4HPgQ8BH1KbYRkOJKIoSrWigxxgiPt7aSKr4mJzMz7f2k4KCyNVLUZTKgQZ6rCBoKBFFUaIVFZMYQkOJKIoSraiYxBAZGTBxog0Z4glnKJGBA4O/lq5lURQlEFRMYohIhRLRtSyKogSKikmMEe5QIq5rWUaPtsJUpcqRtSzTp9vjOkJRFMUV9eZSSqBpcRVF8YZ6c1USQmHn0LUsiqKUBRWTCkKo7ByaFldRlLJQ0WJzVUpCGbPLuZbFVzldy6Ioijs6MqkAhDJml65lURSlLKgBvgLQsKGd0vI1msjOtl5fmzf7Pld2tp0Wcx/lOFmwQCMTK0plRQ3wFZxQ2jk0La6iKGVBxaQCEOqYXZoWV1GUQNFprgqArg1RFCUS+JrmUm+uCkBGhrVzdO3q3c7x6qt2VKEoihIOdJqrAqB2DkVRyhsVkwqC2jkURSlP1GaiKIqi+IW6BisRQ/OgKErlRMVECRmaB0VRKi9RKSbGmAxjzCJjzCFjzBullB1ijNlsjNljjHndGFM9QtUMOZ6e6m+91W7R/qTvKw9K377QuTNccQXExYW2DToSUpToICrFBNgEjAJe91XIGHMpMBToAjQGTgQeC3vtwoCnp/oxY+C99+w2Zkx0P+k744Od3i6Xv/f+TdaOLH7f9jsvf/ArbS9fSXL6Om4ZkMNdd+eGrA2u92zKFOsiffgwvPACNG0KV12loqIokSKqDfDGmFHAcSLSx8vxTGCdiAwvet8FeEdE/uHrvNFmgPcUD8t1HxwdD8tbjKzsbNuxZ2ba8CkpKTYoY0ZGaF2DcwtyWZmzkuwd2fy5+0/+9eyfXODYQnzCweIy+/fD3LlwVjtITLLv582D67rWo1HtRuTvOJZJz57Eh6+cTKfWTagaX9Xv67ven1277KioXz87CmrcGKZNg969ISEB3nlHvdkUJRRU5EWLLYGPXd4vBRoaY5JFZLtrQWNMf6A/QFpaWuRq6Aeeov6673NG/XWuYHeNBOzcN2MG3HQTnHSSfUIH+zpvHkyefKRTLYvgiAhrd61l3vp5LN2ylDU71lAohcXHDyZAXDWoGleVxIREalerzdxlVTnz+HjaNy1kf/5+9lU/QF7ubnYfshv8RoPLZ9PvAzh9TTVOTT2VMxudSdtj2tKoTiO/7llqqhVV98CU110HP/8Ma9daoXGKbqTEVlEqG7E+MskGBorIF0XvqwJ5wAkiss7beaNtZOIp6q/7Pk9Rf133ZWfDGWdYm8Sddx55Qv/zTxtS/uWXobAQnnsOHnqo5FO8s8zEiXZxo+tT/MY9G5n1xyzmrJ/Dlv1bivfHmTiaJjXl5OSTaVK/CTdc3phvph/DqSfXwhjjtQ0dOgor/tjO33v/Zv3u9SxYvYZxk1fT5eq/StyTpklNubDJhXRq3Il6CfW83rPx40uGknEXC7D35dRTrcC4j2C8tV1FR1GOxtfIJNbFZCnwhIhMLXqfDGwDUtxHJq5Em5jEx1t7SJUq3vfl59tOs6DgSBnXfbfeam0rs2d7D6nSubMVlG+/9R1eft78w2ypNp8vsr5g2dZlxccTExJpf1x72h3bjhYNWlCzas3iY57ig7m3wVN8MGcbtu3bxc9//8yiTYtYtGkRBwvsdFm8iafD8R248uQrOTnl5KPOfeyxRwRrxoySYnHMMVC7NvTvD//7H9SpA198UXpo/dWr/RcdRalM+BITRCRqN6wR/g0fxzOxYuJ83xnYXNp5zzzzTIkmUlNFsrJ878vKEmnYsGQZ1301aogMGOD5/FlZIoMHi8THi4A99+DBR1/zUMEhuW7459JqeF9xZDrEkemQa967Rl5Y+IIs37JcDhce9tqGrCyRlBSR+fM9t2H+fHvc/Zqe2nWo4JB8v+57eezbx6TblG7Fdbnvy/tk4V8LpbCwsPjccXEi+fmer+88d16eiDH2Hrlf35WhQ0VuvfXo87jirR2KUhkAFomXfjUqvbmMMVWMMQlAPBBvjEkwxniy77wJ9DXGtDDG1AceAd6IYFVDgqfshu77PGU3dN138CAMHnz0uV09nmrXtvvcvakKCgv4fM3n3D79drad8BKrN23hmNrHMODMAUzuPpm7z76bVqmtiDPefy6e4oP16AHPPHN0fDBXd96mTWHv3pLuvNXiq3Fe4/P41/n/4rVur3Ft82upXa02v2//nVFzRnHvl/fS5ZZFvPqaFIff92R3ct6f9evBGDva8JVt8vbb7eguVFkrFaVS4U1lynMDRgLito0E0oB9QJpL2XuBLcAeYBJQvbTzR9vIxNtTtXOfp6dh930g8ttvvs8LdnMyb16h1Gs1T258547ip/+MT++WuCZzfI5CSmvLkCF2RBAfb0cEPXocqefnn9s6DRsmMnWqSHKyyOzZ9n1Kij3uiYP5B2X6b9Pllg9vEUemQy6c6JAaXR+SK27OlmHDjh7Jud6foUOP3B/3UZAreXm2XGmjDk+jKUWpDOBjZFLuwlEeW7SJiciRTnboUNtZ5eWJvPaanZqpUcP+nZd3pHN073hr1jx6mmvwYNtJO0lIsOcSEflr918yfNZwSR/hkBMedkj/6f1l3vp5smZNYUg7Std2zZ5t/5461XMb/JlCys3PlWmrpknPD3rKWc87pGqvrlK1w3ghYZccOHD0/XGeMylJZNUqK3DeyMqy/xH5+b7blJfn+zyKUlHxJSZROc1VGfEU9Xf4cLjhBrsNH+47EvD111v33wULjuzLzLQGZLD7CwogvVkeby19i0EzBrFs6zJapNdh51d3MuGKCXQ4vgOvvWaOmk4LVbu6doXt22HQIM9t8GcKqXqV6nQ/pTuvOF6h33lXcn6nOJLO/gKuuJOaLb7lrLOFAwfsdNX33x+ZXuvVyy5m9JVt8tVXoWbN0GatVJTKQlR7c4WLaPPmCpbsbBg1yoqJiO0Qr7/edqIrV8Ibb1jX4Py6qylo9xztOm8gMQkuOfESbmrZh9T6dSgo8L4QMlR4coH21BZ3F2hfbNizgf8u+i+Tv1zK3r1QsL4tB78bSINaKfTsCQMHHrHTnHqqbd+77x59HmfbHQ5o1EizViqKJzRqcAXGaWBv1Mg+WScm2s7zvfesG3DbtjDr6wLym79N6/se4LROG1j8zfH8Y/lT/DNlEFs31CExMTJJtLZts262vkhLO7I+xB+Oq3scfdIep8Wee9idU4v99RZR9coMzrv5+2IhAdi6FapWha++8p1A7JFHrPuv6wjPFWfWyoED/a+jolQGVExiGNfgin37wvLl1mtp+XL75Aywj80cPPdBzu77Hg1ShTs6XcXPo5+nTm5zOnaEFi3gwIHgkmi5B1tMTrYLBVNSSgZfrF8/9FNIM2ZA+/aGMxMv4oVLXqLG9rNpePx+fqr6NG0GjOPlVw8Ui8W778JPPx2dQGzLFjsi6dMHmjWzItOlCwwYoFkrFcVfVExiGKc77K5dR9x/f/zRdn6rV8OVg+bCZffwd94aTmiQypMXPcltp99G82bVGDfOxq9KToZly+yUTVk6SPcAlR9/bKfaEhLsyOjjj4+4Ih88aG0/vvDkAu0N90jFd/dLYtn4hznz8EA2b6rG/tRvuGv6YLK2/1EslOnptq2bN1sb0qRJ8MkndmTnDLC5eLGN6zV5MrRrp1krFcUvvFnmK/IWjd5cZSE19YiHlKtbcf7hfPnfov+JI9MhJ41wCOeOltsG7C32EvPmERYo7q7H7u/dvbOmTrWuwlOnej5foAsC3b3VXPlr919y9+d3S/oIh5z08NUyK3tWqfUPtj6KUtHBhzeXGuBjmPh4GyuqVq0jBuPdubt5cu6T/LDuV9avrcK2mX3Z/dMVxMcbqle3T9fOOFOuNoWy4B5CxVNIFXdj9fXX25HA4MHWcystzU5tvfqq3QIJVVKaQT/vcB6jvvgvY9+fySWXwGXpl3FH2zuoElfFa33dUWO7ohwhZsOphGsrz5GJM7RJaqoNBeIttIk/pKbaRX/Oz2Ztz5I+H/WRs553SPUbekn/4b/J7Nl2jYVzoWCwoxH367vW21NYGPcFfs7RgOvCxoYN7ftA74EzlIov8vJE4k76Sq569ypxZDpk6Myhsjt3t9f6uqMLFBXlCOg6k+gg1Glte/aEHTush9QPG37goVkP8WfONpZ/fQof3/kc/3viZGbOhJtvhp077RP49OnWzhCKpFHu3lmevLXcvbPS0mxdXO0WmzeXzWbjDKXii/XrocHeixl70ViSaiTxa86v3PvlvazfvT4s3mWKUllRMYkQvtLalrWTz8gAQRj+5sc8MecJDh0+xOGsLgxqOZpLOyUVu7FeeeURD6lQxpZy78w9de7u3lmhXPDnKaaZO06DftPkpoy7ZBxNk5qyZf8WHpj5APVPXqYLFBUlRKiYRAhPgQhdKUsnf+KJwmm3TeTZb15l5Srh4oa3sGrSPVx6UdUSbqwzZ5b0kLr9drs6PljcO3NPnbu7d1Yg3lqlkZER2JqQ5JrJjOkyhnOPP5cD+QdIvOpRHn39u+LynvLJ33AD/POfoamvolRovM1/VeStPGwmoZ6fzz+cL2PnjpULJzqk6s3d5cI+30nDhlIcYt5pg/DkkRRobClvdh53T7LSvLnC4R3lKaZZad5qhYWFMnHxRHvvejvkiWnvy2efFRbblbKyrC1m6lQbyywxMXR2JkWJZdBAj+UvJv4ai/3p5A/kHZBHZj8ijkyHXDf1Opnw/rLiDtUZ0NBXhxqIaLlG+XV2sllZR4z5jz1WsjOfPt12vuecY1+nTw+dK7I33CMV+2vQ/2jVR0XBIh1Std0kee+9Qo9i5K8IhtK5QlGiERWTKBCTUI1M9uTukXu/uFccmQ656YObJHtHdvFnhwyx0YON8d2hDh1qj5WGv+swZs8u2ZknJ4uccYY9Foy3ViT4bt13cvyDV0q9Ox1Sp/N4iYs/XFzf2bOPiIMx9t56E4fSRFdHNkpFQMUkCsTE1wI7J6V18rsO7pJBnw8SR6ZDbvvoNtm4Z+NRZUK5EC8UdY4FElv9KJe+drU4Mh3y1NynJP9w/lHi8NtvViQ9iYMuflQqCyom5SQm7tMe7omiXCmtw9l+YLvc+emd4si0uUdy9ud4vW5Z7AieqCzrMOLiRH7ZuFyunXqtTQ/80WhJbpBfQhxcpyDdv6vKIrqK4ktM1JsrTHhaU/Lqq9YF+NRT4fXX/Q8guOPgDobPHs5fe/4irW4aT170JCk1vfuresqNUpbYUpVlHUZKCtQ52IonOj9Braq1+PCn+bS8Yyztzi4oLuPqIuzueZeZCRdddLQnmGsq4lB50ClKtKJiEga+/hquuQYOH4axY6FDB3jgATj/fBvRt1s327kkJJTeyTuFZOPejZxY/0TGXDSGxBqJgGdXVmcH5h7QsCwLA/1dFBjr6zCcLs3Nkpvx+IWPs2ldLeIbL+TJuU9SUGgFxd2l2VUccnJsvntfi1Ergugqii9UTELMjBlwxRVw2mlHwp27diyrV9tQ6A89BPfc47uT33lwJw/PfpiNezdyQv0TGNV5FHWr1y2+TihX03sikEWBrvgSufLEW726dj2yXqVpclPyvniC1Hq1+WHjDzw17ynmzCs4KoeJUxycbXrpJd+LUefMiX3RVRSfeJv/Ks8NSAKmAfuBP4GeXsqNBPKBfS7biaWdP1w2E6chNjGxdLtIabaGXQd3FdtIMj7LKI4n5XqdcBt8y3KdaPVqCsTFOTlZ5KufsuWKST0kfYRDal48Vj797HCJ8zm/v8GDRdq1820zGTrUllGbiRLrEGsGeGAK8B5QGzgX2A209FBuJPB2oOcPl5g4DbG+1pQ4DbG+1pTsO7RP7v78bnFkOuSuT++SXQd3ebyOL0Jl8A3EmB+tXk2Bujg73auTm/0uTYdfJxdOdMi4+eOksLCw+DPO++stDYArztD76s2lRIJwrneKKTEBagF5QDOXfW8BT3ooG1Vi4vR+8uUF5Xyi9TYyOZh/UO7/8n5xZDqk3/R+suPADq/X8UUovaz8XRQYrV5NgdbLVXxWbl1Z7OX10o8vSWFhYQlRdD44+BLd5GRbTlHCTbhnBmJNTE4HDrjtux/4xEPZkUWjlh3ACuBOH+ftDywCFqWlpQV3R73g7Fh8dV7OEYmnTvVQwSF5ePbD4sh0yK0f3Spb9231eR1f5OXZcpFckR2trsRlqZerOHz601K5MvMquXCiQxzD35TERJHu3e15wYqFM7yMJ9GdPTv23aeV6CcSMwO+xCQaDfC1gT1u+3YDdTyUnQo0BxoA/YB/GWNu9HRSEXlFRNqKSNsGDRqEsr7FOL2ffAUgXL8e6tblKINuoRTy7PxnWbplKYkJidx+4ihGP9LAoxHbHy+rt96y6XPDaaB3J1pdictSL1f36r6O1nzy0EPMnxvHTwenkpc+jebN7f28+2647jp7X3v0gIsvtk4Vv/8ON94I77xj3Yb37o0OJwSl4hKOYLIB4U1lymvD88jkPjyMTDx8dijwQWnlwm0zEfE+7XHOOTZ4oOtws7CwUF5c+KI4Mh3S4/96yKRpa30OVa+80ve0TVaWvUaPHp6Ph8t2Ecsjk9mzra3E1yjurblfS9XeDuk43iEzs2eKSMmnQed9fe21I1MNU6fakcvs2eXvhKBUbCLx/0eMTXM5bSZNXfa9iQebiYfPPgR8WFq5cHtzuUbRdZ32SEoSSUiwHYsrby55UxyZDrn6vatlxqIVpQ5VExPt5q1Mjx5WTHz9sMJhu4hVm8nnn9v71a6d73nmwYNFrn14ujgyHdJtSjf5ccOPxZ93PjjceKP9jqdOPdpJISvLfjfGaCBIJfSEMpisN2JKTGx9eRfr0VUL6Ih3b64rgUTAAGcBG4HepZ0/nOFUAg1lMv23kp2Tvx1y9+7er2OMfTr2RThGCLHmzZWVJdKnj71frrYP1/q51tn55Ocq/qtyVhWfa8gQKyRwtJOCq2H0jjtE7rknOlymlYqDjkw8i0kS8BF2ncl6itaZAOcB+1zKTQG2Y9eX/Abc7c/5wyEm7u54NWvaLS7Ou/fTvPXzpGtmV3FkOmRW9iwRCewHkZUlcuut9jrWQmKfsI2xgQl9EewTirMe7gb+K6+0o6Zg44KFGneRd4bK/8c/RKpXF/nkE++du3M05XzyKywslBcWviCOTIfc8P4Nsn7X+uKyDRocfV89jVhd/6E1EKQSCiIxMxBzYhLuLdRiUhZ3vF+3/CpXvXuVODId8t6v7xXvD2So6u26NWuWntAp2CcUX212ejsFml8knDhHIa7CW72652lH987dea9chb7gcIH8+9t/F0dw3n5gu4gcCf/vivs/uSch10CQSrCUtzdXuXfs5bGFUkzK8gWu37VeevxfjxJrF5z4OzJJSfF+3cGDRXr18v3DCabzitbpLG+4C9/dd4sMGGDfgx2VuON6f5ydv7so5Obnyn1f3ieOTIfcM+MeOZh/UGrWtKNFV9y/U09CXhGiLyvlT6gihntDxSSMYhLo0HLnwZ1y20e3iSPTIY9/97gcLiwZpsPf851+uvdyzs7+lls8C0awnX20Gto94Un4XDt3pzOD+71w7dxdpxXdz7Xr4C7pN72fHaG8+agk1CyQAQNKnst9tOnp3oRi2lFRRPxfZFwWVEzCKCaB2Dhy83NlyBdDxJHpkHu/uFdy83M9lvXnqT8pyfd1P//cdpI1a4b+CSVaXYA94Un4XDv3wYNF2rf33bm7dv6envzmL9skLUb0lKq9HdLz+f9IckqhV/HyJuTRcr8UxRcqJmEUE79Xo8cfllHfjRJHpkP6ftxXdh7c6bW8UwjatrWiYYx9bdv2iC3En+uuWnVkDj+UTyiRcEEMFZ6EzzYV7fcAACAASURBVHVfVpa9p0lJJcs4O3dPnb+nJ7/e9/1WnK3x4bc/KCE4zmk1X0IeLSM5RfGFLzGJxhXwMYW/OT9qnzuJhRsXUqtqLUZeMJL6CfVLPbcxEBd3ZDMmsOtWrQqpqcHlNPFELOU58bT63TW0fno6TJ4MO3bYJGXZ2TZp2dNP21XxnpKWecoV88YzJzPionsBWBo3iZc+mV+cnGz8ePjf/2DtWnjvPfjqq5Jh8Hv0sMddIyIoSqyhYhIk/uT8GPb6DJLP/YgqcVV4+LyHOa7ucV7LZmfb/BeffQY//mgTLxUU2Ncff7T7e/WCf/6zbLlGQkFZ85yUB56Ezz3cTYsWtpyz809IgFdegVatAstM2TGtI73b9Abg3fXPcsfw39m82SZJ++wzG8LG4YDcXPjuO1i5Eq66Cj75xJZZvTqEDVeUSONtyFKRt0h6c0364mep1qubXDjxyFoSX/hr3L711vLzqIolby5v99PV9uG6iDBYm5JraJybP7y5OFinczqtRw/P047RdM8UxRuozSR8YiLi3R3vrqHrpdpN18tZzztk8pLJfp0rEON2uN0AfVGe1w4EX8LnHt4kVDal/MP5xdGfMz7LkAN5B6LKA87TYlMN66L4gy8x0WmuEOAaYbZjRxtBtv0Fe/g+7jE6XnCAq87syC2tb/HrXIFEuPV03dJyyoeK8rx2IKSnW5tHt24lbSLZ2XYqbvZsOwV1+HDobEpV4qow9NyhHFvnWNbtXscz85/hncxC+vb1/TnXvPLhIhLpnpXKibFiU7lo27atLFq0KGznLygs4JGvH2FFzgpOSjyJJy96kupVqvv12YYN7T+3rw4tO9t23Js3h6jClYDsbBt6OzPTCnFKirXpDBwYvHh4Y9PeTdz31X1s2bmPb166itSNt5W4dkZGyWvn59tOvaAgPPXJzraCMX265zDlCxZY0V24MHz3RIltjDGLRaStp2M6MgkxIsKEHyewImcFSTWSeKTTI34LCcSWcTuW8OSBFYpRiC+OqXMM5zGceXPiiWs5jQdemu1zJBBuD7hyz3ehVGj8FhNjTJoxrs6pxfuNMSYttNWKXT7+/WNmrZ1FtfhqPHLeIyTXTA7o874Sa4Hd755YS4lOsrPh0QGn8uS1d9K4MUxcOp7VO1eSng6jR9sRQq9eRxJmheohITvbJuJyT6z25ptExVSbEnq8feeRTMYWyMhkLTajoTtJRccqPYs3Leb1X14HYMg5Q2ia3DTgc/ia4x82zPO6ByU6cY4EhnS9lL7ndmPt+gKGTBvN1v1bgZIjgVA9JPiyiezYAatW+f58eWTCVIIjauxg3izz7htQCDTwsL8xsN/f80TDFo4Q9Ot3rZfr/+96cWQ6JHNZZtDnC2d8HSUyuEcZ7vnKv6Rqb4ec8UiGrPj9oOTl2YjFtWqFxgOuNJdtb3HI3M+hYV1ih0i76ePDm6tUA7wx5sWiPwcCk4ADLofjsUmp8kSkY2hlLnyE2gC/99Be7vvqPv7e9zcdj+/IQx0fwsOMoFLJiI+3T4lVqtj3+/P20+/9+1i4ciNbfz6Hg18NJyXZkJMDa9YEP9ocMsQ+jY4e7f34Dz/Yp9Vx4zyXGTbM1tnbcSW6KO07h9B+p74M8P6IyTdFf54PLMCm1HWSB6wDnhGRNcFXNTKEUkwOFx7m0W8fZemWpZxY/0TGXjyWhCoJITm3Ett48szbuGcj9311H/vz93NDyxs4p9ZNIfPMK80TMDsb2rWzYXm2bz/6uHpzxR6R9v4MyptLRC4UkQuBycDlzvdF26UickcsCUmoee2X11i6ZSn1qtfjkU6PqJAoxXjyzDu27rF25Irh3RXv8uhrc0PmmVfaGiVvccjUHhe7BLIuLdz4bYAXkVtFZE84KxNrzMyeySerPymOudWglif/BKWy4s0z7/RGp9P39L7s3AH/t+E5uvb6IyTX8ycAp3scsmhdbKr4RzQFXdV1JmVkVc4qXlr0EgB3tb2L5g2al3ONlGjDl2fegsndWDKtC6edmceE5Y9z5727gnbr9HeN0i23RH7NjRIeomldmopJGcjZn8MTc56goLCAbs26cXH6xeVdJSVK8RZ2Ju+QYdHLAzk+4RSmz97GL7VG8/3cgqDcOnWNUuUjqr5zb25e5blh165MA/YDfwI9vZQzwFhge9E2liKnAl9bMK7Bufm5cvfnd4sj0yGPzH5ECg4XlPlcSuUmK0sk6dgd8s9Xeosj0yH/+eE/UlhYWHy8LG6dsRKAUwkdkfzOicFAjxOwnmINgZuAl40xLT2U6w90B9oArYGuwB3hqpSI8OIPL/LHrj9oVLsRD3Z8kPi4+HBdTqngjB8Pd/RK5PmrH6ZafDW+zP6Sz9d8Xny8LOFNYiUApxI6AvnO9+XtC1s9oi7QozGmFrATaCUiq4v2vQVsFJGhbmXnA2+IyCtF7/sC/UTkHF/XKKtr8Psr32fy0skkVEng2UueJa2eRpFRyo6rW+d3677jmQXPEGfiGHXhKE5teCqgQT2V0PHX7r94YOYD9GjZg+6ndC/TWrhYC/TYDChwCkkRSwFPI5OWRcdKK4cxpr8xZpExZlFOTk7AlTpceJifNv4EwP3t71chUYLG1a3z/Cbnc03zayiUQsbMHcOWfVsADW+ihIZ9eft4/PvH2Z+/nzU7wrOSIxrFpDbg7oK8G6jjpexut3K1PQWkFJFXRKStiLRt0CBwF974uHhGdR7Fvzr9i7OPOzvgzyuxQ6SC5rm7dfZq04u2jdqyN28vo74fRW5BbsTcOpWKy+HCw4ydO5a/9/3NifVP5J6z7wlLhI5oFJN9QF23fXWBvX6UrQvskzDN3VWNr0q7Y9uF49RKlBDJoHnubp1xJo77O9xfnFTruQXPMfFV0XQDSlBMWjKJJVuWFC+sDiQlRiBEo5isBqoYY1xD7rYBVngou6LoWGnlFKVUsrNtSPjp022so/R0G1fLW8j4YPHk1lmrWi1GdBpBraq1+GzZfF767l115VXKzOw/ZvPx7x9TJa4Kw84dFtaF1VEnJiKyH/gQ+LcxppYxpiNwJfCWh+JvAvcaY441xhwD3Ae8EYl6RkP+ACW0RDp5lLdFjbk5x1Lj5wf46SdD8x6ZbKk2PzQXLEJ/u5WD37f9zvifxgNwx5l30DLVozk5ZESdmBRxF1AD2ApMAe4UkRXGmPOMMa6+bf8DPgGWA78CnxXtCytRkz9ACSmZmZFPHuXNrTMx90ye730rqakwbsE41u4MTcog/e1WDrYd2Fa8sPqKpldw2UmXhf2aUecaHAmCiRqsebQrLu4h4z0R7jztrogIzy18jm/WfUNqzVTGXTqOegn1ynw+/e1WDvIO5/HQzIfI2plF69TWPHbhY1SJ8/GjDoBYcw2OajSPdsUlmoLmARhjyDgrg2ZJzdh6YCtj5o6hoLDsKqa/3YqPiPCvz17g47lZzPqoIWMcD3FsoyoRmcZUMQmQ8pgKUSJDNAXNc1ItvhoPd3qYpBpJrMhZwcs/vewMJRSw7UN/uxWf4e/8H+M+/J7qcQl8MfwR8vbWjdg0popJgERT/gAltERV0DwXkmok8ch5j1Atvhpf/fEVn67+1G/bh6vgbN0KHTr4Fhz97cYuH/zwA8/OeouzzzK80f9+zm/TJKzeiO6omARItE2FKKHDV8j48k4e1TS5KfecfQ8Az8+ZSM/7fi7Vhfn110sKToMGMGWK76dU/e1GF/6OPtftWsew6c/QOA3uuaCXx4XV4Z7GVDEJkGicClFCRzQHSuzUuBM9WvZg7Tqh0bVPcXzLDR7LtW8P11xjR1qugnPTTTBrlu+nVP3tRg/+jj535+7m8e8e589NuVzb7gKuaX6N13OGdRrTWzjhirwFE4I+K8uGdZ4/3/PxsoQNVxR/KSwslFqXPSEXTnRIv+n9ZE/uHo/l+vQRqVmz5D733+7QoSJDhhw5rr/d6MHffmbl74fkga8eEEemQ7j0Ptmfe8jnefPyROLjy14vYjAEfdQSzVMhSsXHGMOB2ffS8pgT+Xvf3zw590mPHl6ffw65uSX3uf92L77YPqXqbzf68Mfzru/twoDXx7Nq2ypSaqaQsvJh/t5Qzed5wzmNqWJSBqJ5KkSp+DRITKBX4xEkJiSybOuyEh5eTrZt89xpuP52b7wRtmwJ/LerK+jDjz+ed8de/D4Lt3xD9fjqjOg0gpuvSSzXKXgVkzKSnq55tJXyoWdPmPZOCiM6jSj28Pr4949LlElIgCuu8Px55293/nwrBIH8dnUFfWQozWt03vp5fPH3m+TlGe7vcD8nJp5Y/t6I3ua/KvIWjM1EUcob1/n0OX/OEUemQ7pmdpWFfy2UrCyRHj1EQMQYkdRUkcGDPdtB3G0mgVzXE2pzCR2pqd7v4+/bfper37taLpzokHrnfFDiWLhT+KI2E0WpOLjaPj57+VwuTr2Zw4XCwMynaXleFtOnw5gxkJRk5949jRrK8pSqK+gjhzev0Zz9OTz+/ePkHc4jft0l3Nr+qhLHy3MKXmNzKUqMkp1tO+53MoWcE19AmsymScMkpt7+LO1apjBjhnX/vf12OOMMuOsueO89mDnTCsmbbwbWubimGfZVJ00zHDye4qjtz9vP0FlDWbd7HUn5rfns/sf4YUGViE6ta2wuRamAOG0fWzYbBp2dQfv0U2l11g7e/PMxDuQfKPGUOmgQbN8OXbuW/SlVoz9EDnfPu9/XFPDE90+yYtM6Nq48js+GDuOtyZEVktJQMVGUCsC7mVV4+aZhxVkaH5w+hnuGFNChA7zwAohA79522qOsjiIVMfpDsJ5p4fRscz4M5B4S2t31EmMmLWHBN/Vod2AkP8ypzeWXR5dnnYqJolQAtm2DlifVYeQFI9m/vR6vfLKEFbXHM2+eFHtcpaba0UlZPa4qWvSHYD3TIuHZlp4O59w2lfNvm8nVV1Zj2YuP8r9nG5KeHoWedd4s8xV5U28upaLh9P7JyhJJPOl3ueClq8WR6ZB3lr1TXCYrSyQxseweVxXJmyvYtkTqXszKnlXCWy/S13cH9eZSlIqNc9QwfjwMuK4ZT3Z9EINhyq9TmJk9E7Cjhj59yu5xVZGiPwTrmRYJz7bFmxbz4o8vAtD/zP4lgjeWdv3UVGjSBFq3juD0lzeVqcibjkyUiobzSTUx8cjT6GerPyt+qn11xg/FT6pZWSINGwZ3rSFD7Dni4+3rkCGxMSJx4msdhxNf9ynYz5fGmu1r5Nqp14oj0yFv/PJGQNd3rjUZMEAkOVkkP9+WHTYs+LUm+BiZqGuwolQQZsyAf/4THnwQ+ve3nlUvfvs2L895j7/WVePF7qO44+rmEU09HK0Em6I5nCmeN+3dxIMzH2T3od1c2ORChpwzBGOMX9d3dSlu2/bo6webmlldgxWlEnD55ZCcDDk5RxasPXXLTRyfdwkdO+Xxrfybv3b/FTUeV+XpiRSsZ1q4PNt2HNzBiK9HsPvQbk7/x+ncffbdRwmJr+u7Tn95un44F5ZGnZgYY5KMMdOMMfuNMX8aY7z6hhhjRhpj8o0x+1y2EyNZX0UpLzx1xmlpUL36kZhxWzYbZo2+i87NzmZf3j5GfDOCF17dWu4eV+XtiRSsZ1o4PNv25e3jX9/8i60HtnJy8skMP284VeI8D328Xd81QKS364ctp4m3+a/y2oApwHtAbeBcYDfQ0kvZkcDbgV5DbSZKrOOcFx82zM6HO+fFBwywMbmefbZk+UMFh2TozKHScbxDql/XTxav2FE+FZfo8AqLNm+ug/kHi/OS3PnpnV7z1JR2/bg4+1vwdf1gcprgw2ZS7uJRojJQC8gDmrnsewt40kt5FROl0lFaR/bss1ZQ7rijZLC/+4bul4Tu98hZzzsk47MM2Xtob0DXHDzYGn7j4nwHkCyNwYOtCPoi0CCUZSHYoIiBfN7X/TtUcEgenv2wODId0uejPpKzP6fM9U9Ott+7r/oH4xjgS0yibZqrGVAgIqtd9i0FWvr4TFdjzA5jzApjzJ3eChlj+htjFhljFuXk5ISqvooScUpzC733XmuA/+mnksH+Cg/VZMHoxzgj/TjW7V7HyG9HciD/QKnXC/WUlD+5OsKaXraIYIMi+vt5X/fv7PYF9J00lqVbllI/oT5PdH6ClJr+GVo8Xf/gQfj5Z9/1D9vCUm8qUx4bcB6w2W1fP+BbL+VbAMcA8UAH4G/gxtKuoyMTJZYJ1i01Z3+O3PbRbeLIdMiDXz0oB/MP+jxPqKeknFMxvvB3KiaUI6Zw4Ov+FRwukLveGStVezuk6xs3ytqda8N6PZHgpxCJlpGJMeZbY4x42eYC+4C6bh+rC+z1dD4RWSkim0TksIjMB14Arg1vKxSlfAk24GJKzRRGdxlNSs0UVm5byePf2ZDmrjiN+61b2xAs3bt79rQqi3dQIJ5Qvjy+ytuI7w/eRpGFUsjzC59nvZlDeloNkpc/RpP6TYK+XnkuLI2omIjIBSJivGznAquBKsaYpi4fawOs8PcSwNF+dIpSgQiFW2rD2g15ovMTJNVIYtnWZSUExbWTTkiAVat8d9KBTkn56wnVvr13sWjXDm64wa6nGD3ado5VqtjX0aPt/l69Ihvw0JPwvfIKXHRRyWNx8YXUu/R5Jnz+LQW5CTzf/TFmvN209Av4SbnlNPE2ZCmvDXgX69FVC+iIb2+uK4FErICcBWwEepd2DZ3mUmKZUBqw1+9aLzd/eLM4Mh1y1/sPyy19csWYI1kajRH57bcj5T1Nk3iakvI1/eTPVExiot28lenRQ6RGDd/TNeE24ru20XnP2rUTmT37iHedMSJ16ojUrWu/s9VrDssz856TCyc65JQR10r9pitk+vSye1dFGmLFm8vWlSTgI2A/sB7o6XLsPGCfy/spwHbs9NhvwN3+XEPFRIllQj0vvn7Xeun8n1ukam+HNLp5mNw+4GBxZ1izpu3UXT2D3Dtpd/uMN7dl13AepXlCde/uWzBTU63Xki+xCDZsjC9c2zh7tv176tSjQ5YkJ4vUri1Sv77Ib6vz5am5T4kj0yHXvHeNLN+yvFg4U1LCU89QE1NiEolNxUSJdUKZ6zsrSySx8Qa5/H+9pFofh/Sf+qDsO7RPROyTd69eJcXJvZN2FZdAhM5XjK/SnAzi4uyIyZdYBLOewhfubXQfKbq28fTTRdq3F3lwaL50HPqEODIdcu3Ua2X5luXF5c85R+SMM0Jfz3CgYqJiolRAQhVw0dkZbti9QejeW6542yH3zLhHdh3cVdxx3nLLEcFw7aTdR0GhmoIrzeMrNVVk1SrfYhGukYl7Gz0Jn7ONSUki9VNy5boJI6VaH4f0+L8e8lvOkXlD58gkOTn09QwHKiYqJoriFdfOMKXxFunxdj9xZDpkwCcDJGd/jnz+ue3wata05VatOjIqch8FhSqabmnnGTzYrvb3dZ5w2Uzc6+ZJ+JxtNNX3yjUvPyhVezuEa26Ur3/JPmoUWVFsJtG2aFFRlAjj6mp881WpNFoxlsb1GrNh7wYemPkArTr8xfz5dkFcx47QogUcOODZOyhUeeJL8/jKyIDJk+GCCzwfX7DAeoQNHOj7OmXBvY2evOvS0iBn3w4SrhjGjior6XxOCglzx3LjZSce5V3VokV0BN4MFhUTRankuHaGGRnw9sREutcZwynJp7DtwDYemPkAa/evIDUVpk2zkYmXLfOcSz5U0XQzMmDiRCsKnti6FapWha++ivx6Cvc2ehK+H1ZtoJrjAf7RfB271h/LyX8+xZ09jy8OwLl585H7F0upjn2hYqIolRzXztC56O3Gq+vAN6NoVuMc9h7az13vjyCxzbxSO+lQRdP1Z/Hdu+/akDGRXk/h3kZ34Vu+ZTn937+fRidt5YJTm5L936d4878NPI6SwjmCijje5r8q8qY2E0U5gicPLKdxP7XhYTFnvSTc6JDjH3TIC7Pfk8LCwoDO5UqgbsvRmNXRUxud3nU3DJ8l7Z7tLlV7O+TWSY/L/UMPSt26dp1JKDzvyhvUAK9ioii+8OVqnJxSKA+//YF0zewqjkyHPDX3KTlUcKhM54q1ztMb7m3MPXRYHnjndal3p0O40SHmjFclteHhYuGLRlEsC77ERNP2KooC2CmkCRNsaJRt26xtoGdPOwWTng4/bvyRp+c/TW5BLk2TmjLs3GE0qNWgTOeqCDjb+Pb/7WXbyU9R9fglHHdsHMMvuYO+nf5Z3tULC77S9qqYKEolJDvbBiF07+wzMnx39ut2rePx7x5n64Gt1KlWh/s73M8Zjc6IXMXDRFnvxx87/2DMnDFs3r+ZutXrMuzcYbRKbRW5ikcYzQGvKEoxwUTbPby9CQnfPs/X75zJlA/30n7YSC5/cAprsgoj14AQU5b7ISJ8tvoz7vvqPjbv30x6YjrPX/p8mYTEV2TkmMLb/FdF3tRmolRWgjGQu8ajWrOmUN5eMkU6T+wqJ41wSI2uD0rmx1vCW/kwUJb7sffQXhn9/WhxZDrEkemQCT9O8GlD8oU/ccyiCdRmUhKd5lIqK0OG2Cfu0aO9lxk2zD6djxt3ZF92tn1Knz69ZG6OJZuX8NzC58jeuIOfFtTk9TvvpGeH8zEmNjJBBHo/Fm9azIs/vsiOgzuoWbUmg84axLlp55bp2t7uqZMFC6wL9MKF0WNnUpuJGyomSmWlYUM7heOrc8rOtus1Nm8+ss9Xp7vn0B7+88N/mDRrIYWFcNslZ3Nn2ztJrpkckjqX1Z7hD/7ejw7nH2DEh6/zZfaXAJySfAr3tr+XA1sblbluZRX28kRtJoqiAGUPd+Irb3vOhrrkTBvO1k/uZu3vNXl80g+0f+IuXp/zBcE+rIY7m2Jp90NEWM8ccs66ky+zv6RKXBVuPe1Wxl48liVzGwVVN1/31EmgicfKFW/zXxV5U5uJUlkpayBGb1F8Xef8V60SMbW2yZBpj8tJIxxStbdDur88WFZuXVmmuoY7n7mI7/vx1+6/ZMTXI+TCiQ6p1sch9315n6zbuS5kdSstMrJI+MLolxU00KOiKFD2cCeeYm5lZ9vUuM7UuVWrQmrtZJ698mEm3f4QF3dM5rOFWQz8+EGenvc0f+/9O6C6esuf7qQs+efdcb0fTq+qBmk7iDt7Aul3D+Tt2b+wPqs2VzTI4OmLn6Zx/cYhq1uo4phFC2ozUZRKRFmNvp7m9933uc/v5xbkcv2jH7DKfMApLfOJM3F0OaELN7S6gdRaqaXWtaz2nUBw3o9hw+CJcTtp32c6BemfUKXaIXJz49j7yyX8+NrNvDahHrfdFtq6VTSbiYqJolQyZsywI4rbb7dbWpp9An71Vbu9+ebRQRI9iZBrh+pNhLKzof1FW3l4yhS+Xvc1hVJIvInn3LRzufLkK2ma3NRrPePjbUdapYr3tuTn2w65oKDs92PcxL+5b+JHHH/+TE46OZ8aNeCUWueQ+0Nvpr56HMOGwZgxJdsWirqpN1cFQMVEqeyUJdyJuwg1awYrV8Ibb3gXIdcOddPeTbz767t8u+5bBNvvnJJ8CpekX0LHtI7UrFqzxGfDOTIpKCxg4YaFfJH1BZlfL6XwMMTFw9afz+HgoutoENesxP1wHyGEqm5lEfbyJGbExBiTAfQBTgWmiEifUsoPAR4CagLvA3eKyKHSrqNioihlw1WEtmyxuU169fIuQp461K37t/Lp6k/5Kvsr9ufvB6BafDXOPvZs2h/XntMbnU7tarVLnQbKzoYbboAVK2xHX5pbbt7hPJZsXsK89fP4YeMPxdeeOaMa919/Hv3Pu5q0emler+XajlBOUcVSHLNYEpOrgULgUqCGLzExxlwKvAl0BjYB04CFIjK0tOuomChK8ATboeYW5DLnzzl8vfZrfs35tXi/wXBKyik0kJY8ce/JvPdSMy45L7HEQsgZM+DGGyEvDz79FDp1ssbs116zuUUmTxbaX7iLP3f/yaqcVSzfupzftv1GfmF+8Tma1GvCpSddSrfWF5K3r1ZAU1axOEUVCmJGTJwYY0YBx5UiJpnAOhEZXvS+C/COiPyjtPOrmCiK/3hbNNi1K/ToEZoOdev+rcxbP49FmxaxImcFh+Ww3b8VfvkF0tNqcnaLf9C0USNy/qrHSy/WpKrUYvBgaN3mMAWFBezN28vOgztZ8cdO3v9qAx0u3EutWiWvk56YTofjO9Dx+I4cW/dYoOxTVrE2RRUKfImJDy2OeloCH7u8Xwo0NMYki8j2cqqTolQonB1mv362w23c+MgIoEcPGDTICoavDtWfJ/PUWqlc1fwqrmp+FQfyD/Dr1l/5fdvvrG64mgb11rAyaz9TvvyDvLw/iIuDOm3hjDNgWRwsW+52smpwXDpsXFuL7hc0Jj0pnVNTT6VVaivqVK9z1LWd7sG+Rlie3KUvv9wK5YQJVmhchbaijUj8IZZHJtnAQBH5ouh9VSAPOEFE1nko3x/oD5CWlnbmn6U5eCtKJcffqZz33rNTTeGa8xcR9hzaw+Z9m9m8bzNXXr+Xfz95gDqJ1uZRJa4K8XHx1Kpai6QaSSTWSORgzj9wdElmy+bSY4RV1imrshAV01zGmG+B870cnici57qU9UdMlgJPiMjUovfJwDYgpbSRiU5zKUrpROs6iHC4DFfGKauyEBWxuUTkAhExXrayhN1cAbRxed8G2KJTXIoSGqI1dlQ4Vo47p6wOHbJTVjVq2NdDh+x+FZLSiapwKsaYKsaYBCAeiDfGJBhjvD1/vAn0Nca0MMbUBx4B3ohQVRWlwlPWoJDhpqwhYUojPd2OsDZvtiOazZuPjLicyavi4qBmTahVK8YTWYWBqBITrCAcBIYCNxf9/QiAMSbNGLPPGJMGUGQreQr4BlgP/Ak8Wh6VVpSKSLTGjsrIsO6/CxZ4Pr5ggRWTgQODv5Zr1OIxY+y6mt697ZRYYqLdF6oIxrFOVBrgw43aTBSldKLVZgKRsXE4DfMvwSzUWQAACidJREFUvQSffw6TJ9v9DRrYUU+HDnDXXXYabOvWymGkjwqbiaIosUUkRwCBEgkbx/jx0LmzFYwVK+COO+wiSWfOkrvuggsvtK7BoYhgHPN4i01fkTfNZ6Io/uHMVzJ0qM3NkZdnX4cO9T9HeVaWyODBNndIXJx9HTw4uDwkkSA5WSQx0eYm8ZT3ZP58ezwlxb73lAemooHmM1EUxR1n/o6GDb0bk4MdAYQ7U2I42bED+ve3ow5PzgjO0cj2Iv/R8nBGiCZUTBSlEhJIJ+/Ny6k024B78qz0dLs2JD3dvp8+3R6PZk+oSy6xr96cES67DJwhw2IpkVU4UDFRlEpGpDr5SGRKDCci8NVX9m9v7shffGHLQdnckSsS6s2lKJWMSHlpRSJTYjhJTrZC8dlnkJp6dMiVBQvgiivs2pNPPlFvLhUTRalkRKqTj1SmxHAxZAhs2gRff21HUE2bwkMPwTXX2OPvvw9nnWXtJGvXVo6QK+oarChKMZFa2R6tix79JSPDCslLL1lRHD7cGuXfestuO3faaa5WrTTkCqiYKEqlI1KdfLjCnkSK9HQ72rjrLqheHebNg9xcO6XVooWdAjPGToONHx/djgSRQMVEUSoZkerko3nRo7+4u0YnJMBFF1kRmTUrtlydw463BSgVedNFi0plJivLLrSbP9/z8fnz7fFQLCoMxaLHaCGS9y1aQRctKorixDl9062b9drKzraG8Oxs+75bN/8zJJZGRQrtHuuuzuFGvbkUpZKSnW07vnBlSCxPvOWtz8goe9ti3dU5FKg3l6IoR1HWle3RTrhCuERrfpdowYcHuKIoSmzhurrfdTrKubq/a9eyLy50esH5+lw0uzqHGx2ZKIpSYQinXSPWXZ3DjdpMFEWpMITTruFMluU+6nGyYEHFD6miNhNFUcqMP6Hqo4Vw2jUi6QUXi6iYKIrilVjLRxLu1f0VydU51Og0l6IoHonFaZ1ozltfEYiZaS5jTIYxZpEx5pAx5o1SyvYxxhw2xuxz2S6ITE0VpeITi4v0KkIIl1glqsQE2ASMAl73s/wCEantsn0bvqopSuUiMxP69vVd5vbbbbloQe0a5UdUiYmIfCgiHwHby7suilLZidVFemrXKB9ifdHi6caYbcAO4C1gjIh4TLNjjOkP9AdIS0uLXA0VJUaJ5UV6ztX9aheJHFE1MgmQ74FWQCpwDXAj8IC3wiLyioi0FZG2DRo0iFAVFSV20UV6SiBETEyMMd8aY8TLNjfQ84nIHyKyVkQKRWQ58G/g2tDXXFEqJ2rMVgIhYmIiIheIiPGynRuKSwAmBOdRFAU1ZiuBEVXTXMaYKsaYBCAeiDfGJBhjPNp1jDGXG2MaFv19CjAC+DhytVWUio8asxV/iapFi8aYkcCjbrsfE5GRxpg0YCXQQkTWG2OeAW4BagNbgLeBx0Ukv7Tr6KJFRVGUwPG1aDGqxCRSqJgoiqIETsysgFcURVFiExUTRVEUJWgq5TSXMSYHKCW2qFdSgChb8xt2tM2VA21z5SCYNjcWEY8L9SqlmASDMWaRtznDioq2uXKgba4chKvNOs2lKIqiBI2KiaIoihI0KiaB80p5V6Ac0DZXDrTNlYOwtFltJoqiKErQ6MhEURRFCRoVE0VRFCVoVEwURVGUoFEx8YAxJskYM80Ys98Y86cxxmP6H2MZa4zZXrSNNcbEXBj8ANr7gDHmV2PMXmPMWmOM12Rk0Y6/bXYpX80Ys8oYsyFSdQw1gbTZGHOGMeZ7Y8w+Y8wWY8w9kaxrqAjgt13dGPPforbuMMZ8Yow5NtL1DQXGmAxjzCJjzCFjzBullB1ijNlsjNljjHndGFO9rNdVMfHMBCAPaAjcBLxsjGnpoVx/oDvQBmgNdAXuiFQlQ4i/7TVALyARuAzIMMbcELFahhZ/2+zkASAnEhULI3612RiTAnwB/A9IBk4CvopgPUOJv9/zPUB77P/xMcBO4D+RqmSI2QSMAl73VcgYcykwFOgCNAZOBB4r81VFRDeXDaiF/fE1c9n3FvCkh7Lzgf4u7/sCC8u7DeFqr4fPvgj8p7zbEO42AycAq4DLgQ3lXf9wtxkYDbxV3nWOcJtfBp5yeX8F8Ht5tyHI9o8C3vBxPBMY7fK+C7C5rNfTkcnRNAMKRGS1y76lgKenmZZFx0orF80E0t5iiqbzzgNWhLFu4SLQNv8HGA4cDHfFwkggbT4H2GGMmW+M2Vo05ZMWkVqGlkDa/BrQ0RhzjDGmJnYUMyMCdSxPPPVfDY0xyWU5mYrJ0dQG9rjt2w3U8VJ2t1u52jFmNwmkva6MxP5+JoWhTuHG7zYbY64C4kVkWiQqFkYC+Z6PA3pjp37SgLXAlLDWLjwE0uY1wF/AxqLPNAf+HdbalT+e+i8o/X/fIyomR7MPqOu2ry6w14+ydYF9UjRmjBECaS9gDXxY28kVInIojHULF3612RhTC3gKuDtC9QongXzPB4FpIvKTiORi59E7GGPqhbmOoSaQNk8AqmNtRLWAD6n4IxNP/Rf4+N/3hYrJ0awGqhhjmrrsa4Pn6ZwVRcdKKxfNBNJejDG3UWS0E5FY9Wzyt81NgSbAHGPMZmwH06jI+6VJBOoZSgL5npcBrg9EsfRw5EogbT4Na1/YUfSA9B/grCJnhIqKp/5ri4hsL9PZyttIFI0b8C52WF8L6Igd/rX0UG4A1jB7LNYDZAUwoLzrH8b23gRsBpqXd50j0WagCvAPl+1qrKfMP7BTX+XejjB9z52x3kynAVWB54A55V3/MLd5EvABUK+ozcOBjeVd/zK2uQqQAIzBOhwkAFU8lLus6P+5BVAf+Bo/HG+8Xre8Gx6NG5AEfATsB9YDPYv2n4edxnKWM9hpkB1F21MUxTuLpS2A9q4F8rHDY+f23/Kufzjb7PaZC4hRb65A2wzcibUf7AQ+AY4v7/qHs83Y6a13gK3ALmAucFZ517+MbR6JHU26biOx9q99QJpL2XuBLVg70SSgelmvq4EeFUVRlKBRm4miKIoSNComiqIoStComCiKoihBo2KiKIqiBI2KiaIoihI0KiaKoihK0KiYKIqiKEGjYqIoiqIEjYqJoiiKEjQqJopSzhhjGhhj/jbGPOqyr7UxJtcYc1151k1R/EXDqShKFFCUQvUT4HxgCbAI+FFEbi3XiimKn6iYKEqUYIx5HugGfIcNRHiaiOwr31opin+omChKlGCMqY5NndoU6CAiP5RzlRTFb9RmoijRQxPgeGzI8BPLtyqKEhg6MlGUKMAYUxVYiM0O+APwKNBGRNaXa8UUxU9UTBQlCjDGPAn0BFpjswHOwGbI6ywiheVZN0XxB53mUpRyxhhzPnAf0EtEdol9wuuDTaf6UHnWTVH8RUcmiqIoStDoyERRFEUJGhUTRVEUJWhUTBRFUZSgUTFRFEVRgkbFRFEURQkaFRNFURQlaFRMFEVRlKBRMVEURVGC5v8BYD0FCpVmYx0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(123456)\n",
    "\n",
    "# Nuevo set de datos (esta vez, hago dos copias)\n",
    "x_ = np.random.rand(100, 1)\n",
    "\n",
    "# Si quieren elegir a mano los valores de x\n",
    "# x_ = np.array([0.2, 0.8]).reshape((-1, 1))\n",
    "\n",
    "t_ = np.sin(2*np.pi*x_) + np.random.randn(len(x_), 1) * 0.3\n",
    "\n",
    "plot_data_sine(x_, t_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 418,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fijemos los hiperparámetros\n",
    "beta = 1/0.3\n",
    "gamma = 1e-2 / 0.3 # gamma/beta = lambda \n",
    "\n",
    "# Number of points to use\n",
    "n = 2\n",
    "ii = np.arange(len(x_))\n",
    "np.random.shuffle(ii)\n",
    "ii = ii[:n]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Elección de las funciones de base**\n",
    "\n",
    "Una parte importante del modelado es elegir las funciones de base con las que vamos a intentar describir los datos.\n",
    "\n",
    "Acá dejo dos opciones:\n",
    "\n",
    "* funciones de base polinomiales, como las que ya vimos varias veces.\n",
    "* funciones de base normales (ver ec. 3.4 del Bishop)\n",
    "\n",
    "***\n",
    "**Ejercicio**. Pueden divertirse codeando nuevas funciones de base (¡les prometo que es divertido!), como la sigmoidal (ver ec. 3.5 de Bishop), o la propuesta de Manu de funciones de base periódicas. Piensen y evalúen cuántos hiperparámetros tiene cada elección de funciones de base y qué característica del modelo controla.\n",
    "***\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 412,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dimensiones de la matriz de diseño: (2, 5)\n"
     ]
    }
   ],
   "source": [
    "## Features polinomiales\n",
    "\n",
    "# Escribamos la matriz de diseño para polinomios de grado 5\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "\n",
    "polfeat = PolynomialFeatures(degree=4)\n",
    "phi = polfeat.fit_transform(x_[ii])\n",
    "\n",
    "# Necesito también la \"matriz de diseño\" del vector xx\n",
    "phi_ = polfeat.transform(xx)\n",
    "\n",
    "print('Dimensiones de la matriz de diseño: {}'.format(phi.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 421,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dimensiones de la matriz de diseño: (2, 11)\n"
     ]
    }
   ],
   "source": [
    "## Gaussian features\n",
    "\n",
    "mu = np.linspace(0.1, 0.9, 10)\n",
    "s2 = 0.05\n",
    "phi = np.ones((len(x_[ii]), len(mu)+1))\n",
    "phi_ = np.ones((len(xx), len(mu)+1))\n",
    "\n",
    "for j, m in enumerate(mu):\n",
    "    phi[:, j+1] = np.exp(-(x_[ii].flatten() - m)**2 / s2)\n",
    "    phi_[:, j+1] = np.exp(-(xx.flatten() - m)**2 / s2)\n",
    "    \n",
    "print('Dimensiones de la matriz de diseño: {}'.format(phi.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recordamos los parámetros para la posterior de los parámetros $\\om$:\n",
    "\n",
    "$$\n",
    "\\boxed{\\begin{array}{lll}\n",
    "\\mathbf{m}_N &=& \\beta\\;\\mathbf{S}_N \\boldsymbol{\\Phi}\\T \\b{t}\\\\\n",
    "\\mathbf{S}^{-1}_N &=& \\gamma\\;\\mathbf{I} + \\beta\\;\\boldsymbol{\\Phi}\\T \\boldsymbol{\\Phi}\\;\\;.\n",
    "\\end{array}}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 422,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Los valores medios de la posterior de los parámetros son: \n",
      "[[ 0.08289836  0.21117203  0.18733331  0.10130547 -0.00553241 -0.09060352\n",
      "  -0.12698749 -0.11196624 -0.06991266 -0.03162368 -0.01041347]]\n"
     ]
    }
   ],
   "source": [
    "# Escribamos la ecuación para la inversa de la covarianza del posterior\n",
    "sn_inv = gamma*np.eye(phi.shape[1]) + beta * np.dot(phi.T, phi)\n",
    "\n",
    "# Calculemos phiT t\n",
    "tt_ = np.dot(phi.T, t_[ii])\n",
    "\n",
    "# Y ahora usemos solve para resolver la ecuación sn_inv * x = tt_ --> x = sn * tt_\n",
    "m_n = beta * np.linalg.solve(sn_inv, tt_)\n",
    "\n",
    "print('Los valores medios de la posterior de los parámetros son: \\n{}'.format(m_n.T))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 423,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-3.0, 3.0)"
      ]
     },
     "execution_count": 423,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAESCAYAAAAWtRmOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXSV1b3/8ffOHEIGE8JMGGSUAAEiowrFCRwQ9daqOOBVubrEWutw7V13rcu9q7/+qqL92Ypaq0WtqLfVoojaOiAqQy0JIBJGgQTCmIGEzCfD/v2xkxA4BJKQk3OSfF5rPSvk5DnnfDeB53P23s/zbGOtRUREpKEgfxcgIiKBR+EgIiJeFA4iIuJF4SAiIl4UDiIi4kXhICIiXhQOIiLiJSDDwRjzpjHmkDHmuDFmpzHmHn/XJCLSmZhAvAjOGDMS+MFaW2GMGQ6sAq621qb7tzIRkc4hIHsO1toMa21F3be12/l+LElEpFMJ8XcBjTHGvADMAyKBjcDHp9lnPjAfICoqavzw4cPbskQRkXYvPT0911qbeOrjATmsVMcYEwxMBqYDT1prKxvbNzU11aalpbVVaSIiHYIxJt1am3rq4wE5rFTHWlttrV0N9AXu93c9IiKdRUCHQwMhaM5BRKTNBFw4GGO6G2NuNsZ0NcYEG2OuBG4BvvB3bSIinUUgTkhb3BDSS7jwygJ+Zq1d7teqRKRZKisryc7Opry83N+lCBAREUHfvn0JDQ1t0v4BFw7W2hxgmr/rEJFzk52dTXR0NAMGDMAY4+9yOjVrLXl5eWRnZzNw4MAmPSfghpVEpGMoLy8nISFBwRAAjDEkJCQ0qxencBARn1EwBI7m/i4UDiIi4kXhICId0pEjR7j11lsZNGgQ48ePZ/LkySxbtqxNa8jMzCQ5Ofmkx77//ntSUlJISUkhPj6egQMHkpKSwmWXXdbk13zrrbfqv3/ttddYsGBBq9YNCgcR6YCstcyZM4dLLrmEPXv2kJ6ezjvvvEN2drbXvlVVVW1a26hRo9i0aRObNm1i9uzZPP3002zatInPP/+8STWdGg6+onAQkQ5n5cqVhIWFcd9999U/1r9/fx588EHAfdqePXs2M2bM4NJLLyU/P585c+YwevRoJk2axObNmwFYuHAhixYtqn+N5ORkMjMzyczMZMSIEdx7772MHDmSK664grKyMgDS09MZM2YMY8aMYfHixU2uefr06fzsZz8jNTWV5557jnnz5vHuu+/W/7xr164APPHEE3zzzTekpKTwm9/8BoCDBw8yc+ZMhgwZwuOPP97Cv7WTBdyprCLS8Vz79rU+ed0Pb/nwtI9nZGQwbty4Mz53w4YNbN68mfj4eB588EHGjh3L+++/z8qVK7njjjvYtGnTGZ+/a9cu3n77bf7whz9w00038d5773Hbbbdx11138fzzz3PJJZfw2GOPNas9Ho+HunvEzZs377T7/PrXv2bRokWsWLECcEG3adMmNm7cSHh4OMOGDePBBx+kX79+zXrvU6nnICId3gMPPMCYMWO48MIL6x+7/PLLiY+PB2D16tXcfvvtAMyYMYO8vDyOHz9+xtesmysAGD9+PJmZmRQUFFBQUMAll1wCUP+aTfWTn/ykWfvXufTSS4mNjSUiIoILLriArKysFr1OQ+o5iIjPNfYJ31dGjhzJe++9V//94sWLyc3NJTX1xM1Ho6Kizvo6ISEh1NTU1H/f8DqB8PDw+j8HBwfXDyudi4Y1NXzvmpoaPB5Po887tZbWmEdRz0FEOpwZM2ZQXl7Oiy++WP9YaWlpo/tffPHFLF26FIBVq1bRrVs3YmJiGDBgABs2bADcMNTevXvP+L5xcXHExcWxevVqgPrXbIkBAwaQnu4Wv1y+fDmVlW7FgujoaIqKilr8uk2lcBCRDscYw/vvv89XX33FwIEDmTBhAnfeeSdPPvnkafdfuHAh6enpjB49mieeeILXX38dgBtvvJH8/HxGjhzJ888/z9ChQ8/63kuWLOGBBx4gJSWFc1kv59577+Wrr75izJgxrFu3rr5XMXr0aIKDgxkzZkz9hLQvBPRiP82hxX5EAsu2bdsYMWKEv8uQBk73O2mXi/2IiIh/KBxERMSLwkFEfKajDFt3BM39XSgcRMQnIiIiyMvLU0AEgLr1HCIiIpr8HF3nICI+0bdvX7Kzs8nJyfF3KcKJleCaSuEgIj4RGhra5FXHJPBoWElERLwoHERExIvCQUREvCgcRETEi8JBRES8KBxERMSLwkFERLwoHERExEvAhYMxJtwY86oxJssYU2SM2WSMmeXvukREOpOACwfcVdv7gWlALPCfwJ+NMQP8WJOISKcScLfPsNaWAAsbPLTCGLMXGA9k+qMmEZHOJhB7DicxxvQAhgIZ/q5FRKSzCOhwMMaEAkuB162120/z8/nGmDRjTJru/Cgi0noCNhyMMUHAnwAPsOB0+1hrX7bWplprUxMTE9u0PhGRjizg5hwAjDEGeBXoAVxlra30c0kiIp1KQIYD8CIwArjMWlvm72JERDqbgBtWMsb0B/4NSAEOG2OKa7e5fi5NRKTTCLieg7U2CzD+rkNEpDMLuJ6DiIj4n8JBRES8KBxERMSLwkFERLwoHERExIvCQUREvCgcRETEi8JBRES8KBxERMSLwkFERLwoHERExIvCQUREvCgcRETEi8JBRES8KBxERMSLwkFERLwoHERExIvCQUREvCgcRETEi8JBRES8KBxERMSLwkFERLwoHERExIvCQUREvCgcRETEi8JBRES8KBxERMRLQIaDMWaBMSbNGFNhjHnN3/WIiHQ2If4uoBEHgV8CVwKRfq6lU7PWUlFdQWllKZXVldTYGqptNcEmmPCQcMKCw4gMiSQ4KNjfpYpIKwrIcLDW/hXAGJMK9PVzOR2atZbDxYfJPp7NgaIDHDh+gNzSXPLK8sgvy6fIU0SNrTnjaxgMXcO6EhMeQ0JkAj269qBHVA96R/emf1x/ekf3JiQoIP+piUgj2vX/WGPMfGA+QFJSkp+raR+KKorYmrOVjJwMfsj/gT3H9lBSWXLG54QFhxEVGkVoUCjBQcEEmSCqaqrwVHvwVHsorSylyFNEkaeIA0UH4OjJzw8JCiEpJolh3YYxvNtwRnQbQc+uPTHG+LClInIujLXW3zU0yhjzS6CvtXbe2fZNTU21aWlpvi+qnamqqWJbzjbSDqax4dAGsgqzsJz8O4+LiCMpJom+MX3pG9OX7lHdSeiSQEJkAtHh0Wf91F9jayiqKKKwopC80jwOFx/mSMkRso9nk1mQyZGSI17P6d6lO2N6jmFsz7GM6zWOqLCoVm23iDSNMSbdWpt66uPtuucgp+ep9pB+MJ3V+1aTdiiN0srS+p+FBoUyLGEYI7uPZGjCUAbHDyY+Mv6c3i/IBBEbEUtsRCxJsd49uLLKMnYf28323O3syN1BRk4GR0uP8tmez/hsz2cEm2BGdR/F5H6TmdJvCnERcedUj4icO/UcOghrLd8f/Z7P93zOuux1lFeV1/+sX0w/Ununkto7leHdhhMWHObHSl2tewv2svHQRtIOppGRk1HfmwkyQaT0SGHagGlM6TeFiJAIv9Yq0tE11nMIyHAwxoTgejX/hZuQvheostZWNfaczhoO+WX5/P2Hv/PZns/IKc2pf3xI/BCm9pvK1KSp9Oza048Vnl1RRRHrD65n9b7VbDi0gWpbDUBkSCTT+k/jysFXMjh+sJ+rFOmY2ls4LMQFQ0P/ba1d2NhzOlM4WGvJyMlgxc4VrMteV382Ufcu3ZkxcAYzBs6gV3QvP1fZMkUVRazZv4Yv9nzB9rzt9Y8PSxjG7GGzmdJvis58EmlF7SocWqIzhEN1TTVr969l2fZl7MrfBbhhmMl9JzNr8CxG9xjdoc4A2le4j093f8oXe7+g2FMMQHxkPNcOvZarhlxFl9Aufq5QpP1TOLRjVTVVrNy7kr9k/IXDJYcBiAmP4arBVzFz8EwSuiT4uULfqqiq4MvML1m+Yzn7j+8HICo0iquHXM3sYbOJjYj1c4Ui7ZfCoR2qqqni8z2f8+eMP9fPJ/Tq2ovrh1/PjIEzCA8J93OFbctay4ZDG3h367tsydkCQERIBNcMuYYbRtxAdHi0nysUaX8UDu1Ija3h66yvWbp5aX1PoV9MP25OvpmLki4iyATkLbHa1Pbc7fw548+sP7gecCExZ9gcbhhxA5GhuuOKSFMpHNqJDYc2sGTjEjILMwHoE92H20bfxtR+UzvUfEJr2Zm3k7e+f4v0Q+kAxIbHcnPyzcwcPFMT1yJNoHAIcFkFWfxx4x/ZcHgDAIldEpk7ai7TB0zXTe2aYFvONpZsWsK23G2AG367e+zdTOgzQaEqcgYKhwBVVFHE0u+X8vGuj7FYokKjuGnkTVwz9Bq/X6zW3lhr+eeBf7Jk0xJ3jycgpUcK94y7h/5x/f1cnUhgUjgEmBpbw6e7P+WN796gyFNEkAli1uBZTIq6lddfjuGttyA3F7p1g1tvhQUL4Pzz/V11+1BVU8Unuz7hrS1vUewpJsgEce3Qa7l11K06/VXkFI2Fg2Y2/WDPsT089uljLF6/mCJPEaO6j+K5mc/RP+c+Lr8khshIWLsWKirc18hImDQJPvnE35W3DyFBIVw77FpevuZlrhp8FdZaPtjxAfd/dD9fZ31NR/lAJOJL6jm0obLKMpZ+v5TlO5ZjscRHxnPvuHuZ2m8qe/YYJk2C5cth8mTv565bB7Nnwz/+oR5Ec+3O382LaS+yI28HAON6juP+C+8P+NuKiLQFDSv5WfrBdBavX0xOaQ4Gw+xhs5k7am79aZcPPwyREZZf/UcxFBTA8eNQUuK2sjLweHjrNQ81nipuu632RY2B0FAICYHwcIiKclvXrhAXB+edB2GatwA3H/HZns9YsmkJxZ5iwoLDuCX5Fq4ffr0m/KVTUzj4yfGK47zyjxfYuvFTEvPKSC6LYWbUGBIKKuDgQTh6FHJyyN2eQ0LQMUx1desWEBUF3bu7rUcP6Nv3xDZwIAwaBD17QlDnGGEsKC/glQ2v8FXWVwAMPm8wD016iAFxA/xbmIifKBx8rbgYdu6EXbvc9sMPFG7bRNXuXZyXX0rQqX/NMTHQq5c7YCcm8tJ7icz/RQJB8XHuU39MzImeQEQEhIdTGRTOmHHBbN1We2pmTQ1UVkJVlZugqOtpFBVBYSEcOwZ5eZCTA0eOwOHDkJ3tHm8oIgKGDIHhw2HYMBg5EkaNgqFDXc+kA9pwaAPP//N5ckpzCDbB/PiCH/OT5J/o2gjpdBQOrSU/HzIyYOvWE9v27e6g20BRYgxZCSEcSewCAwYwfvINxA0dDf36QZ8+7qDfQI8ebvL5TPMJu3fD1KnuGH9OSkpg/37IzIQ9e9wL79zp2rFnjwsdcMGQnAzjxrktNRXGjHFDWB1AWWUZr3/3Oh/t+giAgXED+fnkn6sXIZ2KwqG5PB7Ytg02bz6xbdnihoLqREXBiBFuq/3U/X1MOc/lruAIxUSERHBXyl3MGjzrrBdiPfywOyvpV79qfJ9f/MJ1EJ59tpXaeDoVFS4kvv/ebRs3ui031/08LAxSUtys+dSpbuvd24cF+d6Wo1t47h/PcbjkMCFBIdySfAs3jrhRcxHSKSgcziQvD777DjZtOrFt2+aGa8B9Ur7gAjfUMmqUG3YZOdKN29eO1ZdWlvLKhlf4bM9nACQnJvPQpIeafEbM7t0E7tlK1rqexvr18O23blu/3k2Ug5u3mDbNbT/6ESR5LxUa6MqrylmycQkf//Ax4NaP+Pnkn9M7un0Hn8jZKBwaM2cOfPDBie9793afjMeMgdGj3dchQ9wZQY3YlrONZ9Y9w5GSI4QGhXLHmDu4bth1zb5twyefwB13wD33uC0pCfbtg1decdsbb8CsWc1vok94PC5E16yBr7+Gr746MZcxeDBcdpnbLr3UzaG0E5sOb+K5b58jtzSX8OBw7h57NzMHz9QtOKTDUjg05rXX3JBJXSAkJjb5qVU1Vbyz5R3+nPFnLJZBcYN4dMqj9Ivt1/w6au3eDYsX43WF9AMPBPj1DTU1bhjqyy/h889dWBQXu57VxIkwc6ZLtvHjA/7MqBJPCS+lvcSqrFUAXNj7Qh6a+JDWjZAOSeHQyg4VHeLptU+zK38XBsONI25k7ui5OtulTmWlG3769FP4+9/dMJS1LnxnzoRrroErr4TYwD3gfpP1DYvXL6aksoS4iDh+NvFnjO893t9libQqhUMrsdaycu9KXkp/ifKqchK7JPLzyT8nuXuyz9+7XcvNdSHxySduy893Q3UXX+wmU667zl13EWByS3N5Zu0z9YsLzR46m3kp8wgN7pin+Erno3BoBSWeEl5Y/wJf7/sagIv6XcSCCQuICos6yzPlJNXVbmZ9xQo3A791q3s8OdnNAc2Z406dDZBx/hpbw1+3/ZU3N79Jta1mYNxAHp/6OH1j+vq7NJFzpnA4RzvzdvLUmqc4UnKEiJAI7ht/HzMGztBEZWv44Qf48EN3YsA337j5i759XUjccIPrXZzhhIC2sjNvJ4vWLuJQ8SHCgsO4b/x9XDboMv0bkHZN4dBC1lqWbV/GG9+9QbWtZlDcIB6f+jh9Yvq0+nsJbvhpxQpYtszNV5SXQ0KCG3q64QZ3BlREhN/KK6ss48W0F/ky80sALkm6hAcmPKBbgUu7pXBogcLyQp5d92z96mzXDbuOO8fcqfHmtlJSAn/7mwuKFSvcLUG6doWrr3ZBMWsWREf7pbSVe1fyYtqLlFeV0zOqJ49PfZwhCUP8UovIuVA4NNPmI5t5Zt0z5JflEx0WzcOTHubCPhe22utLM3k87jTZ995zw09Hj7qLE6+80gXFtddCfHyblnTg+AGeWvMUewr2EBIUwrwx85g9bLaGmaRdUTg0UY2t4Z0t7/DOlnewWEYmjuTRKY/SrUu3VqhSWkV1tbsR1XvvuV7Fvn0QHAzTp7ugmDOnzW7pUVldyZJNS/hw54cATOwzkYcmPkR0uH96NCLNpXBogvyyfJ5Z+wybj27GYLhp5E3cknyL7rETyKyFDRtOBMX27e7xiRNPnPk0fLjPy1i7fy2//fa3lFSW0K1LNx6f8jgjEkf4/H1FzlW7CgdjTDzwKnAFkAv8wlr71pmec67hsOnwJhatXURhRSFxEXE8MvkRUnqmtPj1xE+2bXMh8f777sI7cLcenzPHXUsxcaLrZfjA0ZKjPLXmKXbk7SDIBHHnmDu5fvj1GmaSgNbewuFt3PrWdwMpwEfAFGttRmPPaWk4VNdU89b3b/GXrX/BYhnTYwyPTH6E8yLPa2n5Eij273fXUXzwgZuvqKpyV2hffbWbo7jiCjfB3Yqqaqp447s3WLZ9GQCpvVJ5ePLDxITHtOr7SAdUXe1WgczLcxeJ1m3HjrmtoODkrbAQCgupKSwg6Nt/uhtgtsA5h4MxJgnYb095gnEfi/pZa/e1qDLv94kCjgHJ1tqdtY/9CThgrX2isee1JBwKywv59epfsyVnCwbDraNu5aaRNxFkAvveP9IChYXuyuzly93XggJ3+/Fp01xYXHWVu8FiK1l/YD2/+cdvKPIU0a1LNx6b8hgXJF7Qaq8v7UB1tTs1+8gRtx09Wr/yI0ePup/l5LivubkuAM50PK5b/jcuDmJjsbEx7LOFbKs8yMBn/siw0T9qUZmtEQ7VQC9r7dFTHk8AjlprW6WvbowZC6yx1nZp8NijwDRr7bWn7DsfmA+QlJQ0Pisrq1nvVVZZxsN/f5jSylIem/IYo3qMOvcGSOCrrHR3k12xAj766MQ8xfnnu/s+zZzpJrfPsVeRW5rLU2ueYlvuNoJMELePvp0bR9yoYab2zuOBQ4fc2i5126FDJ7bDh92Wk3Ni4ayGQkJcD7ZbtxNfu3Vz1/PUbfHxbjvvPPc1Lu6kVRmPlR1j0dpF9fOjd4+9m+uGX9ei5rRGONQAPay1Oac83h/Yaq1tlXtIGGMuBv5ire3Z4LF7gbnW2umNPa+lw0rZx7PpGtaVuIj2c1tpaWV79rjexN/+BitXQmmp+484ZQpcfrm77Xhqaouu0q6qqeLNzW/y3rb3ABjfazwPT3pYd3gNVGVlblXH/fvd19NtOTnezwsJcWux1229ermvPXqc2OrWco+LO6dbw3x3+DsWrVtEQXkBseGxPDrl0XOaH21xOBhjflv7xweAJUBpgx8HAxMAj7V2aourO/n9TtdzeASYfmrPoSG/ryEtHUNFBaxeDZ995q7Q3rjRPR4dDZdc4noU06bB2LHNCou0g2k8u+5ZijxFxEfG89iUx3SzxrZmrRvO2bcPsrLctm+f2/bvd19Pd+CPj3e3c+nb1y3xW7f17n3ia0KCz29FX11TzTtb3uF/M/4Xi2V099E8MuUR4iPP7fqecwmHL2v/OA1YB3ga/NgDZAKLrLW7zqnCE+9XN+cwsu41jTFvAAdbe85B5KxycmDVKvjiC9er2FX7zzw62i3dd9FFbqnUCRPOerV2w2Emg2HuqLn8eOSP6+e4du+G55/3XstjwYIAX8sjUFRXw4EDJw78mZkn/lwXBOXlJz+na1fo39+t7d6vn1thq+7Pdeu9d/H/rVHyy/J5es3T9fOjNyffzM3JN7fK/GhrDCstAR6y1h4/52rO/l7vABa4B3e20sf46GwlkWY5dOjEyndr1rgFjqx1wwTJye5U2QkT3KJGyclu0ruB6ppqln6/lHe3vnvSp79vV8Vzxx1w771w993ueJWVBa++Cn/4Q4CtAugvFRUnPuGfeuDPzHRDPnVL+9bp3t39ZTa2xcYGzN1/G5N+MJ1n//EsxyuOc17EeTw65VFG9xjdaq/f3k5ljQf+CFwO5AFP+Po6B5EWKShwtx//9tsTX+uWSw0LcwHRcNnZkSMhMZGNhzbyzLpnKKwoJMgTy9dP/pyPXx0XeOuHtxVr3WmbdQf/hltdABw+fPLZPMa4IZ0BA9xWd8Cv+3NSEkRG+qlB5+7U06JTeqTwyJRHWn1+tF2FQ0soHCQgWOsmuNPTIS3NrbO9adPJY9mJiTBiBOWDB7IyLJu/5h3nQEIXbr3zVm6Z8K+nXU3wF79wH5yffbYN29JaKivdqZyHDrlhn4MH3deGk7z797vJ4IbCwtwBvm5reODv39/NAZzSM+soDhcf5uk1T7MzfydBJojbRt3Gv1zwLz45003hIOIv1roD45YtkJHhvm7fDjt2uAueGiiMjyLy/OGEDTzfHRB794bevTlQ04ufPNid1Ttrz3bx5zrclZWux5Sf7+rPy3PhV3f+/tGjLgzqTunMzfU+fz8k5MSEbp8+J4/114VBYmLArzfuC19nfc3i9YsprSxtk1uxKBxEAlFeHpMTf+B/n/+atNVv0/VQLj3zKhhcFEqXw3neE6jgbv/R8Pz3mBg3GR4d7SZPu3RxwylhYW4LDXUH2brNWnf+fU2NG6P3eNwBv6LCvV9ZmTudt6QEiouhqAiOH3cXEhYUuMcbExnpxvkbnsbZu7c7tbPh2T3du3fKA/+ZlFeV8/u03/P53s8BmNJ3Cg9OfJCuYa17Ff+pGgsH/y+vJdKZJSSwJzGByisncsU99/O7b3/H2uy1APyo/3TuG3wzXXIKOJh+iP95MIeX/qf2k3ndLRXy893B+/BhdwCvO7CXlp75atvTCQ93CymFh0NU1IktJsYd1GNiTlyhWxdOdRdtdevmDvhRWjK3JXbl7WLR2kUcLD5IWHAY94y9h5mDZ/r1gkn1HET87OGH3QfuX/3KrTz46e5PeXnDy3iqPfSM6smjUx7ltWeGNX/Oobra9Qg8nhO9hepq94ndGPc1JMT1LkJCAv6snY6oxtawbNsy/rT5T1TbagbEDuCxqY+RFJvUZjVoWEkkQO3e7S6ZWL6c+rOV9hfuZ9HaRewp2EPBsSC+e/MW0l7/MUOH6PbxHcXRkqP8Zt1v2JKzBYBrh17LvJR5hAW37SR7Y+GgQT8RPzv/fHcdw+zZ7qyk3buhZ5d+3D94ETUZN/DtP2sY8ZOlvLz33zlYdNDf5co5stayKnMVP/3kp2zJ2UJcRBwLpy1k/vj5bR4MZ6JwEAkAs2a56xgqKtwF15GRMP3iUIaV3sWyBf+HCwZ0Y0feDn76yU/5ZNcndJQef2dzvOI4T655kmfWPUNJZQkT+0zk+VnPM773eH+X5kXDSiLtQImnhJfSXmJV1ioAxvYcy08n/lTL17Yj6w+s57f//C0F5QVEhERw77h7uXzQ5X6/S6/mHEQ6gNX7VvPC+hco8hTRJbQL88fNZ8bAGX4/wEjjij3FvLLhFb7Y+wUAyYnJ/GzSz+jRtYefK3MUDiIdxLGyYyxev5hvD3wLuNuAP3DhAyRGJfq5MjnV+gPreX798+SX5RMWHMbto2/numHXBVSYKxxEOpC6Sc3fp/+eksoSIkIiuCvlLmYNnhVQB57OqrC8kJfTX+brfV8DMKLbCB6a+BB9Yvr4uTJvCgeRDii/LJ+X0l5iXfY6AC7odgELJiygX2w/P1fWOVlr+TLzS17Z8ApFniLCg8O5ffTtXDvs2oBdfljhINJBWWtZs38Nv0//PQXlBYQEhXDjiBu5aeRNAXVqZEd34PgBXlj/ApuPbgbcXVQXTFgQMHMLjVE4iHRwxZ5iXt/0On/b/TcAekb1ZP74+VzY50I/V9axeao9/CXjL7y77V2qaqqIDovm7rF3t5sTBRQOIp3E1pytvLD+BbIKswCY2Gci94y7h55de57lmdIc1lrWZa/j1Q2vcrT0KABXDLqCeSnziA4/86qAgUThINKJVNVU8dHOj1j6/VLKqsoICQphzrA53DTyJiJD2+8COIEisyCTVza8wndHvgNgYNxA7k+936e31vYVhYNIJ5Rfls/rm15nZeZKAOIi4pg7ai6XD7qc4CDdp6m58svyeXPzm3y+53MsluiwaG4bfRszB88M2Anns1E4iHRiO/N28of0P7A9bzsAfaP7cmfKnUzsM7FdjIv7W4mnhGXbl/HBjg8oryon2ARz1ZCruCX5lnY1hHQ6CgeRTq7urKY3vnuDQ8WHABgaP5S5o+cytudYhcRplFeVs2LnCt7b9h7FnsSKj0gAAAyuSURBVGIAJvWZxLyUeQF5zUJLKBxEBHDzEX//4e+8veVtCisKAXeR1s3JNyskapVVlvHRro94f/v79X9HyYnJ3DHmjnY5r3AmCgcROUl5VTkf7/qYd7e+S5GnCIDzzzuff7ngX5jSb0q7HUM/F4XlhXy06yNW7FxR/3cyNH4ot42+jZSeKR0yOBUOInJadZ+SP9jxAQXlBYC7RuKaoddw2aDLiArr+Et/7i/cz/Idy1mZuRJPtQdwvalbkm/psKFQR+EgImfkqfbwxZ4v+Ou2v3K45DAAESERzBgwgysHX8mg8wb5ucLWVVVTxbfZ3/Lxro/rr2oGmNB7AtePuJ6RiSM7dCjUUTiISJPU2BrWH1jP8h3LTzpoDokfwuWDLueipIva9Rk6mQWZfL7nc1ZlrqqfT4gIiWB6/+lcN/w6+sb09XOFbUvhICLNllmQyae7P2Xl3pWUVJYAEGyCGddrHJf0v4QLe1/YLoadso9ns3rfatbsW0NmYWb94/1i+nHVkKv40YAftYt2+ILCQURazFPtYc2+NazKXMXGwxuxuONGsAlmVPdRTOgzgbG9xtInuk9ADMVU1VSxNWcraQfTSDuYxv7j++t/FhUaxbT+07hs0GUMjh8cEPX6U7sIB2PMAmAeMAp421o7r6nPVTiItI2C8gK+yfqGddnryMjJoMbW1P8sITKB0T1GM6LbCIZ3G05SbFKbXIld4ilh97HdbM3ZypajW9ieu52K6or6n0eFRjGp7ySm9ptKSs8UQoNDfV5Te9FewuEGoAa4EohUOIgEtqKKItYfXE/6wXS+O/Jd/Rh+nfDgcJJikxgQN4D+sf3p2bUnPbr2oEdUDyJCIpr1qb3G1nCs7Bi5pbkcLDrIgaIDZB/PZs+xPfUX9TXUP7Y/qb1TSe2dyvBuwwkJCjnn9nZE7SIc6hhjfgn0VTiItB/WWjILMsnIyWB77na2527nSMmRRvcPCw4jJjyGmLAYIkIiCAsOIyw4DIulxtZQXVNNWVUZpZWllFSWUFBecFIvpaGQoBAGxg1kWMIwkrsnk9w9mdiIWF81tUNpLBzadZQaY+YD8wGSkpL8XI1I52aMYeB5Axl43kCuGXoN4HoWWYVZZBZksr9wP0dKjnCk+AhHS4/iqfaQW5pLbmluk98jNjyWhMgEekX3ok90H/rG9KV/XH+SYpPUM2hl7fpv01r7MvAyuJ6Dn8sRkVNEh0fXf5I/VXlVOYXlhRyvOE5FdQWeag+eag8GQ5AJIjgomMiQSKLCougS2oW4iDgFQBtqs79pY8wqYFojP15jrb2orWoREf+LCIkgomtEwC+j2Vm1WThYa6e31XuJiMi5Cag+mjEmBFdTMBBsjIkAqqy1Vf6tTESkcwm02y7+J1AGPAHcVvvn//RrRSIinVBA9RystQuBhX4uQ0Sk0wu0noOIiAQAhYOIiHhROIiIiBeFg4iIeFE4iIiIF4WDiIh4UTiIiIgXhYOIiHhROIiIiBeFg4iIeFE4iIiIF4WDiIh4UTiIiIgXhYOIiHhROIiIiBeFg4iIeFE4iIiIF4WDiIh4UTiIiIgXhYOIiHhROIiIiBeFg4iIeFE4iIiIF4WDiIh4UTiIiIgXhYOIiHgJmHAwxoQbY141xmQZY4qMMZuMMbP8XZeISGcUMOEAhAD7gWlALPCfwJ+NMQP8WJOISKcU4u8C6lhrS4CFDR5aYYzZC4wHMv1Rk4hIZxVIPYeTGGN6AEOBDH/XIiLS2QRkOBhjQoGlwOvW2u1n2G++MSbNGJOWk5PTdgWKiHRwbRYOxphVxhjbyLa6wX5BwJ8AD7DgTK9prX3ZWptqrU1NTEz0cQtERDqPNptzsNZOP9s+xhgDvAr0AK6y1lb6ui4REfEWMBPStV4ERgCXWWvL/F2MiEhnFTBzDsaY/sC/ASnAYWNMce0218+liYh0OgHTc7DWZgHG33WIiEgA9RxERCRwKBxERMSLwkFERLwoHERExIvCQUREvCgcRETEi8JBRES8KBxERMSLwkFERLwoHERExIvCQUREvCgcRETEi8JBRES8KBxERMSLwkFERLwoHERExIvCQUREvCgcRETEi8JBRES8KBxERMSLwkFERLwoHERExIvCQUREvCgcRETEi8JBRES8KBxERMSLwkFERLwEVDgYY940xhwyxhw3xuw0xtzj75pERDqjgAoH4P8CA6y1McBs4JfGmPF+rklEpNMJqHCw1mZYayvqvq3dzvdjSSIinVKIvws4lTHmBWAeEAlsBD4+w77zgfm13xYbY3a08G27AbktfG57pTZ3Dmpzx3eu7e1/ugeNtfYcXtM3jDHBwGRgOvCktbbSx++XZq1N9eV7BBq1uXNQmzs+X7W3zYaVjDGrjDG2kW11w32ttdXW2tVAX+D+tqpRREScNhtWstZOb8HTQtCcg4hImwuYCWljTHdjzM3GmK7GmGBjzJXALcAXbfD2L7fBewQatblzUJs7Pp+0N2DmHIwxicC7wBhcaGUBv7XW/sGvhYmIdEIBEw4iIhI4AmZYSUREAofCQUREvHSKcDDGxBtjlhljSowxWcaYWxvZzxhjnjTG5NVuTxpjTFvX2xqa0ebHjDFbjDFFxpi9xpjH2rrW1tLUNjfYP8wYs80Yk91WNba25rTZGDPOGPO1MabYGHPEGPNQW9baWprxbzvcGPNSbVvzjTEfGmP6tHW958oYs8AYk2aMqTDGvHaWfR82xhyuvT/dH40x4S19304RDsBiwAP0AOYCLxpjRp5mv/nAHNyk+GjgWuDf2qrIVtbUNhvgDuA8YCawwBhzc5tV2bqa2uY6jwE5bVGYDzWpzcaYbsDfgN8DCcBg4NM2rLM1NfX3/BDuYtrRQG/gGPC7tiqyFR0Efgn88Uw71Z7h+QRwKe6q50HAf7f4Xa21HXoDonD/kIY2eOxPwK9Ps+9aYH6D7+8G/uHvNviyzad57m+B3/m7Db5uMzAQ2AbMArL9Xb+v2wz8CviTv2tu4za/CDzV4PurgR3+bsM5tP2XwGtn+PlbwK8afH8pcLil79cZeg5DgSpr7c4Gj30HnO6Txsjan51tv0DXnDbXqx1CuxjI8GFtvtLcNv8O+A+gzNeF+VBz2jwJyDfGrDXGHK0dYklqkypbV3Pa/Cow1RjT2xjTBdfL+KQNavSX0x2/ehhjElryYp0hHLoCx095rBCIbmTfwlP269oO5x2a0+aGFuL+TSzxQU2+1uQ2G2OuB4KttcvaojAfas7vuS9wJ26oJQnYC7zt0+p8ozlt3gXsBw7UPmcE8D8+rc6/Tnf8grP/vz+tzhAOxUDMKY/FAEVN2DcGKLa1fbR2pDltBtykF27u4Wp74rbp7UmT2myMiQKeAn7aRnX5UnN+z2XAMmvtemttOW4seooxJtbHNba25rR5MRCOm2OJAv5Kx+45nO74BWf4f38mnSEcdgIhxpghDR4bw+mHTjJqf3a2/QJdc9qMMeZfqZ3Ista21zN3mtrmIcAA4BtjzGHcAaNX7RkeA9qgztbUnN/zZtz6KHXa2weeOs1pcwpujD6/9gPP74AJtZPzHdHpjl9HrLV5LXo1f0+ytNFEzju4LnQUMBXX3Rp5mv3uw01S9sGd3ZAB3Ofv+n3c5rnAYWCEv2tuizbjbubYs8F2A+5skJ64oSa/t8NHv+cZuLN1UoBQ4DfAN/6u38dtXgK8B8TWtvk/gAP+rr8F7Q0BInArZf6p9s8hp9lvZu3/5QuAOGAlTTgJpdH39XfD2+gvNx54HygB9gG31j5+MW7YqG4/gxtyyK/dnqL2FiPtbWtGm/cClbguad32kr/r92WbT3nOdNrp2UrNbTPu9vcHakPiQ6Cfv+v3ZZtxw0lLgaNAAbAamODv+lvQ3oWcWBmzbluImzsqBpIa7Ptz4AhujmUJEN7S99W9lURExEtnmHMQEZFmUjiIiIgXhYOIiHhROIiIiBeFg4iIeFE4iIiIF4WDiIh4UTiIiIgXhYOIiHhROIj4gDEm0RhzyBjzXw0eG22MKTfG/NiftYk0hW6fIeIjtcs2fghMAzYBacA/rbV3+bUwkSZQOIj4kDHm/wGzga9wN4ZLsdYW+7cqkbNTOIj4kDEmHLdc4xBgirX2Wz+XJNIkmnMQ8a0BQD/cbZYH+bcUkaZTz0HER4wxocA/cKuXfQv8FzDGWrvPr4WJNIHCQcRHjDG/Bm4FRuNWK/sEt4rXDGttjT9rEzkbDSuJ+IAxZhrwCHCHtbbAuk9h83BLOP67P2sTaQr1HERExIt6DiIi4kXhICIiXhQOIiLiReEgIiJeFA4iIuJF4SAiIl4UDiIi4kXhICIiXv4/vxkp/zjf/JYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Podemos plotear la predicción semi-bayesiana\n",
    "y = np.dot(phi_, m_n)\n",
    "\n",
    "plot_data_sine(x_[ii], t_[ii])\n",
    "plt.plot(xx, y, '-r')\n",
    "plt.ylim(-3, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Encontremos ahora los parámetros de la posterior predictiva\n",
    "\n",
    "$$\n",
    "p(t^\\prime | \\b{t}, \\beta, \\gamma) = \\norm(t^\\prime | \\mathbf{m}_N\\T \\boldsymbol{\\phi}(x), \\sigma_N^2)\\;\\;,\n",
    "$$\n",
    "\n",
    "donde la variaza es\n",
    "$$\n",
    "\\sigma_N^2 = \\frac{1}{\\beta} + \\boldsymbol{\\phi}(x)\\T \\mathbf{S}_N \\boldsymbol{\\phi}(x)\\;\\;,\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 424,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Busco la media en los 100 puntos que voy a usar para plotear\n",
    "# El vactor phi es la transpuesta de la matriz de diseño\n",
    "mean = np.dot(m_n.T, phi_.T)\n",
    "sigma2 = np.dot(np.dot(phi_, np.linalg.inv(sn_inv)), phi_.T) + 1/beta\n",
    "\n",
    "# Me quedo solo con la diagonal para los errores puntuales\n",
    "err = np.sqrt(np.diag(sigma2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 425,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PolyCollection at 0x12dba99d0>"
      ]
     },
     "execution_count": 425,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAENCAYAAADKcIhSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deZxcZ3no+d9zTu1rr+rW1mrJWm3Zkm3Z7GDikEBuHLhAJmyJDQGTzJCQIZCQOyHJMAwx+SSXJA6XZYCw3ACTxbnEhNzP5M4AMeAY5FUWshZLaqkltdR7LV3bOeedP051uyV1t7rVtXX38/18ClxVb1W9R1V9nvM+7ybGGJRSSqnlsJpdAaWUUiufBhOllFLLpsFEKaXUsmkwUUoptWwaTJRSSi1boNkVaIauri7T39/f7GoopdSK8vjjj48YY7rnem5NBpP+/n4OHjzY7GoopdSKIiID8z2naS6llFLLpsFEKaXUsmkwUUoptWxrss9EKdVaKpUKg4ODFIvFZldFAZFIhE2bNhEMBhf9Gg0mSqmmGxwcJJlM0t/fj4g0uzprmjGG0dFRBgcH2bp166Jfp2kupVTTFYtFOjs7NZC0ABGhs7Nzya1EDSZKqZaggaR1XM93ocFEKaXUsmkwUUop4OLFi7ztbW9j27Zt3H777bzkJS/hoYcewjMGx/Moux5Fx6VQcZmqOOTL1VvFYariUKi4lByXsuvheB6eMSx1v6jTp0+zd+/eyx47dOgQ+/fvZ//+/XR0dLB161b279/PT//0Ty/6Pb/2ta/N3P/Sl77E+973viXVazG0A14pteYZY3jDG97Avffey1f/69/gGsPJU6d5+OF/Ild2EGA6LDiOQzAw69Rppv/vhcAxXV4A2xICloUtgiVLTyHdfPPNPPXUUwDcd999/PzP/zxvfvObLyvjOA6BwNyn8+lg8ra3vW1Jn7tU2jJRSq15//o//l8CwSDveNe7yVccio7Lxs2b+bX/5X1YIvzNV7/CW970Ru557c/wC6/7WcbHx3nrL76Zlxy4jZ965cs5/OwhLBEe+Nj/wYOf/CQigiXCi267lZMnT3H0xPPceOMe3vnud3PjTTfxmte8hqmpKQAef/xx9u3bx759+/jUpz616Drfdddd/NZv/RYHDhzgL/7iL7jvvvv4+7//+5nnE4kEAB/+8Id55JFH2L9/P5/85CcBOH/+PK997WvZsWMHv/M7v1OTf0NtmSilWso9X7+nLu/78Fsfvuy+MQbXGMqux1OHDnHL/lsxxmDN03J4+qknefTgE3R0dPDB33o/+/bv5xt//w987zvf4f53vZMf/vjxOV8nIgjw/IkTfPEr/5W//C+f4d63v5Wv/e3f8cvveAfvfOc7+au/+ite+cpX8qEPfWhJx1Qul2fWGbzvvvvmLPPAAw/wp3/6p3zrW98C/DTXU089xZNPPkk4HGbXrl38xm/8Bps3b17SZ19JWyZKqTXFGEPF9chXXKYqLq5nEPyU1HQK6gO/+Ru85MBtvOqlL5553avvvpuOjg4AHv3hD3jL294OwKte/WrGxsbIZDILfm5//1b27d+PJcKtt93OwOkBLo6OMT4+wYGXvAzH83jHO96xpGP5pV/6pSWVn3b33XeTTqeJRCLceOONDAzMu37jomnLRCnVUq5sQdSKMQbHM5Rcv3NcYKYVsufGm/jmf/vHmbL/+S8fZGRk5LJgEo/Fr/kZgUAAz/Nm7pdmzdUIhcMz/21bNkW34H++gON5OB4UHXemrovpW4nHX6jT7M/2PI9yuTzv68Kz62LbOI5zzc+6Fm2ZKKVWPdfzmKq4FBx3JpU1+2T9qle/mlKxyOc/+5mZxwrVPo25vPRlL+dvv/F1AB753vfo7OwklUrRt2ULTz/1JABPPfkEp0+fWrBebW1tpNNp/v2HP8AS4Rtf/xqeMUxVXBzPW9JosP7+fh5/3E+1/dM//ROVSgWAZDJJNptd9Ptcr5YLJiKSu+LmisiD85S9r/r87PJ3NbjKSqkW5RlDoeKSr7h4cwSRaSLC1//uH/j+I4+wd+cO7nrZS3jvu9/FR//Pj8/5vr/3kT/gySee4MW338of/P5/4rNf+CIAr/+Pb2R8bIw79u/js5/+L2zfsfOadfz0//V5fvv9v8lL77gdP3bITECZcvw03GK85z3v4Xvf+x779u3j0UcfnWm13HLLLdi2zb59+2Y64OtBljoOupFEJAEMAT9njPm3OZ6/D3i3MeblS3nfAwcOGN0cS6nWceTIEfbs2VPT96y4HkXXW7BTvdV51fNz0LIIB6yGHsdc34mIPG6MOTBX+VbvM3kTcAl4pNkVUUqtDMYYio5HxfMu6xdZiabrXvE8nLJHKGARsqyWXHqm5dJcV7gX+IpZuPl0q4iMiMgxEfmIiMwZIEXkfhE5KCIHh4eH61NbpVRTuZ4hX3GpeN68Ka2VyO+oF0qOPwrNmdXJ3ypaNpiIyBbgVcCXFyj2b8BeYB1+K+atwJwDtY0xnzPGHDDGHOju7q51dZVSy7SclLupzhfJV5wVndZayHQry1T7UwrVfqB6uJ7vomWDCfDLwPeNMfMOhzDGnDTGnDLGeMaYQ8BHgTfPV14p1ZoikQijo6PXdRKbTmsVHRdh9bRG5jM9u77i+a2Uiru0UV/XMr2fSSQSWdLrWrnP5FeAB5b4munlcJRSK8imTZsYHBxkqSloYwwVz8zMG1lrf/7T64FZIgSt2gXS6Z0Wl6Ilg4mIvBTYCPzdNcq9DnjCGHNRRHYDH7nWa5RSrScYDC5pVz+AfNnhqYsZyq5HyF79LZL5GGMouYaAJezpTNAdD1/7RXXQqmmue4GHjDGXzbQRkb7qXJK+6kN3A8+ISB74NvAQMPfAcKXUqjFZqnBwaJKK5xEOtObopkYRESIBCzA8M5zhJ8NZKk3ooG/peSb1ovNMlFq5RqZKHBrOYok//0K9wG+leIRtixu7krRHQzV9/4Xmmeg3oZRaMYZyRZ65lMEW0UAyB7+VYuN4hicvTvL8eK5uI76upN+GUmpFOJ8t8JORLEHLImCt3bTWYgRti5BtMTBZ4MfnJ8hXlr+Q47VoMFFKtbxz2QJHRnMEbQtbA8miWCKEbYspx+VH5yc4ly3UdAjxVZ9Xt3dWSqkaOJct8NxojpDtb32rFk+qASUgwnOjOZ655I9+qwcNJkqplnUhV+ToaI6wBpJlsS0hYluMFMqcHM/X5TM0mCilWtKlfJEj1T6S1bg8SqOJCLYITp1SXRpMlFItZ2SqxLMjOYKW9pGsFBpMlFItJVOq8OxwloCIBpIVRIOJUqplTFX8JVIEdPjvCqPBRCnVEkqOx5NDk3jGELT11LTS6DemlGo6xzM8fWmSsucR0kCyIum3ppRqKmMMPxnOkCs7hHSJlBVLvzmlVNMYYzgxnme4UCZsr+3Vf1c6DSZKqaY5nytyJlPQQLIKaDBRSjXFRLHCsdE8IVsnJa4GGkyUUg1XdFwOXcpgWegyKauEBhOlVEO5nuGZSxkc4+meJKuIfpNKqYYxxnB0NKcjt1Yh/TaVUg1zPlfkQr6oHe6rUEsGExH5rogURSRXvR2dp5yIyCdEZLR6+4ToL1SplpQpvdDhrn+mq09LBpOq9xljEtXbrnnK3A+8AdgH3ALcA7y3URVUSi1O2fV45lIGS7TDfbVq5WCyGPcCf2aMGTTGnAP+DLivuVVSSs1mjOHwSIayq2turWaBZldgAX8sIg8AR4H/zRjz3TnK3AQ8Pev+09XHriIi9+O3ZOjr66ttTVXDVFyPXMUhW3LIlB2KjkfZ9ah4L2xFaosQDthEAxaJYIBUOEAyFNATWZMMTBYYK1SI6L//qtaqweR3gZ8AZeAtwMMist8Y8/wV5RLA5Kz7k0BCRMSYy7cTM8Z8DvgcwIEDB+qz1ZiqOWMM+YrLaKHMhVyRqYqHiP+4JYIIWMhVqZNCxSFfhkumhIhgDESDFr3xCOtiIWJBW/P2DTBeKHNyIq8d7mtASwYTY8xjs+5+WUTeCvwc8OAVRXNAatb9FJC7MpColcfxPC7mS5zJFCg6HsYYApYQtmVRJ6UrZ1QbY6i4Hqcm8pyamCIWtNiSirIuHtENmOqk7HocHslii+gM9zWgJYPJHAww16/xMH7n+4+q9/dVH1MrVMnxOJOZ4ly2iFcNICFLEFleikRECIgQsPzAUnI9jozmOD4+RV8qyqZUhIDOe6gZYww/GclSdg2RgP67rgUtF0xEpA14EfA9wAF+CXgl8P45in8F+ICIfBs/4Pw2V7de1ApQdj0GMlMMZooYAyFbsJYZQOYjIgRFCFr+bOyTE3nOZArc0B5jfSKiV9E1MJgtMlZdCVitDS0XTIAg8DFgN+ACzwFvMMYcE5FXAP9ijElUy34W2AYcqt7/fPUxtUJ4xnAuW+T5iTyeZxq+6J9tCbZl43r+zOwzkwX2dCVpiwQbVofVJlt2ODGu80nWmpYLJsaYYeCOeZ57BL/Tffq+AX6nelMrzGSxwpHRHFMVl6AlhJqYDpkOKmXX44mhCTYlo2xrj+s+5EvkeIZnL/l7uGsLb21puWCiVr/p1NLZTAFLFt+p3ghB2yJgDIPZAsOFMnu7k6TD2kpZrONjOQqOSyRgN7sqqsE0mKiGypQqHB7OUnC8lt3HQkSIBGwqrscTFybZ3hFnUzLSMgGvVQ3nS5zPFXU+yRqlwUQ1hDGGwWyRE2N5RFgRI3yCtoVnDMfGcowXK9zYldARX/MoOS5HRnMELe0nWav0L0PVXcXzODSc5fhYjqAthFbQlaslQsS2GJkq8fiFSYqO2+wqtRxjDEdGc7iep31Ma9jK+atWK1K+7PDj8xOMFEqEWzStdS3Taa8px+HHFybIlpxmV6mlnMsWGS2UV9RFgqo9/fZV3YxMlfjxhQnKrkfEXvnLl4RtG88zPD40wchUqdnVaQn5isNxHQas0GCi6sAYw8DkVHXJ8ZWV1rqWYLV19cxwlkv5YrOr01SeMfxkOAvosvJKg4mqMc8Yjo7lZiatrcYcesDyZ9A/O5zlQm7tBpQzkwUyZYfQKvyO1dLpaC5VM45nODycYaRQJrLK0x7+4pAWPxnJ4hnDxmS02VVqqGzJ4dTElK4GrGZoMFE1UXY9nr44SbbsrPpAMs22hLBYHB3NYYmwPhFpdpUawvX8za5EdJa7eoEGE7VsRcflyYuTFCvemrtStUQIWhZHRrJYAj3x1R9QTk3mmaroLHd1OQ0malmmKi5PDk1Q9jzCa/TkMp3yOjycxULojoebXaW6mSxVODNZWFWDKlRt6C9CXbdc2eHxCxOUPUPYXpuBZJptCQHL4tmRLBPFcrOrUxeuZzg8rJtdqblpMFHXJVt2eGJoEtcY3bOiKmAJlsDTlzLky6tvYuPz43mKjktQv281B/1VqCXLlCo8MTSBMUbTHVcIWhbG4PchraKlVyaKFQazmt5S89NfhlqSyVKFJ4YmwaBXqPMI2RYVz+Ppixkcz2t2dZbN8fwteDW9pRaiZwO1aJPFCk8OTSKigeRawrZNvuJweDiLv4fbynVS01tqEfTXoRZloljmyYvVQKLLsC9K2LYYKZQ5MZ5fsQFF01tqsfQXskQnxnIcHs6QKVVW7AliqSaKZZ66mNFAskQiQti2OJMpcH4FLrui6S21FDrPZInyFZdLUyUu5kskQgH601G6YuFV+8c2WapoIFmG6YUuj43mSQQDpCMrZwvg6fSWTk5Ui9FyZwcRCYvIF0RkQESyIvKUiLxunrL3iYgrIrlZt7vqXceQZRG2LQoVl2eHszw6OMa5bAHXW10tlUypwlNDkwgaSJbDv7KHZ4YzK2aEl6a31FK14i8lAJwFXgWkgd8H/lZE+ucp/6gxJjHr9t1GVFJECNoWkYCNawxHR3P8YHCMs5nVEVQypQpPXpwEtLO9FoK2heN5HLqUbfnfh6vpLXUdWu4sYYzJG2P+yBhz2hjjGWO+BZwCbm923eYTsPygIsDxsZUfVLIlxw8kOvy3pkKWRbZc4fhYrqX7205O6OgttXQt/2sRkR5gJ3B4niK3isiIiBwTkY+IyJz9QCJyv4gcFJGDw8PDdamrbcllQeXRc2MM5YotfeK4Urbs8IQGkrqY7pA/lytyIdeaOzVOFiuczWh6Sy1dS/9iRCQI/A3wZWPMc3MU+TdgL7AOeBPwVuBDc72XMeZzxpgDxpgD3d3d9aoy8EJQ8YyfLnjs/ATjxUpdP7MWcmWHJ4cmwRgNJHUiIoQsi6NjObIttuSKv7S8prfU9WnZM4aIWMBXgTLwvrnKGGNOGmNOVdNhh4CPAm9uYDUXFKh21Jccf2Xdpy9OMlVpzQ7Y6bW2jAaSurMtQYBDlzJU3NaZIe+ntzz9/tV1aclfjfgbYnwB6AHeZIxZ7GW9AVrqkmq6oz5sW4wVyzx2fpyT4/mW6k/RQNJ4Idui6LocGWmNGfKTpen0Vkv9+agVpFXPHJ8G9gD3GGMK8xUSkddV+1QQkd3AR4BvNqaKS+Pny22ClnB6ssC/nxtjeKrU9BPJ9KKNmtpqvLBlMVwoM5id9yfeELq0vKqFljt7iMgW4L3AfmBo1vyRt4tIX/W/+6rF7waeEZE88G3gIeDjzan54lgiRAIWrjE8cynDoUvNm3swXijroo1NNN0hf2JsislS8/rUdPSWqoWWmwFvjBlg4VRVYlbZDwIfrHul6iBgWdhiGCmWGTs3zvb2OBuTkYZteTs6VeKZYX+rWZ2Q2DyWCJYYDl3KcOeG9oaPoprQ0VuqRvQX1EQiQsS2sUU4Npbj4IUJcg0Y4XMpX+SZalpDA0nzBW2Limsa3n/i6OgtVUN6JmkBtuWnO3IVlx+fn6hbB70xhtMTUzw7nCUgQsDSE0irCNnCaKHM2Uzj+k+eH89T0vSWqhH9FbWI6fx50BZOT07x2Plxxgq120vcM4bnRnM8P5EnZFvYGkhailQXhHx+vDH9J2OFMueyBd1yWdWM/pJajN9Bb1NxPZ66OMlPhjOUnOXNRSg4Lk8OTXI+VyRiW5rSaFFWdUHIes8/qbiev/aWJQ3ro1OrnwaTFjU9N2UoX+bRc2MMTEwtOfVljOFirsiPzo2TKVeI2JaePFpc0LYo17H/xBjDc2M5yq6n/WWqplpuNJd6gYgQCQiuMTw/kedstkB/OkZvIkzgGieCbNnh1HiekUKZoGXpiWMFCdvCcLX/pC8dq+l7X8yXGM6XNL2lak6DyQpgi2AHbBzPcGwsx4nxPBsSETqiIeIhm4htYYCi4zJV8TibmWK86ACGsLZGVpzp/rPnJ6ZIhYO01WhDrYLj8txojqClvwlVexpMVpCAJQQsfwHJwVyBc7ki0yvIGGMuO0GEbcFf3kytRJYIFoZDw/78k+W2JDxjePZSBoPB1laqqgP9Va1AVnV+Sti2ZpZoCVf7WKZveuW58gVtC8c1HB7O4C2z/+TEWJ5s2SGkgUTVibZMVoF6j85yPIeiM0WhkmOqkqfkFik7RUpugbJbxvHKlN0yrnHwPNf/f2MwvHACtMRCEGyxCVhBbCtAwAoQsiME7RBhO0I4ECUSiBIJxIgG4sSCiTUfFEO2MF6scGIsz46O+HX9e1zMFzmbLegADFVXGkzWsKIzxURhlInSKJnSOJPFMbKlCbLlCbKlSfKVLPlylqIz1ZwKihALxImHkiRCaZKhNMlwG6lwO6lwG+lIJ22RTtLhDoJ2qDl1rDN/lQSLs9kCsZDNpmR0Sa/PVxyOjOQIaT+JqjMNJivY4GmLf/hShH/9ZojJMSHdYXjN68u86b4im/o9jDFMFEcZnrrAyNQQY4VL/m3qEmOF4cUHiepJPRqMEwnEiASiMy2JkB0maIUI2EFsCWBbAWyxsar9NVJdZs3Dr49rXFyvguNVcDyHsluk7JYpuUVKToGSU6Dg5P0WkFNgqpJjqpJjOH9hwSomw2nao+toj3TRFeuhM9ZDV6yXrth64sHkij6RTk9oPDaaI2pbdMbCi3pd0XF5aigDoJNUVd1pMFmhHv1OkI/9Vpx73lriU/8wRrDzAj85Nci/fuci9350iNtfd4pAx3kcd/7Z1AE7SFuki7ZwB+lIJ+lwO6lwO8lwG4lQikQoRSyYbFq6yfEcCpU8+XKGXCVDrjRJtjzBZLUVlSmNMVEcrbaoJsmWJjnD8aveJxKIsS6+ge74BtbFN9Cb2ERPYhPtke4VE2RsEYwlHBrOcmuvRTq88Aivsuvx5NAkZc8lbNsNqqVayzSYrDC5coYnj57how9e4uceOMpU22k+c/I87onqMvY3wY5NwmNP29xyh0N3OkV3fD1dsV46oj10RtfREVtHR3Rdy1+xB6wAyXCaZDi9YDnPeGRK436Lq3iJ0alLjEwNMTp1kZGpIYrOFGcmT3Bm8sRlrwvZYXoTm1mf7GN9cgsbkn2sT/QRDiwtldQoAcui4no8MTTJ3q4k3fG5WyiVaiApuhpIVONoMGlhmdI4g5mTDGZOMZg5xbnMKbKlCU4es2n/acMF24OsX7Y92k1vYjM9iY30xDfyP05vI/j0Jj7wB809hkawxKKt2n+yjT2XPWeMIVfOMJw/z6Wp81zMDTKUO8vF3CC5cubqICNCV6yXjcl+Nqb62ZTaxqbUViKB2k4evF5B28Lx/CHDO9wEm2ZtW2CMYXiqxPHxPGXHEA5oIFGNo8GkRZScAmczJzkzeYKzkyc4mzlJpjh+VbmQHWbsyE5++a3rubG/j97EZnoTmwkHIpeV6/lFi19/Y4oP/MFEow6hJYnITOtmW8flgSZXzjCUO8uF7ADns2e4kB3gYn6QkfwFRvIXeHro0ZmyXfH19KVuYHPav21I9hOwmvPnE7AESyyOj+U4OZGnPRykPRLkQr5EruwSsCAc0CHAqrE0mDSBfwV5ntMTxzgzeYKBieNczJ+DK+YShANRNqW2sjG1lU3JrWxM9dMVW88r39XBL//FOIEFvr3ejR6TY62bwmoFiVCK7R03sb3jppnHHM/hYm6Qc9lTDE76rcILuYGZAPPEhe8DYFsBNqW20Zfezpa2HWxJ7yAd6WhY3acXBPWMYaxYYaRQxhKpTlbV7101ngaTBii7JQYnT3Jq4igDk8c5M3GcqUrusjK2ZbM+uYXN6Rvoq179dsc2zHliSHcYhgYtNvXPv7Ls0DmLdEdz95dfiQJWgI0pP8V158ZXA36AuZAd4GzmJGcnn2dg8jgj+QsMTBxjYOIYjwz4r22LdLKlbSf9bTvpb9tFb2IztlXfVJMlQsjW4KGaT4NJHeTLWU5PHOXUxFFOjx/lXPYUrnf5Pu/JcBv9bTvZkt5BX9t2Nia3LnquxGteX+Zb3wjzax+efyOlh78R5jWvr91+KGtZwArMpLfY/BoApio5zk4+77csJ48zMHGcieIoE0OPzqTHwoEIfekdbG3bRX/7LvrS2wnZixvWq9RK05LBREQ6gC8APwOMAL9njPnaHOUEeAB4d/WhzwMfNo3c+xSYKI5ycvwIp8af49T4c1zKn7+yoqxP9tHftpOtbbvZ0raDtkjXdacj3nRfkfe+PsXLX1Nh7+1Xb/P77OMBvvX1MJ/9Zua63l9dWyyYYFfXPnZ17QP81OVQ7iwDk8c5Pe5fSIwXhjk+eojjo4cAv/W5MbmVre272da+my1tO4kFE808DKVqpiWDCfApoAz0APuBfxaRp40xh68odz/wBmAf/oqH/wqcAj5Tr4oZYxiZGuJs5ignx5/j1PgRxgsjl5UJ2EE2p26YuSLdkt5BNBivWR029Xv8/p/n+d13Jfj5t5a45y0lejd6DJ2zePgbYb719TC//+f5BdNgqrakesGwPtnHizfdDfij8U6NH/VbqePPcT53Zmb02PdOfwtE6E1sYlv7Hra27WZr+25S4bZFf+a1Jq0qNZd6XWtLgy/ir0lE4sA4sNcYc6z62FeBc8aYD19R9ofAl4wxn6ve/1XgPcaYFy/0GQcOHDAHDx5cct0+c/Az/PPx7zI6NX7ZeliRQMxvdbTvZmv7LjalthGwarNs+EIGT1s89OWrTyZvvFdPJq2oUMkzMHl8pgV7NvP8VenPrvh6trbtYlv7Hra176Y92j3ne82etPrzbynRu8ljaNDiW98I83D1YuIlr67/9r+q9WVKE/5vbuI5Toz+hDs2vojfe/l7ruu9RORxY8yBOZ9bbDARkT7g7JUppGqqabMx5sx11e7qz7kV+IExJjbrsQ8CrzLG3HNF2UngZ4wxj1XvHwC+Y4xJzvG+9+O3ZOjr67t9YGBgyXX7xPc/wbdPfIeQneSGjt1sbdvNtvY9rE/2zSwfotRiVdwyZyef5+TEc5wcP8KZieOU3dJlZdoindWLFP/3ti6+gXMDNu99fYpPfDE3b5rzd9+V4LPfzOhFxRpjjGG8OMyp8ef8zMnEUUZmLUXkGcON3Tfx+V/4z9f1/rUKJi6w3hhz6YrHO4FLxpiaDFsRkVcAf2eM6Z312HuAtxtj7pqjTjcZY56r3t8BHAOshfpNrrdlcnbyrL9ont1FSCeEqRpzPIfz2dPV/jc/NXbl+mmxUJILj++lU3bz67+6jQ3J/jlHjH36gSiVkvCbf9ikRTpVQ0z31Z2qplJPTxxlsjh2WZmgHWJL2062te9hc2on+3v3cGtv13V93kLBZCl9JgLMdYJOAMXrqdg8ckDqisdSzMz1XrBsCsjVqwN+c3ozY8VJJoqaPlC1F7AC9KW305fezl3992CM4ULuzExa7NTEc2RLkxwd/zE33/EoDz7mT2LtS2+nv30X/W272JLeTjgQ5Z63lPj1N6Y0mKwyFbfM2cxJBiaO+VMNJo5RqOQvKxMNxulv28W2aot29gTbsuvVbYXtawYTEfnL6n8a4I9FZPav0wbuBJ6qYZ2OAQER2WGMmV61bx9wZec71cf2AT+6RrlVx1+B1+B4/g6Lfvg0CELQlrrvcaLqT0TYkNzChuQWXtb3sxhjGC1c5A0fP8fL/qeDDGSeY3TqIifGDnNi7PD0i9iQ6GNTcgcTydsYmeqlM9qjExlXqIniKAMTx2fmp801zaAt0jlzMbGtbTc9iU1N+b4X0zK5ufr/AuzBH2U1rQw8AfxprSpkjMmLyEPAR9s4J4EAABxeSURBVEXk3fijuV4PvHSO4l8BPiAi38YPdr8NPFirurQizxjKrgEM8WCA9Ql/j3DPQMlxyVVcLuVLeMYjaFm69PgqItV1w9omd/Py1B285RaPbGlyZk7TwMQxzmVOcT47wPOXBrBf8V3+5PsO8VCSvvQO+tLb/fkyqW01HV2oaqPsljiXOTUz4u/M5ImrUlaI0Jvc7I8UbdtFf9vOeQdpNNo1g4kx5tUAIvLXwPuNMY2YvPA/A18ELgGjwK8bYw5X+1P+xRgzPTj/s8A24FD1/uerj606ZjqICGxMRtiSjhKZp99mR3ucC7kiA5kCRcfVbXxXmdmTVpPhNDf33MnNPXcCL6y28LkvDbCr4xix0GHy5SxHhp/gyPATM+8xvdbYxtRWNqe3sT6x5ar13VT9OF6FC9kznMuenlmL72Ju8Kphu5FAjM3pG9jStoP+tl30pW9omUVHr9RyQ4Mb4Xo74AGevuj3mQTtxo3eqrgeroENiTD9bbF5g8iVHM/jxFie87kiAcsioK2UVWHwtLXo0Vwbt7iMFS7NzNQ/O/k857Onr0qVIMK62Ho2pbaxIbllZkn+eOjK7ku1VIVKngu5s5zPnuZCdoBzmdNczA9e9R2ICL2JzTMtyL70DnriG2t6IVh2PbpiIfZ2X9/3WqsOeNVgfkrLIxqw2dOVpC2ytLkrActid1eSzmiII6M5Srq/xaqwtEmrQmd158lb178MeOGq+GzmJOcyJ2euii/lz3Mpf35mMUvwl/2Z3vOlN7GJnvgm1sU3tOyeL81Ucctcyp/nYn6wutXBIEO5M1dNagZAhO74+pktDqZXol7Jy+1oy2SJGtUy8Vsjhi3pGP3p2LL7PoqOy9MXM+Qrjqa9VolaTlqtuGWGcmf9VZKzA5zPnWEoe+aqeS/T0pEOumPrqxuvra9ukdxDe3Rd05bmbwTXc5kojjBa8DdgG5kaYiR/gYv5c4wXR65a+Rv8FaZ7E5tZn+hjQ2oLG5P9bEhuaUpArmfLRIPJEtU7mBhjKLkeYdtib3eK9BJbIwupeB7PDmcYK1SIaEBR1zA9euxibpAL1Q3FLuYHGc6fvzpNNk2EdLidjmg37ZFu2iKdpCOdtEU6SIU7SIXbSYRSLfnbM8YwVcmRKY0zWRojUxpnojjKeGGEieII44URxovD8y5HIiJ0xnpnWnA9iU2sT/TRFeut++rRi6VprjXC9Qxlz2N9PMLOzjgBq7YBK2hZ7FuX5rnRHEO5orZQ1IKmR491xXq5ad0L5w/XcxkvDjOcP89w/oJ/dV64yOjUEOPFUSaLY0wWxzjF0XnfNxFKEQ+lSITSxIMJYsEE0UCcSDBGJBAjYkcJByIE7TBBK0TQDhGwAtgSwBIbSyxkZtUJ/+TuGQ/Xc/CMi+M5lN0Sjleh7BYpuSVKTpGSM0XByVNwppiq5MiXs+QrWfLlDLny5PxBcpZ0pIOO6Dq6Yj10xdZXW2gb6Iz1rOpW2bWs3SNvIabaNyIi7O1Ksi4erttJ3hJhT2eCgAhnswXCtqVzUlYBY0zDLgxsy54JMnuuGJXqeA6TxVHGCsOMF4eZLI4xURxhojhGtjxBpjjOVCVHtjRJtjQJnG1InRcrEoiRDKdJhztIRzpJh9tpi3bRHumiPdpFe6S7bpP+VjoNJk023cmeDAXY250iGqx/c1hE2NERRwTOZDSgrETG+K1Yf/rXCyvBNntuUcAKzHT4z8fxnGpLwG8NFCp5pio5pio5im6BQiVPySlQdkuU3CIVt4xjKjieg+tV8IxXvfmtCMHfXVLEwhYbWwLYVqDaogkSsiOE7QjhQISwHSUajPu3QJx4KEk8mCRRbSVpoLh+GkyaqOJ5uJ7fyb61LdbQE7qIsL09jiXC6ckpDSgriFftV+uIBNmcihENWkRsm/O5Is+P53Ecj1ALpzADVoB0pKOh2xyr+tNg0gR+J7shZAs396ToiDbnakhE2NYWw2AYmCxop/wK4HgGx/PoT0fZ2ha/7AJgcypKdyzEsbEcw1Nl/T5VQ+m66Q3mev5VZXcsxIs2tDctkEwTEW5oi7M5GaXoenXbOEctn+P567Hd3J3ihvbEnC3JSMBmb3eKrmiIkn6fqoE0mDSI3xpxcY3hxq4ke7uTDZ1Fv5DpPpSNyYiegFqU6xkcY7ilO0l3fOGJbZYIe7tTpEJBSq5+l6oxWuNstsq51Rx3OhzkRRvb6U1EWi79ICLs6kjQEw9rC6XFeMZQ8Txu6kzQGVvcDGnbEm7pSRENWJRc3SBL1Z8GkzoyxlByXFzPsLMjwa09aaItvKmWiLCnK6kpkhYy3dm+syNBT2JpCzGGbIt9PWkEv2WjVD1pMKmT6b6RVCTInRva2ZSKtlxrZC7TKZK2SFADSpNNr4bQl4qyKXV9S2/EgjZ7OhOUPf0uVX1pMKmxmdaIMezqTHBbT5pYA+aO1JJtCbesS5EIBap7p6hmKLkeXdEQ29uXt/fIuniYDYmIprtUXWkwqaGK5/kngFiYF29sZ2NyZbRG5hKwLPb3pIkELErutZeYULVVcj1iwQA3dSeX/RsSEXZ2JIgFA5Q1oKg60WBSA54xFB2XgGWxb12Km9elFr3nSCsL2Ra39qYJWpaehBqo4nnYIuzrSdVsfTbbEvZ2JzH4v1elak2DyTJMp7Qcz7CtLc6LN7QverTNShEJ2Nzam8YSoaIBpe5cz+B5cMu6ZM0HayRCAba1xShrX5iqAw0m18HgpyFKrqE77qe0+tuWv+dIq4oHA+zvSWHwd29U9TG93tauzjhtkfpMZt2cipIKB7WlqWpOg8l1qHj+woy3r0+zt3t1pLSuJRUOcsu6FK6nw0zrwRhD0fXYnIyyfolDgJfCEuHGLr8fRr9HVUstE0xEJCwiXxCRARHJishTIvK6BcrfJyKuiORm3e6qdz17E2H2r0tze2+adLh2G1etBB3REHu7kzMLVKraKbmGjkiQ7R3xug/aiAVtdnTEqehwYVVDLRNM8BedPAu8CkgDvw/8rYj0L/CaR40xiVm379a7kj3xCN113G+k1XXHw+zu9AOKduTWRsn1iAT8nTUbtXLzhkSE9kiwuoy9UsvXMsHEGJM3xvyRMea0McYzxnwLOAXc3uy6qcttSEbY3h6n5GpAWS7H87AE9vekG7pWm4iwuyuJoOkuVRstE0yuJCI9wE7g8ALFbhWRERE5JiIfERFdUr9BNqeibE3HdJb8MrjG4Hpwc3eqKRNbowFNd6naaclgIiJB4G+ALxtjnpun2L8Be4F1wJuAtwIfWuA97xeRgyJycHh4uNZVXnNEhK1tMV26/jpN77C5szPe1G0INN2laqVhwUREvisiZp7b92eVs4CvAmXgffO9nzHmpDHmVDUldgj4KPDmBcp/zhhzwBhzoLu7e75iaglmlq5PRDSgLMH0mltb0lE2Jq9vza1a0XSXqpWGpYWMMXddq4z4vdpfAHqAnzPGVJbyEUxviK0aRkTY1ZnANYahfEl397uG6UCyLhbmhrblrblVK9PprqOjOSzR709dn1ZLc30a2APcY4wpLFRQRF5X7VdBRHYDHwG+Wf8qqitNL13fHdO9UK6l5BpS4eDMXI9WsSERIa3pLrUMLRNMRGQL8F5gPzA0a+7I26vP91Xv91VfcjfwjIjkgW8DDwEfb0bd1fTS9X5A0U75uZVcl1jQZt+6VMutliAi7Omsprv0u1PXoWVGPxljBlggTWWMOQMkZt3/IPDBBlRNLdJ0QHl22DAyVSasKa8ZZdcjaFns70m1zHbNV4oFbW5oj3NsLKfpSrVkrfmrVivW9OZaXbGQpryqKq6HJcKtvemWX3pnUzJCOhzUfWzUkmkwUTU3HVB64pryml5peX9PiniwZRIB85Lq2l2ILlWvlkaDiaqL6QUFe9fwsOGK52GA/b1pUitoHbdY0GZ7uy5Vr5ZGg4mqG0uEPZ0JNiWja27pFcfzMMZvkazEBUE3JqMkQwHKOvdELZIGE1VX/paxcba1+Wt5rYWRQhXXwzOwb12qbvuS1Jslwo3dSUDTXWpxNJiouhMR+tti7O5MUHE9nFV8tVt2PRDhtt407U1cJqUW4kF/Z8a13u+lFkeDiWqYjckot6xLzaxLtdqUXH/v9ttXWB/JQjanoqQ03aUWQYOJaqiuWJgD69sIWELRdVfFFa8xhqLjEbYtbl/fRiLU+qO2FkvTXWqxNJiohkuEAtyxvp1kMLjiUyjTa22lwwEOrG9rylLy9Tad7tLRXWohGkxUU4QDFrevT7OxOtJrJa5Y61X3bV8XD7O/J02oRWe214Kf7tLJjGp+q/fXr1qeVV1x+KauJG71Cn+lXPlWXI+KZ9jZ4de/1dbaqrXpdJcIa2JEnlo6DSaq6XoSEe7Y0EY8aFNs8fko0/0jtiXc1pNmcyq6Ztaw8iczxqmsoKCvGkeDiWoJ8aDf57CtLU7FNS2Zn3c8P63VFQtx54Z20pHVMWJrKTYmI7TpUvVqDhpMVMuwqlsBT3dkt0pfynRrxAB7u5Pc3J1c1f0jC5neu8bSnRnVFdbmX4RqaclwgDvWt7GrM4EHFB23Kamv6ZFaJdejNx7ixRva6YlH1kxaaz7RgM2uzoS/9liLtR5V86yeAfFqVRERNiajrIuFOZOZ4mymiGc8QraFVeeTuTGGimfwjL8r4s6O+KqZhFgrPfEwl6bKjEyViQTWdnBVPg0mqqUFbYsb2hNsTsU4k5liMFPEw59pHhCpaSvhhZn5QjoSYGs6RnskuOZbInMREXZ3JnisNI7jeQQsTXKsdRpM1IoQsi22tyfoT8cYzpc4kymSr7iAP7LqegOLV22FGPycb288zKZUVFsiixCyLfZ0JnjmUgZbjAbdNU6DiVpRApbF+mSU3kSEqYrLaLHMxVyJbNlBRDDGIAgiYAlIdSdoU/1fz0wvC+I/Y1nQFQ3REw/RGQ2v+vkitdYVC7MxGeVctqjprjVOg4lakUSEeChAPBSgLxXD8TzyFZdc2SFbcim5LqXqCsUC1eBiEQ3YJEI20aBNMhQgFrD1inqZtrfHGS9WKDrumh3lplowmIjId4EXA071oXPGmF3zlBXgAeDd1Yc+D3zY6BCTNSdgWaTDlr8RVbLZtVlbbEvY253kxxcmcI3B1uC8JrXqZcT7jDGJ6m3OQFJ1P/AGYB9wC3AP8N5GVFAp9YJEKMAOnR2/prVqMFmse4E/M8YMGmPOAX8G3NfcKim1Nm1MRuiMhijpYpBrUqsGkz8WkRER+YGI3LVAuZuAp2fdf7r62FVE5H4ROSgiB4eHh2tYVaUUvDA7PmQLlVW4+ZlaWCsGk98FtgEbgc8BD4vIDfOUTQCTs+5PAgmZo0fVGPM5Y8wBY8yB7u7uWtdZKYU/XHhvdwrXmJZesFPVXkODiYh8V0TMPLfvAxhjHjPGZI0xJWPMl4EfAD83z1vmgNSs+ykgpx3wSjVPWyTItrZ4Sy7WqeqnoaO5jDF3Xc/LgPmGhxzG73z/UfX+vupjSqkm2pKOMl4qM1FwCOv8kzWhpdJcItImIj8rIhERCYjI24FXAv99npd8BfiAiGwUkQ3AbwNfalB1lVLzEBFu6koR1P6TNaOlggkQBD4GDAMjwG8AbzDGHAMQkVeISG5W+c8CDwOHgGeBf64+ppRqspBtcfO61KxVB9Rq1lKTFo0xw8AdCzz/CH6n+/R9A/xO9aaUajHpcJDtHXGOjeWI2JauNrCKtVrLRCm1ymxKRuiJhylpumtV02CilKorf7n6JPFgQAPKKqbBRClVdwFLuGVdClsER/ePX5U0mCilGiIatNnbncTxdELjaqTBRCnVMB3REDs6EjqhcRXSYKKUaqhNyQgbk1GKGlBWFQ0mSqmGEhF2dMTpjAYpaf/JqqHBRCnVcJYIN3WniAV0hNdqocFEKdUUQctiX0+KgCWUNaCseBpMlFJNEw3Y3NqTRkCHDK9wGkyUUk2VCAW4eV0K1wPX0w75lUqDiVKq6TqiIW7qTlDxPFwd4VUXjufhGUNbOFiX99dgopRqCeviEXZ3Jim7nk5qrCFjDCXXBfxVCDalonX5nJZaNVgptbZtSEbwjOHoWI6wbWHpKsPL4hlDyfVIh4Ps7U4SCdh1+ywNJkqplrIpFcX1DCcm8hpQlsHxPBwP+tNRtrbF6/7vqMFEKdVy+tJREDgxrgFlqYwxlDyPgFjsW5egMxZuyOdqMFFKtRwRYUs6hgDHNaAs2nRaqz0S5Mau+qa1rqTBRCnVsvrSMSwRjo3lCFoWtqUBZT4V1x8Jd0NbnC3paMN3tdRgopRqaZtSUYKWcHg0h/EMAUsHoc5mqq2RkG2xvztNOlKfob/X0lLfiojkrri5IvLgPGXvqz4/u/xdDa6yUqoBehIR9nUn8Yx/Ba58rucHknXxMC/a2N60QAIt1jIxxiSm/1tEEsAQ8HcLvORRY8zL614xpVTTdcbC3NZr8fSlDCXXJWRZDU/ltApjDGXPw0K4sStJTzzc9H+LlmqZXOFNwCXgkWZXRCnVGlLhIHesb5tZbXgt7ofizpo7cufGdnoTkaYHEmjtYHIv8BWz8K/lVhEZEZFjIvIREZm3pSUi94vIQRE5ODw8XPvaKqUaIhKwuX19G12xMMU1NFt+eia76xl2diS4tSdNtIGjta5FWjGyi8gW4CSw3Rhzap4y2wADDAA3Af838FVjzB9f6/0PHDhgDh48WMMaK6UazRjDmckCz0/ksUUI2q18bbw8rmeoeB5tkSC7O5PEgs0JIiLyuDHmwFzPNexfX0S+KyJmntv3ryj+y8D35wskAMaYk8aYU8YYzxhzCPgo8OZ6HoNSqnWICFvaYtzW24ZtCUVn9aW9jDEUHQ8Xw65OvzXSrEByLQ3rgDfG3LWE4r8CPLDUjwCanzhUSjVUWyTInRvaOTKSZbRQXhXzUYwxOJ7BNYaeeJjtHQnCLd7yarnaichLgY0sPIoLEXmdiPRU/3s38BHgm/WvoVKq1YRsi1vWpbixK4mHoei6K7aV4lSH+4YDFrf1tnFTd6rlAwm0YDDB73h/yBiTnf2giPRV55L0VR+6G3hGRPLAt4GHgI83tqpKqVYhIvQmIrxoQzudkRBF11tRc1I8Yyg6LgB7OhPcuaGdtibOG1mqluyArzftgFdqdTPGMFYoc3QsT9FxWzr15RlD2fWwLaEvFWNzKtKys/wX6oBvqUmLSilVCyJCZyzMi6MhzmWLnJqYarmg4nh+v4gl0N8WY3MyuqJHpGkwUUqtWpYIm1NR1icinM8VOD1RoOi42JYQEGn4ZD9jDBXP4BlDyLbY0RGnNxEm2KItkaXQYKKUWvUC1RTSxkSUS/kSA5kppip+/0TQtrDrGFSMMTjG4HoGBNrDITYlI3TGQqtqWX0NJkqpNcO2hPXJCL2JMJmyw8VciaF8iZLnYgw1abEYY/AMOMbMjChLhAL0JsKsi4UbusdII2kwUUqtOSJCOhwkHQ6yvSNOpuQwVigzUiiTKzszwcQYgyWCJSAIs2OMAYx5IXgYzMzrApawLhqiIxKiIxoiHFj5aaxr0WCilFrTLBHaIkHaIkG2tcdxPI+pistUxSVbdihUXEqeoey4uGb26yAUsAhZFtGARSIUIBa0iQVtwvbaW9FYg4lSSs0SsCxSYYtUOEhvsyuzgqz+tpdSSqm602CilFJq2TSYKKWUWjYNJkoppZZNg4lSSqll02CilFJq2TSYKKWUWjYNJkoppZZtTe5nIiLDwMB1vrwLGKlhdVYCPea1QY95bVjOMW8xxnTP9cSaDCbLISIH59scZrXSY14b9JjXhnods6a5lFJKLZsGE6WUUsumwWTpPtfsCjSBHvPaoMe8NtTlmLXPRCml1LJpy0QppdSyaTBRSim1bBpMlFJKLZsGkzmISIeI/KOI5EVkQETeNk85EZFPiMho9fYJWYF7dS7heD8kIs+KSFZETonIhxpd11pZ7DHPKh8SkSMiMtioOtbaUo5ZRG4TkX8TkZyIXBSR9zeyrrWyhN92WEQ+Uz3WMRF5WEQ2Nrq+tSAi7xORgyJSEpEvXaPs/yoiQyKSEZEvikj4ej9Xg8ncPgWUgR7g7cCnReSmOcrdD7wB2AfcAtwDvLdRlayhxR6vAL8CtAOvBd4nIm9pWC1ra7HHPO1DwHAjKlZHizpmEekC/jvwWaAT2A78Pw2sZy0t9nt+P/AS/L/jDcA48GCjKllj54GPAV9cqJCI/CzwYeBuYAuwDfjfr/tTjTF6m3UD4vg/vp2zHvsq8MAcZX8I3D/r/q8C/97sY6jX8c7x2r8EHmz2MdT7mIGtwBHgdcBgs+tf72MGPg58tdl1bvAxfxr4k1n3/wNwtNnHsMzj/xjwpQWe/xrw8Vn37waGrvfztGVytZ2AY4w5Nuuxp4G5rmZuqj53rXKtbCnHO6OaznsFcLiOdauXpR7zg8B/Agr1rlgdLeWYXwyMicgPReRSNeXT15Ba1tZSjvkLwMtEZIOIxPBbMf/SgDo201znrx4R6byeN9NgcrUEkLnisUkgOU/ZySvKJVZYv8lSjne2P8L//fx1HepUb4s+ZhH5j4BtjPnHRlSsjpbyPW8C7sVP/fQBp4Cv17V29bGUYz4OnAXOVV+zB/hoXWvXfHOdv+Daf/tz0mBytRyQuuKxFJBdRNkUkDPVNuMKsZTjBfwOPvy+k/9gjCnVsW71sqhjFpE48CfAbzaoXvW0lO+5APyjMebHxpgifh79pSKSrnMda20px/wpIIzfRxQHHmL1t0zmOn/BAn/7C9FgcrVjQEBEdsx6bB9zp3MOV5+7VrlWtpTjRUTeRbXTzhizUkc2LfaYdwD9wCMiMoR/gllfHf3S34B61tJSvudngNkXRCvp4mi2pRzzfvz+hbHqBdKDwJ3VwQir1Vznr4vGmNHrerdmdxK14g34Bn6zPg68DL/5d9Mc5X4Nv2N2I/4IkMPArzW7/nU83rcDQ8CeZte5EccMBIDeWbc34o+U6cVPfTX9OOr0Pf8U/mim/UAQ+CTwSLPrX+dj/mvgH4B09Zj/E3Cu2fW/zmMOABHgj/EHHESAwBzlXlv9e74RaAP+PxYx8Gbez232gbfiDegA/huQB84Ab6s+/gr8NNZ0OcFPg4xVb39Cdb2zlXRbwvGeAir4zePp22eaXf96HvMVr7mLFTqaa6nHDPw6fv/BOPAwsLnZ9a/nMeOnt/4GuARMAN8H7mx2/a/zmP8IvzU5+/ZH+P1fOaBvVtkPABfx+4n+Gghf7+fqQo9KKaWWTftMlFJKLZsGE6WUUsumwUQppdSyaTBRSim1bBpMlFJKLZsGE6WUUsumwUQppdSyaTBRSim1bBpMlFJKLZsGE6WaTES6ReSCiPzhrMduEZGiiPxiM+um1GLpcipKtYDqFqoPA68CngIOAj8yxryzqRVTapE0mCjVIkTkz4FfAL6HvxDhfmNMrrm1UmpxNJgo1SJEJIy/deoO4KXGmMeaXCWlFk37TJRqHf3AZvwlw7c1typKLY22TJRqASISBP4df3fAx4A/BPYZY840tWJKLZIGE6VagIg8ALwNuAV/N8B/wd8h76eMMV4z66bUYmiaS6kmE5FXAb8N/IoxZsL4V3j34W+n+rvNrJtSi6UtE6WUUsumLROllFLLpsFEKaXUsmkwUUoptWwaTJRSSi2bBhOllFLLpsFEKaXUsmkwUUoptWwaTJRSSi3b/w9MH9jUrp4Z2AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_data_sine(x_[ii], t_[ii])\n",
    "\n",
    "ax = plt.gca()\n",
    "plt.fill_between(xx.flatten(), mean.flatten() + err, mean.flatten() - err, color='LightBlue', alpha=0.7)\n",
    "# ax.set_ylim(-5, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 426,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAENCAYAAAD34uk0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOy9eZCk6V0e+Hx5X5WZVVlH1t3HdPdMT8+oRzMgZISQBZJFBKyFxRqBMcjyosCLFA7HmjAOe3fHsexyLOHdMLDekA8EGLCDQBKYc9fGInaEdcwM0sz0TE93T1fXfWVl5X1nfvvH079+3y8rqyqrKs/q74n4oqqy8njz+97v9/zun2GaJmzYsGHDho2TwNHvBdiwYcOGjeGDTR42bNiwYePEsMnDhg0bNmycGDZ52LBhw4aNE8MmDxs2bNiwcWK4+r2AXmB8fNy8cOFCv5dhw4YNG0OFV155JWGa5kSr/z0W5HHhwgW8/PLL/V6GDRs2bAwVDMNYPux/ttvKhg0bNmycGDZ52LBhw4aNE8MmDxs2bNiwcWI8FjGPVqhWq1hbW0OpVOr3UmwA8Pl8mJubg9vt7vdSbNiw0QYeW/JYW1vDyMgILly4AMMw+r2cxxqmaWJvbw9ra2u4ePFiv5djw4aNNvDYuq1KpRJisZhNHAMAwzAQi8VsK9CGjSHCY2t5ALCJY4BgXwsbPYNpApUK0Giow+UC3G7+tNEW7DNlw4aN841aDUilgGwWKBaBUokE0gqGAfj9QCjEY2TEJpRD8Ni6rQYB29vb+OEf/mFcunQJzz//PN773vfiC1/4Qk/X8ODBA9y4ccPy2Ouvv46bN2/i5s2bGBsbw8WLF3Hz5k1893d/d9vv+Vu/9VuP/v7c5z6HT3/60x1dtw0bR6LRABIJ4O23gW9+E1heBpJJksdRM4xMEygUgJ0d4P594LXXgHv3+NpGo3frHwLYlNonmKaJj370o/ixH/uxR4J2eXkZv//7v3/gubVaDa4eaj/PPPMMvvGNbwAAPvGJT+B7v/d78QM/8ANtr0nI44d/+Ie7vlYbNiyo1YDdXQr/Wu3w54mLyuGgtVGrAdUqUK9bn2eaQDrNw+kEJid52NaITR79wp/92Z/B4/HgJ37iJx49tri4iM985jMAqK1//vOfRy6XQ71exxe+8AV88pOfxP379xEIBPDZz34Wzz77LF588UWEQiH8w3/4DwEAN27cwB/8wR8AAL7ne74H73vf+/AXf/EXmJ2dxe/93u/B7/fjlVdewSc/+UkAwIc//OG21/yBD3wAN2/exEsvvYQf+qEfwuuvv24hllAohFwuh5/+6Z/GW2+9hZs3b+LHfuzHMDo6io2NDXzkIx/BO++8g+///u/HL/zCL3TkPNqwAYBCfncX2Ng4SAAAXVDRKBAM0i3ldLZ+n3odyOeBXI5urlzO+r/NTWB7G5iYAKanD3+fVusrFvnepRJjLkJWEu9zOACPh4fXCwQCPAY0HmiTB4Dv++3v68r7/scf+o+H/u/WrVt497vffeTrX331Vbz22msYGxvDZz7zGTz33HP44he/iD/7sz/Dj/7ojz6yDg7D3bt38du//dv4V//qX+Fv/s2/id/93d/Fj/zIj+Dv/J2/g1/+5V/G+9//fvzUT/3Uib5TpVJ51CfsE5/4RMvn/NzP/Rx+8Rd/8RGJfe5zn8M3vvEN/OVf/iW8Xi+uXbuGz3zmM5ifnz/RZ9uw0RLZLLC6SuGsw+OhlTA2RkujHTidQDjMA6CQTybpAiuX+VijQQLZ2wNmZ4FYrLWAl1jL/j5JqB23Vz5v/dswSHiRCMnP52vve/QANnkMCH7yJ38SL730EjweD77+9a8DAD70oQ9hbGwMAPDSSy/hd3/3dwEAH/zgB7G3t4dMJnPke0qsAgCef/55PHjwAKlUCqlUCu9///sBAH/7b/9t/PEf/3Hb6/zBH/zBE383APiu7/ouRCIRAMD169exvLxsk4eNs8E0aWlsbVkf93qBmRlgdPTsWrvHA8TjwNQUiWBzU5FUrcZYyu4ucOECLRqAZLa9DWQyR8dX2oFpknhyOWB9nZ8xNgaMj/fddWaTB462ELqFp59++hEZAMCv/MqvIJFI4IUXXnj0WDAYPPZ9XC4XGppGo9dKeL3eR787nU4UmzWzU0Bfk/7ZjUYDlUrl0Nc1r6V2lD/aho3jUC4DS0tWTd3hoCtpaqrzrh7DIBmNjtKSWFujVQIwwP7WW7QQGg3+3QoeD58TDNIS8niU28s0+dpKhYe4uMTaERSLJJGNDa5lcpLv1wfY2VZ9wgc/+EGUSiX8y3/5Lx89Vjhs0wH4ju/4Dvzmb/4mAOBLX/oSxsfHEQ6HceHCBbz66qsA6OZaWlo68nOj0Sii0SheeuklAHj0nqfBhQsX8MorrwAAfv/3fx/VahUAMDIygmw2e+r3tWHjSGQyFNY6cYTDwNNP00rodoxgdJSfNT3NzyqVSGRf/zrw5puMZQiCQWBuDnjmGR6XLpHcxsYYh/H7eQQC/HtsjN/h4kXgxg3g2Wf5ezRKchSYJt1pt28Dd+9aYzM9gm159AmGYeCLX/wi/sE/+Af4hV/4BUxMTCAYDOLnf/7nWz7/xRdfxCc/+Uk8++yzCAQC+LVf+zUAwMc+9jH8+q//Op5++mm85z3vwdWrV4/97F/91V/FJz/5SRiGcaKAeTN+/Md/HH/9r/91vOtd78JHPvKRR1bJs88+C6fTiXe96134xCc+gdHR0VN/hg0bFuzuMr4h7iDDYNxhaqq363A4+JnZLFN6xeIvFoEHDyj0L106e4zC7SahjI3RMtnf5znQiTOT4REOk6jEfdZlGOZZfXJDgBdeeMFsHgb11ltv4amnnurTimy0gn1NbByJtTXGEgQeD3D5MrX2XiObJUlUKiSyvT0G1UdHGUB3uSjIu0VqhQLPRTJ58H/j44z5dKDJqGEYr5im+UKr/9mWhw0bNgYbpsnA9N6eeiwQAJ54oiMC8sRrWV+3kphhkMSef56xCHFbra3REllc7LwrLRCgO2tmhkH8ZFJZY4kELZSZGaYUd8mNZ5OHDRs2BhetiCMapeB09DhkW6sB77xjjS+4XMDCAi0OgCm19++r5+zt0Tq5fLn9mpCTwOtlplc8TrJKp/l4vU733t4eyasL1tljHTB/HFx2wwL7Wtg4ANOka0gnjliMsYReE4dkVOnEEQ4D168r4gBoCV29SteRIJtlYPuIbMQzw+ejJXblijXOIuteX+/4Rz625OHz+bC3t2cLrQGAzPPwDVABlI0BwMqK1ac/Pk4tu9cV16kUe2Tpwn92loK6ldvMMKjtz86qx0olEki3xw4Ioc3MWM9TF8j2sXVbzc3NYW1tDbu7u/1eig2oSYI2bACgppxIqL8nJuge6jX29mj9CJxOusweFrweiXicbqWlJVpR1SpJ6MqV7gb5DYNpxGNjdPlVq1xLh/HYkofb7ban1tmwMYjY2bFWjcdi/SGOnR3GDQReL11DJ7GQR0cZF7l3j6m2tRpw5w7fJxTq/Jp1eL10oVWrXbHWHlu3lQ0bNgYQ+/tWgR2J0AXUa2xuWtfh9wPXrp2ubmNkhEJc2onU670t7OtSRppNHjZs2BgMFApWF1EwyOB4r2Mc29tMuRWEQiSOswjhYJAEIu/RaPStMrxTsMnDhg0b/Ue1yjRY6dMmLqJeZ1Xt7DDlVRAOM0bRiTRbv/9cEYhNHjZs2OgvTJPEIdlMTieJo9ddYxMJq6sqFGJ9RicJzOdrTSDNrdiHADZ52LBho79YWbEKz070hDop0mmuQxAMds/yOYxAOtD1upewycOGDRv9g/SEEszNqUFMvUKhwKpwqfkKBDrnqjoMPh8/Qw+i37nT/TqQDsImDxs2bPQHxaJV2x8b63133HKZWr/EWjweWhzdJA6B328lqVqNa+lmJXoHYZOHDRs2eo9Gg9q+CG2fr/cpufU66y9kMJnLdXjVeLcgDR7FPVapkECGYFiaTR42bNjoPVZWlIvG4eh9vyrTZOW3rEE64/ajRY4E5iUluVRSRYUDjIEkD8MwPm0YxsuGYZQNw/hc0/++yzCM24ZhFAzD+C+GYfShgsiGDRunRjJpbXa4sNCzAUaPsLGhOtAC7JnV7YrvoxAOs+2JIJ9nBtoA994bSPIAsAHgZwD8W/1BwzDGAXwewP8IYAzAywD+Q89XZ8OGjdOhUrHGOWIxHr1EMmltfxKPM97Sb4yOWtuwZDLWoskBw0CSh2manzdN84sA9pr+9TcA3DJN83dM0ywBeBHAuwzDeLLXa7Rhw8YJIa6iep1/e73A/Hxv11AsslmgIBKxdr/tNyYm2BFXkEx2pZ16JzCQ5HEEngbwTfnDNM08gHcePm6BYRifeuj6etnunGvDxgBge9taTX3hQm+ymgT1urWK3eezuooGBdPT1nkgW1usfB8wDBt5hACkmx5LAxhpfqJpmp81TfMF0zRfmJiY6MnibNiwcQiKRWu/qOnp3scYHjxgai7A4Hy3pvt1AgsLnJgoWF1l08gBwrCRRw5AcwVRGEC2D2uxYcNGO5CJgBL8DQZJHr3E9jaHOgkuXOhPZlW7MAxaRcGgemxpaaD6YA0bedwC8C75wzCMIIDLDx+3YcPGIGJri1XcAIVir6cB5vPWuMHkpHV07KDC4WANiNfLv6UH2IBUoQ8keRiG4TIMwwfACcBpGIbPMAwXgC8AuGEYxsce/v9/AvCaaZq3+7leGzZsHIJikbMxBLOzvdX463Vr65FgkC1QhgVSuChtTKQKvVrt77owoOQB4J8CKAL4aQA/8vD3f2qa5i6AjwH4XwHsA3gPgI/3a5E2bNg4Aq3cVZOTvV3D8rK1W28/5oOcFc3t6SuVgSgiHEjyME3zRdM0jabjxYf/+0+maT5pmqbfNM0PmKb5oL+rtWHDRkvs7Ch3lcPRe3dVImENMi8usnfVMEIGYwmamzn2AQNJHjZs2BhylMsHs6t66a4qlayzOSYmhiPOcRQiEWsRYTpt/Y49Ro+nrdiwYeOxwPKycqv4/Z3vlttokKAqFR5SeAjQullaYnzA7eYM8WGKcxyFiQl+X6mQ392lNRWP93wpNnnYsGGjs9jbA7Ja9nwn3FWNBtt1ZLNMVy0WD3fZ7Oyo3lmGwcFLy8t0/UQiKntpWDE7SwJJJvn3+joJpMctVmzysGHDRudQq1lngE9Nse34aZHNUkju71uti8NQKFibLk5MMN6STPJYXeV6RkcpbIc1BnLhAglE6j4ePFBWVo9gk4cNGzY6h7U1NYvC47H2aToJUim6Zo6a7e3z8TM8Hqay1uvA22+TGKpVPj4xcTArqVDgsb7OKu7JyZ4K3Y7AMJiBdfs24ztSA/Lkkz2LLdnkYcOGjc4glzvYav2kMzoyGRJQq3neXi+JIRTi0dxaZHlZuW6cTuD6dWrjxSLXlsnw0N1dqRSPQIBEF4mcbL39hNPJGpDbt0mW9TprQJ58sicDrWzysGHDxtlhmtZW69HoyQRxpULSaO7fZBhsEhiLWVt1NCOdts5CX1hQLqlAgMfkJAVsOk2Sy2TU8wsF1k4EAowp9HqO+mkhY3PffpsWlkwivHat6327bPKwYcPG2bGzo6wFh+NkrdZ3d0kcunvJ4aCwn5w8Xouu1axt1iWe0QpOJ/83NkZ3jwTX5bMLBQrfaJQZWsMQXA8E2OTx3j2SeLHIGpAnnuhqXY1d52HDho2zoVKx1nTMzLQXiK7V6KdfWbESRywG3LhBC6Ad98vKimrX4XZbayGOgs/H5964QZLSBW0qBbz5JlurDPA0v0cIh60z4HswSMq2PGwML2o1ZuMUiwfz/Q2D2qvXSyHh99NP3ss52Y8LdKvB72+vBUkuR+1Y79Hk91OYn6RV+/7+wSpy1wnFmttNSykeZxBd4jaNBklxf5/ZTWfJGusFYjErkSeT/G5dqnGxycPGcKFYVP7qVkHVZuj1BobBrJpIhG6LkwoZGweRzVqF98LC8a6SvT26mXSNfnKSlsZJyL1Ws8ZZxsfPFvB2u0kSExN8X2mtUiwyKB2Ps1J+kHtjTU+TkGUA3vY2v1enizRhk4eNYUCjQYGTSKgb+jQwTZVxs7ZGApmcHHyNclBhmtb2GGNjR1sNpknNfntbPeZyUWCfRuivrFjTgjulYQeDzFja3eV6Gw2ufXOTe+fixcGOhczPk0BkfsnaGgmkw0WENnnYGFyYJgljc7N1C2rD4I0eClF4eL0URqLR1moMipbLqipZf++9PR6RCAXPIA8HGkTs7lqD5EcJb5lfrlspfj+Duqcp1GvlrupkdpFhULGIRBg7kGK8fB546y1+3qD2ypJBUnfvqnXv79vkYeMxQSpFrVbaaQsMg5kwsRhdUMe5OXSNtlJRaZp68Vk6TY1yYoJmv+3OOh7V6sEg+WHB7UaDgXE9NTYapYA7TQyqlbuqW6m1Xi/TXre2+H1NU80ImZwkYQ6iG0sGSb39Ni1rPZjeIdh3iY3BQqVCwZBuGlUvfttY7PTCXSqOJyZIHjs7qj+QafJvCY4OS55/v7CxoZITfL7Dg+T1OlNI9fGpZxW63XJXHYV4nMrK0pKagy4t5y9d6klR3onhdHa13sNOPbExOEgkgFu3rMThclE43LhB8uiUVRAMUvN96ilra4pqleb+6mrfh+0MLAoFa0He/HxrImg0DhLH9PThz28H3XZXHYVgkPslGlWP5XJ0Yx3VRqWf6OK5scnDRv8hbgC9jTdAd4SQRrdSbAMBdl29fNmqPe7sMMOm2W1mwxokj0RaW2mtiGNu7vS9roDeuqsOg9PJvTI7qx6rVukeEiv2MYFNHjb6i0KBmpuuTfp8NLd7qVVGo+yFpGuVxSLXpgvAxx3JpDofhtG6ktw0qQzoadLz82dPF11dVe6qLtYvtIV4nEqHWMKSEKDHgc45bPKw0T+kUtTYxIcMUJt86qmTFYp1Ci4XtUq9VqFWA+7ceey0ypZoNJi6KpicbJ2y+uCB1fU4O3v22eWplPUa9FKxOAwjIwe72G5ukkSGoSr9jLDJw0Z/sLXFDBxxUzmdDDwuLva/CnxiorVWubPT33X1G1tbyo3ndjN+0Yz1dauQj8fPPuWuXre6q2Kxwel+6/WSQHT3WTLJuFk780eGGDZ52OgtTJOxDV2DlRtwkPLmQyGuye9Xj62uWgvcHidUKtbvPjNzUPNPJNR4VIAkrMcGTovVVWvvqkEbKet0Mi12YkI9ls2e+5iZTR42egfxheuZOq1M/0GB5PjrrcDX1uiaeNwgldYAkwxiMev/MxmrdRCJnKyz7mHIZA7OCBnEOhzD4Np0siyVSCDttNEZQtjkYaM3aDRoykvLBIAC6MqVwRQGAhm4o6fzbmw8XhZIPm91RTXXaJRKVArEzx8I0AV51uK5ev1gq3U9oWEQEY9bv7tkYp3DpAubPGx0H/U6g8569s3UFIvxBrE6txniltD92mtrVgvqPENPzY1GrURarzN2Jf59GU7UibjV2ppy+7hc7bda7zdGR6lwiFtPJvw1F74OOWzysNFdSL6/XkQ1MzN4fuvj4HAwE0vPAltetlpS5xHJpLp2hnHwui0t0fIA1DnqRLV1NnuwEHGQLdRmjIxYky6kRUvzpMQhhk0eNroHcVXpJvv8fOssnWGA9AvSu/Dev38uXRIAjk/N3diwatOLi53pUFyvWwcZRaMdb+rXEwQCjOdJ48dWMb8hhk0eNroD02xdYXzWfP9+Q2IgIkRNkxrlecyq2dmxuo100k+nrYkDU1OdE/Dr68PprmoFySTUE0KWl89FzMwmj6Owvc0MkvOqWXYLrSqM5+a6MpCmL3C5rIH+Wo1EeZ56YdVq1rRbPTW3UqG7ShAOdyYlF+CekUFGAC3VQWw6eBK43cza062yc5C1Z5PHUdjd5fH228Drr1MjEv+ujcOxsmKNBczMHE0cjQazUsplpjWWy/x7kIWx10v/vgT8i0WrQB12NHfNHR/n76IY6AHyixc7k/jQaFizq2Ti43mAy8UYiB4z29iwugWHDEMUgeox8nlr24xKhZrY1hY39eSk3ba7FTY2rD7dyUnl7mg0eF7zefa0Kpd5HFWJ63RSUMss8lCIdRf9bk0BcC2Li8o/n0pRmxzWmI6gWLRq/3pq7tqaNYB+6VLnAtlra+qeczq7MoOirxCXpz7bZGuL90UnamJ6DJs8DoN0W00mmSGhC7h0moffT6160HPPewWZ+ieQMa87OxSsudzJe/7U6ySa5vGzfj9JPBIhmfQr5TcWo7AVH/bGBtczzIqFrg2Hw6oVSCplbdEyO2stoDwLzqO7qhUk6eL+fWWd7+yQQIaMLG3yOAyGwXS7kREG7NJpCkc9u6RYpBYRDPJG0vPfHzfoFcbibiqVgDfeOP61hkHt1eHg7zKtrV4/nGyKRR5bW3zt2BgFeT/mkc/OktwkxrO0xOaOpxmv2m9kMtY9Lqm5lYrVpRSNdi6G1ZxdFYmoCvZ6nZ9dqdCVWatxT8i+cDio0TudPN8eD0lnkOuHxGLTx/ImEvxOi4unW7ucp2pVHbWa+jk93fFmo0NJHoZhfAnAtwF42J8Z66ZpXuviB/JmiUZpVu/s8GKLkMznWQQ3Nsab7TxqTEehWKQmVa3SUisWqTk2WwsALYZgkIfPR3fUUedLYiHlMt9PZpHrpFKr8Zrs7PD9Jyd5LXrVYFFmRr/1lrpZl5ZouQ6yEGuGadJ1JBgf5/k0TQp3fXpfJ7Vkya4yTZ4/j4dKmcS/TgrD4LoDAR4jI4PX/kb2jMOh2q/s7fEcHFY8W61SISuX1c9K5XjXL8DCRZs8HuHTpmn+655/qterahU2N2lqiyBLJqm1zcwMf0pqu6hWKTS3tqhFuVzc/CK4DYOuj2iUGuVJidXt5hEKKW200VAacjqtmuYBFDjSeHF8vLPTB49b56VLVCJMkyS3sdG5LKReYG9P9WFyONTgpu1ta+bcxYudO6f7+8xUy2aphMXjVvfVaWCaB12dbjdJRFydgxAzMwxFFBInTCZJ0vE4CUIs7FJJkfdpcJbXHoJhJo/+wuVSA270NtT1Ots5pNPcGOfZCqnXga98he4q06TAmZvjufF6KbzHxzsvvB0OZQmaJgXP3h59yGINSqrpzg67ncbj3SeRUIgCV2IGW1skzmFwZzYa1kFG8Tj3bj5vfbxT7o9slufn1VcV+Ytw12EY3EvijtLdmwD3YKPB6y2aeCtBKVZxMqlc0mNj3EP9JJJ6nVbB7q7K5qxWub7Z2fYtV4eD50fOkxwuF392wfIyzCEcWvLQbfU0AAPA2wD+iWmaX2p6zqcAfAoAFhYWnl/W/bUnQb3e3ubKZilE9VReyRgZpFbjnUI6DfzX/6o0JmldIdlV/UgiqNe5Hr24TeBwkOinprovLO7eVdk0Hg8nFA6CpnsUNjZUsoPbzfG/APDmm8p1FAqdzRXXaFB47+xQm15fV+dJ5rmIm2lkhL/7fCf/vFqN7y9uzmz2cLeOKCLj490n+UaD65KMw+aMTrHeBaEQCUSseIeD50MOyUIUwugCDMN4xTTNF1r+b0jJ4z0A3gRQAfBxAL8M4KZpmu+0ev4LL7xgvvzyyyf/oESCN9X8fHsEYJp8vl5cBVBgnUSLGGTUarSs7tyxuhcWF4FnnhmMzDPTVGmzze2w3W5aB7FY965HtUqhKxrw6CgF46CiWmVig1htFy7w/KysqGvsdDIJoNXkwONgmioTT6yMTEZZaF4v986FC92LTRQK/Mz9/daxOICfLfGyTpB9uWwlikLh+GzD7W2u0ePheZmY4HkPBPqSgHHuyKMZhmH8CYA/NE3zl1r9/1TkUavxhhKNJRxm1lU7N08ux4Cprv2OjHQ2J74fSKUoUPb3VadVh4PVs+9+92CSYypFQm8mEb+f17Nb425TKQZ9BRcvDm7B24MHKmjr91NYZTKMRQiEUE6KvT2ef/1ekISCkREqG7OzfP9eoVLhHtZjPDqcTgrtycn2NfpajeSgk0U7cQY9uB8M8ufenjUlOhxmQWofJmweRR5DLMksMEEXVudQLvNiCXlkMsCtW3TJxONHC8pQiK6KBw9ULnc2y8DyE09Yp9MNAyQLZ2eH50U0xkiEhHj9+mASB6BiI3t7XLdovsUiOwfEYipO0+nPHR9Xbr2VFQrLQYuBFYvWYUvz8wfnaESjJyeOfJ4Kht5NGeD3L5e5byS9ttcFch6PcmEWCrxGe3vK8qrX6T3Y3ub3jsetSqNeeySE0W5WmNersg2DQcqCZlIIBLgfJdYkRN6pVvcdwtBZHoZhRAG8B8Cfg6m6PwjgswCeM03zTqvXnNptVa/zAjbPrg4EqCm1QwKbm9aAo8yG6JbG22mUy0zDLRRUPn6jwRsqFqOWOizWVKNBgSBVvQKXiwRyGs36uM/TYwaRCK/9IEGPz8j67t9Xvne3m8pBu9e4Xm8968Tt5p6Re0pw7dpg3Av1uorH6HFLSYX1+2kB1OvtE4XTSYKQrghCCu2iWXaEQrw+PYyfnTfLww3gZwA8CaAO4DaAjx5GHGeC00mtKBajJia+0kKBVkQ7Vsj0NDfN/fsUJjIY5uLFwYgPHIV0mu4Fsb42N6k5xePU3i5fHh7iAKi1TU/zeq6uKquwViMpJpOM3XTKt+xw8P3uPNya6TQ13E6T1GmRTiviAEig0lFBsLjY/jUWt6aeOm0Y3C/xOIXuW2+p/8Xjg0Ec1SpJwumka3Fnh1bq/v7BGEU4TJdW8x4xDFVXotcxnQXT03xfsfRzOcoOfdBUHzF0lsdpcGrLQ4dpqk2ln7NQiERwnMApFGh66jfWaf3IvcD2trVgbG+PG1n89oO89naRTlPY6f54p5NCVBoBdgLNgeenn+6/+8o0aRWJlj0xQWF165ZSFsbH2ysGrNf5HfVRtQCVo7k5KhyNBolDPk9mXXTS3Wma/Bw59DRe/ZAKbKlcP0wGikur2fVmGHR5LS7yOwppdMt1u7NjneYYCPRsfPO5D5gfh46Qh6BUopaqb6h2U3LLZWoOutk7aELYNGll6X7wapU3omhSk5ND2citJQ5zTUYivKadEPLN7qtolFZbP7G7q9rJOBxMzX3w4OQpxvk8rWqdgDEz/tIAACAASURBVCUtOhBQQnp5mZ8pxHTxIs+ttBo5jRzSX3fa9zgOLpeqJapUVIqsYfCYnOxNDZF+vQC60fRJhV2CTR6dJA+Am2lri24c/fxNTlo7kLZCrUY3hp7lsbjYWU33tJB51Ho1scdjFQxnzfUfVORyFHC6v9vlUtrlWZHNKvcVQPLol9uyXmcmoWQDzc7yu+pB8qtXj6972NqiQCsUVCW0309lSCcdPS0X6F8d0GHQOzd7PNZaCl04S9cA/f6Q10sAvpsB7UTCeo18Pl6nLlqxNnl0mjwErVJyg0FmkhzlxqrVaIHo+eb9tkCqVa5JJzXpGCvrdLsZIO+3y6VbkLGrzVbI+DiVgrP6mZeXVSDZ7ab7qh++6/V1VYvk8dAFcvu2sgqOsiwbDcY2bt2iIBU3rMSTmrsJy+AoSVAYGene/HqxBqRZovx0udRPqbiWw+M5+TXIZnkOW2WSTU9zv3RLuUomrXNjvF4SSJdqQGzy6BZ5ACoDSR9+5HKRQI7S3JoJxDCojTa3Z+gFSiWuRSfB2Vn+Lb56wzg4zOa8IpvlNdXPh9dLV8tZWpA3a/z9cP9VKhT8IswvXqSLUtxVPh8VBF2DFrdNMknSWV62xu78fu4XUSr0liIrK/y+Tiefd/06n+dwqEMEbTsC1zStzxPCkKOXSKVaD4jz+ViI2q3OEvv7JBCR3aIAdKMFiU0eXSQPwfa2NZiut+s4DM0uLIeDm6CXArpYJHGIMDAM1RZa13DO0xjZdtAqCGwYFJJnOQ/NmuOTT3ZuJkY7WFpS3ykQoJas+9L11Nnmdi+5HJModJkxPk6rWdJR/X4VE1hdVVacYfC9e/ldewHTVIWQOqEC3R3VkEox1iTXwu2m7OhwDZlNHr0gD4A3l7QmF4yPs5L5MK2oWmWxmj5B7dq13hQSFgokDtGEHQ5aTF4vM2NEOx309hrdRDJJ4ar3RgqHz9b08t49NTNDKrp7oTXn83RPCS5c4HeT6zw1RSWhWqWFoY8dSCbVwCu3mxbyU0/Rcmrl59/f570gOO/KR6NBotzaOthHS2a8d3rWTCbDGKVcI5eLdSAdJOijyGNwyhXPA0IhapL6JkkkKKAPa8wmGoMIonr9YEpvNyAzSIQ4ZETmyIh1Q3q9QzfhrKMYG6OQ1G/ITIbkqtdInAQLC0rgNo987Sb01Gupupfr7POpDtFvvKGm20lySCJBJeLCBbqevvM7uS9aEYdkJOqfdZ6JA+B5iMeZtTY1ZVUGZL/cv3/QxXUWhMPWmg/xZDQH9LsEmzw6DY+HloPexyibpcZ3WGWq12utHK1UrAK80xCLQwhNiCMUsnYGdjgYhxmAgqS+wuvlNY3H1WOSYNBc99MOPB7rnPNWLo9OY3+fljGgYhK6kAmFmE6sV983GiQNn4/7Mx6nJf3UU4dbxpKxpysfvexb1W9It4IbNw4mwOzvM9704MHphly1gmQ/SlZYo8F9qcdguwSbPI5CJsPjpELc4WAgUh8EVCqRQJozNASBgNU1lM9btbdOQWIcQhwuFzdfMKh6/AgWFoavD1e3IPGO5uKsrS1qeycV/lNTKsApLT26BdO0pspGIsraKZW41xIJayM/j4ffeWKCloPDQYXouNqCBw+syof0sHrc4PEoK605LXlvjySyvHxwdMBpEAhQuRHvhWnSytHv5S7AJo+jsL5OQfvNb/Ln9vbJNIZ4nDePmLBiVuozonVI517B/r61t81ZUSy2dlUFArRG9MBpLDZYxYuDgnCYAkEPguZy1NoPu66tYBjWTKtkUlkGnYY0tAQo+CsVFQzf2LAqCF6vSkvWXS9TU2ps6mHY3LRqvAsL/ZkpP0jw+2m9P/mkNY1Z2tS/8cbB+qLTwOfjZ0gDRxkdLHGqLsAmj8Mg/Y7W13lD7O9TO3zjDQqKVm2+W2F0lFqBblbeu3ewcZxAWkELmm/I00Kq23XiuHqVN3e9bs3ckHblNlpD4lQyphXgeb1372RurHDYms4pExk7iVpNDXkCuHZpy5JI0H0m6bLxON1Tu7vWvT03d3xthrS+F0xO2sqHjmCQe6Y53V1I5NYt1YD0tBCXua4MrK11zaodoq52PYZpUsBms3RdSeMz6axZLPKmlIrasbHDs2+CQV7Ue/eUBii58rrvWzA3R01EArIPHlCrOG0et/jn9YIusThkLbIucTUMUOvngYRhqJGsS0vq3G5t0YJop98ZwGudTqspc4kEFYhOYWNDuSgdDr7/+jofi8W4f4NBFZfQ631kxvZxc0gKBWv6cTcLAYcdIyOUBdksr41ube7v8wiFaOmdpgrf7VayRt57e5tKRIdjT7aEOAzVKm+uy5dVsNDhoHC4e5dsnsvxhl9bA15/nYHCwzJwfD5eVN2M39iwNjwTGIZV+IhlcJoAunTxFXIwDGs63+7uwS6q3Zrmdh4xMkI3lu6SEDdWOxZjq+B5O0OE2kFzJtfOjko7lil109Pcl6ZJl6ZOHJcuHU8c1SoFlexNj8fqqrXRGkIi164dLAzO5ShL3niD8uak+0Hc0Tr5dMF9aNd5HIVCgVphKqXMyUaDAcZMhtqDy0XXQzSqAoMyzjIWO6jBt+ofFYupwrzmz799W7kyxsZIKu1CMi/0LJtLl9Sman7/iQnbXXUWbG1R+J+031lz48ROXQeZ1dFokDTEFSUKxNNPU8ttjoVJll1zq5FW6377bXVvOJ1ns5AfZxQKagRts0w2DDWQKxxun5hNk9fd5bIm75wAdpFgJ4oEZXRlMqlulnqdJJBK0c00OsoLLPENl4vCY3LSmnFimjTzdY0/Gm2tse3tWbOu2m2iKBkXuvar98+q15l7LgJLMjZsd9XZ0KrfmWTSHTXCuHls7fXrZ8t0S6dpEYjVWiwqt+rCAvDe9/LvVsTRTpcD0+R6JUlACOk4whkGSMt2vY27tHeXo1luStxI+mnJ4XarflrtdMCVlkDN2W8Ct5tyZnS0J50obPLodIW5jO5MJJQ/uVwmGWQyvIHGx9VmcTrp9pqYUCQiWoEeOD9sVrHeUM/hoHZ3nGDRZ0gAByt833lHEYvTydz9duaz2zgetRqvmU7cMhjqKDeQPtVvZITB1dPANBmAzee5D/b2lKCZngb+6l/lekolWg56EkW7Uy71uefA4HSGbgeNBu/XUolHuawmBlar3WntDvCc6517vV7+9PsPEkujQXmyu3t4er/bTaUzEuF+6YLiZ5NHt9qTNBpqdKW4BOp1Co1UimQwNqYIw+XizTsxoSyMtTVrOl2rUZONBt1L8hmtmtfp2Nqy5vVL2wlB83CZS5e618TtccbOzsFeULEYU3Rb1T6USnRfyfNP27Z9e5vWjygdoRA/b3IS+I7voLAql0kcEugXP3k7rS2a92w8fmq3SNdRqahZ49I2vlMFep2E283rok8iFEtRlNVk8vB6IoeD1zkUIpEEAh0hE5s8etHbKpWi0BYtwTRVvCQatQoBn483mzzWPKtY0vp0AVMqWftNHTblrbmnUHOcJJdjIE60rEiEN3+tpsxx6VwqprfLRW3J41FzD2y0h0KB10MXWEd16NUtRq+XcYmTBJ+rVeAb3+Bn5nLcN5EIlZann+a1rlRIHOJaczhUoehxaFZM2p022AvUarz/5CgUTpd8oLuZ5JB7obkTMGAdaCUTDOt1NbWwVuO5PmnCi8ejZqCPjFBuZLO8x1Opo7+bYfD5Qkby+hPCJo9eNUYE6HZYW1NWgpBIPk+tU7+AUhTo9R4c+9pq1GTzMJhmzTSXo/9armkoRItDbqRslpqtaC8+H+MgJ82McTrV6M2Rka6ZzOcG9TotPd3NI6m+8bj1/NdqJHdxh87PH92ZuRl37wKvvkqySiRUwD4eZ0yrVrM24jxJJ+fmPXpYnK4XME3eY/k8930+fzKLQlxG4j4Sd5LH0729LCQiLrNymd+hVGqPWJxOdb+NjPD90mkexxUZntI6tMmjl+QhSCapoYl2Jz5M06Q1IBtUirPicQqX5lGTetNEgBqlBNpdLjUfoVxW7U9yOW6s5mZ0q6sq80paqHRiiIxhULuRQN55HRZ1VrTq0BsKkcB1a04X0k4n+yS1E2xNpYA//VMqB8mkSg+Xqnink8QhguYkQe5m4hgZ4d7sFXFIlmMup452Ba4oOn4/D59vsJQd0+T9K+61QoHf9bjv5/Hw2oXD/E7FIhVEOT+6q+7bv93am61N2ORxWvKQJnEyq+CkPXoaDb7H1payBmo1Wider9VN4PPR/C+XrdlVzaMm63VaD0JK4TBjKC+/TPKpVrnOCxesxJBIKHeIy0WrZWKCzxHTXDfHTdNqfusa03GugFCI7ozR0cG6SQcBcn314jCHgxaGBJwl4C2adDtDo2o14E/+hARSLHKPvfACr8XCAhWWO3dON3ysOUYWCpE4unltazUrURQKxweyRYkJBpW7ZlhdrGJZyffPZg+/7+p17hWHQw3aajRUQkCpBHzrt9rkcRqcmjxef92acun3KyIJBtv3IZbL1Dj1AsJSiYI+FLJqb5OT/By9VUXzqMlsln7tTIYma6WiNFMZ5iTZWG43H9vaUqb6zMzZKoAlCJnLcQ2HtWlxOumqm5wc3pu4G5A255ubVoEYDvPaeTzW1F3DYLzisHPYaABf+5pKzd3YAJ5/XtUFXL7M/+m1Re0mSTTH47pFHOWylSza6fXkdqsgsRDGeS5OLJVUgk4iQetE5EgzPB5ah6EQz8spp1YeRR5ttycxDGMBwKrZxDaGYRgA5k3TXGn9yiFFtXqw46WYgKLBi0ksGzcQsDYmk3zxapX+4UqFGly1qvLF19dVMM40eZP7fLzwyaQar7m6qtqWZLN0I2SzJI9Cgdql280NMj1tTd976y2VIhoKnT0zRnzDEm+pVrmOZNIqoOp1bvSdHX5+PG536QVUvCMSYVaUCMpMhlbl/DwFfyhEQSpdcVsN5DJNxjnEWt3b4/Ok3mhxkZ+hX5fFxfaIQ58ECHSOOExTKR/iimqnK7Eob3J0aW73QKBaVVaDyB3d6o9GKWvk/DWTbbVqVeq6QKon6W21BGAawE7T42MP/3e++i5LTEA2eKuGZbWa6kwq2Uv1uspUkkCc16tcXqOj1OZ0t4U+RRDg++RyJBHJHBFrIx4nSTQavLkTCa61UgHe9z5aFT4f11Ct0k8tN6bMVu/0RnK76XIZH1f+9t1dawAzmeQRjZK87CpkKhtPPUViECFdr5MIkklqi7JPZB5Hc2B7ZYVxsFqNeyQYVNl1i4vcm3q9yezs8fUY0pFVH8F7WA1SO6hUrFlQ+Xx7LqhAQGUaBYPtxX2GCbUa7xE5JIheKh0+PE4g50esCpdLuZhrNeVxEHTh3J3kHQ0Ara54CEAHx2MNCJxOasuisYvFIUJwf58XuV23n8ulMjrEVZVKKasD4GPyfqUSiUGC33LTbm2RQBoN1RZcUvpcLpV9AVAg7e1xI/l8DJqWSqoK9hQwTRN1s45ao4ZqvYqG2Xh0mLI9RpwwRqZhZLNw7ybhyObgdDjhNJxw7e/DSKUowKanz7f22A4k3jE6ah0SlMnw2tfrKha1tsYCUcHODoknmVRt1q9e5bUdH6eQ1jO8pqaO93vXanSX6crN6CgJqR2lo1q1Bn0LhfasCqdTpaWKG2rY42XivWg+hCxO06vO6VSBfz0RQL82jYbqfCEpvaepFzoGx5KHYRj/4uGvJoCfNQxDV8GdAL4VwDc6vrJ+o1aja0huBD1YFYnwENNQD0zpGoPenkDPFTcMvt7vVxqnx8Pn1uu86ep1FVt5mMq3nfLgP/0XP/7zV/zYz7oRCY/he17YwQe/JYOJxXG+lwQJs1klOB66zar7SeR2VpCv5JB11JBx1ZFxN5D2NJAxSyhUC8hX8ijWiihUCyjXyijWiijXy6jUK4+Ok8JbrmM0VUIwTyHidLjgcrjgdnlQGI+gMj4Kn9sPv8uPgDuAgDuAoCeIoDuIkCeEkCeEEe8Iwt4wwt4wIt4IfC4fjPPk3w6FSO66FSI1A6urKnNuf5/CPJ3m49vbfE6hQDKORlX6qV6PEYsdH+cqlaydnwEmVczPtyYOyRAqFtV90u5QLEkYEbIQa3lYIO7oSuUgScjfZ4knO50qlViq0P3+9pQth0PJqIUFXpcuWPrtWB7PPPxpAHgKgC49KgBeBfCLHV5X/yFB5qPgdqvAudywDofKUtLzuA/bSKOjdD2sr1PgNxpqTGi5/Ciu8to7AfzSvw3iW79vCv/b/7yB8ckqNn0X8YU/uYS/8dkq/vFPb+Hqk3vY3/5L7AYBx5tvoZRJolTIIukoYzUCFGuHzx+puRwo+Vwo+F0o+l2ouQ7X+pyGE26nG07DCYfhgMNwwOmgJWOAAsCE+chKqbvrSAWDyBRLiCRyCBRrqDdqKNdKwFoG5vY6VmN+FALtp/h6nB5EvVFEfVGM+kcx6htFLBDDmH8MMX8MsUAM44FxBN3B4SEZsULGxljPUyxSWPj9tEpGR5UFK3EM0e7DYVUzEo0enBx4XCFfKsXP0JWf2Vm+p64kif+9WGxfc9Zjg83V04OI40ihUy1MnE5rqxJJaPF6O3d+JAOtC2g728owjF8F8PdN0zyk5/jg4tTZVlKJDShzUc8X9/vbM60lj1u3UEol1R9LZoYkk2pTGgZ9mR4PVnfq+Jn/PY2/9fFdTJXvIGNUkasVkDWL+PoTASzlK0h9tYRrzxpwBwBHAygEufnqDgNbkwF4aiYCFRMx+DHW8CHoCSLgCSLoCsDvDsDv9iHgDsDn8sHn8sHjH4F7NAZ3NAbPaAxeTwBelxduh/tMwtg0TVTTSdSWl1DNZVCpV1F5aNnkR7xIT4RRMCvIVXLIV/PIVXLIlrPIVrLIlrNIl9NIl9NtW0A+lw8TgQlMBCYwGZzEVGiKP4NTiIfiCHvDg0kukpG1tXVwpr0QxdIS3UsuF91U8TiJR+/MelyQW4LxW1tW4SgZcsXi8f53HQ6HcqcIYQxSfEsSWaSflbiRdJLoVAaqkIPenUH/fQhiOHaq7mnJQ7Kd9Cyqs6K5H9bDFgb5fArJrQcorNzDfj2PvUYeqXIKW0YOd+/mEWw0cCVcQSRTQSBfRS7kxn7Uh0LAheTECNK3/RjPefCBRSDccMGYn0MgHIPz6RsYmZhHODKJYDgGw+/npj5JwZVoL2IKdyJjyjQZQ9KHFQG8oebmjp1CV6qVsF/cR6qUwn5pH8liEvvFfewV97BX2MNecQ+7hV2UakeH43wuH+LBOOKhOKZHpjEzMvPoiPlj/SeWcpnuqfv3uWe2t1UHZ0DtzcuXSRQS+wBUkaAupOp1pbyk06pBpghNt5vnvx2BL/2YdKVqENxPoqy1CkZ3ihx0YnC7W/8+7DEb2OTRnwrzZlQqqG1tILF8GzuZTWznd7Bb2MF2fgeJ/C52XWVkwl6UPQ7EdwrwluswTMBVb+DPv+bB990cwQuZKsLOAKJ1D0bqTnieuIZQZBJBXwjrzgX880/dxj//qYeuNr8feO651oEy8afKze7z8YaSNODjyMTtVkQyMnLq4DsAkufamjWwCyhXyxnMd9M0UagWsFvYxU5+59GxndvGVm4L2/lt5KuHdCwF4HV6MTMyg9mRWcyGZzE7Mou58BzmwnPwu3uccpxMAp//PDP1xI0k7WGeeorxDsOw9idbXFTXVQ6J3e3vk4j0+z8UYrZe8/XU94vsGb+//66nWu2gNV8qdSbeoFsMrQjiHBBDO+hInYeN9lGqlbCWWcNqehUbG28jtXIH6e1lJAoJNEwllE2HgWzIg9SUB87AOGZCD7XeZ+K4lHZgqurBeGAcf/GiF//4EztwRULU1oNBauZLS9RAg0FMpfaxm33oL63VqG0epgFKUL651bPLRYE0Pk4hJJXlzWnK1SrdbYkEPyMUUmRyUheFy6VGnS4vq9qadJo1DwsLp+74axgGA++eIC5EL7R8Tq6Sw2Z2E1u5LWxkN7CZ28RGdgMb2Q2ky2kspZawlFo68LrxwDjmRuYwH5nHfHge85F5LEQWEPZ2aZ6F1GxIUDuf57UqlXjeJFkC4PVdWGg9pbJSsTbwBHgNJybo9vJ6rUqFz9f/jDixlvR4y2HFce1A4ka6G0l3Kz0mxHBWDCV5GIYxBuDfAPgwgASAf2ya5m/1eh2lWgmr6VWspFcsx25uGyM5upg8FeWSMWAg5h9HbHQakbknEFu4htnoAmbDsxj1jR50kaytAevruB5+gMRKFfG5Kt0TPh9vqMuXH/VK2lspIhR+WJvicKjMmoUFq6viqBxyaZ2iV8KL1mmaqkuofnOZJq2VbJbr9XoVkYRC7d+I4TCrqPVso1qN7pqxMX6Ps1g4hyDkCeFK7AquxK4c+F++ksd6dh3rmXWsZdawnuXPjewGEoUEEoUEvrFtTTSMeCNYiCxgPjyPxegiFiILZyeVh63a1zcNfPU/lLH39Qrc5TIC3vsIP7uIb/m2HcQmHp6bkRHgmWcOkri4S5NJCs+REdUb6epV6+iAfkHcTUIQksnVXKzbDoQM9GC0HDY5dARD6bYyDOO3wfnrfxfATQB/COCvmKZ5q9Xzz+q2qtarWMusYTm9jOXUMlbSK1hOL2Mnv6NqGwC4K3VEshVE83VM+cYRD8UfBWanQlMYj1+CZ3qu/VGSpgl85Sv4jX/2Djxu4Ad/6GEH1GvXGCvI53lzffWr+IM/BFKeCfzIjz4cOSnFZLHYwcH3Ur2qa3LtZs/oPa+kBsHrbf19HA4KKb15WzvIZpn5owsNj4fE2IPpacehYTawndvGamYVq+lVrGaoQKxmVg+NsUS8ESxGFh8RymKEP4MeLROmWXiKAL1zB9/4ShGf/b/r+G9v3MbNC0mMhB3IVnz42r1RvPKKgY99DLh2FbQgQiFed2n7L20tAKuraWqKbqp+CNNGw0oQ8vMktQ+y95otJZsgOoZzFfMwDCMIYB/ADdM07zx87DcArJum+dOtXnNa8vjTe3+KL97+IjZyGxZ3k8BpODEbmsEV5wSu1KOYMyKIh+KIBWJwOR4adQ6H6u90UpfOw/YQ67ez+F9+Yh2f/rk53Pi2h4N9Ll1inOCNN3D3Kwn8+/9zCz/xk05MfOQF/r9eV03vrl5VBYVHQe/sKe2uj2uCKEVh0m5F/OGtbl69C+jIyNHZJq3amAP07Yt/f8BgmiYShcQBS3Qls3KAVBz1BryVBuLOCC5545h3j2PePY54cApToSl4XQ8TNHZ3sXUrgRdfBP7Rj27iYuxhS5p6nedwehpvp6fwf/18Dv/o5yOYuTmtXC+iXPh8VtINBGjJdSmF8wCkDkWvBzlJga3MptDjLVLzMID7oF8wTRP7pX0qt6llLKep6Lodbvzsd//sqd7zvMU8rgKoCXE8xDcBfKf+JMMwPgXgUwCwsLBwqg+qNWpYy67BgIHZkdlHGuNidBELvinMlD1wJQ8ZyuLzUQuMxU7nDtjbe+S+mX1yBB//P96Dz/xUDn/jo3V89KN1xCvvYCMfwX/+d2587Y/c+HufmsbEEw/nH8/PM4Om0WC8YHmZBWjHaWNi1usxBol56K0ldO1QgucAhUE+z1hIo2Fth20YfC+JlQAUADKbIBSykol0Bo5EuH5xtW1u0jLpVDv5DsIwDEwEJzARnMDzM8/zQdOEWSwisbeK9e172N5dwk5iGdupNWznt1Gtr+MB1vFAfx8YGPPHMO8aw+WsC69+eRzf9v4w5qZKwNgkY1Kl0qO2FNeen8P05jz+xVtj+LkPP2xkmEhYext5vXzd9et8XbeEbq1mbS1eKJxszsagZnANGNKl9CPlREhiJb2CbCV74LluhxsNswGH0VlrbBgtj+8A8Dumaca1x34cwN8yTfMDrV5zWstjv7iP/dI+5sJz8Dg9FIiplKrNOLg4CruJifZmJByGfJ5zF+TajI4Cly7h/hsFfP7n7+L/+aMaKvs5PDWyhqc/PIfv/29qmB0t8nNXVngDT09TiIyN8fHpaboozormITy53OE+aelpJJaJ36/y/lsJA6/X2rVY2i5UKnRj6efc6Wy/wV8voNfy6AV1R2jYDbOBvcLeo2D9Vm4La8VtrFX3UHKamNwromEAL30N+FgsB2PEiYgvirHQBMaC45hwRxELxDB68TrWLn4Y3/u+FN768x26qPQ9ahg8T7EYz7HsibOmXDcTRT5/sviEWBF6m41+Z3ANEMSSWMus0S0q8dXMCjLl1uV2QXfQouSKezTia6P1fgucN7fVcwC+bJpmQHvsfwDwAdM0v6/Va86cqpvL0RLY328dbJbGgBMTZ9/81Sq74Eomid/PfkZiNZRKLF68c0dVoz/xBF1TEidYW+PzQiEKkNFRZtI8/XR3CrYqFTVz4LB22jrplErWxm6HrckwrAWZUkipWyjj47S0euXjlgKz5kZ2J3HDAFyvjAnVv6PLhVqjhp3X/iu21+9gK7OB33rxNXzwIwns19Jw1GrYGQ+g4nEivl1A3QEU/S44Ryfxe1+YwT/5+yyCnAxOYdIbRdg9AkNa4zQjEFD9247bt/os8JO2ItFHoupE0e8A/YCg1qhhO7eNtczao2M1s4q1zNqhqeR+l/9RMoYci5FFjPnHOlqbdN7cVncAuAzDuGKa5t2Hj70LQMtg+alRq9EFtLd3uNmtWxmduGCmyewiuSmdzoOdTN1u1SNLXEPy2VevMo1zYYEEIk0V9/f5Mxi0NtbrFDweaxPJalVlYGWzPH86WQBce6FAQpDOwPJ/cUdJ6249VbhapRvLMHgutrZImpcv8zo09xFrB3oCgHQllTRlvYdRuXy6+gFpMaITxWFJBgBc+2nM1AOYid8EyqP4fzem8D0LAcQiJSRdFWxEndjbXkKlcAc71RT2ixmkcltwxTP48+XbMB0GMiMepCJeePwhzAWmcbkRwWI1iClXhBX3wQkECuC5XVuj63BsjHtaYhR6nKLdWeBCxG8ENgAAIABJREFU+PqogkGb3NcHNMwGEoUENrMqFXw9u/7I4qybrTMgQ54QU8HDTAWfC89hMbo4EAWsQ2d5AIBhGP8ebNT434HZVn+ETmdblcvU8JshvuNYrPMm9sqKmhUCMLNKd3/JvI9MRrVbn55WbqALF3iz3r3L9W9s8P0qFQoHv59t25vH03YblYqVTFq5NqSYrVBQWVwuV+s+SPW6in0IpJ2LEJj+uD4hUf88KajrFMRf31xQdxINu1xmfUujwWu3vIzf+EIQfqOMH/hrWVpZxSLfc24OWFlBrV7Dv/udNJZCPrzwEwbuOzNYLWy21Fz9xRrC2QqChSqC8GLaM4pxRxgTjhHEnCGMOgIIRyYwEpuFOxI9Oq7kcFjdTq06vD4mME0T+Wr+URHqdp6FqHJs57dRa7QmYAMGJgITmA3PYj48j7nwHGbDjLFGvJG+ksR5szwA4L8H8G/B2SJ7AP7eYcRxakjXz3yeN6r4jLuVKqqPiQWYZtkcN1ldVTUYbjfw/vfzNeIyWVpiIPnJJ0kgs7MUwHfvqpkOL70EfOQjvcu0ASiAYjHVcqRc5vcQMqnVrC4qgT7vwOm0WhRzc/xOMuLXNFkxnc+TUMW1JanFp4W8vlZTs1pcLp5/qTqWRAMZ5Sst7xsNEqV0Vm7ns5aW+LrtbbYOqdfxkRt7+He/ksHd567gyniRa5if57mKxXDvm2V8/v+7gF/5rTHMv//GI+FtmiYy5Qw2MuvY3HuAncQyEnurSGId6coGjP0ssoUkypUG1lssJ+QNIRCMwheJwT8Why8+i/DYNEbHZhEei2NsdAYeyQo7xxBiSBaTSBaT2CvsParz2S3sYje/i93C7pGNRwFgzD+G6dA04qE45sJzllY4HudgJX+0g6G0PE6KU8c8MhmV8tpNszuXYwyjKUBuQfMcaQmAV6t8rR5nuHiRxCOVyIkE8OqrJIyxMVpO3/7t7aXv9gKFgtUyOS7XX7oOywCcjQ3lWhEBPzWlXGTHQXdZScqxEI6QhRDGaf30Ml9aCtfkCATUFMn793nNVldZLCluyWIRt/an8Uu/OYYPfAB438fnMHVzBmv5Ufyb3wnjy//6LfzMP6vj2/+KyX0hbfz1auym+9w0TeQqOSQKCeylt5FNrCOf3EQunUCqlsOemUPFZaDidqLidqDuovVW8jpR9LtQ9DpR8rkQ8IbY2dg3iqgvypb5vsij9vkhTwgjnpFHLfaDnqBKY+8DTNNEuV5GoVpArpJ7dGTLWWTKGWTKGaTLaaRKKaRLaeyXmDRzmNWgw+fyYTLA5ptTQcaepkdIFvFQHD7XADWIbBPnKmB+GgxEb6vDcFyAHGBe/7176u+xMTUtDqDQu3NHpWYahiKQd96hQN7c5CTCUIjkcfEiYyQTE93/jieBxDmESI7rsyXIZNTcAikUm5hQlfZSIa8LVAl4n2Yoz1khMR9pjlkuK8upWOR+8Hi4xmAQWFjAZjaEX/vqNfzKH11Ceq+GmVgZH/9oCX/3r61ivrFCK6dVnOww6PUTerZTLof6fhKZxDpSheSjxpOZcgbp0kPBWk4jW82h4DJR9rpQ8jpR8jpR9RxPrm6HG/6H81v8Lj+8Li+8Ti/cTjc8Tg/cDjfcTjcchgMuhwtOwwnDMA6kmjbMBuqNOn9qA8qqjSoq9QrKNXZrLtVKKNVLKFY5p8ZsOdPuaPhcPrb697Pt/3hg/NExEWTX5pAn1Pc4RKdhk8egkodpMiVX+gw5nSQOPZhaKPA5IuCCQQr9ZuFQq/F5YoEYBq2XcJiukFSK7qv79ykk5uZIIJOThw/7GQToZCKpwYe5oXI5kqRYIACF4dSUigudprBMLA85xDXVHEvR4yh68F0ytPQgtB58l7bo0mm5XleusJERZtNJ/7DR0YPB60aDSoI8PjFxcNSszI04Sf1Eo8Fzmk7z/Betbhlx52QrOeSkbX6tgLSjin1HGSlHBSmjjH0Uka8WkK/mka/mWxbc9hIepwcBdwAhd+jAsLERzwiivugjK2rUT4tqGK2GTuA8xjyGB42GtdunnsGzssJsLiGG+Xnglha6kcwiQMUErl6lW0P87Xq3z6tXlQtL3CCXL5NEHjxgXUS5zDjB6ipbVxgGBdqlSwNXdAdAtYOXGI0E1mW2vPyUc+zxUAjLKF6AglWyiSR2IudOXElSqChZQnoH1dMQq96hOJdTxKVndgmZ7OyoQUwS4Je5EzJ+dmqKpNAq68nhIFlsbfG7lUoqPieEcRp3m8OhOgIAXI+Wkm0Ui4+EL0JHjLfVLBzT50PV60LRBRSMGsr1Msq1Mkq1Eir1CqqNKqr1KmqNGmqNGoeJNeqPhos1zAYMw3g0dEzGGzsMB9xONydUOmjBeF1eeJwe+F3+R3Nq/G5/X91m5wn2WewkGg1VhX1cda3M9BBMTVmD2LUaiUPcWdLmRIRRK4hWLIF3t5sZY9evMxPL4aB7xzAYB3njDeC97+V633pLuboGGc2xCcAq3J1OnkuPx0rMojlPT1OD1xvlSTBbF/RCMhL8bmVtyE+9nb1M96tW1RqlQ3GzlR8KcZ2rqySGep1rdLn4fLdbNZd0OqkMSAcA6WElbrpr16zWh9QedRJuNz9bCjP17sxytCI3qfEpFmEA8Dw8IvqoVX8UCPnaH7Vqo++wyeMskHYckjmUz7eX+pnPW4kjElEppjKTYW1NEYdk1xw3kEpu3GjUSjwPHtCyGBujwHnwgERVKDAO8vzzfO3duxS8s7P9d2Md1iSwVZqv06lmYeuvr9d5nvX24+JGisUOF1ISHzmqtkGq50VROC6jS69zkayyVEolMsja5udpNVWrvF4+H4lf4hFiDQkZyZAuw2B9z/37fK+dHboku1mx7XRaLRPgYDHhUV1xDxsNIAWUzfNDbFIZKNjkcVI0GiSL/X3+PK54SjKDJIgrgVJpMS7DfKTYzzSpQcrvjQbdTeGwtYBNHxnaPCHN7eb7Ly8rIbO0xJvV76dQun2bz00m+T0uXuQaEwlqwk88wc/s5uAbSWWVArzmyu2TNM6T1if66FMhwGRSkbH+mT4ftWjTVOdSj5c0r1UEXS53fHW19PWSGIMUyklm1d4e1+T3k/QkQ6xYpMUglfNPPMHHZL/pwf10mofbTUtkYoKfJe1gNjaOn13eaYi7Tx9CJhZV8zyOw+4duUea58jIKOjmgsvHvACxX7AD5u3ANFVrDGk4eBikg6kIMr1oql6n0JagtttN4tC1w6Ul1T4boACZnGxvnc2tMzIZfl6hwO/gdJJUfD7gtdcovABFYnpKsmGoYkhJU22u3hZ3jvjpm8+ZHjzW02FFSLdbtaxDr2AWominMK1eZ1Bar6WR7z45+ajJ4KO1S9+mvT1e82yW76EHuXXIjAw5gkFFFvq5ktG7uRwtQEBV6O/t8dzIWNmnnrLWvdTrXMve3uG91TweEoqQ5/XrgzVDXIc+GkD/eZJ90dz6RPaDTSgdgZ1tdVryKBZ5oyaTh2ua0lX2uDbjYlFIINcw6KfW4xxra0zXFHSimWGxqNJ4RZufneX3+fKX1ShScZE0ZyP5fKqKvdc4TYbQcSgUaFk1x40k6CzuvHS6dY8u/flCFOHw0U0GTZNCf3OT10EswXKZlk8gwH0mPb8uXgSeffbogtRyWXUobha2y8skq1iM+6e5ZmjQ0ZxSLb+3W+wpLkJR4kIhu+HiKWFnW50G+bxy7TRD3B3RaPuFaBsb1gygxUUrcWxuWoljfLwzXXD9frY5uXNH3UClEokrFmPMY2ODN+bMjHKXSdFduUwNORrlmtqplG4X0p9Kr9LWXXzdaJwXCPC7p1K0RAoFFSi/fZvnIRjk922ehOj3K0VBXE1HQfqKbW9bEycSCZ7jhQW+/8aGigtI2/TjOhl4vVQCZmb4GXpsZ3JSdSHe2eF7tWu9DgJ0K05HpXJwyuBhTTibYyniEZD3tcnkzLDJ4zBIuqbc1G43he3Y2MlbWSeTTKMUTE2pVh0AhcvGhvo7GqVg6RQCAfrO795VPv87d5jaOz9PIb21xcfn5nijRSJ0sUimkLicJNtGqrCbe0TJT92dpbtvdPdXv25gcZtJWm8iYRXuInjEpTUzw2Nk5Gh3iGRd5fMkp0zmoHurXOZ7XLrE8/fOO8oN6vcD7363NV5wHAxDNaXM59U+GxkheRQKbEnzzDMkm3aVnUGExFMiWntxiafoM2daZThKKrc+S0YfTNbvBJEhhE0eR2F8XAUwT7vBdN82wI0/O6v+3t1VsQeAm/nSpc5v5lCIfvR791Qm0717FIrFh/2SNjet7pwrV0hq2ayyOKRLbjze3aFCnUapRIGeSlk1UhEguRxJvly21nw4HBRM9+6p5IfmzsZCrKXS4fEwp5MB7WRSZc3JtQD4Od/yLWer+A8GeY2lHf/LLysrcmODZDY6yv13XObesKBVpl2tZp030yoLUiyY7W1VzxKN8v7spHV9jmHHPLqJSoX1E+KT9vlYQS7umERCFQECvAGuXOlusC+VYjqnXHdpgZ5K0a22u0vykvnQV6+qtt3NGp3bTQLptDurEzBNEp1kJB01zc7vp+CIRvk9JM51VMyjXUg/sViMVo5YmNJ3TNJPb94Ebtw4++fpuHOHXQckgC4tbQxDDQgbtOvWDUj2VjZLAj0upT4UUhb2Y+7esgPm/SCPep03rrR0cLlU6xGAAkq3SIJBEkcvBuQkkwzYCrxeleqbyfB3KTCTynVx8WxtHQzQiutEqpr7ZY1IhpnU3RyVFSfZZdHo4Vq4dP/NZFTCwXHweFSwdnRUkUOpxFbrpkkiXl9X9REXL7JRZadRrbIQtFikUiDfWeB00vKcmBgeC7ITaDQUkRynWIRC3NeiWDxmsMmj1+TRKrPq6lVlWjcTRyBA4ujl5my2emo1FZ/I5ZR7BuDPK1eoocucie3t1hloUoUciRwMOHcS+mRCvbX7YXA4uCY5TnOupdWMxIBkkqOkL0sBX6u1Sg+znR3+HoupdOgPfah7SsP6uoqD1OskrOaiPJ+PMbZB6bLca5RKvFdTqcO7N8iI6VhMFWU+BrDJo9fk0TzU6cIFFSBvFtr9IA6B3ubdNLm28XFr3YZo7y4Xg+6SIdZoqCyf5mIugWFYa17ERXaSG096POnzwSXT5rhuuF6vIot+BkW3t2lt7O3R+hgZUf20PvSh7gax63VaH0Ks8/O8Bq3ckGNjTJh4nF011SpJZH+/dS0NwHshFuO9Mqg1NB2CTR69JA8RFAK9VqN5Jkc/iUOwuan88JKNIu1JJicp8CS/Xtp9N2uo+bya8d5OgZdkzejFc4Aa6qR3pNUr54+Dy0WykiD4IASFZTJgIsE4R72uUrzf977epNDqe9LlYmzF4eDjm5tWEnY6ef0HrVV/P1Ctck8nkwetNcHICM+VNBk9Z7DJo1fkkUrRXSXQ525sbdGFIBgE4hDoro2tLQqTmRlqoBcvMsAupCBpprrvXKAHqTOZzgScj4LHo4rAZMzuoOHOHVqaKyskkJkZXvvnnutd8V6jwW7NknauKzQyzljvagDwvC4uDuY57QdKJSpI0gWgGdIiZnz8XFluNnn0gjxyOVVHAVgzp3ThDPQ2ON4uVldpGdXrJItQiEJGsqnu3rU2alxctNaqtEKlwvOiu5qO6wnVCi6XKiDUexoN+k26u6uKMDc2eC6jUca/nn66ty009Dibw0HrQz9/2SwJTid8w2BK9vT0udSqTwVJKkkk1GhnHZI8MjV1LojXrjDvNkola7GX9CYyDGqdUpgEUEO+fHmwiAOgL7zR4FonJynsRGDEYqzKvntXtTh58IDWyNTU4e8pPZt0SBPCSkUFnfXW6tIOXZ+7MYx9iioV9g/b2qJ7SKrWL1xoPcyr24jFuA6JFTU3TZQGnVtb1rnwm5sUkhcuDHeBYacggfNIhNdYWsSIUiRNL/f26DqNx89tIoJNHmdFtUpftrh1XC5lVSwt0WcqiERUDcUgYnFRxRfSabX2QIBpxkIgkn4snWpl1Gs7kEK7QYhHdBNf+xoFbzKp3IDz89wb/fruc3O8fgCFm0xYFDgcXOfYGJUeyTwqFlmvFI+rFjY2qNjMzNAya24RA6g072CQ5+4knQOGAAMqxYYE9TqJQ7JWHA5mJLlcvEl14hgba3+2dD+xuMi1xuMUEvv7JMHdXbo5rl2zVvNub9PN1Y854IOK116jZSaxn4UFHtJav1/QZ2/I6NtW8PlU6xp9v25tkUSaxtE+9hBX1ZNP8pBhWYJ8np6JN9+kMnFOQgUDLskGGDI3Wk9TvXSJAvb2bWua3+QkA8/DoLEZBl0UU1OqUDCZ5NCiSoUW1ZUrVi1qf59keZoW6+cJpskA+RtvcF8kk9RKr1xhMPUoF1+voLfGkVbzrSDZdtevW90uYoWIa8uGFcEg5cCNGweLL4tFKmLnhERs8jgNpKW2fuMtLiri0IOOc3PU4IYJ0hb8iSdUhXQiAXz1q/xdMq70NFPpStvtDKtBhVihb75JIbGzQzfls89S+PZ6KNNhCASsiQ56WnkrSIsa3QoRq+XOnfaq7h9HeL20Np95hla8bsGVSopEdO/EkMEmj9NgZcWaaTE7S4387betGUmXLg2GtnkayPqfeko9trwMvP66+v/8vDXeUS6TQDKZ3q6136hWKUhlX2xv07X3rnepjsaD5K6cmVHrEQvpOExOci/oYwRyOQrAvb3urPM8wO2mfHjmGVqheqJMqUSX75tvWsc1DAkGaEcPCdbWrNlTU1N0Yel+f6eT2lqz73PYYBg0vy9cUI/dumVtrDg1ZU0CEA1cn9F+nlEo0I2TTjODaWuLcYUrV2h5SAxskODxWK3G9fX2YlY+H2NeetC80WB85/799oc1PY5wuXjeWpFIsch7ptndPeCwyeMk2NiwDmyKRqltb26qx3w+amjHDfMZFhgG8G3fplJuazUSyNKSIpDRUQoVqRswTdaNLC2d70D6/r6yNtfXuQ8iERXfuHx5cNtXxOOK1CqV9sneMCj8rl2zfrf9fWrQh/WGskFIM8obNw66s/J5WrD37g1FUoJNHu1CWjkIAgFeYN19FQ5bO+eeF7hcwHvfq6rK02mSg17bEggcdGskk9Smzptf3DSpSIi1mUxS+x4dZXxjZoYWxyDn94sQE2xtnayAMxjk9ZakCoAk9PbbPDdDHgzuOlwuurNu3Dg4FyedJhE/eKC6AgwgbPJoB839qgD6K3WhODlJgTFoxX+dwugoq6LFFbe1RaGpZ1m53XTX6QJFsnPa8asPA6pVfmdRJEolfr/RUVYUj41RqA5DTr/e2K9et06zbAcOBxMBLl+2uuY2N6lBD7DgGxi43Ywd3rhxsGPD3h4z99bXB9IlaJPHcdjZUcQhvZv0qmiHg5lJ8/PDkYp7FiwsMEA+Pk4hurNDN8XbbytBIQJlcVGdj3qdLqylpYG8CdpGNkuiEL+0zGyJRFSzx+ees5LnIEOSHgSJxOEdko9CNHowpVeC6a1aeNg4CI+HscXr161940yTitobb/B+GyCLziaPo7C7q7rgVqv8OxhUfkqPh77f5hYc5xUuFwlEfPoy/6BUontK99OOjx904SWTFCjDlo3VaHAf3LmjXDuNBm9qGUlrGMDzzzMeMEwIh63C6rjU3cMgM190V5jMaF9dHSihN9Dw++nBuHrV6gKu1Xgeb90aGEK2yeMwpFJMvQTog9zaolkpxBGN0j3xuPX7GR3ldx8bo892a4tCololgegphxIH0c3xSoVun+Xl4bBC8nlaG3pA2TDossxmlXX19NN03wwj5ubU98hmTy+c9GC61AcBPHfnMfbVTYyMUPm6dMmqgJXLJOQ7d05nJXYQQ0UehmF8yTCMkmEYuYfH2137sHCYzL+xwUySmRmlYc7NHfTzPk5YWOB3D4cpLGTwVaPBTBE9I83ppDl+8aL1fCUS1KIGtUagViPBNRc+BgLcA3pQeGGBNR3DCp/POr9jbe1sWXKh0MG4j6Q0n5fYV68gsca5OWs8VVyoy8un61TdAQwVeTzEp03TDD08rnXtU0ol1Y5jfp4/vV5qVcNa+NcpuN2qYlq6xeouq7U1ZoroAmhsjP5cvfalWuXzbt8enBRP0yQZ3rplredxOHjdazUGMCXGE48D73nP8Me79NqDcvnsdTouFxUsPRYosa/l5fOdwt1pGAb3XqvMrESC8ZDmoV49wDCSR2/gdHKzT03x91jsYCrq44xoVLmjfD4euqtib4/BZN1V4XbTDJceYIJ8ns+9e/fwiW29QDJJ0lhZsfbpikQoBKX9diqlXDTPPXc+UrMldVSwudmZbKnJyYOxr0Ti8W5lc1q4XNyHzUF1abF/61ZPLbthJI+fNQwjYRjGlw3D+EDXPkV60ziddLlcuHB+03BPC5mHDVCYNs/vEFdFsw99dJRaVPOQoUyGQuXuXb6mF0HWRoNa9htvUCvWyc7jUWN3l5cp7DY3levy0qXjB2INE8bH1QCjRuPwrrsnhcS+9L0hKdyD6rYcZPh8DKpfuWIdOFWpcA/fvt0TJWyoJgkahvEeAG8CqAD4OIBfBnDTNM13Wjz3UwA+BQALCwvPLy8vn+5Da7XHN7bRDiRVVzAzw/PVnGEzMUGB29zjqVym1tRKY/J4KNBGRztfqZ3L8TOTyYOBe5eL7qixMX6P/X1+l+Vl7oe5OcZ7rl8/fwpFNstgrKC5Bf9ZkUjwnOoulliMitog9f8aFpgmz+nGxsGu1pLUonsEToihGENrGMaXAHznIf/+smma72vxmj8B8Iemaf7SUe/dsxnmjys2NqzV9zIp7/59q+vD66UV18r1J1r9YWa3309TfWSErz+p0K7VKBhlxnorl4zTSTfL1JRqWifP29riY3Nzaq7JeXVh3r+vur36/bQaOhnTKRb5GbrbyuejJXcORrf2BfU675/mWhCJ1TW3QmkTQ0Eep4FhGH8M4I9N0/wXRz3PJo8uwzTpapLiObebWjnA+EFz2+nJSZW91oxKhQHrROLo+SA+H8nI5+PnORw8DIM3kqQPl0oq+eGo95qcpAZsGLwJ9XkV6TRdcDIga27ufCdNVCr0n4t10I3v22hwb+huK4eDn6Vnftk4GcplJqw0u4pHRqjUnRDngjwMw4gCeA+APwdQA/CDAD4L4DnTNO8c9VqbPHqAapUFgCLww2H6ZQ2DAmJlxeqq8HgoKA7rPNxo8AbY32cspNOZJP9/e+ca29Z5n/HnT+ou2tTVcmTVkp0LnARIgsEounZFiuZDOxQbug4DhnRod0PWAcGKbQg2FMOSbkHXZh82LCvaBVgToLtiWDIgGDrsw27dimbzgLVN2thxIjtxZMmyZd0vpKh3Hx6+PSRFSTwUzzmk+PyAA1jksfi+pHie8793dPC1fT8qgK4sH9so3dfWVuC6GRzkHfJRp7QlTyrFdNFDuD/2pNrfxuAgs/mOmkswTlZW+Pn5WpC77ioPstfIfuLRSs78TgBPAzgHoADgdQAfP0g4REx0dtIl5WdkLy8z4DoxwTv6TIYXCV9dnsvRdZHJ0C9b6VdPpeizHRrihWVpiRf31dX6iqPM6GY6doxHJhO4YnI5ut4qg7c9PVyHzwzr6WmeoU5Rc+IErb/NTb7/165FI5rDw/xc3norSPe+fZsB37Nnj65rMGp8keHCAoWkDuE4iJaxPA6DLI8YqYx/nDlTnmWzsMCAaaVLKpulW6iW4GyhQGvAH/l8eb+xdDo4fBpxd/duv30+T/fU/Hy5n9h3nL19O6g/SaX4ZWwnn3xlMsTdd0c3g90LlC84Bfh5jY/TZdbqdTQtylGxPEQrMD5Oy8C3Kbl6lRdv38ZlaIgXoNnZ8uDe0hKPTIZ3vQMDe18w0mn+vnpbw6yv87WrzZEeGGAK8vXr5YWLU1PtJRwAP4vh4cAie/ttxrKiyIpKpZhx5dOiC4Vg3O3yMt//KNxmom4kHqLxnDkTFIH5liXnzgVf/o6OIDBamabrXVPpNIVmcLDcxVQvuRwtCe8SqaS/n2vKZCgspdXl4+OtPxWyXiYmGHvy1t7MTPno4UYzOMjPYno6EO+VFcbTJifb93NoQuS2EtHgO+36GoreXqa3VguCbm4yQHvrVvXCwFQqiFX09vIorVCvZGeHF7qNjUCM9prM1t9Pd5nvw7SwwAuXZ2iIYtjO3LxJa8ATR0NQ52j9lbpAAVpCvl2QiJwjkW11GCQeCbG8TKvD/41lMvSb7+X28G3vb906uDVGKkULJp3mv52jaGxv75/iC9CKGRyke6w0ILu0xI6lfr39/UHNSrtz6VKQit3XR0syjjjE6irFvLJeaGrq6Ix6bmIkHhKP5Ki8kx8YYBbNQRceXwG+V0FfWMxoufiW8pVdA1ZXmSnmg+69vRQOdRcgW1t0Hfn359QpWmxxUCgw3lJZQDo2xnUomB4ZCpiL5BgaokXhawYWF5mWeZCAZDLBneXWFq2YjY3g2G8WiO+z1d0d/J7SIV6VrK7SQvIXxq4uWkgSjoDubvYi8/2uZmYowo1uG1MN318um6WI+M9+bo43F1NTSulNAH07RPSMjVFA/JyPWgXE0929u+rYV5Fvb9PNZMYjnWY8pNa70eVluqq8cHR00OLYL6bSroyNMeFgfZ3v+ZUrjGPFdec/NMQbgatXg3ohH1s7eZLiJhdjbOidFvFQ2eJicZEX7XqnCabTtBD6+njX2ddHV1NXV+0Xs8XFcoujs5PCcRRarEeBGe/y/fu7tlY++CsOvFVY2UhxdpZdeptlLkwbIPEQ8TExUe4nX1piEVojYhphmZ0tD477efTtVssRlt7e8jnlMzN7Z7JFyegoa058axmAVsjFi+WuLREZEg8RL6dO0b3g2diIbf4AAFoZ09Plsyr8hEhZHLUxNhak6nr3VRKJN93dtBT93B2PnwSpkbeRIvEQ8TM+Xu7+yOd5x1jayTYK1tcpVKUXlUxVE9/2AAAOWUlEQVSmvIBRHIwZA9j+81tfb9zgqHrwVkhp/6Z8njcJly4lYxm1ARIPkQzDw+UZTb4VReXo2kbgHN0rr79efiEZGVE6br309JRXms/NBUHsJOjqYufYyhHHvjq9crSwODQSD5EcvvNnaZrl2hpdDtWaJ9bD7du8eFy/Hlg1vo/S5KRqBA7DiRPljRKvXEn+Aj04yPbxlc0U5+c5anhurvHt/dsUiYdIFh9vGB8PvuzOsb/U977H+pDS+Rq1UCjQNfX97++eWJfJsL2GBg41hqmpwHLL5ykgSZNO0yq6995ycSsU+Pf06qtsudIGBdJRIntdJI8Zg+i+CMwHz3d2eKc4N8cA7cBAeUqup1CgQGxsMIOr2vCodJqvceKErI1G0tlJAbl8mT8vLTF2FVf1+X709tI1urRUfhOSz7NWZHaWfxNDQ/qbqAOJh2gefM+kpSXGP0rjE+vru4dAlVoqe5FKBXPJFduIhmyW76+v+Xj33WDwVjOQzdICuXmT7st8no9vbdFSmpmh2A0Pq8gwBPo2ieYjm+WxuEj30+JidYHYTzR6e3lHOTIi0YiDU6doMfoivelpuo2apVLfjK7K4eHAmvW1ILkcLd6ZGZ4zOto8624EvgNDg9G3SjQvAwM8CgUKiLc+KntbmTH7x7drj6vnkggwY6aTn2OfzzPedM89zeUSSqUC9+X8PEXEB/m3t2mZzM7yb2hkJLrJiVGyvU0R92K+vh7JFEyJh2h+0mneMQ4Plz9eank00wWqXenspIBcusSfV1cZW5iaSnRZVUmn6arys9rn5oJOB84Fg8O6umjBDg01b/eBnR3G+VZWeFSra1ldlXgI8UMkGM3HsWPMdPJdlG/d4kWrtK9ZM+FjYqOjtG5v3Cjvj5XL0RKZneU+sllaJUl38c3luN6lJQpG6Y2Uc8EwNH9ksw3PMJR4CCEay9gYL1h+9vm1a3QjllaANxt+QNjgINd+8ybXX+oe9Rfi2VnG0Y4do1srk4nHTZrPMwZYOkrZC8XmZnBsbTEFvreXazxxIpLxvRIPIUTjmZzkRczfxfv4R9J37LXQ28tRtxMTvLP3Q8lK07+3twPXFkA3WF9fkEre08ML+GGTNQoFvsbsLNeRz9PqyOX4/ubzFJDubr63IyNcQ2XW2NrabrfvIZF4CCEajxlw551sk57L8cJ7+TILQlslmcEsSNrwcQXvKqqspC8UgphDKakU4yadncHY5HSav9tf4J0LjkIheC0/ktnPrKlcW19fMOOkMjusszMYguZHFjQYiYcQIho6Olikd/FiMFv+jTcoIK3WiDKVCoQEoPvKB6nX1vZuy7KzE7iTDmJnh8J0+/be/d36+uj+y2QCq8a7qLzl09cXS6qxxEMIER09PWxYeOkSL465HAWk1ac1+rRwnwiQy1FENjbKYw+19NHK5aq7xjxdXQx2j4zQiujpCdxivb2JFTZKPIQQ0dLfTxfW5ct0v2xuUkxaXUBK6eriURmYLhQoDvl8MDZ5Z4fHykqQ3dXRwZhEOk0x6OigYJw8SWunCTMLJR5CiOg5fpwzQKanAwG5eJEC0mourDCk04GVAgQ1JHNzLN4rfQ6gReEr4UsHXDUhEg8hRDwMDvIO+q23ghTTixcZF2mVIHq9FApM/71xo/rY5Wx2d4v7JkfiIYSIj4EBurD8/PhcjkO67ryzeRopNpJ8noIxP797rnoqRQvjxImWFE+JhxAiXrJZBtHffJO+/0KBQfTJyYbXIiTG5iZdU7du7U6z7egIqtpbuGln665cCNG6HD/OlN3Ll4NCtytXmLE0MdG6rdFXVykai4u7n+vuZgB8aKh191eCxEMIkQx+fsvly0Ezv/l5CsjZs7zYtgLOUSzm5oK2IaX09zOlt0mzpuqlqeTPzB43swtmtmVmL1R5/hEze93M1s3sX81sMoFlCiEaRVcXLZDSFNf1dbZ2v3GjuUfFFgoUjFdfZRJApXBks8wmO3cuSBY4QjSb5TED4GkAHwFQ1j/YzEYAvAjglwG8DOD3AfwtgPfFvEYhRCNJp2lpzM8D77xDwdjZ4b8XFhgLaaZ26OvrXOvCwu6iPrMgCN5Ma46AphIP59yLAGBm5wFMVDz9CQCvOef+rnjOUwBumtk559zrsS5UCNF4Rkfp4pmeDtp5rK2xP9bwMIc4JVUTUihQLG7e3D0OGQiK+o7aFMJ9aCrxOID7AXzH/+CcWzOzN4uP7xIPM3sMwGMAcPr06bjWKIQ4DH19wH33sYvs9etBw0DfIn1kJL7UVt9ryrcOqeZC6+ujYByRIHgYWkk8MgDmKx5bAlA1Odw59xyA5wDg/PnzTew4FUKUYUYrY3CQrqvlZT7uHN1F8/OsCRkZYVyhkZXY+TyFYmmJr1ut15Sf/TE6ygaFbUps4mFm/wbg4T2e/i/n3I8d8CtWAVSWXx4HsFLlXCFEq9PTw+rzlRVgZqZ8wp9vf24WDGXq72ecoVYx8R1v19eDmd/7db/t66NgDQ01feuQOIhNPJxzHzrkr3gNwKf9D2bWD+DO4uNCiKPKsWPMyPIzLkprKJzj4946AZji29kZzNAoPXd7m9aFH6p0ED09FIvBwZasAo+SpnJbmVkHuKY0gLSZ9QDYds5tA3gJwB+a2U8D+EcAvwvguwqWC9EmHD/OI5dj/GNxsXrwemtr73kYB2FGV1Q2y0OCsSdNJR4AfgfAkyU//xyAzwN4yjk3XxSOPwXwFwBeAfCz8S9RCJEoXV2MidxxB4VkeZlup/V1up1qrQ0xo5XS00OXVyZTfYSrqIq5Zi7CaRDnz593Fy5cSHoZQoio8QOnvGuqshlhR0fgzuruPnKFe43GzP7XOXe+2nPNZnkIIUT9pFLBpD0RKbLPhBBChEbiIYQQIjQSDyGEEKGReAghhAiNxEMIIURoJB5CCCFCI/EQQggRmrYoEjSzeQBXD/ErRgDcbNByWoF22y+gPbcL2nM4Jp1zo9WeaAvxOCxmdmGvKsujSLvtF9Ce2wXtuXHIbSWEECI0Eg8hhBChkXjUxnNJLyBm2m2/gPbcLmjPDUIxDyGEEKGR5SGEECI0Eg8hhBChkXgIIYQIjcQDgJkNmdlLZrZmZlfN7NE9zjMz+5KZ3SoeXzJrzVFkIfb8hJm9amYrZjZtZk/EvdZGUeueS87vMrMfmNm1uNbYSMLs18x+xMz+w8xWzWzOzD4b51obRYi/624z+2pxrwtm9rKZnYp7vY3AzB43swtmtmVmLxxw7q+b2ayZLZvZ18ysu97XlXiQLwPIARgD8EkAXzGz+6uc9xiAjwN4EMADAH4CwK/EtcgGU+ueDcCnAAwC+CiAx82sVWfH17pnzxMA5uNYWETUtF8zGwHwTwD+DMAwgLsA/HOM62wktX7GnwXwo+D3eBzAbQDPxrXIBjMD4GkAX9vvJDP7CIDfBvAIgEkAZwF8vu5Xdc619QGgH/xju6fksa8D+GKVc78F4LGSn38JwLeT3kOUe67yf/8EwLNJ7yHqPQM4A+AHAH4cwLWk1x/lfgF8AcDXk15zzHv+CoBnSn7+GICLSe/hkPt/GsAL+zz/VwC+UPLzIwBm6309WR7APQC2nXOXSh77DoBqdyv3F5876LxmJ8yef0jRRfdBAK9FuLaoCLvnZwF8DsBG1AuLiDD7fR+ABTP7lpndKLpwTseyysYSZs9/DuADZjZuZn2glfKNGNaYJNWuX2NmNlzPL5N4ABkAyxWPLQE4tse5SxXnZVow7hFmz6U8Bf7NPB/BmqKm5j2b2U8BSDvnXopjYRER5jOeAPBp0JVzGsA0gL+OdHXREGbPbwB4B8C7xf9zL4Dfi3R1yVPt+gUc/L2visQDWAVwvOKx4wBWajj3OIBVV7QBW4gwewbAoBwY+/iYc24rwrVFRU17NrN+AM8A+LWY1hUVYT7jDQAvOef+xzm3CfrB329m2YjX2GjC7PnLALrBGE8/gBdx9C2PatcvYJ/v/X5IPIBLADrM7O6Sxx5EddfMa8XnDjqv2QmzZ5jZL6IYaHPOtWTmEWrf890ApgB808xmwYvKHcUMlakY1tkownzG3wVQegPUajdDnjB7fgiMDywUb4aeBfDeYvLAUaXa9WvOOXerrt+WdJCnGQ4AfwOa6f0APgCac/dXOe8zYBD1FJih8RqAzyS9/oj3/EkAswDuTXrNcewZQAeAkyXHJ8BslpOgKyvxfUTwGX8YzDZ6CEAngD8C8M2k1x/xnp8H8PcAssU9fw7Au0mvv849dwDoAfAHYIJAD4COKud9tPhdvg/AAIB/QQ1JMnu+btIbb4YDwBCAfwCwBuBtAI8WH/8g6Jby5xno0lgoHs+g2B+s1Y4Qe54GkAdNXn98Nen1R7nniv/zIbRgtlXY/QL4VdD/fxvAywDek/T6o9wz6K76SwA3ACwC+E8A7016/XXu+SnQWiw9ngLjV6sATpec+xsA5sA4z/MAuut9XTVGFEIIERrFPIQQQoRG4iGEECI0Eg8hhBChkXgIIYQIjcRDCCFEaCQeQgghQiPxEEIIERqJhxBCiNBIPIQQQoRG4iFEzJjZqJldN7MnSx57wMw2zexnklybELWi9iRCJEBxJOjLAB4G8H8ALgD4b+fcLyS6MCFqROIhREKY2R8D+EkA/w427nvIObea7KqEqA2JhxAJYWbd4CjQuwG83zn3SsJLEqJmFPMQIjmmALwHbKF9NtmlCBEOWR5CJICZdQL4Njj97hUATwJ40Dn3dqILE6JGJB5CJICZfRHAowAeAKfdfQOcAPdh59xOkmsTohbkthIiZszsYQC/CeBTzrlFxzu4nwfHg/5WkmsTolZkeQghhAiNLA8hhBChkXgIIYQIjcRDCCFEaCQeQgghQiPxEEIIERqJhxBCiNBIPIQQQoRG4iGEECI0/w+/8qxbq7uB1gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# También podemos obtener muestras de la multinormal\n",
    "import scipy.stats as st\n",
    "X = np.random.multivariate_normal(mean.flatten(), sigma2, 10)\n",
    "\n",
    "plot_data_sine(x_[ii], t_[ii])\n",
    "\n",
    "for pp in X:\n",
    "    plt.plot(xx.flatten(), pp, color='r', alpha=0.2, lw=3)\n",
    "    \n",
    "ax = plt.gca()\n",
    "# ax.set_ylim(-5, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "**Ejercicio**. Exploren el efecto de cambiar las funciones de base, sus hiperparámetros, la cantidad de puntos usados para el entrenamiento, su ubicación en el eje $x$, la constante de regularización (es decir, el prior de los parámetros), etc.\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  },
  "nav_menu": {
   "height": "279px",
   "width": "309px"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Índice",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "239px"
   },
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}