04.00_intro_estadística_frecuentista.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from fast_kde import fast_kde\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import interact\n",
    "import seaborn as sns\n",
    "plt.style.use('seaborn-colorblind')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estadística Inferencial\n",
    "\n",
    "En muchos casos describir los datos y graficarlos puede ser suficiente para nuestros propósitos. Pero en otros casos esto es solo el primer paso de un análisis que intenta no solo describir los datos si no comprender los procesos o mecanismos que dieron origen a los datos observados. Esto se conoce como estadística inferencial, es una forma de modelado y será el tema central de los próximos capítulos. \n",
    "\n",
    "Los modelos son una descripción, simplificada, de un sistema bajo estudio y no solo una descripción de los datos. Los modelos se construyen intentando capturar lo mejor posible la información relevante, por lo que un modelo más complejo no necesariamente es mejor que uno simple. Decenas, cientos o miles de datos pueden pasar a ser descriptos por unos pocos parámetros que definen al modelo (en cierto modo estamos comprimiendo información).\n",
    "\n",
    "* Los modelos, al ser más simple que el fenómeno estudiado, permiten pensar el problema con mayor claridad.\n",
    "* Crear un modelo implica pensar acerca del problema y evaluar cuales son los factores que se consideran más relevantes y cuales pueden ser despreciados. \n",
    "* Los modelos estadísticos (como otros modelos formales) pueden ser estudiados analíticamente y/o numéricamente lo que puede contribuir enormemente a comprender el fenómeno subyacente a los datos.\n",
    "\n",
    "Dentro de la estadística Inferencial se suelen distinguir dos grande paradigmas:\n",
    "\n",
    "* El paradigma Frecuentista (tema de este capítulo)\n",
    "* El paradigma Bayesiano\n",
    "\n",
    "En los últimos años se ha desarrollado, además, un conjunto de técnicas agrupadas bajo el rotulo de _machine learning_. El _machine learning_ es algo así como una re-invención de la estadística por parte de informáticos más preocupados en resolver problemas que en demostrar teoremas.\n",
    "\n",
    "Una de las diferencias, algo caricaturizada, entre la estadística inferencial y el _machine learning_ es el énfasis en el modelo. En _machine learning_ se usan modelos que en general son de caja negra es decir los parámetros del modelo no son de interés principal, lo importante es predecir. En cambio la estadística inferencial pone énfasis en entender los parámetros del modelo como medio para entender el problema y eventualmente predecir el comportamiento del sistema. Esta linea que separa ambas aproximaciones se está tornando cada vez más difusa, el flujo de ideas entre estas dos disciplinas (que por mucho tiempo se desarrollaron en paralelo), es cada vez mayor. Los métodos de _Machine Leanring_ usan muchas ideas de la estadística. Además es cada vez más grande el interés por parte de la comunidad de _Machine Learning_ en desarrollar métodos Bayesianos. En el fondo estas tres grandes grandes ramas de la inferencia tienen un mismo sustento teórico que es la teoría de probabilidad."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Métodos frecuentistas\n",
    "\n",
    "Al analizar datos, se puede pensar en dos grandes alternativas:\n",
    "\n",
    "* Pruebas de hipótesis\n",
    "* Estimación del tamaño de un efecto junto con la determinación de la incertidumbre\n",
    "\n",
    "Ambas aproximaciones pueden encararse desde una posición frecuentista o Bayesiana. En este capítulo veremos la forma frecuentista.\n",
    "\n",
    "Dentro del paradigma frecuentista la primera de estas estrategias es la que tradicionalmente se ha enseñado y divulgado con más énfasis. Mientras que la segunda viene ganando adeptos y se la ha llamado \"[nueva estadística](https://www.amazon.com/Understanding-New-Statistics-Meta-Analysis-Multivariate/dp/041587968X)\" (aunque las ideas disten bastante de ser realmente nuevas).\n",
    "\n",
    "Empecemos viendo las pruebas de hipótesis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prueba de hipótesis\n",
    "\n",
    "La prueba de hipótesis,también llamada test de hipótesis o prueba de significación, es un procedimiento para establecer si una propiedad dada de una población es compatible con lo observado en una muestra de dicha población. \n",
    "\n",
    "En la práctica, como la población es desconocida, lo que se hace es pensar una _hipótesis_ la cual nos permite generar una **distribución de muestreo**, la cual nos permitirá calcular probabilidades, ya sea de forma analítica o realizando simulaciones. Luego se compara lo observado con lo calculado a partir de la distribución de muestreo. Por último se interpreta esa diferencia en términos de la _significancia estadística_. \n",
    "\n",
    "Si alguna vez calcularon algo en estadística que luego buscaban en una tabla llena de números, entonces ustedes hicieron la versión analítica. Nosotros veremos ambas formas.\n",
    "\n",
    "\n",
    "De forma general, resumida y algo más formal, tenemos:\n",
    "\n",
    "1. Calculamos un _estadístico_ $T(\\mathcal{D})$ es decir una medida cuantitativa que describe nuestros datos $\\mathcal{D}$\n",
    "\n",
    "2. Generaremos una hipótesis nula $H_{0}$, es decir la hipótesis del no-efecto (la no-diferencia) y asumiendo $H_{0}$ calculamos $T(\\mathcal{\\tilde D})$\n",
    "\n",
    "3. Calculamos $valor p \\triangleq p(T(\\mathcal{\\tilde D}) \\ge T(\\mathcal{D}) ) | H_{0})$. Es decir la probabilidad de obtener un estadístico al menos tan extremo como el observado, asumiendo  $H_{0}$ como cierto.\n",
    "\n",
    "4. Evaluamos si el _valor p_ es _pequeño o grande_. Para esto usamos un valor predefinido, el cual suele ser $\\alpha = 0.05$. Si $p < \\alpha$ entonces se dice que \"rechazamos $H_0$\" si en cambio $p \\ge \\alpha$ se dice que \"fallamos en rechazar $H_0$\".\n",
    "\n",
    "La lógica de este procedimiento es similar a la prueba por contradicción usada en lógica y matemática. La cual consiste en algo como lo que sigue. Para probar X, asumimos por un momento que X es falso, Si ese supuesto nos guía a un resultado contradictorio o sin sentido entonces concluimos que X es cierto.\n",
    "\n",
    "De forma parecida en una prueba de hipótesis, para evaluar si el efecto que observamos/medios existe, asumimos una hipótesis nula que dice que no lo es, y basados en ese supuesto calculamos la probabilidad de observar nuestros datos. Ese es el _valor p_. Si el _valor p_ es _bajo_ (según un criterio predefinido) concluimos que los datos no parecen estar de acuerdo con la hipótesis nula. Ojo que NO podemos establecer que la hipótesis nula sea cierta o falsa."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### El problema de la moneda y la prueba de hipótesis\n",
    "\n",
    "Para hacer esto  más concreto veamos un ejemplo sencillo, que nos acompañará por un buen rato. Supongamos que arrojamos al aire 200 veces una misma moneda y observamos 90 caras y 110 cecas.\n",
    "\n",
    "¿Es compatible la idea de una moneda que cae la mitad de las veces cara y la mitad ceca con estos datos?\n",
    "\n",
    "En el ejemplo de la moneda un estadístico podría ser la cantidad de caras $z$ para $n$ tiradas. La hipótesis nula sería que tenemos el mismo número de caras que de cecas, $H_{0} = 0.5$. \n",
    "\n",
    "Un par de cosas a notar son que:\n",
    "\n",
    "* Bajo el paradigma frecuentista no es posible aceptar una hipótesis nula. Por eso los frecuentistas dicen \"fallamos en rechazar $H_{0}$\". \n",
    "* El computo del *valor p* asume que $H_{0}$ es cierta. Por lo tanto el *valor p* NO es la probabilidad que $H_{0}$ sea cierta, tampoco es la probabilidad de $no H_{0}$ (la alternativa) sea cierta, ni siquiera es $p(\\mathcal{D} | H_{0})$.\n",
    "* Bajo el paradigma frecuentista no es posible preguntarse sobre la probabilidad de un parámetro tome tal o cual valor. Los parámetros tienen valores fijos (aunque desconocidos). Por lo tanto NO podemos averiguar $p(\\theta = 0.5 | \\mathcal{D})$, es decir ¿Cual es la probabilidad, dado los  datos, que nuestra moneda NO esté sesgada?\n",
    "\n",
    "Ahora que conocemos el modo de pensar frecuentista veamos como aplicarlo para resolver el problema de la moneda. Primero lo haremos mediante una simulación y luego analíticamente.\n",
    "\n",
    "\n",
    "A continuación vamos a generar los datos con Python, _como si_ hubiéramos arrojado la moneda 200 veces y hubiéramos observado 90 caras. Vamos a usar el número 1 para representar caras y el 0 para representar cecas. Manipular números es mucho más simple que manipular _strings_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "monedas_obs = [1] * 90 + [0] * 110\n",
    "#monedas_obs = [1] * 32 + [0] * 44\n",
    "# una forma alternativa a la linea de arriba, pero usando NumPy\n",
    "# monedas_obs = np.repeat((1, 0), (90, 110)) \n",
    "n = len(monedas_obs)\n",
    "caras_obs = sum(monedas_obs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Valores p por simulación\n",
    "\n",
    "Supongamos que son un poco escépticos sobre el uso de la distribución binomial, para resolver este problema. Una alternativa al cálculo analítico que acabamos de hacer es usar Python para simular que arrojamos monedas al aire!\n",
    "\n",
    "1. Vamos a asumir que la moneda tiene la misma posibilidad de caer cara (1) o ceca (0)\n",
    "2. Vamos a lanzar una moneda 200 veces (igual número de veces que en nuestros datos)\n",
    "3. Vamos a repetir el punto 2, una gran cantidad de veces (por ej 10000 más o menos), en cada repetición vamos a contar la fracción de  caras obtenidas.\n",
    "4. Luego, vamos a contar todas las veces que obtengamos una fracción de  90/200 (o menos) caras."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAd8AAAFXCAYAAADj40TtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XecVOWh//HPmZntvQy7wMLC0hEE0ViiSGIvIIoIiCBe\nDV69ptzEGPUXY4tXzc3NjYnGfq8aK4II2LArqEjvLCBLWdhle58t087vjwWuhbK7zM6Z8n2/XiTg\n7sx+9+Ew3z3PPOc5hmmaJiIiIhI0NqsDiIiIRBuVr4iISJCpfEVERIJM5SsiIhJkKl8REZEgU/mK\niIgEmSMYX6SysjHgz5mRkUhtbXPAnzcSaaw6TmPVcRqrjtNYdVwkjZXTmXLEj4Xtma/DYbc6QtjQ\nWHWcxqrjNFYdp7HquGgZq7AtXxERkXCl8hUREQkyla+IiEiQqXxFRESCTOUrIiISZCpfERGRIAvK\ndb4iIhJ+TNMEd4CvuY1NxDCMwD5nGFL5iojI4bmbad30NoY9NiBPZ/rcxI8YD3FJAXm+I2lra+X+\n+/9AbW0tiYmJ/P7395GRkfGDz6utreXmm2/ghRdeJS4ujqamJu655//R0tJMTEwsd999P1lZ2fz8\n5zceekxx8R4uvng8N9/8i+PKqGlnERE5IsMeixETH5hfASrxY3nzzXkUFAzk8cef5aKLLuWFF/7n\nB5+zfPkyfvObW6ipqT7039599y0GDBjA448/y7nnns8rr7wIwGOPPc1jjz3NnXfejdPZg1mzbjju\njDrzFRGRkDF//nzefXcxzc3N1NXV8S//8jN+8pNzD3183769PPzwH7/zmPPPv4iJEycd+vOGDeuZ\nPv1aAE4//Uyef/6H5WuzGTzyyOPccMPMQ/9twICBFBfvBsDlcuFwfLci//73v3Dzzb8gMTHxuL9P\nla+IHJHb52f53no2lDdiMwxi7QY9U+IYm5+B0+pwErFaWlr461//QV1dLbNnz+Kss8YdKsK8vD48\n9tjTR328y+UiOTkZgMTERFyuph98zo9+dPoP/ltqahorVnzNjBlX0dDQwD/+8cyhj+3Y8Q0ul4tT\nTjn1eL61Q1S+IvID26tc/OcXu/l4Zw0ut+8HH4+xGYwtyGLS0GwmDe+Bw6Z3sCRwRo8eg81mIzMz\ni5SUVOrq6sjOzgY6duablJREc7MLgObm5kNFfCzPPfcM06dfy+WXX8mOHd9w112/44UXXgPggw/e\n5bLLrgjEtweofEXkW1o8Pv62rJhHvy7G4zfpn5HAeSMzOb1POnYbeHwmhZUuPt5ZzSc7qvhkRxV/\n/WoPt57Zj0nDe2DTKlYJgG3btgJQU1ONy+X6zmKpjpz5jhw5imXLvmT48BF8/fWXjBp1Uoe+bkpK\nyqGizsjIwOVyHfrYqlUrueaaWZ39Vo5I5SsiAJQ0tDJ1zga2VzfTOzWOB88bxMWDs3/weROHwR1n\n98dlt3P3O1t4dUMZ//ZWIS+uK+WRS4bSPyPBgvTSXUyfO+jPVVNTza9+dTNNTU3ceuvt2O2du9PR\nFVdM5oEH7uHmm28gJiaGe+55AIDXXnuJvLw+nHXWuMM+bvbsm3n44T/y5pvz8Hq93H7777+TKS0t\nvVM5jsYwTdMM2LMdQXfcz9fpTOmW541EGquOi9ax2lXbwuRX17G3oY3rTurF3T8tIDn26D+bHxyr\nPXUt3PNJEe9uryIxxsZd4wq44eTeupbzW8L1uLLiOt+lSz9k06atx30pTyg42v18deYrEuW+qXZx\n5avrKWtyc/vYfvzmx/mdKs789ASeu+IEFhRWcOeH3/D/PtrBV3vr+PslQ0mO00tMODMMo9uvyY1W\n+pchEsVqWjxc/fpGyprc3HfOAG4+tU+XnscwDK4YnsOP+6Zz48ItvL2tiu1Va3h+0ggGZh3/ZRkS\nPSZNmsTYseE3S9BZWqIoEqV8fpObFm2huL6VW8/M73LxfltOchzzpo3iX0/JY3t1M5e8uIYV++oD\nkFYksqh8RaLUn5bu4rNdtZw3IJPbzuoXsOeNsdv443kD+fulQ2ls83LVa+v5cEf1sR8oEkVUviJR\n6KOiah5ZVky/9HgenzCsWy4RmjYyl39eORKAa9/YyPwt5QH/GiLhSuUrEmUaWr3cungbMTaD5yaN\nID0+ptu+1vkDs5g7bRRJsXZueauQD3ZUddvXEgknKl+RKHP/Z0Xsb3TzmzPzOaFHx3b+OR6n5qXx\n8lUnEmu38bMFW/iquK7bv6ZIqFP5ikSRL/fU8s91+xnmTOIXp/cN2tc9LS+N5yadgM9vMmPeRrZU\n/HCvXZFoovIViRKtXh+/WbwdmwGPXDKEWHtw//mfU5DFPyYMo8nt4/o3N9PQ6g3q1xcJJSpfkSjx\n9Mp97KptYfYpeZzUM9WSDJcP68Etp/VhZ20Lv3p3K0HYYE8kJKl8RaJAhcvNI8uKyUqI4bYz+1ma\n5ffj+nNGnzTe2V7Fkyv3WZpFxCoqX5Eo8Kelu2hy+/jd2H6kxlu7sZ3DZuPpicPpkRTLHz/bqfd/\nJSqpfEUi3JaKJl5ev58h2YnMHN3T6jhA+05Yj1wyBK/f5Lfvb8ev6WeJMipfkQh336dF+E2475wB\nIXXT+/MGZDFxqJNVJQ28sLbU6jgiQRU6/xJFJOC+3lvHp7tqObtfBucUZFkd5wceOG8gqXF2Hvh8\nJ2WNbVbHEQkala9IhDJNk4eX7ALgjrH9rA1zBDnJcfzhJwNobPNxzydFVscRCRqVr0iE+mJPHV/t\nree8AZmc0jvN6jhHNHN0T0bnpvBmYQUbyyL/VnIioPIViUimafLw0vaz3tvH9rc4zdHZDIP/N649\n40MHMotEug6V7/r165k5cyYAe/bs4eqrr2b69Oncc889+P3+bg0oIp336a4aVpY0cNGgLEblplgd\n55jG9cvgzL7pfFRUw3Ld/1eiwDHL95lnnuGuu+6ira19McRDDz3Ev//7v/PKK69gmiYff/xxt4cU\nkY4zTZP/+nIPAL87K7TPeg8yDIM7z27P+uDnO7XzlUS8Y5Zv3759efTRRw/9efPmzZx66qkAnH32\n2Xz11Vfdl05EOm3Z3npWlTRw4cAsRuR0/12LAuXUvDTOH5DJsr31fLa71uo4It3qmFvdXHjhhezb\n939bwJmmiXHgxttJSUk0Nh57gURGRiIOh/04Yh6e0xn602mhQmPVceE+Vk+8uRmAey8Z1u3fS6Cf\n/0+XjeDDvy7hqdWlTDk1P6DPbbVwP66CKRrGqtP7zNm+dZG+y+UiNfXYG7TX1jZ39ssck9OZQmWl\nVkZ2hMaq48J9rNaXNfL+tkrO7JvOgERHt34v3TFWeXE2zu6XwSc7qvhkYykjw+D96o4I9+MqmCJp\nrI72Q0SnVzsPHz6c5cuXA7BkyRJOOeWUricTkYD6+7JiAH51RvDu1Rto/3ZqHgCPr9xrcRKR7tPp\n8r399tt59NFHmTp1Kh6PhwsvvLA7colIJ+2obubtbZWMyk1mXL8Mq+N02U/7ZzI0O5GFhZWUNLRa\nHUekW3Ro2jkvL4/XX38dgP79+/PSSy91aygR6bwnV+7FBH55et9D6zLCkWEY3HxqH3717jaeWVXC\nvecMsDqSSMBpkw2RCFDd7Ob1TeX0TYvnksFOq+Mct0nDc+iRFMuL60tpbPNaHUck4FS+IhHghbWl\ntHr93HhKHnZb+J71HhTnsHH9mF40tvmYt7nc6jgiAafyFQlzbV4//7umlJQ4O9NPzLU6TsBcM6on\nDpvB82tLtemGRByVr0iYe7OwggqXm5mjepEc1+mrB0NWTnIclwzOprDSxcqSBqvjiASUylckjJmm\nyZMr9mI34Gcn97Y6TsDNGt0LaJ9WF4kkKl+RMPb13nq2VLoYP8RJXlq81XEC7qz8dAZkJrBoawU1\nLR6r44gEjMpXJIw9t7YEgBsi8KwX2i87mjW6F20+k9c2lFkdRyRgVL4iYaq8qY23t1UxzJnEaXlp\nVsfpNlNH5hJnN/jnOi28ksih8hUJUy+t34/Xb3LdSb3CelONY8lIiGH8UCc7a1tYoYVXEiFUviJh\nyOv38891pSTH2rnqhByr43S7aSPbL6Gas1FTzxIZVL4iYej9b6rZ3+hmyoiciLq86EjG5mfQOzWO\nBYUVNHt8VscROW4qX5Ew9PyBS2+uOykyF1p9n80wmDoilya3j3e2VVodR+S4qXxFwkxxXQtLdtdy\nal4qQ51JVscJmikj26fXX9PUs0QAla9ImJmzqRwTuObEnlZHCaqCjEROy0vjiz117K3XrQYlvKl8\nRcKI3zR5bcN+EmNsTBga/ncv6qyrR+ZiAnM36exXwpvKVySMfLGnjr0NbVw+rAfJsZG/0Or7Lhvq\nJDHGxmsby3TNr4Q1la9IGHllw34Aro6yKeeDkuMcXDrYye66Vpbvq7c6jkiXqXxFwkRdq4d3tlUy\nIDOBU3unWh3HMgev+dXCKwlnKl+RMPHmlgrafCZXn9gzone0OpYz89PpkxrHwq2VuNy65lfCk8pX\nJEy8saUcA6JiR6ujsRkGU0bm4nL7eGe7rvmV8KTyFQkDxXUtrNjXwFn56fRMibM6juWmjDgw9aw7\nHUmYUvmKhIE3CysAuHJ4dJ/1HtQ/I4Ez+qTxRXEdxXUtVscR6TSVr0iIM02TNzaXE2s3uHRIttVx\nQsbBhVdzN5dbnESk81S+IiFuc4WLrVXNnD8gi7T4GKvjhIwJQ5wkOGzM3VSua34l7Kh8RULc/C3t\nZ3aTNOX8HclxDi4enM3O2hbW7G+0Oo5Ip6h8RUKY3zR5s7CClDg75w/MtDpOyDm48lvbTUq4UfmK\nhLAV++opaWjj0sFO4h12q+OEnHH9M3AmxbCgsAK3z291HJEOU/mKhLC3D9y7dmIU3kShIxw2G5OG\n5VDT4uXjohqr44h0mMpXJESZpsk726tIjbMztl+G1XFC1pQRB6aeN2vqWcKHylckRK3d30hJQxsX\nDMwm1q5/qkcyIieZodmJfLCjmrpWj9VxRDpE/6JFQtTBKefxurb3qAzDYPKIXNw+k0Vbtd2khAeV\nr0gIMk2Tt7dVkhhj46f9tcr5WCYP74EBzN2kDTckPKh8RULQlkoXu+taOW9AFgkxWuV8LL1S4zkr\nP53l++rZre0mJQyofEVC0P9NOWuVc0cdvOZ3ns5+JQyofEVC0DvbKomzG5xXoCnnjhp/cLvJzdpu\nUkKfylckxOyobmZrVTM/7Z9JcpzD6jhh4+B2k7tqW1hd2mB1HJGjUvmKhJiDU86Xasq50/7vml9N\nPUtoU/mKhJi3t1XisBlcOCjL6ihh5+x+B7ab3KLtJiW0qXxFQsieuhY2lDdxVn466bp9YKc5bDYm\nDc+htlXbTUpoU/mKhJB3tlUBWuV8PKacoO0mJfSpfEVCyNvbK7EZcPEg7WrVVdpuUsKBylckROxv\nbGNVSQOn56XhTIq1Ok7Y+vZ2kwsLtd2khKYula/H4+HWW29l2rRpTJ8+naKiokDnEok6727XxhqB\ncmi7SU09S4jqUvl+/vnneL1eXnvtNW655RYeeeSRQOcSiToH3+/VJUbH7+B2kyv2NWi7SQlJXSrf\n/v374/P58Pv9NDU14XBoIwCR41HV7OarvXWc3CuVnilxVseJCFeNyAW03aSEpi61ZmJiIiUlJVx8\n8cXU1tby5JNPHvXzMzIScTgCvzm805kS8OeMVBqrjrNirBZ+vQe/CdPG5IXV31UoZ73uxwnc/sF2\n5m+t5OHLR2AYhqV5QnmsQk00jFWXyvf555/nrLPO4tZbb2X//v3MmjWLt956i7i4w//EXlvbfFwh\nD8fpTKGysjHgzxuJNFYdZ9VYvbp6HwA/yQufv6twOK4uHpTN/C0VvLe+hB/1TrMsRziMVaiIpLE6\n2g8RXZp2Tk1NJSWl/UnT0tLwer34fL6upROJcvWtHpburmVkTjL56QlWx4koh7ab1NSzhJgule91\n113H5s2bmT59OrNmzeLXv/41iYmJgc4mEhXe31GNx29qlXM3OLtfBj2SYllYqO0mJbR0ado5KSmJ\nv/3tb4HOIhKV/u/evdpYI9Dat5vswZMr9/FRUTWXDNYPOBIatMmGiIWa3F4+21XL4KxEBmUlWR0n\nIl2lqWcJQSpfEQt9XFRDq9evKeduNKJHMsOcSXxYpO0mJXSofEUstHBrBaBdrbqTYRhMPiFH201K\nSFH5ilikye3lo6IaBmYmcEIPTTl3pysPbDf5+iZtNymhQeUrYpEPd1TT6vUzcVgPyzeAiHQHt5tc\nWdLArlptNynWU/mKWGTh1vYp0IlDe1icJDpMObjd5GYtvBLrqXxFLNDY5uXjomqGZCcy1Kkp52C4\ndEg2iTE25m4qwzRNq+NIlFP5iljg/R3VtPlMLtNZb9Akxzq4eHA2u+taWVXaYHUciXIqXxELLCps\nX+U8cahWOQfTVSe0Tz3rml+xmspXJMgaWr18squGYc4kBmdryjmYzu6XTo+kWBYUVtDq1X70Yh2V\nr0iQvfdNFW6fqbNeCzhsNqaMyKGu1ct726usjiNRTOUrEmSLDmysMXGY3u+1wjWjegLw0vr9FieR\naKbyFQmiulYPn+2qZUSPZAZk6k5gVhiQmcgZfdJYuqeOPXW65lesofIVCaLF26vw+E0mDtOUs5Wm\nn9h+9vvqBu14JdZQ+YoE0YIDG2voEiNrTRjqJCXOzqsb9+P16z6/EnwqX5EgqWnxsGR3LaNyk+mf\nkWB1nKiWGGPnyuE57G908+nOWqvjSBRS+YoEyXvbq/D6tbFGqJihhVdiIZWvSJAcvH3gZbrEKCSc\nmJvCyJxkPthRRXlTm9VxJMqofEWCoKrZzdLdtZzUM4X8dE05h4prRvXEZ8Ic7XglQabyFQmCd7dX\n4TN1B6NQc+XwHsQ7bLyyfr9utiBBpfIVCYKFB/ZynqAp55CSFh/D+CFOdta28PXeeqvjSBRR+Yp0\ns0qXmy+L6zi5Vyp90uKtjiPfo4VXYgWVr0g3e3tbJX4TLtfGGiHpjD5p9M9I4K1tldS3eqyOI1FC\n5SvSzRYd2FhjwhCVbygyDINrRvWk1evnjS0VVseRKKHyFelG5U1tfFVcx6l5qfRK1ZRzqJo6Ige7\nAS9r6lmCROUr0o3e3laFiVY5h7qc5DguGJjNxvImNpQ1Wh1HooDKV6QbLSyswEBTzuHgmlG5gBZe\nSXCofEW6yf7GNpbvq+f0PmnkpsRZHUeO4ZyCTHKTY3ljSznNHp/VcSTCqXxFusnb2yox0R2MwoXD\nZuPqE3NpbPPx1oFFciLdReUr0k0Wbq3AZsD4IdlWR5EOuvrAfX5f2aCpZ+leKl+RblDS0MqKfQ38\nuE86Ocmacg4X/dITGJufzrK99RTVNFsdRyKYylekGxyctrxMG2uEnYM7XumyI+lOKl+RbrBwayU2\nAy4drPINNxcPziYj3sFrG8vw+PxWx5EIpfIVCbC99a2sLm3grPwMnEmxVseRTop32Jl8Qg5VzR4+\n2FFtdRyJUCpfkQBbtLV9i8KJuoNR2Jp+cOpZC6+km6h8RQJs0dZK7AZcMlirnMPVCT2SGdMzhU92\n1lDS0Gp1HIlAKl+RANpT18La/Y2M7ZdBVqKmnMPZNaN64jfhtY1lVkeRCKTyFQmghQdWOV+ujTXC\n3hXDepAYY+PVDWWYpml1HIkwKl+RAFpUWIHDZnCxppzDXnKcg0sHOymub2X5vnqr40iEUfmKBMjO\n2mY2lDcxrl8GGQkxVseRALhqRA4A8zaXW5xEIo3KVyRADm6sMXGYppwjxdj8DHKSY1m0tZI2r675\nlcDpcvk+9dRTTJ06lUmTJjF37txAZhIJSwsKK4ixGVw8KMvqKBIgdpvBpOE9qGv18lGRrvmVwOlS\n+S5fvpy1a9fy6quv8uKLL1JWptWAEt2+qXaxucLFT/tnkhavKedIMvkETT1L4Dm68qAvvviCwYMH\nc8stt9DU1MTvfve7QOcSCStvbmnfWOOK4ZpyjjQjeiQzNDuRD4uqqWv1kK4friQAulS+tbW1lJaW\n8uSTT7Jv3z5uvvlmFi9ejGEYh/38jIxEHA77cQU9HKczJeDPGak0Vh3X2bEyTZO3tleREGNjxhn9\nSI7r0j+rsBQtx9V1p+VzxzuFfLqvkRvPyO/Sc0TLWAVCNIxVl14l0tPTKSgoIDY2loKCAuLi4qip\nqSEr6/DvddXWBv7WXE5nCpWVjQF/3kikseq4rozVxrJGtlW6uGyok5aGFlq6KVuoiabj6sL8NO4A\n/rmimCsGZnb68dE0VscrksbqaD9EdOk935NPPpmlS5dimibl5eW0tLSQnp7e5YAi4ezNwgNTzlrl\nHLF6p8Zzal4qXxXXUd7UZnUciQBdOvP96U9/ysqVK5k8eTKmaXL33Xdjtwd+Wlkk1PlNkwWFFaTE\n2Tl3QOfPiCR8XD60Byv2NfDW1kp+dkqe1XEkzHX5zSktshKBVSUN7GtoY+qIHOK7YV2DhI4JQ538\n/qMdLNhaofKV46ZNNkSOg1Y5R4+c5Dh+3DedFfsadKcjOW4qX5Eu8vr9LNxaQVZCDGPzM6yOI0Fw\n+YH39RcWVlqcRMKdyleki74srqOq2cP4oU5i7PqnFA3GD8nGbsDCrRVWR5Ewp1cMkS5acGDKeZJW\nOUeNrMRYxvbLYO3+RnbXRctFZdIdVL4iXdDm9fP2tip6psRyWp80q+NIEB28V/OiQp39StepfEW6\n4LNdNdS3eZk4tAe2I+zsJpHpkiHZxNgMFuh9XzkOKl+RLji0sYZWOUed9PgYftI/g00VTeyoDvzu\nfRIdVL4ineRy+1j8TRX90uMZnRv5e9DKDx1c9bxAU8/SRSpfkU76sKiaZo+fK4b3OOLNRCSyXTQo\nmzi7oVXP0mUqX5FOen1T+/2rJw3PsTiJWCUlzsG5A7LYVtVMYWWT1XEkDKl8RTqhvKmNT3fWMDo3\nhSHZSVbHEQtpww05HipfkU54Y3MFPhOmjtRZb7Q7f0AWiTE2FhRWYJqm1XEkzKh8RTrINE3mbCoj\nxmZwhaaco15SrJ0LBmaxs7aFjeWaepbOUfmKdNCm8iYKK12cPzCLzIQYq+NICJg4VKuepWtUviId\nNOfAQqupI3ItTiKh4twBmSTH2lmoqWfpJJWvSAd4fH7mb2m/g9G5AzKtjiMhIt5h5+JB2extaGN1\naYPVcSSMqHxFOuCDHdVUNXuYNLwHsbqDkXyLVj1LV+hVRKQDXly/H4AZo3tanERCzbj+GaTHO1i4\ntQK/pp6lg1S+Isewt76VT3fWcErvVIY5k62OIyEm1m7j0sHZlDW5Wb633uo4EiZUviLH8PL6/ZjA\nzFE665XDu/zADTYWaLtJ6SCVr8hReP1+Xt24n5Q4O5cN1R2M5PDO7JtOdmIMb22txOv3Wx1HwoDK\nV+QoPi6qYX+jmyuH55AUa7c6joQoh83G+CFOqpo9fFlcZ3UcCQMqX5GjeHFd+0KrmVpoJcdwhVY9\nSyeofEWOoLiuhY92VnNSzxRG5ui+vXJ0p/VJIzc5lre3VeL2aepZjk7lK3IEL6wrxW/C9WN6Wx1F\nwoDNMJg4tAd1rV6W7K61Oo6EOJWvyGG0eHy8vH4/WQkxTBzmtDqOhImDx8qbW7TqWY5O5StyGAsL\nK6hp8XLNqJ7EO7TQSjrm5F6p9EmN471vqmj1+qyOIyFM5SvyPaZp8uzqEmwGzDqpl9VxJIwYhsHE\nYT1ocvv4uKjG6jgSwlS+It+zurSBDeVNXDgwmz5p8VbHkTBzcNXzG1vKLU4ioUzlK/I9/7umFIAb\nTtZCK+m8ETnJDHMm8cGOampbPFbHkRCl8hX5lvLGNhZtrWBQViJj89OtjiNhyDAMrhqRg9tnslDb\nTcoRqHxFvuXZ5Xtw+0yuH9MbwzCsjiNhavLwHGwGvL5JU89yeCpfkQO8fj9PfLWHpFg7U0bkWB1H\nwlhuShxn98tgVUkDO2uarY4jIUjlK3LAe9urKalvZeqIHFLiHFbHkTA3ZUQuoLNfOTyVr8gB/7um\nBNCOVhIYFw/KJinWzrzN5fj9ptVxJMSofEWArZUuviyu49xB2QzOTrI6jkSApFg7E4Y4Ka5v5fOd\n1VbHkRCj8hXh/856/+00baohgTP9xPap52e/LrY4iYQala9EvYZWL69vKqN3osElBclWx5EIclpe\nGoOyEnlj435d8yvfofKVqDdnUxnNHj8z8lw47Lq8SALHMAyuObEnbV4/c7XwSr5F5StRzW+a/O+a\nEuLsMC3PbXUciUBTRuYQYzd4aX0ppqmFV9JO5StR7fPdtRTVtDChl0lWnM56JfCyE2O5YkRPtlY1\ns6q0weo4EiKOq3yrq6sZN24cRUVFgcojElTPrW5faDUrr8niJBLJZp/eF4CX1u23OImEii6Xr8fj\n4e677yY+Xnd9kfC0t76V93dUMzrTYHS2NtWQ7nPOwGzy0+NZUFhBXasWXslxlO+f/vQnpk2bRo8e\nPQKZRyRoXlpfignMzNP2f9K9bDaD607qRYvXz6vrdfYr0KUf9+fPn09mZiZjx47l6aefPubnZ2Qk\n4nDYu/KljsrpTAn4c0YqjdV3eXx+XtlYTnqcjWuGxZMYY2B6249RjVXHaaw67pfnDOI/l+7k+TV7\n+f0lw7HbtMbgSKLhuOpS+b7xxhsYhsGyZcsoLCzk9ttv54knnsDpdB7282trA39m4XSmUFnZGPDn\njUQaqx96a2sF5Y1t3FDgw9Psph4wvW4SQWPVQTquOs7pTMHnauOKfAevFHmYs3wP5w/MsjpWSIqk\n4+poP0R0qXxffvnlQ7+fOXMm99577xGLVyQUPb+2FIBr+rQAgZ+VETmc6wbH8EqRl2dW7FT5Rjld\naiRRZ0d1M0v31HGG02BQqopXgueEDDunZvr4rNhFkW41GNWOu3xffPFFBgwYEIgsIkHxz3UHznp7\nuyxOItHouv4+AJ5dscfiJGIlnflKVGnx+JizsYysOIOLe2vBiwTfRbk+esb7eXVTOfW67ChqqXwl\nqry1rZLaVi/T+niI1dEvFoixwfX9fTR74Z8H1h5I9NHLj0SVF9aWYgDT+rRaHUWi2NX5XpLsJs+u\n3IPb57dPKXs+AAAYxElEQVQ6jlhA5StRY3NFEytLGhiXY9AvWQutxDppMTCtr4/9zX4WFlZYHUcs\noPKVqHFwodUM7eMsIeCGAi82TJ5YtlN3O4pCKl+JCk1uL3M3lZObaHBurg57sV6fRJNLevrYVO3m\niz11VseRINOrkESF+VsqaHL7mJ7nxqGjXkLEjQO8ADzx9U6Lk0iw6WVIIp5pmjy3pgS7AVPz2qyO\nI3LISRkmP8rw8dHuRrZX6brzaKLylYi3oqSBzRUuLuwFvZK00EpCy+wDZ79Pfr3b2iASVCpfiXjP\nrSkBYJYWWkkIuiDXT36in7lbKql0ua2OI0Gi8pWIVuFy89bWSgan2TjDqbNeCT12A35W4KXND8+t\n3mt1HAkSla9EtJfWleLxm8zs04Kh3SQlRE3p4yMtxuS5Nfto8fisjiNBoPKViOX1+/nnuv0kxRhM\n6uW1Oo7IESU6YEa+l+pWk7mbyq2OI0Gg8pWItfibakob25jcx09qnA51CW3X9fcSY5g8uXwnfm26\nEfH0iiQR6+BCqxm9tdBKQl9uPEzs7WNHnZePi2qsjiPdTOUrEembahdL99RxRg8bQ9MdVscR6ZCD\nm248vqzI4iTS3VS+EpGeW9O+j/O1OuuVMDIs1WRsto8vS5rZWNZodRzpRipfiThNbi9zNpWRm2Bw\nQU+r04h0zuxDW07usjiJdCeVr0SceZvLaWzzcU1fN7F2XV8k4eUnTj+DU/ws2FZDaYPuOx2pVL4S\nUdr3cS7FYcA07eMsYcgwYHaBF68Jz6wstjqOdBOVr0SUz3bXUljp4tLekJuoHa0kPF3R24czzuTF\ndftpatM16pFI5SsR5R/L27fnm91Xi1UkfMXZYVY/Lw0ek5fX77c6jnQDla9EjI3ljSzZXcuPexiM\nytLlRRLeZuZ7ibeZPL1iN16/3+o4EmAqX4kYj69oP+u9KV+XF0n4y4yDyX187G3y8d72KqvjSICp\nfCUilDS0smBLBUPTbPwkR4e1RIbZBV4MTB5fttPqKBJgepWSiPDUyn34TJid79LdiyRiFCSbnJfj\nZ3V5Kyv21VsdRwJI5Sthr6rZzT/XlZKbaDAxTxvSS2S5UZtuRCSVr4S9p1buo9nj5+Z+rcRpUw2J\nMKdl+jkxzc+7O+rYVdtidRwJEJWvhLXaFg/Pri7BGW9wdb5uQi6RxzDaz35N4KkVe6yOIwGi8pWw\n9vSqfbjcPv61oI0Eh856JTJd0tNHr3g/r24so67VY3UcCQCVr4SthlYvz6zaR2acwYy+ekGSyBVj\ng+sLfLR44YUD96mW8KbylbD19Kp9NLT5mN2/jaQYHcoS2a7u6yXZYfI/q/bi9mnTjXCnVywJS7Ut\nHp5YuZfMOINZ/XTWK5EvNaa9gMuafby5pcLqOHKcVL4Slh5bXkxjm49bBrSSorNeiRI39PdhN0ye\n+HonpqnL6sKZXrUk7JQ3tfHsqhJ6JtqYqRXOEkV6J5pc2tPHlmo3S/bUWh1HjoPKV8LO35YV0+L1\n84sBLVrhLFFn9oFNNx7/SltOhjOVr4SV4roW/rmulPxkG1P76KxXos/odJPTMn18WtzExnLdOjNc\nqXwlrDy4ZBdun8lvBrqI1W5WEqV+Pqj97Pe/lxZZnES6SuUrYWN9WSPzt1QwMsPg8jyr04hYZ5zT\nz6h0P+/sqGNrpcvqONIFKl8JC6Zpct8n7T/l/35wIzad9EoUMwz41aD2S+z++oXOfsORylfCwsc7\na/iiuI5zcg3O7GG3Oo6I5c7L8XNCqp+F22soqmm2Oo50kspXQp7H5+e+T4uwGXDnoAar44iEBMOA\nXw7y4DfhL3rvN+x0qXw9Hg+33XYb06dPZ/LkyXz88ceBziVyyPNrS9lW1cy0fD9D0x1WxxEJGRf1\n9DM0xc8bhdUUVjZZHUc6oUvlu2jRItLT03nllVd49tln+eMf/xjoXCIAVDW7+c+lu0mNNbhtoKbW\nRL7NZsAdwzyYwIOffmN1HOmELp1GXHTRRVx44YVA+0IYu13vwUn3eGjJLurbvNx3gpvsBB1nIt93\nTg8/p2b6eH9nPcv31XNaXprVkaQDDPM4Nghtamri5ptvZsqUKUyYMOGIn+f1+nA49MIpnbN2Xz0n\nP7KEoRkOll7o6fbrek2vm8RhP8WRnNmtX0eiV3PRCnwN5QF/3uWVcP77cGbfFJb+chyGocsBQl2X\n30Dbv38/t9xyC9OnTz9q8QLU1gZ+utDpTKGyUru7dEQ4jpVpmtw0dx2mCXcNrKelyUZLd39Nr5tE\nCLuxsko4HldWOThWnjoX/sbAH8lDY+H8nFg+LG7k5WW7uXBQdsC/RrBE0nHldKYc8WNdes+3qqqK\n66+/nttuu43Jkyd3OZjIkczfUsGKfQ1c3Nvg7Bwtyhc5ljuGebAbJn/4cButXm29Guq69Kr25JNP\n0tDQwOOPP87MmTOZOXMmra2tgc4mUarJ7eX+z4qIsxv8flC91XFEwsLgFJPr+vnY3eDhiRV7rY4j\nx9Claee77rqLu+66K9BZRAD4+7Ji9je6+dVgH/kpMVbHEQkbvxniYWGJnUe+2sOUEbn0To23OpIc\ngebzJKTsqm3hiRV76Zlo4+YCzaaIdEZqDNw53EOL1+Sej7ZbHUeOQuUrIcM0TX73/nbafCZ3DXGR\nFKPDU6SzJuf5GJPhZ9H2Gj7ZWWN1HDkCvbpJyJi3uZzPd9dyTi5M6G11GpHwZDPgoZFuHIbJb97d\nQmOb1+pIchgqXwkJ1c1u7v64iASHwR+HNqDLFEW6bniayc8HeSlt8nL/pzusjiOHofKVkHDvJ0VU\nt3j47WA3fbXISuS4/WKQlyEpfl5YV8bS3bVWx5HvUfmK5T7cUc2cTeWMzLDxL/08VscRiQixNviv\n0W5smPz7O1toaNX0cyhR+Yqlalo8/Pq9bcTaDP5yQgMx3byFpEg0GZ3ePv28t9HDrYu3chy7CUuA\nqXzFUnd+8A0VLje3DnEzLEP7f4sE2q8Hezk5w8fCrVW8uqHM6jhygMpXLLOwsII3Cys4OcvgxgJN\nN4t0B4cNHh3jIdVhcueH29le5bI6kqDyFYvsqWvh1sXbiHcY/PeIRhw2TTeLdJc+iSb/OcpNi9fk\nX+Zv1OVHIUDlK0HX5vUze8EWGtp8/McJrRSkarpZpLtd2svPDf29fFPTyi1vFeLX+7+WUvlK0N3/\naRHryhqZ0tdkSl+/1XFEosZdwz38OMvH4h3V/NcXu62OE9VUvhJUCwsreGZ1CYPTbNw/PPD3eRaR\nI3PY4ImT3fRJ8PNfX+5hYWGF1ZGilspXgmZNaQO/eGcryTEGj5/YoL2bRSyQGQf/c6qbZIfJLW9t\n4cs92oDDCnr1k6DYV9/KzDc24vb5+cfoZoak631eEasMSzV55hQ3pmky642NbKlosjpS1FH5Srdr\naPUyY95GKl0e7hnu5pxcqxOJyFlOP38Z7aHB7WfanHXsqm2xOlJUUflKt3K5fUyft4EtlS5m9vfx\nL/19VkcSkQOuyPNx93A3ZS4vk15Zy+46FXCwqHyl27R6fcyav4kV+xqYmGdy//A23a1IJMTMHuDj\nzmEeShrdXPHyWopVwEGh8pVu0eLxccObm1myu5YLesJ/n9iijTREQtS/DfTyu6HtBTzx5bXsqNaV\nCN1N5SsB19DqZdrrG/iwqIZxOfDYaBexumGCSEj7xSDvoTPgCS+tYUNZo9WRIprKVwKqwuXm8lfW\nsWxvPeN7w7NjXCQ4dJiJhIN/G+jloZFualo8XP7yWr7QZUjdRq+KEjAbyhq56IXVbKpoYkY/H4+O\nbiZexSsSVmb08/HYGA9tXh9T5qznlfX7rY4UkfTKKAExf0s5E15aS0lDG7cN9fDgCLfe4xUJU5f1\n9vHy6W6S7X7+/b1t/PGzIu0FHWAqXzkuzR4ft72/nZsWFWI3TP7nlGZ+OcirVc0iYe6MbD+LznLT\nP9HHo1/v5YY3N9Ps0aWCgaLylS5bX9bIec+t4oW1pQxNs7HojAbO76nWFYkU/ZNNFo51c3qml3e2\nV3H5y2spb2qzOlZEUPlKpzW1ebn74x1c9MJqdtS08LMBPhb9uJFBaQ6ro4lIgGXEwstneLgqz8O6\nsiYueH41q0rqrY4V9lS+0mF+02RBYQVnPruCJ1fuIy/JxiunubhnuFsrmkUiWKwN/jK6/VKk8qY2\nJr68jv9dU4Kp94G7TKcqckymafLprloe/HwnG8qbiLXBrwZ7uWVAm0pXJEoYRvulSCPT/Px8TQx3\nfPANq0oa+K+LBpMYoxuldJbKV47Ib5os/qaKR7/ey+rSBgAm9oFbCxrpn+pAEyci0Wes0897Z7v5\n15UxzNtczuaKJp6bdAIFGYlWRwsrKl/5gbLGNuZuLueVDfspqmnf5/WCnvDrAY2MyHCgw0YkuvVK\nMJl3ppv7Njt4cY+LC55fzV8uGsLEYT2sjhY29CoqQPtNEBZ/U81rG8v4bFcNfhPi7DAt32R2vovB\naXZ0uIjIQXF2ePBELydnmNyxIYbZC7fwwY5qHjp/EKnxeq04Fo1QFDNNk9WlDby2sYyFhZXUt3kB\nGJNtZ3LPZsb39JIRZwP0fo6IHN6VfXyclOHnl2timLu5nGV76/jPCwdz3oAsq6OFNJVvFCptaGXu\n5nJe21h2aFq5Z6KNGYNNrsxpZFB6zIHP1Hu6InJsBckmb57l5m/bHfxjh8n0uRuZMMTJf5w3kNyU\nOKvjhSSVb5Ro9vh4b3sVr20sY8nuWkwg3m5weV8bk3MaObMHOGw2IOZYTyUi8gMxNvjtUC/je/m4\nY4ODt7ZV8lFRNTef2odbTutDSpzq5ts0GhHMNE2W76tn0ac7eW1tCU3u9q3hfpRt48qeLYzv6SEt\nzo7OcEUkUIammsw/08Pre/38eauD//5qDy+sLeGmU/sw66RepMfrB3xQ+UakvfWtzN1UxpxN5eyq\nbZ9W7pVocH0/L1f2bKEg9eBfu97LFZHAsxkwra+Pib18PLvLzhM7TP7j81389as9zBjVk1kn9WJQ\nVpLVMS2l8o0QLrePt7dV8vqmMpbuqQMg3g5X9jW5tsDL6GQv7TcZ0l+5iARHggN+McjHrH4+Xtlj\n4392xfD0qhKeXlXCj3qnMm1kLpcMziYrMdbqqEGnV+Iw5jdNvt5bz5yNZSzaVonrwLTyaU6DybnN\nXNLLJDXWIC0tgfp6r8VpRSRapcbATQP9XF/QxuL9NuYU21ha0sDKkgZ+9/52ftw3nQlDnVw8KBun\nM8XquEGh8g0zpmlSWOli0dZK5m0up7i+FYA+SQY39vMwqVcb/VIOTifrDkMiEjpibXBZbz+X9fZT\n0uzjrRJ4tyyGpXvqWLqnjtvf/4azCjL5aX465xZkMjQ7CSNC70+q8g0TWytdLNxawaKtlXxT3QxA\nosNgSl+TK3NdnN7DfmBaWe/jikjo651octMguGmQm9IWg3dL4b2yGL7YWcPSnTXc/+lOeqXEcU5B\nJucUZHJ2fkZEbd4ROd9JhGnx+Fi2t55Pd9Xwyc6aQ4Ubb4dLesOlPVycm2OQFGOgwhWRcNYrweRn\nA+BnA9y44xJYtKONz8oNPq+Gl9bv56X1+7EZMDInmdP7pPOj3qmMyk2hb1p82J4Zd6l8/X4/9957\nL9u2bSM2NpYHHniA/Pz8QGeLGm1eP99UN1NY2cSG8iZWldSzoawJj7/9dl3xdri4N1zqdHFuLiTH\n2NDlQSISiZzxMLmPn8l9wGe2sL7O4JMy+KomhvUVTawva+Kple2fmx7vYHB2IgMzE+mXnkBOcizO\npAO/EmPITool1h6ar5VdKt+PPvoIt9vNnDlzWLduHQ8//DBPPPFEoLOFDb9p0ur103bgV6vPj/tb\nv2/2+Kht8VLT4qH2wK+aFi+lDa0U17dS2tiG/1u3xXQYMCLD4LQMD2dntnGq0068HVS4IhJN7AaM\nyTAZkwHgpsUH62ptrK/xs6HBwaZGk9UlDazY13DE50iMsZEa5yA1zkFKnJ2Ug7+Pbf99+5/tOJNi\nuWRwNvGO4Mwkdql8V69ezdixYwEYPXo0mzZtCmioo/H4/Mx8YxOlTW58Pv8RP+9oExHHmqU42ofd\nPpM234FiPVCwHn/XbihtADkJBqdkwaAkH0OS2hiW4mdUlp2EQ3//mlIWEQFIsMMZ2X7OyAbwAl7c\nftjjMih2QVWrn8o2qHLbqfI4qGozaPBCg8dNVZObXbXgOXJt8Mxlw5g4PCco30uXyrepqYnk5ORD\nf7bb7Xi9XhyOwz9dRkYijgD9NNHm9dHk9VPT2IzJ4UvPNPlWgxrfKVPzGD1pfut/D/e8sfb2bRmT\nYiE+AeJtBnEOg3hb+10+4m3mt/7fJM5mEm/3kxULGbGQEW+QmWAnK95GbgLE2Q+mswPdd61bWkby\nsT8pypl+H4bNETWXOgSCxqrjnM4UWpsz8DharI4S8jr7euXMglM6+LmtPpMGN9R7oPHA/ze0+fH7\nPEwalUlyenCO6S6Vb3JyMi6X69Cf/X7/EYsXoLa2uStf5ojenj4apzOFysrGgD5vsJlAaxC+TiSM\nVbCkJmqsOkrHVccdGqvEAe2/5IiCcVylHvj1fS0eaAng1z7aD6ddehNxzJgxLFmyBIB169YxePDg\nriUTERGJQl068z3//PP58ssvmTZtGqZp8uCDDwY6l4iISMTqUvnabDbuv//+QGcRERGJCrp2RURE\nJMhUviIiIkGm8hUREQkyla+IiEiQqXxFRESCTOUrIiISZCpfERGRIFP5ioiIBJnKV0REJMgM0zzW\nfX5EREQkkHTmKyIiEmQqXxERkSBT+YqIiASZyldERCTIVL4iIiJBpvIVEREJMofVAb7P7/dz7733\nsm3bNmJjY3nggQfIz88/9PH333+fp59+GsMwmDBhArNmzQLgiiuuIDk5GYC8vDweeughS/IH07HG\n6qA//OEPpKWl8dvf/rbDj4lEXRkv0LF1uLF6/vnnmTt3LpmZmQDcd9999OvXLyqPra6MVUFBgY6r\nw4zVhg0bePjhhzFNE6fTyZ///GdiYmIi87gyQ8z7779v3n777aZpmubatWvNm2666dDHvF6vef75\n55sNDQ2m1+s1L7jgArO6utpsbW01J06caFVkyxxtrA569dVXzSlTpph//vOfO/yYSNWV8dKxdfix\nuvXWW82NGzd26jGRqitjpePqh2Pl9/vNyy67zNy9e7dpmqb5+uuvm0VFRRF7XIXctPPq1asZO3Ys\nAKNHj2bTpk2HPma323n33XdJSUmhrq4Ov99PbGwsW7dupaWlheuvv55rr72WdevWWRU/qI42VgBr\n1qxh/fr1TJ06tcOPiWRdGS8dW4cfq82bN/P0009z9dVX89RTT3XoMZGqK2Ol4+qHY7Vr1y7S09N5\n/vnnmTFjBnV1dRQUFETscRVy085NTU2HpmKgvXC9Xi8OR3tUh8PBBx98wP3338+4ceNISEggPj6e\nG264gauuuordu3cze/ZsFi9efOgxkepoY1VRUcE//vEPHnvsMd57770OPSbSdWW8dGy1+/5xcuml\nlzJ9+nSSk5P5+c9/zqeffhq1x1ZXxqpXr146rvjuWNXW1rJ27Vruvvtu+vbty0033cSIESMi9rgK\nufTJycm4XK5Df/b7/T8Y5AsuuIDzzjuPO+64gwULFjBhwgTy8/MxDIP+/fuTnp5OZWUlPXv2DHb8\noDraWC1evJja2lpuvPFGKisraW1tpaCgoEPjG6m6Ml7jx4/XscV3x8o0TWbNmkVKSgoA48aNY8uW\nLVF7bHVlrM4880wdV3x3rNLT08nPz2fAgAEAjB07lk2bNkXscRVy085jxoxhyZIlAKxbt47Bgwcf\n+lhTUxMzZszA7XZjs9lISEjAZrMxb948Hn74YQDKy8tpamrC6XRakj+YjjZW1157LfPnz+fFF1/k\nxhtvZPz48UyaNOmoj4l0XRkvHVuH/3c4fvx4XC4XpmmyfPlyRowYEbXHVlfGSsfVD8eqT58+uFwu\n9uzZA8CqVasYNGhQxB5XIXdjhYOr4bZv345pmjz44INs2bKF5uZmpk6dypw5c5g3bx4Oh4MhQ4bw\nhz/8AZ/Px5133klpaSmGYfDb3/6WMWPGWP2tdLtjjdVB8+fPZ+fOnd9Z7fztxxz8STPSdWW83G63\njq3DjNWCBQt48cUXiY2N5YwzzuCXv/xl1B5bXRkrHVeHH6tly5bxl7/8BdM0Oemkk7jrrrsi9rgK\nufIVERGJdCE37SwiIhLpVL4iIiJBpvIVEREJMpWviIhIkKl8RUREgiz8r1QWkR+47777WLNmDR6P\nh+Li4kOXZlx77bVceeWVFqcTEV1qJBLB9u3bx7XXXssnn3xidRQR+RZNO4uIiASZyldERCTIVL4i\nIiJBpvIVEREJMpWviIhIkKl8RUREgkyXGomIiASZznxFRESCTOUrIiISZCpfERGRIFP5ioiIBJnK\nV0REJMhUviIiIkGm8hUREQkyla+IiEiQ/X+b7+4fC3tIJwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe510a8beb8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dist_muestreo = []\n",
    "for i in range(10000):\n",
    "    monedas_sim = np.random.choice([0, 1], size=n)\n",
    "    caras_sim = np.sum(monedas_sim)\n",
    "    T = caras_sim / n\n",
    "    dist_muestreo.append(T)\n",
    "y = np.array(dist_muestreo)\n",
    "\n",
    "T_obs = caras_obs / n\n",
    "valor_p = np.sum(y[y <= T_obs]) / np.sum(y) + np.sum(y[y >= 1-T_obs]) / np.sum(y)\n",
    "\n",
    "kde = stats.gaussian_kde(y)\n",
    "x = np.linspace(min(y), max(y), 150)\n",
    "density = kde.pdf(x)\n",
    "\n",
    "plt.plot(x, density)\n",
    "\n",
    "plt.fill_between(x[x <= T_obs], density[x <= T_obs],\n",
    "                 alpha=0.3, color='C2',\n",
    "                 label='p = {:.3f}'.format(valor_p))\n",
    "plt.fill_between(x[x >= 1-T_obs], density[x >= 1-T_obs],\n",
    "                 alpha=0.3, color='C2')\n",
    "\n",
    "plt.legend()\n",
    "plt.xlabel('T');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Interpretando los resultados\n",
    "\n",
    "Recapitulemos lo que acabamos de hacer es calcular la probabilidad de obtener fracciones de  caras respecto del total de las tiradas (curva azul), este número va entre 0 (ninguna cara observada en las 200 tiradas) y 1 (100% de caras en las 200 tiradas). Esta curva azul se llama _distribución de muestreo (sampling distribution)_. Y es la distribución de nuestro _estadístico_ $T$ para una muestra finita de $n=200$. \n",
    "\n",
    "Esta distribución es central en estadística frecuentista. De hecho a partir de ella hemos calculado la probabilidad el _valor p_, la suma de las dos áreas naranja en las colas de la curva azul. ¿Por qué tenemos dos colas? \n",
    "\n",
    "Esto deriva de la pregunta que estamos intentando contestar: \"Son nuestros datos compatibles con una moneda no sesgada, es decir con una moneda que cae la mitad de las veces cara y la mitad ceca\". Dado que es posible obtener una moneda sesgada de 2 formas (o sesgada hacia las caras o sesgada hacia las cecas) el _valor p_ lo calculamos para uno de estos casos y por simetría obtenemos el _valor p_ para ambos casos. Esto se llama prueba de _dos colas_, si hubiéramos argumentado que solo nos interesa un caso sería una prueba de _una cola_. Ambos análisis son totalmente válidos, cual usar depende la pregunta a contestar\n",
    "\n",
    "\n",
    "¿Cómo se interpreta el _valor p_? \n",
    "\n",
    "Bueno, desde el punto de vista frecuentista las probabilidades viven en _asintópia_ un mundo donde las probabilidades son frecuencias que se obtienen en el límite de infinitas repeticiones de un mismo experimento. Entonces la interpretación sería:\n",
    "\n",
    "    Si arrojáramos 200 veces al aire una moneda NO SESGADA y esto lo repetimos infinitas veces entonces el ~8.5 % de las veces obtendremos 90 caras o menos. Y el ~17% obtendremos un resultado tan sesgado (o más) que el observado tanto para caras como para cecas.\n",
    "\n",
    "Es decir hemos encontrado la respuesta al pregunta ¿Cuán _esperable_ es este resultado para una moneda que _debería_ caer la mitad de las veces cara?\n",
    "\n",
    "Otra pregunta _natural_ sería ¿Es realmente _diferente_ este resultado de $H_{0}$? Para dar respuesta a esta pregunta se suele usar el concepto (o \"confupto\" es decir concepto muy confuso) de significancia estadística, que consiste, como ya adelantamos, en comparar el _valor p_ con el valor predeterminado $\\alpha$:\n",
    "\n",
    "* Si p < 0.05, entonces la diferencia es estadísticamente significativa\n",
    "* Si p < 0.001, entonces la diferencia es estadísticamente altamente-significativa\n",
    "\n",
    "En este caso como p > 0.05 diríamos: \"Se falla en rechazar la hipótesis nula con un p=0.154 y un nivel de significancia de 0.05\".\n",
    "\n",
    "¿De donde salen estos valores predeterminados? \n",
    "\n",
    "La verdad, de una confusión histórico-didactica que llevó a la fusión de dos métodos \"rivales\" de inferencia, junto con un accidente histórico que dejo cuasi-congelado el valor de $\\alpha$ en 0.05. Fuera de ese detalle, la intención de fijar un valor de $\\alpha = 0.05$ (o cualquier otro) es la de controlar el nivel de errores de tipo I. Es decir las falsas alarmas o más formalmente controlar la cantidad de veces que estamos dispuestos a rechazar una hipótesis nula verdadera (es decir una que NO deberíamos haber rechazado). Lamentablemente se ha vuelto muy común que estos límites se utilicen para definir si un resultado deber ser publicado. El límite entonces los fija un criterio editorial (auto-impuesto por los propios científicos) que dependerá de la revista/disciplina, siendo 0.05 común por ej en psicología y en varias ramas de la biología."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Valores p analíticos\n",
    "\n",
    "Los _valores p_, para este problema, pueden ser calculados de forma analítica, basta aplicar la distribución binomial.\n",
    "\n",
    "¿Por qué esta distribución?\n",
    "\n",
    "Precisamente por que modela eventos que tienen dos posibles resultados, cara-ceca, si-no, apagado-prendido, sano-enfermo, etc. La distribución binomial la vimos en el capítulo anterior y tiene la siguiente forma:  \n",
    "\n",
    "$$p(x|n,p) = {\\binom {n}{k}}p^x(1-p)^{n-x}$$\n",
    "\n",
    "\n",
    "Hay que tener cuidado que $p$ en la distribución binomial NO es el _valor p_ del que estábamos hablando! Si no un parámetro de esta distribución. Este es un buen momento para  volver al capítulo anterior y revisar _como es_ esta distribución.\n",
    "\n",
    "Según el paradigma frecuentista el procedimiento sería:\n",
    "\n",
    "1. Vamos a asumir que la moneda tiene la misma posibilidad de caer cara (1) o ceca (0), p = 0.5\n",
    "2. Vamos a lanzar una moneda 200 veces (igual cantidad que en nuestros datos)\n",
    "3. Dados los puntos 1 y 2 y asumiendo que la moneda sigue una distribución binomial (con n=200, p =0.5) vamos a calcular la probabilidad de obtener 90 o menos caras (y luego multiplicar por dos). Esto es el _valor p_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAecAAAFJCAYAAAChG+XKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGiJJREFUeJzt3X9Q3PW97/EX7DckwK4ET5DpnEjGcENnztiRH+2pHsuY\ng0VTa3saV7OQip4jR9vYuR570k7QtjT1RMCrczzTmDATb61ezrSSYTwdsa09RVJzhsEjoBuHmZqe\n4ep6bXstMeyY3c3Jgvs9fzjZHpKwm2QX9r3wfPzl8tldPt/3rD75rl+WAtd1XQEAADMKc70BAAAw\nH3EGAMAY4gwAgDHEGQAAY4gzAADGEGcAAIxxcr2B06anT2T1+crLSzQzE8vqc640zDBzzDA7mGPm\nmGHmsj3DigrfgmvL9szZcTy53kLeY4aZY4bZwRwzxwwzt5QzXLZxBgAgXxFnAACMIc4AABhDnAEA\nMIY4AwBgDHEGAMAY4gwAgDFmPoQEAIClcOrUf+qhh76jmZkZlZSU6Fvf+p7Ky8vPut/MzIx27GjX\nM8/8WKtXr9aBAwf00kuHJEmRSETHj7+v55//hV599RX19n5fa9YU69OfvkZ//dd/m/EeOXMGAKwo\n//IvA9q48X9o//7/rS1bPq9nnvnBWff5938f1d///dd0/Pj7ya/dc889euKJA3riiQO67LLL9O1v\nf0+JREI9Pf+gPXv+l3p7f6B33gnpyJFgxnvkzBkAcNHef/jKrD7fn3xrMuX6z342qH/7t18pFosp\nHA7rb/7mb7V58/XJ9Xff/X/q6fmHeY9pbt6iv/qrW5K333jjiLZvv0OSdPXV1+rpp8+Oc2Fhgf7p\nn/arvb3trLWXXx6Wz+fTn//51ZqZOS6f7xL96Z+ulyR94hNX6Y03grrqqtrzP+hzIM4AgLxy8uRJ\nPf74PoXDM7r77jv1mc9cJ8f5KGfr11+uJ544kPLx0WhUXq9XklRSUqJoNHLWfT71qasXfHxf39Pa\nvfthSdLateU6deo/FQq9rfXrL9fo6Ig2baq52ENLIs7ACtawf1SSNHHvNTneCfJVujPdxVBbW6/C\nwkJdeumfyOe7ROFwWOvWrZN0fmfOpaWlisWikqRYLJYM9fl4663/K6/Xq/XrL5ckFRQU6NvffkiP\nPdatVauKtHFjtcrK1mZ6iMQZAJBfjh59U5J0/Pj7ikaj8y7mOp8z50984iqNjo7oz/7sSr3yyoiu\nuqruvL/3+Piruvrqv5j3tVdfHdU//uMTchxHDz74Td100xcu4GjOjTgDAPLK8ePv6+/+bocikYh2\n7twlj+fC/lrU1q23as+e72rHjnatWrVK3/3uHknSs8/+s9avv1yf+cx1Cz72nXdC+tSnPj3va+vW\nVejuu+/U6tWrdcMNW7RxY/WFH9QZClzXdTN+lizI9t9zrqjwZf05VxpmmDmrMzz9dvaZrL69bXWO\n+WS5zPBnPxtUKPS2duz4n0v+vbM9wxX595wBAMhXvK0NAMgb2fj/ufmAM2cAAIwhzgAAGEOcAQAw\nhjgDAGAMcQYAwBiu1gZWkIV+v/nMdau/7wysFJw5AwBgDHEGAMAY4gwAgDHEGQAAY4gzAADGEGcA\nAIwhzgAAGEOcAQAwhjgDAGAMcQYAwBjiDACAMcQZAABj+MMXwAqQ7g9eLHR//gAGkBucOQMAYAxx\nBgDAGOIMAIAxaf+fcyKR0O7du3X06FEVFRVpz5492rBhQ3J9eHhY+/btk+M48vv92rZtmyRp69at\n8nq9kqT169eru7t7kQ4BAIDlJW2ch4aGFI/H1d/fr2AwqJ6eHvX29kqSZmdn1d3drYGBARUXF6u1\ntVVNTU3y+XxyXVd9fX2LfgAAACw3ad/WnpiYUGNjoySptrZWk5OTybWpqSlVVVWprKxMRUVFamho\n0NjYmN58802dPHlSd911l+644w4Fg8HFOwIAAJaZtGfOkUgk+fa0JHk8Hs3NzclxHEUiEfl8vuRa\naWmpIpGI1qxZo/b2dt122216++23dffdd+vFF1+U4yz87crLS+Q4ngwPZ76KCl/6OyElZpg5CzMs\n9Fzc5SUW9n6apb3kK2aYuaWaYdo4e71eRaPR5O1EIpGM7Jlr0WhUPp9PV1xxhTZs2KCCggJdccUV\nWrt2raanp/Wxj31swe8zMxPL5DjOUlHh0/T0iaw+50rDDDNnZYaJDxMX9TgLe5fszDGfMcPMZXuG\nqUKf9sfp+vp6HT58WJIUDAZVU1OTXKuurlYoFFI4HFY8Htf4+Ljq6uo0MDCgnp4eSdJ7772nSCSi\nioqKTI8DAIAVIe2Zc3Nzs0ZGRtTS0iLXddXV1aXBwUHFYjEFAgF1dHSovb1druvK7/ersrJSt956\nqx544AG1traqoKBAXV1dKd/SBgAAf1Tguq6b601I2X/7jLdwMscMM2dlhhf68Z2nWfn4TitzzGfM\nMHOm3tYGAABLizgDAGAMcQYAwBjiDACAMVxCDSxjF3sh2JmPt3JhGLBScOYMAIAxxBkAAGOIMwAA\nxhBnAACMIc4AABhDnAEAMIY4AwBgDHEGAMAY4gwAgDHEGQAAY4gzAADGEGcAAIwhzgAAGEOcAQAw\nhjgDAGAMcQYAwBjiDACAMcQZAABjiDMAAMY4ud4AgOxr2D+6KM83ce81WX1eAOfGmTMAAMYQZwAA\njCHOAAAYQ5wBADCGOAMAYAxxBgDAGOIMAIAxxBkAAGOIMwAAxhBnAACMIc4AABhDnAEAMIY4AwBg\nDHEGAMAY4gwAgDFp45xIJNTZ2alAIKC2tjaFQqF568PDw/L7/QoEAjp48OC8tffff1/XXXedpqam\nsrtrAACWsbRxHhoaUjweV39/v3bu3Kmenp7k2uzsrLq7u/XUU0+pr69P/f39OnbsWHKts7NTa9as\nWbzdAwCwDKWN88TEhBobGyVJtbW1mpycTK5NTU2pqqpKZWVlKioqUkNDg8bGxiRJjzzyiFpaWnTZ\nZZct0tYBAFienHR3iEQi8nq9ydsej0dzc3NyHEeRSEQ+ny+5Vlpaqkgkoueee06XXnqpGhsbdeDA\ngfPaSHl5iRzHcxGHsLCKCl/6OyElZpi5XMyw0LM4l5Pk8vXAazFzzDBzSzXDtHH2er2KRqPJ24lE\nQo7jnHMtGo3K5/Opr69PBQUFGh0d1a9//Wvt2rVLvb29qqioWPD7zMzEMjmOs1RU+DQ9fSKrz7nS\nMMPM5WqGiQ8Ti/K8uXo98FrMHDPMXLZnmCr0aeNcX1+vQ4cO6aabblIwGFRNTU1yrbq6WqFQSOFw\nWCUlJRofH1d7e7u2bNmSvE9bW5t2796dMswAAOCP0sa5ublZIyMjamlpkeu66urq0uDgoGKxmAKB\ngDo6OtTe3i7XdeX3+1VZWbkU+wZwDg37R5fk+SfuvWZRvw+w0qWNc2FhoR566KF5X6uurk7+c1NT\nk5qamhZ8fF9fXwbbAwBg5eFDSAAAMIY4AwBgDHEGAMAY4gwAgDHEGQAAY4gzAADGEGcAAIwhzgAA\nGEOcAQAwhjgDAGAMcQYAwBjiDACAMcQZAABjiDMAAMYQZwAAjCHOAAAYQ5wBADCGOAMAYAxxBgDA\nGOIMAIAxxBkAAGOcXG8AQOYa9o/m5PtN3HvNkn5fYKXgzBkAAGOIMwAAxhBnAACMIc4AABhDnAEA\nMIY4AwBgDHEGAMAY4gwAgDHEGQAAY4gzAADGEGcAAIwhzgAAGEOcAQAwhjgDAGAMcQYAwBjiDACA\nMcQZAABjiDMAAMakjXMikVBnZ6cCgYDa2toUCoXmrQ8PD8vv9ysQCOjgwYOSpA8//FAPPPCAWlpa\n1Nraqt/85jeLs3sAAJahtHEeGhpSPB5Xf3+/du7cqZ6enuTa7Oysuru79dRTT6mvr0/9/f06duyY\nDh06JEl69tlndf/99+vxxx9fvCMAAGCZcdLdYWJiQo2NjZKk2tpaTU5OJtempqZUVVWlsrIySVJD\nQ4PGxsb0uc99Tps3b5Yk/e53v9Mll1yyCFsHAGB5ShvnSCQir9ebvO3xeDQ3NyfHcRSJROTz+ZJr\npaWlikQiHz2x42jXrl365S9/qe9///uLsHUAAJantHH2er2KRqPJ24lEQo7jnHMtGo3Oi/Ujjzyi\nb3zjG9q2bZt++tOfqqSkZMHvU15eIsfxXNRBLKSiwpf+TkiJGWZuKWZY6MnNtZ1L+frgtZg5Zpi5\npZph2jjX19fr0KFDuummmxQMBlVTU5Ncq66uVigUUjgcVklJicbHx9Xe3q6f/OQneu+99/SVr3xF\nxcXFKigoUGFh6v94zMzEMj+a/6aiwqfp6RNZfc6VhhlmbrFn2LB/dNGe+3xUfe9fJUkT916zqN+H\n12LmmGHmsj3DVKFPG+fm5maNjIyopaVFruuqq6tLg4ODisViCgQC6ujoUHt7u1zXld/vV2VlpW64\n4QY98MAD+vKXv6y5uTk9+OCDWrNmTdYOCACA5azAdV0315uQlPWf6PgpMXPMMHPL/cz5NM6c7WOG\nmVvKM2c+hAQAAGOIMwAAxhBnAACMIc4AABhDnAEAMIY4AwBgDHEGAMAY4gwAgDHEGQAAY4gzAADG\nEGcAAIwhzgAAGEOcAQAwhjgDAGAMcQYAwBjiDACAMcQZAABjiDMAAMYQZwAAjHFyvQEAF65h/2iu\ntzDP6f1M3HtNjncCLA+cOQMAYAxxBgDAGOIMAIAxxBkAAGOIMwAAxhBnAACMIc4AABhDnAEAMIY4\nAwBgDHEGAMAY4gwAgDHEGQAAY4gzAADGEGcAAIwhzgAAGEOcAQAwhjgDAGAMcQYAwBjiDACAMcQZ\nAABjnFxvAMD5m+zY9NE/VP2f3G5kAaf3d2XPf+R4J0B+SxvnRCKh3bt36+jRoyoqKtKePXu0YcOG\n5Prw8LD27dsnx3Hk9/u1bds2zc7O6sEHH9Rvf/tbxeNx7dixQ9dff/2iHggAAMtF2jgPDQ0pHo+r\nv79fwWBQPT096u3tlSTNzs6qu7tbAwMDKi4uVmtrq5qamvTyyy9r7dq1evTRRxUOh/WlL32JOAMA\ncJ7SxnliYkKNjY2SpNraWk1OTibXpqamVFVVpbKyMklSQ0ODxsbGtGXLFt14442SJNd15fF4FmPv\nAAAsS2njHIlE5PV6k7c9Ho/m5ubkOI4ikYh8Pl9yrbS0VJFIRKWlpcnH3nfffbr//vvTbqS8vESO\nk92IV1T40t8JKTHDzC3GDAs9tq/lXIxj5rWYOWaYuaWaYdo4e71eRaPR5O1EIiHHcc65Fo1Gk7H+\n/e9/r6997Wvavn27vvCFL6TdyMxM7II3n0pFhU/T0yey+pwrDTPM3GLNMPFhIuvPmU3ZPmZei5lj\nhpnL9gxThT7tj9/19fU6fPiwJCkYDKqmpia5Vl1drVAopHA4rHg8rvHxcdXV1enYsWO666679M1v\nflO33nprFg4BAICVI+2Zc3Nzs0ZGRtTS0iLXddXV1aXBwUHFYjEFAgF1dHSovb1druvK7/ersrJS\ne/bs0QcffKD9+/dr//79kqQnn3xSa9asWfQDAgAg3xW4ruvmehMSb4NZxAwzl+0Znv494juN/p7z\nM+/cISn7v+fMazFzzDBzpt7WBgAAS4s4AwBgDHEGAMAY4gwAgDHEGQAAY4gzAADGEGcAAIwhzgAA\nGEOcAQAwhjgDAGBM2s/WBpB7pz+2M1+c3m+2P8YTWCk4cwYAwBjiDACAMcQZAABjiDMAAMYQZwAA\njCHOAAAYQ5wBADCGOAMAYAxxBgDAGOIMAIAxxBkAAGOIMwAAxhBnAACMIc4AABhDnAEAMIY4AwBg\njJPrDQBY2GTHplxvISOn939lz3/keCdAfuHMGQAAY4gzAADGEGcAAIwhzgAAGEOcAQAwhjgDAGAM\ncQYAwBjiDACAMcQZAABjiDMAAMYQZwAAjCHOAAAYkzbOiURCnZ2dCgQCamtrUygUmrc+PDwsv9+v\nQCCggwcPzls7cuSI2trasrtjAACWubR/lWpoaEjxeFz9/f0KBoPq6elRb2+vJGl2dlbd3d0aGBhQ\ncXGxWltb1dTUpHXr1unJJ5/U888/r+Li4kU/CGC5yfe/RnUm/joVcGHSnjlPTEyosbFRklRbW6vJ\nycnk2tTUlKqqqlRWVqaioiI1NDRobGxMklRVVaW9e/cu0rYBAFi+0p45RyIReb3e5G2Px6O5uTk5\njqNIJCKfz5dcKy0tVSQSkSTdeOONevfdd897I+XlJXIcz4XsPa2KCl/6OyElZpi5xZhhoSc/LxfJ\nZBa8FjPHDDO3VDNMG2ev16toNJq8nUgk5DjOOdei0ei8WF+ImZnYRT1uIRUVPk1Pn8jqc640zDBz\nizXDxIeJrD/nUrjYWfBazBwzzFy2Z5gq9Gl//K6vr9fhw4clScFgUDU1Ncm16upqhUIhhcNhxeNx\njY+Pq66uLgtbBgBg5Up75tzc3KyRkRG1tLTIdV11dXVpcHBQsVhMgUBAHR0dam9vl+u68vv9qqys\nXIp9AwCwbKWNc2FhoR566KF5X6uurk7+c1NTk5qams752PXr15/161UAACC1/LyqBACAZYw4AwBg\nDHEGAMAY4gwAgDHEGQAAY9JerQ1g6Sy3z9Q+E5+xDZwfzpwBADCGOAMAYAxxBgDAGOIMAIAxxBkA\nAGOIMwAAxhBnAACMIc4AABjDh5AABiz3Dx85Ex9GAqTGmTMAAMYQZwAAjCHOAAAYQ5wBADCGOAMA\nYAxxBgDAGOIMAIAx/J4zkEMr7febz8TvOwPnxpkzAADGEGcAAIwhzgAAGEOcAQAwhjgDAGAMV2sD\nObDSr9I+F67cBv6IM2cAAIwhzgAAGEOcAQAwhjgDAGAMF4QBS4gLwdLjwjCAM2cAAMwhzgAAGEOc\ngUU22bGJt7MvAnPDSkacAQAwhgvCgEXCWV92cIEYViLOnAEAMCbtmXMikdDu3bt19OhRFRUVac+e\nPdqwYUNyfXh4WPv27ZPjOPL7/dq2bVvaxwDLGWfMi4MzaKwkaeM8NDSkeDyu/v5+BYNB9fT0qLe3\nV5I0Ozur7u5uDQwMqLi4WK2trWpqatJrr7224GOA5YooLw0ijZUgbZwnJibU2NgoSaqtrdXk5GRy\nbWpqSlVVVSorK5MkNTQ0aGxsTMFgcMHHAMsFMc6tM+dPrLGcpI1zJBKR1+tN3vZ4PJqbm5PjOIpE\nIvL5fMm10tJSRSKRlI9ZSHl5iRzHc7HHcU4VFb70d0JKzHBhf/mD/78kjzmXd7LyLIvh4o8vW7PB\nwvj3OXNLNcO0cfZ6vYpGo8nbiUQiGdkz16LRqHw+X8rHLGRmJnbBm0+losKn6ekTWX3OlYYZZo4Z\nZgdzzBwzzFy2Z5gq9Gmv1q6vr9fhw4clScFgUDU1Ncm16upqhUIhhcNhxeNxjY+Pq66uLuVjAABA\namnPnJubmzUyMqKWlha5rquuri4NDg4qFospEAioo6ND7e3tcl1Xfr9flZWV53wMAAA4PwWu67q5\n3oSkrL/dwls4mWOGmWOG2cEcM8cMM2fqbW0AALC0iDMAAMYQZwAAjCHOAAAYQ5wBADCGOAMAYAxx\nBgDAGOIMAIAxxBkAAGPMfEIYAAD4CGfOAAAYQ5wBADCGOAMAYAxxBgDAGOIMAIAxxBkAAGPyOs6J\nREKdnZ0KBAJqa2tTKBSatz48PCy/369AIKCDBw/maJf2pZvjCy+8oNtuu00tLS3q7OxUIpHI0U7t\nSjfD077zne/oscceW+Ld5Yd0M3zjjTe0fft2tba26r777tOpU6dytFO70s3w+eef19atW+X3+/Wj\nH/0oR7vMD0eOHFFbW9tZX1+yrrh57Be/+IW7a9cu13Vd9/XXX3e/+tWvJtfi8bj72c9+1g2Hw+6p\nU6fcW265xZ2ens7VVk1LNceTJ0+6119/vRuLxVzXdd2vf/3r7tDQUE72aVmqGZ724x//2N22bZv7\n6KOPLvX28kKqGSYSCfeLX/yi+/bbb7uu67oHDx50p6amcrJPy9K9Dq+99lp3ZmbGPXXqVPK/jzjb\ngQMH3Jtvvtm97bbb5n19KbuS12fOExMTamxslCTV1tZqcnIyuTY1NaWqqiqVlZWpqKhIDQ0NGhsb\ny9VWTUs1x6KiIj377LMqLi6WJM3NzWn16tU52adlqWYoSa+99pqOHDmiQCCQi+3lhVQzfOutt7R2\n7Vo9/fTTuv322xUOh7Vx48ZcbdWsdK/Dj3/84zpx4oTi8bhc11VBQUEutmleVVWV9u7de9bXl7Ir\neR3nSCQir9ebvO3xeDQ3N5dc8/l8ybXS0lJFIpEl32M+SDXHwsJCrVu3TpLU19enWCyma6+9Nif7\ntCzVDP/whz9o37596uzszNX28kKqGc7MzOj111/X7bffrh/+8Id65ZVXNDo6mqutmpVqhpK0adMm\n+f1+ff7zn9fmzZt1ySWX5GKb5t14441yHOesry9lV/I6zl6vV9FoNHk7kUgkB3rmWjQanTdU/FGq\nOZ6+/cgjj2hkZER79+7lp+1zSDXDF198UTMzM7rnnnt04MABvfDCC3ruuedytVWzUs1w7dq12rBh\ng6qrq7Vq1So1NjaedVaI1DN888039atf/UovvfSShoeHdfz4cf385z/P1Vbz0lJ2Ja/jXF9fr8OH\nD0uSgsGgampqkmvV1dUKhUIKh8OKx+MaHx9XXV1drrZqWqo5SlJnZ6dOnTql/fv3J9/exnypZnjH\nHXfoueeeU19fn+655x7dfPPNuuWWW3K1VbNSzfDyyy9XNBpNXuA0Pj6uTZs25WSflqWaoc/n05o1\na7R69Wp5PB5deuml+uCDD3K11by0lF05+7w9jzQ3N2tkZEQtLS1yXVddXV0aHBxULBZTIBBQR0eH\n2tvb5bqu/H6/Kisrc71lk1LN8corr9TAwIA++clP6s4775T0UWyam5tzvGtb0r0WkV66GT788MPa\nuXOnXNdVXV2dNm/enOstm5NuhoFAQNu3b9eqVatUVVWlrVu35nrLeSEXXeGvUgEAYExev60NAMBy\nRJwBADCGOAMAYAxxBgDAGOIMAIAxxBkAAGOIMwAAxhBnAACM+S/+DIkBCkm/jwAAAABJRU5ErkJg\ngg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe50da94898>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "p = 0.5\n",
    "\n",
    "dist_monedas = stats.binom(n, p);\n",
    "x = np.arange(0, n + 1)\n",
    "y = dist_monedas.pmf(x)\n",
    "valor_p = dist_monedas.cdf(90) * 2 # por simetría\n",
    "\n",
    "T = x/n\n",
    "T_obs = caras_obs/n\n",
    "\n",
    "valor_p = np.sum(y[T <= T_obs]) * 2\n",
    "plt.vlines(T, 0, y, 'C0')\n",
    "plt.vlines(T[T <= T_obs], 0, y[T <= T_obs], 'C2',\n",
    "           label='p = {:.3f}'.format(valor_p))\n",
    "plt.vlines(T[T >= 1-T_obs], 0, y[T >= 1-T_obs], 'C2');\n",
    "\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Problemas con las prueba de hipótesis\n",
    "\n",
    "Esta aproximación tiene varios problemas\n",
    "\n",
    "1. Impone una forma de pensar dicotómica, que suele ser inapropiada en la mayoría de los estudios científicos.\n",
    "2. Los *valor p* omiten muchos factores, información previa, nivel de precisión/incerteza de la estimación.\n",
    "3. Un efecto estadísticamente significativo no tiene por que ser un efecto relevante.\n",
    "4. Los _valores p_ dependen de las intenciones del observador! Si el experimento A es tirar $N$ veces una moneda y contar el número de caras y el experimento B es tirar una moneda hasta que obtener $z$ caras. El _valor p_ no necesariamente es el mismo aún si en ambos experimentos observamos que $N=10, z=7$). La razón es que las _distribuciones de muestro_ son diferentes en ambos casos, aún cuando los datos sean idénticos.\n",
    "5. En general los *valores p* son interpretados erroneamente. Un trabajo mostró que solo el 62% de los encuestados fue capaz de contestar de forma correcta sobre la definición de los *valores p*. Lo interesante es que en realidad ninguna de las respuestas provistas en este estudio correspondía a una definición correcta de los _valores p_! Los *valores p* NO son:\n",
    "\n",
    "    * La probabilidad de que $H_{0}$ sea cierta\n",
    "    * La probabilidad de que $\\neg H_{0}$ falsa\n",
    "    * La probabilidad de cometer un error al rechazar $H_{0}$\n",
    "    * La probabilidad de que los datos observados se hayan dado por *azar*\n",
    "    * Una forma de indicar que NO hay efecto si p > 0.5 (o el nivel que sea)\n",
    "    * Una medida de la *relevancia* de un efecto.\n",
    "\n",
    "Gran parte de los problemas derivados de los *valores p*, provienen de usarlos como parte de la maquinaria de \"prueba de hipótesis nula\". La otra parte del problema de usar los _valores p_ proviene quizá de que estos no derivan de un sistema formal de cálculo de probabilidades si no que fueron introducidos de forma totalmente _ad hoc_. Como veremos más adelante la Estadística Bayesiana provee de algo llamado factores de Bayes que curiosamente es todo lo que usted siempre quiso que un _valor p_ fuera, pero que el _valor p_ NO puede ser.\n",
    "\n",
    "Una forma de usar los _valores p_ sin meter la pata es no pedirle que ofrezcan la información que no pueden ofrecer. Un _valor p_ debería ser usado con el fin que originalmente se postuló: Una forma aproximada, cuando se cuenta con poca información sobre un problema, para intentar estimar (de forma objetiva) si vale la pena seguir mirando los datos o haciendo experimentos a fin de replicar o extender resultados. Entonces en general valores \"bajos\" de los _valores p_ indicarían que es posible que tengamos algo interesante. Por ejemplo algunos autores recomiendan de forma _muuuy_ general interpretar los _valores p_ de la siguiente manera \n",
    "\n",
    "* p < 0.01 Hey parece que tenemos algo! :-)\n",
    "* p > 0.1 humm al parecer no hay mucho que decir :-(\n",
    "* 0.01 < p < 0.1 La vida rara vez nos ofrece respuestas claras! :-|"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimaciones puntuales e incerteza\n",
    "\n",
    "En vez de poner a prueba hipótesis una alternativa es realizar inferencias sobre parámetros de un modelo estadístico. Estás estimaciones suelen incluir dos partes:\n",
    "\n",
    "* Estimaciones puntuales\n",
    "* Estimación de la incerteza asociada a la estimación puntual"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimación por máxima verosimilitud\n",
    "\n",
    "Ahora veamos una alternativa a los valores *valores p* y la idea de la prueba de hipótesis, llamada estimación de máxima verosimilitud, o en inglés maximum likelihood estimation (MLE).\n",
    "\n",
    "Para entender este método continuemos con el ejemplo de la moneda. Como dijimos podemos modelar este problema usando la distribución binomial. cuando $n=1$ la distribución binomial se convierte en la distribución de Bernoulli, la cual podemos escribir como:\n",
    "\n",
    "$$f(x | p) = p^x(1-p)^{1-x}$$\n",
    "\n",
    "donde como ya vimos $x \\in [0,1]$.\n",
    "\n",
    "Dado un conjunto de observaciones $x_1, x_2 ... x_n$, las cuales son indepentietes entre si (como las tiradas de una moneda) tendremos que es posible escribir: \n",
    "\n",
    "$$f(x_1, x_2 ... x_n | p) = p^x_1(1-p)^{1-x_1} p^x_2(1-p)^{1-x_2} ... p^x_n(1-p)^{1-x_n}$$\n",
    "\n",
    "Esto que acabamos de escribir es el _likelihood_ (de nuestro modelo) el cual se puede escribir de forma más sintética como:\n",
    "\n",
    "$$\\mathcal{L}(p; x_1, x_2 ... x_n) = \\prod_{i=1}^n p^x_i(1-p)^{1-x_i}$$\n",
    "\n",
    "¿Cómo nos ayuda esta expresión a resolver nuestro problema? Bueno lo primero que hay que notar, y que he omitido hasta este momento, es que desconocemos el valor de $p$, de hecho esto es lo que queremos estimar. Lo cual se resuelve como un problema de optimización, queremos encontrar el valor de $p$ que maximiza el _likelihood_ (de ahí el nombre del método). Analíticamente esto se resulve haciendo\n",
    "\n",
    "$$\\frac{\\partial \\mathcal{L}}{\\partial p} = 0$$\n",
    "\n",
    "Es común usar el símbolo $\\hat{}$ para indicar el valor del parámetro que maximiza el _likelihood_ en nuestro caso sería $\\hat p$, o alternativamente $p_{mle}$ (o incluso ambas notaciones en simultaneo).\n",
    "\n",
    "En la práctica se suele usar el logaritmo del likelihood (log-likelihood), esto facilita los cálculos matemáticos, por ejemplo el producto se convierte en una suma. Esto también provee de ventajas numéricas, básicamente previene el _underflow_ que podría ocurrir al multiplicar sucesivas veces probabilidades (números entre 0 y 1).\n",
    "\n",
    "Resulta que para nuestro modelo, el valor de $p$ que maximiza el likelihood es:\n",
    "\n",
    "$$\\hat p = \\bar x$$\n",
    "\n",
    "donde $\\bar x = \\frac{1}{n}\\sum_i^n x$, es decir la media. Que siguiendo nuestra convención de llamar 0 a las cecas y 1 a las caras es también el número de caras dividido el número total de tiradas.\n",
    "\n",
    "Por lo tanto para nuestro problema $\\hat p$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.45"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "caras_obs / n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Una forma de corroborar que este es el valor de $p$ que maximiza al _likelihood_ (o _log-likelihood_) es usar el método de la grilla. Este método consiste simplemente en evaluar una función en varios puntos de una grilla (en general equidistantes). \n",
    "Entonces lo que vamos a hacer es dibujar el valor del _log-likelihood_ para varios valores de $p$, en todo el rango posible de valores de $p$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgEAAAFdCAYAAACaQr5wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4lOXd9vHvzGSy78lAEkhC2ATCZoiIGkEL7ha1EOBR\nsWgVXEht3Wsf4aWPorQuVbFqqyBSKxpT0Wqr4g4CguwEEAhLAtlJAtmXmXn/iEQpW0ImuWcy5+c4\ncih3ZvnlYsicc60mp9PpRERERLyO2egCRERExBgKASIiIl5KIUBERMRLKQSIiIh4KYUAERERL6UQ\nICIi4qV8jC6gM5WUVLb7MSIiAikvr3FBNd5N7eg6akvXUVu6jtrSddrbljZbyEm/p56ANvLxsRhd\nQpegdnQdtaXrqC1dR23pOh3ZlgoBIiIiXkohQERExEspBIiIiHgphQAREREvpRAgIiLipRQCRERE\nvJRCgIiIiJdSCBAREfFSCgEicpx33/VhzJhAYmODGTMmkHff9arNRUWOkZOzm/HjLyMnZ7fRpbic\nQoCIHOPdd32YMSOA7dst2O0mtm+3MGNGgIKAeK3XX1/Aiy++yuuvLzC6FJczOZ1Op9FFdBZXnB1g\ns4W45HG8ndrRdf67LZ1OJ9UNdiob7NQ02qn+4b81jQ7qmxzU2x3UNTX/f5PDicPppMnR/OUE/prR\nh+L9/sc9T0xSHfe9vB9fixl/HzO+FhP+PmYCrZaWryBfCyG+FgJ9LZhNpk5sBdfQ69J11JbHWr16\nJc8++yQOh4Orr76WqVOnnfS2drudW2+dis3WjT/+8c/YbCGMGXMRgYGBmM0WLBYLr766uNXPfaqz\nAxTtRdxcXZOdwqoGiqsaKKluoKSm+f/Lahspr22iyu6g8HAdh+uaOFLf/OVoT7TPHXjCy4X7fbnv\no52teggTEOJnIcTPh3B/HyL8rYQHNP83KtBKdKAv0UHN/7UFWYkJ9iPc3weTBwYHkdOx2+08/fQ8\nnnnmBbp1686tt95EWtpokpJ6n/D2mZlvkpiYRE1N9THXn3vuZcLDw11am0KAiMHKaxvZV1FL3uE6\ncg/XceBwHQcO15NfWU9BZT2HahtP+xiBVjPh/lZiQ/zoHx1ImJ8PIX4+BP3kE3qA1YyfxYyfT/Mn\neavFjNVswmI24WM2YTGZMJvggfeayMvxPe45eiY18rurB9Bgd1Jvd9DQ1NyjUNtkp6bB0dzr0Gin\nqt7eEkYq65vIPVxHdnH1Cao+lp/FRPdgP2JCfOkR4k+PUL8fvvxJCPMnIdyfYF/9ypLONXv2wzid\nTgoKDlJWVsa99z7E+eentekxtm/PpmfPeHr06AnAuHGXsmLFVycMAcXFRaxa9Q033XQLb731hkt+\nhlPRvyiRTlDbaCenrJbdZTXsLK1mT3kte3/4qqhrOuF9Aq1m4kL8GNQtiJhgP7oH+9ItyBfbD19R\ngVYiA6z0j4+gssJ1R7b+7wN2Zsw4/vojDzi4bnDMGT1mo91BRV0T5bWNlNU2UlLdSGlNA6U1jRRX\nN1BUVU9RVQOFVfWsO3iENc4jJ3ycqAArCeH+JEUEkBQRQO+IAHpHBtIvMpBQf/0682bvvuvDn//s\ny86dZvr3d/Cb3zRw3XUn/rfVFjk5u0hLG8Mf/vA4mzZtZP78p48JAXfeeSs1Ncf/+7vrrrs555xz\nASgpKaZbt+4t37PZurFt29YTPt9zzz3FHXf8+rheAJPJxG9+cydms5lrrvkF11zzi3b/bKAQIOJS\nDqeTveW1bCuuYltJNduKq9leUsX+ijr+u4fe12KiV3gA5/YMIzHcn4SwAOLD/H/48iPUr3Xd4/5W\nC64ceW3+xVnLs8/++Av17rvb9wvVajG3hJfTaXI4KKpq4OCReg4eqePAkXpyD9eRW1HL/oo6sour\n2FBw/E/cPdiX/lGB9IsK5KzoIAbamr/C/K1nXLd4hqOTWY86OpkVatv1uq2vr6eiopybb74NgKSk\nJCorj33t/eUvr5zx4/+3b75ZTnh4JAMGDGT9+u+Oex6brRvl5WX85jd3kZjYi+HDU9r9nAoBImfI\n4XSyp6yW9QVH2FxYyebCKrYUV1HdYD/mdpEBPoyKD6NfVCD9o4LoFx1In4gAeoT6YzG75xj4ddc1\nueRT1JnwMZvpEepPj1B/IOy479sdTg4eqWNPeW3zV1kNuw41fy3fX8Hy/RXH3L5HqB+DbEEM6R7C\n4O7BDO4eTGKYv+YfdCF//vOJw+Wzz/q263W8d28OPXsm4OfnB8D33++gb99+x9ymNT0BNls3iouL\nWr5XUlKMzdbtuPts2bKJb775mtWrv6GhoYHq6ir+8IdHeP75P7fcPiIiktGjL2LbtmyFAJHOVFXf\nxNr8I6w9cJh1+UdYn1/J4foff8GYTdAvKpAh3UNI7hbEoG7BDLIF0S3IV284LmQxm0gIDyAhPICL\nko79XlVDEzlltewoqWZ7SXMvzPaSapbllLEsp6zldmF+PgyLDeHs2BCGx4SQEhdKbIhfJ/8k4io7\nd554tfvJrrfW7t07KSoqpL6+HofDwYIFL3PHHb8+5jat6QkYMGAQeXl55OcfxGbrxqeffsLs2Y8e\nd7vbb5/J7bfPBGD9+u9YsuTvzJr1f9TU1FBTU01gYBC1tbWsXfst06bd2q6f7Si3DAHLli3jo48+\n4qmnnmr587x584iNjQUgIyODkSNHMn/+fL788kt8fHx4+OGHGTp0qJFlSxdTXtvIytwKVuVVsDrv\nMFuLq46Zdd8r3J9xfSNJiQ1lWEwIyd2CCfK1GFewEOzrw7CYEIbFHLskqrSmga1FVWwpqmJrURWb\nCiv5el85X+8rb7lNXIgfab2jGBIdwDk9whjSPRirRVupeIL+/R1s3378v73+/R3tetzdu3cxZszF\nTJ8+Dbu9ialTb2bo0OFtfhwfHx/uued+7rknA4fDzlVXjad37z4t37/vvl/z0EOPEB1tO+H9Dx06\nxB133AE0rzS45JLLGDXq/DP7of6L2+0T8Oijj7JixQoGDhzIM888A8AzzzzDoEGDuOyyy1pul52d\nzbx581i0aBEFBQVkZGSQlZV1ysfWPgHuwx3bsbbRzqq8wyzfX87yfeVsKapqGcf3tZg4OzaUc3uG\ncW7PMFLiQogKPP34dmdwx7b0BBV1jWwqrGRjQSXr8o/w3cEjlNb8uBIj0GpmRFwoo+LDOS8+jBFx\noQRYFfJaqzNfl/89J+Col19u35yAmTOn88ADD5OQ0Ksd1bVfe9vSo/YJSElJYdy4cbz11lst17Kz\ns9m+fTuLFi1i6NCh3Hfffaxbt460tDRMJhNxcXHY7XbKysqIjIw0sHrxNDllNXyWU8bne8tYmVtB\nXVPzJwdfi4nzE8JJSwzngoRwhseG4O+jN4CuJNzfyphekYzp1fw7w+l0UmWx8NHmfNYePMK3Bw4f\nM8fA12LinB5hXJgYTlpiBGfHhqinwE10xGRWgIMHD9CzZ4JrinRThoWAzMxMFi1adMy1uXPncuWV\nV/Ltt98ec/2CCy5g3Lhx9OzZk9mzZ7NkyRKqqqqO2TQhKCiIyspKhQA5pSaHg+8OHuGjXYf4eHcp\nOWW1Ld8baAviZ70jGdMrgpE9wwjUpz6vYjKZ6B0VRPrgGNJ/WApZVtvIt3mHWZlXwcrc5q9vcitg\n+T5C/CykJURwce8ILk6KJDH8+E+i0nk6YjLru+/+26WP544MCwHp6emkp6e36rYTJkwgNDQUgLFj\nx/Lxxx8zYMAAqqt/XEdZXV1NSMjJuzwAIiIC8XHBp7lTda1I63VWOzbaHXy+q5TMTQUs3VrAoR+6\nfAN9LVw7OIarB3XnsrNs9PTgX+J6TbrOT9vSBpyVEMlNP/z5UHUDX+aU8tmuUpbtLOU/u5q/APrb\ngrh6UHeuGtidtKRIfH3US6DXpet0VFu63XDAf3M6nYwfP54lS5YQExPDqlWrSE5OZtiwYfzpT3/i\nV7/6FYWFhTgcjtP2ApSXt39DFY2/ukZHt6Pd4WRFbjnvbivmPztLKf9hQ55uQb7cNDyWy/tFk5YY\n/mMXf2OTx/696jXpOq1py9GxIYyODWHO6CT2VdTyxZ4yvthbxtf7ynn6qz08/dUeQvwsXJwUyRX9\normkT5RXbmSk16XreNWcgP9mMpl49NFHmTlzJv7+/vTp04dJkyZhtVpJTU1l8uTJOBwOZs2aZXSp\nYjCn08nWoioys4t4d3sxRVUNAMQE+3LriB6MH2DjnB5hbrs2XzxPr/AAbk7pwc0pPahvcrAyr4JP\ndx/i492HeH9HCe/vKMFqNpGWGM6V/W1c2T+6VRsmiXQWt1sd0JG0OsB9uLIdD9U08E52Ef/YXMj2\nkuYhonB/H8YPsDExuTsje4Z55Il2raXXpOu4qi2dTiffl9bw750lfLizlC1FVUDzXhIXJIRzzcBu\nXNk/mmg3WWHSEfS6dB2v7gkQORGn08nX+8v5+8YC/rOrlAa7E6vZxFX9o0kf3J2xvaPw05isGMRk\nMjHAFsQAWxD3XNCL3Ipa/r2zlPd2lLSsOHjw452M7hXBxOTuXNE/WocjiSH0qhOPcriukbe2FLFw\nw8GWmf1nRQdy/dBY0gd379KfrMRzJYQHcPvIeG4fGU/e4Tr+taOEpduL+WJvOV/sLSfQaubyftFM\nGhzDmF4RGrKSTqMQIB5h96EaXlqbxzvZRdQ0OvC1mEhP7s60lDhS40K1La94jPgwf+48N547z40n\np6yGd7KLyMou4p/bivnntmLiQvyYPKQ7k4fE0Dsi0OhypYvTnIA20jiXa7SmHZ1OJ6vyDvPimjw+\n3n0IgPhQP355dhzXD4vVp/4f6DXpOka1pdPpZF3+EZZsKeTd7cVU1jcfQnVefBg3DY/j6rNsHje8\npdel62hOgHgVh9PJR7tKeXZVbsuRsSPiQrlzZDxX9o9WV6l0OSaTidQeYaT2COMPY/vy752lvLm5\ngOX7K1iVd5j//XQ3k4fEcNPZseodEJdSCBC3YXc4+eD7Ep5euZ/tJdWYgCv7R3PnyHhG9jz+SFmR\nrijQamFicncmJndnT3kNizcWsGRzIX9Zk8df1uQxtnckt6X24KKkyC696kU6h0KAGM7hdPL+jhL+\ntGIfuw7VYDbBxOTu/Oa8BPpHBxldnohhekcEMvviPjx0YRIf7ixhwfqDfLanjM/2lNEnMoBbR/Rg\n8pAYrSyQM6ZXjhjG6XTy2Z4y5n61l63FVfiYTVw/NIZfn5egLk+Rn/DzMfOLQd35xaDubCqs5JXv\nDvDu9mJ+t2w3T3y9j2kpcdw6ogfdg/2MLlU8jCYGtpEmu7jG7upGfvvuVr49cBgT8IvkbjyQlkRS\nhOfu328UvSZdx5PasqS6gdc25LNw/UFKaxrxtZiYmNydu86Np1+U8T1ontSW7q4jJwYqBLSRXtjt\nc+BwHY9+tYd/bisG4PK+UTw0OolB3YINrsxz6TXpOp7YlrWNdt7eWsSLa/LYU16LCRg/wMZvzk8k\n2cB/V57Ylu5KqwPE41U1NDF/dfPEpromB6nxYcwancSo+PDT31lETirAauGXZ8dx47BY/rOrlD+v\n3M97O0p4b0cJl/WN4t4LEhkeG2p0meKmFAKkQzmdTv69s5Tff7qb/Mp6YoJ9+f2Y3tx5UV8OHaoy\nujyRLsNiNnH1WTau6h/N53vKeHrlfj7+4TCjy/tF8eCFSYb2DIh7UgiQDpNbUcvDy3bzSc4hfC0m\n7jk/kYxRCQT5WjBrrb9IhzCZTIztE8XPekeyYn8F81bs5aNdh/ho1yHGD7DxQFovrbqRFgoB4nJN\nDgcvrjnAkyv2Udvk4MLEcOZd2p++UZrxL9JZTCYTF/aKIC0xnC/2lvHE1/t4f0cJH3xfwg1DY7k/\nrRcxIVpN4O0UAsSlvi+t5tcf7mBDQSXRgVaeuuIsJgzqpr39RQxiMpn4We8oLk6K5KNdh3jsqz0s\n3lTAO9lF3D6yJzPPTSDET28F3kp/8+ISTQ4Hf1lzgD8u30uD3cnE5O48Nq4vEQFWo0sTEZrDwBX9\no7mkbyRLthQyb/k+nlmZy+KNBTw0OokbhsZqS24v5FknUohb2lNew9WLN/Dol3sI97fy+oTB/OXn\nAxUARNyQj9nMjcPiWD39XH43OomaRjv3fbSTS15bx6rcCqPLk06mECDtkrm1kLEL17G+oJJfDOrG\n8lvP4fJ+0UaXJSKnEeRr4bfnJ7J6+rlMHtydrcVVXPOPjdy2NJv8I3VGlyedRCFAzkhVfRN3/Ws7\nd32wAxPw4s8H8tL4Qfr0L+JhYkL8eP7qgXx0Uwoj4kJ4b0cJF7yylpfW5NHkcBhdnnQwhQBpsy1F\nlYx9bR2Z2UWkxIbw+S2pTEjubnRZItIOKXGhfDg1hWevPAs/i4lZn+dwyWvrWHvwsNGlSQdSCJA2\nydxayNWLN7C3vJaMUfH868az6RWu/f5FugKzycT/DI3lm9tGcv3QGLKLq7lq8QYe/GQnVfVNRpcn\nHUAhQFql0e7g98t2cdcHO/CxmFg8YTCPXNQHq0UvIZGuJirQlz9fOYD3bxhO/6hAFq7PZ/Sra/li\nb5nRpYmL6Te4nFZxdQMT3tzE39Yd5KzoQJb9cgSXafKfSJc3Kj6cz25O5bfnJ1BQWc/ktzZz94c7\nOFzXaHRp4iIKAXJK35dWc8Widaw+cJjxA2z856YUekdq5z8Rb+HnY+Z3o3vzyS9HMLhbMG9uKWTM\nq9/x9b5yo0sTF1AIkJNavq+cqxavJ+9IPQ9e2Iu/XTOIYF/tLyXijYbEhPDxL1N4IK0XRVX1TFyy\niUc+201dk93o0qQdFALkhJZsKWTy25upa3Lwl58P5N4LemnrXxEvZ7WYuS+tF/++KYW+kQG8vPYA\nl762ji1FZ37WvRhLIUCO4XQ6eXLFPn794Q6CfS28PXkYE7X8T0R+4uzYUD69OZVbUuLYUVrDFa+v\nZ8H6gzidTqNLkzZSCJAWTqeT2Z/n8McV+0gI8+ffU1M4PyHc6LJExA0FWi08cWl/3kwfQoivDw99\nsotfLc3WpEEPoxAgANgdTu77aCcvrT1A/6hA/nXj2Tr6V0ROa2yfKD6/JZVRPcP44PvS5m3E848Y\nXZa0kluFgMrKSm6//XZuvPFGJk+ezIYNGwDYuHEj6enpTJkyhfnz57fcfv78+UycOJEpU6awefNm\no8r2eI12B3d9sJ3FmwoY0j2YpTcMJ1bnjItIK8WG+PHP64dxz/mJ5B2u4+d/38BfV+03uixpBbea\n6r1w4UJGjRrFtGnT2LNnD/feey/vvvsus2fP5vnnnyc+Pp7p06ezbds2nE4na9asITMzk4KCAjIy\nMsjKyjL6R/A4jXYHt723jX/vLOWcHqH8I30IYf7a/19E2sbHbOah0UmMig9jxnvbmPHOZlYMi2Xu\nJf3w83Grz5vyE24VAqZNm4avry8AdrsdPz8/qqqqaGhoICEhAYC0tDRWrlyJr68vaWlpmEwm4uLi\nsNvtlJWVERkZaeSP4FHsDid3fbCdf+8s5cLEcF6fMIQgX4vRZYmIB7soKZJl00Zw6/vNvYvbSqpY\ncN1g9S66KcNCQGZmJosWLTrm2ty5cxk6dCglJSXcf//9PPzww1RVVREcHNxym6CgIPLy8vDz8yM8\nPPyY65WVlacMARERgfj4tP9NzmYLafdjGM3hcPKrtzexdHsJaUmR/Oe2cwny69yXQ1doR3ehtnQd\ntWX72WwhfJMRyYx3NvP3dQe59PX1LL35HEYlRhhdmsfqqNelYSEgPT2d9PT0465///333HPPPTzw\nwAOMHDmSqqoqqqurW75fXV1NaGgoVqv1uOshIadupPLymnbXbbOFUFLi2WtinU4nDy3bxWvr8zk7\nNoRF1w6i5kgt7W+d1usK7egu1Jauo7Z0HZsthKfG9WVAeACzPt/NRS98w/NXD+Tagd2MLs3jtPd1\neaoA4VYDNbt37+buu+/mqaeeYsyYMQAEBwdjtVrJzc3F6XSyYsUKUlNTSUlJYcWKFTgcDvLz83E4\nHBoKaKX/+3IPC9fnM8gWxJJJQwnp5B4AEfEOJpOJ6ef05I2JQ7BazEx/bxvPrNyv/QTciFv99n/q\nqadoaGjgscceA5oDwIsvvsicOXO47777sNvtpKWlMWzYMABSU1OZPHkyDoeDWbNmGVm6x/jr2gPM\n/zaPvpEBZE4ZRkSAJgGKSMca2yeKD248mxvf2cLjX+8lp6yGp684C1+dQmo4k9OLIpkruvk8ubvw\nPztLmfbPrdiCfPnPTSnEh/kbVosnt6O7UVu6jtrSdU7UlkVV9fwyayvrCyoZ0yuChdclE6yeyNPy\nmuEA6TgbC45wx7+2EWA180b6EEMDgIh4p+7Bfrx7/XAu7xvFV/vKmbBkE6U1DUaX5dUUArxA3uE6\nbnhnC7WNDl4aP4hhMZr9LCLGCLBaWPCLZKYMiWFDQSXj/76BvMN1RpfltRQCurjK+iZuyNxMSXUj\nj47ry+X9oo0uSUS8nI/ZzLNXnsXMc+PZXVbL1X9fz/el1ae/o7icQkAX5nQ6yfhwBztKa7htRA9u\nS+1pdEkiIkDzyoFZF/fh/13ch4LKBq77x0ayi6uMLsvrKAR0YX9Zk8e/d5ZyfnwYc8b2MbocEZHj\n3HluPE9e3p/SmkYmvLmRLUWamNmZFAK6qJW5FTz65R66B/vy8jWD8DHrr1pE3NNNw+P48xVnUV7b\nxIQ3N7GpUEGgs+idoQsqqqrntveyAfjbNYPoHqw9u0XEvV0/LJbnrhrA4bomJry5UccRdxKFgC6m\n0e7g1qXbKKluZPbFfRgVH376O4mIuIHJQ2J44ecDqWqwM/mtzZoj0AkUArqYJ7/Zx7cHDjN+gI0Z\n52gioIh4lonJ3Zt7BOqbmPTWJvaUdeapJt5HIaAL+e7gYZ5dlUtCmD/PXHEWJpPJ6JJERNps0uAY\nHr+kHyXVjUxcsokD2kegwygEdBHVDXZmfrADpxOev2qADgUSEY/2qxE9+P2YJA4cqSf9rU0UV2tn\nwY6gENBF/N+XOewpr+WOkfGcl6B5ACLi+e4+L5GMUfHklNUy5a3NVNU3GV1Sl6MQ0AV8sbeMBevz\nGRAdyEOjexldjoiIy/zvmN5MHR7L1uIqbn1vG412h9EldSkKAR6uoq6R3/x7Bz5mEy9cPRB/H4vR\nJYmIuIzJZGLepf0Y1yeSz/eU8dAnu/Ciw287nEKAh5v1WQ4FlQ3cn9aLIToYSES6IB+zmb9eM4gh\n3YNZvKmA51bnGl1Sl6EQ4MFW51WwZEshQ7oHkzEq3uhyREQ6TLCvD29MHEKPUD8e+2ov/9xWZHRJ\nXYJCgIdqcjh48JNdAMy7tJ+2BRaRLi8mxI9/pA8hxM/Crz/cwTrtKthueufwUK+uO8j2kmpuHBZL\nao8wo8sREekUA23BvHJtMk0OJzf/cytFVfVGl+TRFAI8UGFlPfOW7yPC34ffj0kyuhwRkU51cVIk\nvx/Tm8KqBm5duo0GrRg4YwoBHuj/fZFDVYOd31/Um6hAX6PLERHpdDPPjefagTa+PXCYRz7bbXQ5\nHkshwMMs31fOP7cVkxIbwo3DYo0uR0TEECaTiWeuGMBAWxAL1+fzj00FRpfkkRQCPIjd4eR/P9uN\nCZh3WX/MOhtARLxYkK+FRRMGE+7vwwOf7GRTYaXRJXkchQAPkrWtiO0l1UwZEsMw7QkgIkKv8ABe\nGj+IBruT6e9t09bCbaQQ4CEa7A7+uHwfvhYT96f1MrocERG38bPekcw8N5695bU88MPSaWkdhQAP\n8fdNBeQermPa2XH0DPM3uhwREbfyu9FJpMSG8E52EW9tKTS6HI+hEOABahrtPP3NfgKtZu4+L9Ho\nckRE3I7VYualawYR4mfhwU92klNWY3RJHkEhwAO8uu4gxdUN3H5OPLYgLQkUETmRXuEBPHlZf2oa\nHUx/bxv1Tdo/4HQUAtzc4bpGnl+dS7i/D3eM7Gl0OSIibu26Qd25fmgMW4qqmLd8r9HluD2FADf3\n4poDVNQ1MXNUAmH+VqPLERFxe4+N60diuD9/WZOn8wVOw61CQGVlJbfffjs33ngjkydPZsOGDQAs\nW7aMcePGMXXqVKZOncqaNWsAmD9/PhMnTmTKlCls3rzZyNI7REVdIy+tzaNbkC+3juhhdDkiIh4h\nyNfCs1cOwOGEuz/cQV2T3eiS3JaP0QX81MKFCxk1ahTTpk1jz5493Hvvvbz77rts3bqV+++/n8su\nu6zlttnZ2axZs4bMzEwKCgrIyMggKyvLwOpd7/UN+dQ0Org/rSeBVovR5YiIeIzzE8L5VUoPXl1/\nkCdX7Od/L+ptdEluya16AqZNm8aUKVMAsNvt+Pn5Ac1v+FlZWVx//fU88cQTNDU1sW7dOtLS0jCZ\nTMTFxWG32ykrKzOyfJdqsDt4Zd1Bgn0tTB0WZ3Q5IiIe5/cXJZEQ5s/8b3PZUKBhgRMxrCcgMzOT\nRYsWHXNt7ty5DB06lJKSEu6//34efvhhAC644ALGjRtHz549mT17NkuWLKGqqorw8PCW+wYFBVFZ\nWUlkZORJnzMiIhAfn/Z/orbZOn63vsXf5VFY1cBvR/emT3xEhz+fETqjHb2F2tJ11JauY3Rb2oCF\n/3M2Y19axW8/2sn6e0bj54L3ACN0VFsaFgLS09NJT08/7vr333/PPffcwwMPPMDIkSMBmDBhAqGh\noQCMHTuWjz/+mAEDBlBdXd1yv+rqakJCTt1I5eXtXzdqs4VQUtKx+1M7nU7mfbYLswluGGTr8Ocz\nQme0o7dQW7qO2tJ13KUth4T7Me3sOF7bkM/D72Xz0GjPO369vW15qgDhVsMBu3fv5u677+app55i\nzJgxQPMb4vjx4yksbN4BatWqVSQnJ5OSksKKFStwOBzk5+fjcDhO2QvgSVbsryC7uJqfn2UjITzA\n6HJERDzarIt6Exfix/xvc9mjTYSO4VYTA5966ikaGhp47LHHAAgODubFF1/k0UcfZebMmfj7+9On\nTx8mTZqE1WolNTWVyZMn43A4mDVrlsHVu86La/MAuH1kvMGViIh4vmA/H/5vbB9+tXQbv1u2iyWT\nhmLSKayUihTwAAAgAElEQVQAmJxOp9PoIjqLK7qmOrqLa2dpNWmvrGVkz1A+uDGlw57HaO7SVdgV\nqC1dR23pOu7Wlk6nk0lvbearfeUsvC6Zq86yGV1Sq3nNcIDAy2sPAHD7OeoFEBFxFZPJxOOX9MNq\nNvHIZ7upadTeAaAQ4FYO1TTw9tZCEsP9uaJftNHliIh0KX2jArljZDwHjtTz7Kr9RpfjFhQC3EhW\ndjH1die/SumBxazxKhERV/vt+YnEhfjxwrd5miSIQoBbWbKlEB+ziQnJ3Y0uRUSkSwrytfB/Y/vQ\nYHfyv5/tNrocwykEuIns4iq2FlcxrnekjgsWEelAV59lIy0hnE9zyvhmf7nR5RhKIcBNvLWleR+E\nSUNiDK5ERKRrM5lMPHJx81kC//flHrxokdxxFALcQKPdwTvZRUQG+HBp3yijyxER6fLOjg3l52fZ\nWF9QyYc7S40uxzAKAW7gi71llNY0ct3A7vha9FciItIZHh6ThMUEc7/aQ5PDYXQ5htA7jht4a0sR\nAFOGaihARKSz9IkM5IZhsewuq+XNzYVGl2MIhQCDldc28vHuUgZEBzK0e7DR5YiIeJX7LuhFgI+Z\nP63Y55UbCCkEGOzd7cU02J1MGhKjvaxFRDpZTIgf08/pSWFVA6+sO2h0OZ1OIcBgb20pxGyCdO0N\nICJiiJnnxhPh78Pzq3KprG8yupxOpRBgoJ2l1WwoqOTipEi6B/sZXY6IiFcK87dy57nxHK5vYuGG\nfKPL6VQKAQZ6f0cJAOmD1QsgImKkm8/uQYifhZfW5FHrRXMDFAIM9MnuQ/iYTYzrrb0BRESMFOrv\nwy0pPSitaeTNLd6zUsDnTO5UX1/P119/zffff09TUxM9evTg7LPPpm/fvq6ur8sqrKxnY2ElFyaG\nE+p/Rn8NIiLiQtNTe/Ly2gO8sDqXqcNisXrBvi1tfvdZsGABL774IjU1NQQHB+Pn50dlZSW1tbWc\nd955zJkzh4SEhI6otUtZlnMIQDsEioi4CVuQL9cPjWHB+nze3V7MpMFdf++WVoeAuro6/vCHP5CQ\nkMDixYs566yzjlnSduTIEVatWsUTTzzBrbfeSkpKSocU3FUs2300BEQbXImIiBx117kJvL6xgOdW\n5TIxuTvmLr50u9Uh4B//+Ae//vWviYk5cTIKDQ3lsssu49JLL+WVV16hV69eREZGuqzQrqS20c5X\n+8rpHxVIUkSA0eWIiMgP4sP8+cWgbry9tYiPdpVyZX+b0SV1qFYNeJSVlXH55ZefNAD8lMlk4rbb\nbqOgoKDdxXVVK/ZXUNvk4BINBYiIuJ1fj0rABDy7KrfLnzDYqhAQGRlJXFwcTzzxBI2NjSe8zaFD\nh5gxY0bLn5OTk11TYRf08e7mE6suUwgQEXE7/aODuKJ/NBsKKlmdd9jocjpUm6Y+vvbaa0yePJl9\n+/Ydc3358uVcc801fP31166srUtyOp0syzlEhL8PqT1CjS5HREROYEZqTwAWbujaWwm3ef3D9u3b\nue6663jnnXdoaGjgscceY8aMGZSWlmI2d/3lFO21taiKgsoGxvaJwkftJSLilkbFhzEgOpAPvi+l\nqKre6HI6TJvehR5//HGCg4Opra3lkUceYcyYMfz973/H4XAQHx/P4sWLO6rOLuPjH1YFaChARMR9\nmUwmbk7pQZPDyRubuu4ctzaFgOuuu44PPviA1NRUnE4n5eXlAFxzzTW89957WhbYCkd3Cbw4SSsn\nRETcWXpyd4J8Lby+sYAmh8PocjpEm/ujV69ezc6dO4/ZI2D9+vXs2LHDpYV1RUd3CTwvPky7BIqI\nuLlgPx8mDe5OfmU9n/zQi9vVtCkE3HXXXTz00EMcOXKEgIAArrvuOsxmMwcOHGDq1Kk8+eSTHVVn\nl/DpHu0SKCLiSaadHQfAwvVd83TBNoWAzz77DKfTydChQ3nvvfd4/PHHWbRoEbGxsdjtdl599dWO\nqrNLWL6vAoCxOjBIRMQjDLQFc158GF/tKyenrMboclyuTSHAbDZz++238+abbxIfHw9Aamoq7733\nHpdeemmHFNiVrDl4mOhAK30itUugiIinuDmlBwCvbeh6vQFtGphetGgR55xzznHXQ0JCeO6558jM\nzGxXMTU1Ndx7770cOXIEq9XKvHnz6N69Oxs3buSxxx7DYrGQlpbGzJkzAZg/fz5ffvklPj4+PPzw\nwwwdOrRdz9+RDh6p4+CRei7vF3XMfAoREXFvV/aPxhZkZcnmQn43OolAq8XoklymTT0BJwoAP5We\nnt6uYt5++22Sk5N54403GD9+PH/7298AmD17Nk899RRvvvkmmzZtYtu2bWRnZ7NmzRoyMzN5+umn\nmTNnTrueu6OtPXgEgJE9wwyuRERE2sLXYmbqsDgO1zfx/o4So8txKbeaoj5t2jTsdjsA+fn5hIaG\nUlVVRUNDQ8vxxGlpaaxcuRJfX1/S0tIwmUzExcVht9spKytz20OL1hxo3npyZA+FABERTzN5SAxP\nr9xPVnYRU4Z0nSOGDQsBmZmZLFq06Jhrc+fOZejQodx0003s3LmThQsXUlVVRXBwcMttgoKCyMvL\nw8/Pj/Dw8GOuV1ZWnjIEREQE4uPT/m4cmy2kzfdZX1SFr8XMuCGx+Lmghq7gTNpRTkxt6TpqS9fp\nSm1ps4UwKjGC5fvLsftZiQn17/Tn7wguCwF1dXW88847xMbGMnbs2NPePj09/aTDB6+//jo5OTnM\nmDGDpUuXUl1d3fK96upqQkNDsVqtx10PCTl1I5WXt39mp80WQklJZZvuU9XQxKaDh0mJC+WIC2ro\nCs6kHeXE1Jauo7Z0na7Yltf0j2L1/nJeWbGHGefEd9rztrctTxUgzmjz+muvvZYXXnjhmGsLFizg\nsccew8/Pj9dee+1MHpaXX36ZpUuXAs2f7C0WC8HBwVitVnJzm490XLFiBampqaSkpLBixQocDgf5\n+fk4HA63HQrYkF+J3an5ACIinmz8gG5YTJCVXWx0KS5zRj0BI0aMICkp6Zhr55xzDjfccANpaWn0\n6NEDp9PZ5lnwEyZM4MEHHyQrKwu73c7cuXMBmDNnDvfddx92u520tDSGDRsGNC9PnDx5Mg6Hg1mz\nZp3Jj9Ip1hzUfAAREU9nC/LloqRIPttTRk5ZDX0iA40uqd1MTqfT2dY71dXVsWTJEg4cOEC/fv24\n+uqrCQoK6oj6XMoVXVNn0i0z+a1NfLG3nOyM87EF+ba7hq6gK3YVGkVt6TpqS9fpqm35TnYRd/5r\nO/dekMiDFyad/g4u4HbDAQ888ABPPvkk27Zt4+WXX+aKK65g69atZ1xgV+ZwOvku/wi9IwIUAERE\nPNzl/aIItJrJyi7iDD5Du50zCgHLly9n3rx5/OMf/+Dzzz/nt7/9LTNmzODAgQOurs/j7SipprLe\nrvkAIiJdQLCvD5f3i2ZfRR3rCzy/p+OMQkC3bt3o169fy5+vu+46/ud//oenn37aZYV1FUfnA5zT\nI9TgSkRExBUmDOoOQFZ2kcGVtN8Zrw548803j7k2duxYli9f7pKiupK1B7RToIhIV3JRUgRRAVaW\nbi+myeEwupx2OaPVAUuXLiU3Nxez2cz48eMZOnQo33//PX5+fq6uz+OtOXiYcH8f+kV5/ixSEREB\nq8XM+IE2Fq7P5+t9Ffyst3suT2+NMwoBzz77LJ999hmrVq3ihhtuICgoiOrqakaNGsWePXvo3bu3\nq+v0SEVV9eyvqGNcn0jMOjRIRKTL+MWg7ixcn8/7O4q9LwQMGDCAAQMGcNddd1FbW8uaNWtYvXo1\nq1ev5qqrriIqKorzzz+fP/7xj66u16OsOToUoP0BRES6lHN6hBIdaGVZziEcTqfHftBr97bBAQEB\njBkzhjFjxgBQXl7OqlWrWL16NQcOHKBnz57tLtJTrT26SZDmA4iIdClmk4lxfaJYsqWQTYWVnB3r\nmZO/WzUxsKKigm+++aZVDxgREcGFF17IxRdf7NUBAGBDQSVmEwyP7TqHaIiISLNL+kQB8MnuQwZX\ncuZaFQLCw8MpLS3lnXfeOe1t9+3bx9y5cxk9enS7i/NkTqeTnaXV9AoPINCqUwNFRLqai5MisJpN\nLPPgENDq4YBrrrmGhQsXMnnyZC6++GIGDhxIWFgYFouFw4cPk5eXx6pVq6ipqeFPf/oTFot3v/GV\n1jRSXtfEuRoKEBHpkoL9fDg/IZyv9pVTUFlPbIjnrZBr05yAm2++mSFDhvDKK6/w7LPPtlx3Op0k\nJiZy4403cv3113t9AADYdaj5yOB+0VoaKCLSVV3aN4qv9pXzac4hpg6PM7qcNmvzxMDU1FRSU1Op\nq6vjwIEDVFZWEhsbS0xMTEfU57F2HqoGoF+U+x+sJCIiZ2Zcnyh+/+lulu32khBwlL+/P3379gWa\nVwSUlZURGem5ayVdbVdpc09Af20SJCLSZSVFBNA/KpCv95dT22gnwMPmgLVqYqDT6eSbb77h22+/\nPeH3w8LCeO2113B4+PaJrrTzkEKAiIg3uKRvFDWNDlbmVhhdSpu1KgQ88cQT3HbbbUybNo2vvvrq\n+Acxm/nZz3523HkC3mzXoRriQvwI9mv3VgwiIuLGLj26VDDH81YJtCoE7N69m2XLlnHttdfSq1ev\nE95m0KBBvPXWW66szWNV1jeRX1mv8wJERLzAOT1DCfPzYdnuQzidTqPLaZNWhYARI0Zgt9t5/PHH\nSUxMPOFttmzZQk5ODna73aUFeqJdGgoQEfEaPmYzY/tEcuBIPdtLqo0up01aFQKmTp3KAw88wJo1\na056mzfffBOTyYTZfEanE3cpO7U8UETEq1zat3lIYJmHDQm0asA6JCSEu+66ixkzZhAUFERKSgqp\nqalccMEF9OnTB4CioiLi4uIweeghCq6064flgf21PFBExCv8rHckZhN8llPG3eeduMfcHbV61tqF\nF17IP//5T5555hm++OILPvnkE0wmE7GxsUyaNInIyEgCA/XJF35cHqg5ASIi3iHc38rgbsFsKDhC\nXZMdfx/PWCrYpqnrSUlJPPfcc1RXV7Nu3Tq+++471q5dywsvvIDdbueiiy7qoDI9y85DNUT4+xAd\naDW6FBER6STnxYezuaiKjQWVjIoPN7qcVjmj9WtBQUGMHj265ZCg+vp6vvnmG60OAOqbHOyrqCW1\nR6iGRkREvMi58WG8/N0BVuUd7toh4L/5+fnxs5/9zOuPDgbYU16DwwlnaT6AiIhXGRXffGDc6rwK\nwDPmBbRpKv9rr73G2rVrqa4+8RKI/v37u6QoT7ZT8wFERLxSdKAv/aICWXPwCE0esoNum3oCnn76\naRobGzGZTMTHxzN48GAGDRpEcnIyycnJhISEdFSdHqNlu2AtDxQR8Tqj4sNYvLGA7OJqhsW4/3ti\nm0LAhx9+yM0330xYWBiJiYnk5+fz+eefU1dXh8lkYsCAAUyaNIlJkyZ57X4Bu3R6oIiI1zovPpzF\nGwtYlVvR9ULA7NmzueSSS3jwwQdbrlVUVPDCCy/wr3/9i7PPPptnnnmGjz/+mL/97W/4+Hjfvvk7\nS2sItJrpEepndCkiItLJRvVsnhewKu8wt4+MN7ia02vTx/V169a1rAg4Kjw8nN///vdceumlhIWF\nkZWVxa5du3j11VddWqgnsDuc5JTV0DcyELNWBoiIeJ2eYf7Eh/rx7YEKjzhHoE0hICoqiu3bt5/w\ne1dccQVLly6lZ8+eTJs2jX/9619tLqampoY77riDG264gWnTplFUVATAsmXLGDduHFOnTmXq1Kkt\n2xfPnz+fiRMnMmXKFDZv3tzm53O13MN11Nud2i5YRMSLnRsfTlltU8scMXfWphBwzTXX8OKLL7Jt\n27bjvldUVER5eTkAycnJ5ObmtrmYt99+m+TkZN544w3Gjx/P3/72NwC2bt3K/fffz+LFi1m8eDEj\nR44kOzubNWvWkJmZydNPP82cOXPa/Hyupu2CRUTkvPijQwIVBldyem0atL/zzjvZsWMHkyZN4sor\nr2TMmDFER0ezd+9e5s+fz+DBgwFoaGg4o/kA06ZNazmFMD8/n9DQUACys7PZvn07ixYtYujQodx3\n332sW7eOtLQ0TCYTcXFx2O12ysrKiIyMbPPzukrLwUFaHigi4rV+3C/gMNPO7mFwNafWpndqq9XK\niy++yFtvvcWrr77K+++/j8lkwul00q9fP/7whz8AsHHjRhISEk75WJmZmSxatOiYa3PnzmXo0KHc\ndNNN7Ny5k4ULFwJwwQUXMG7cOHr27Mns2bNZsmQJVVVVhIf/uCNTUFAQlZWVpwwBERGB+LhgP2eb\n7cQzPnOrGgE4r5/tpLeRH6mNXEdt6TpqS9fx1raMjg7GFuzLmoNHiI4OdsnusR3VliZnO2Yu5OXl\nUVxcTGRkJL169Wr5QTdu3EhtbS3nnXfeGReWk5PDjBkz+PTTTzly5EhLr8BXX33Fxx9/zIABA6iv\nr+e2224D4Nprr2XBggWnDAElJZVnXM9RNlvISR/nitfXsamwiv33XojV4p1LJFvrVO0obaO2dB21\npet4e1ve8u5WPvi+lO9uP5eE8IB2PVZ72/JUAeKM36kaGhqIj49nxIgRJCUlHZN0hg8ffkYB4OWX\nX2bp0qVA8yd7i8WC0+lk/PjxFBYWArBq1SqSk5NJSUlhxYoVOBwO8vPzcTgchg4FAOwpq6VXuL8C\ngIiIlxvVs7mnelXeYYMrObU2DQc0NTXx/PPP8/bbb1NRUUFYWBjnnXceU6dOJSUlpd3FTJgwgQcf\nfJCsrCzsdjtz587FZDLx6KOPMnPmTPz9/enTpw+TJk3CarWSmprK5MmTcTgczJo1q93P3x71TQ7K\n65oY3D3Y0DpERMR45/3kHIHJQ2IMrubk2hQCXnrpJRYuXMi1115LYmIiBQUFfPvtt9xwww389re/\nZfr06e0qJjo6+oT7C6SlpZGWlnbc9YyMDDIyMtr1nK5SWtMAgC3I1+BKRETEaIO6BRPiZ2H1gS7U\nE7B06VLuuecepk2bdsz1d955hzlz5tC7d2/GjRvnyvo8RnF1cwjophAgIuL1LGYTw2NCWL6/gqr6\nJoL93HMH3TYNXpeUlDBo0KDjrk+cOJFf/vKXLFiwwGWFeZriKvUEiIjIj5K7NQ8PZ5ec+ORdd9Cm\nEJCYmHjSHQMvvPBCduzY4ZKiPJF6AkRE5KdaQkBxlcGVnFybQsDEiRP5y1/+csIgUFhYSGCg926S\nU3I0BAQrBIiIiGeEgDYNUkydOpW1a9cyadIkrrjiCkaPHk10dDS5ubnMnz+fCy+8sKPqdHvqCRAR\nkZ/qHx2I1Wwiu6iLhACTycSzzz7La6+9xuuvv87777/f8r0LLriAhx56yOUFeoqjIUBzAkREBMDX\nYqZ/dCDbS6qxO5xYzO53uuxpQ8DKlStJTk4mLKx5zaPZbOaWW27hlltuYc+ePZSXl9OjRw9iYtx3\nHWRnKKluxGyCqACr0aWIiIibSO4WTHZxNXvLa+nrhufKnDYE3HLLLZhMJmJjYxk8eDDJyckMGjSI\n5ORkevfu3Rk1eoTi6gaiA33dMumJiIgxBncL5m2K2Fpc5Zkh4MMPPyQ7O5vs7Gy2bdvGX//6V6qr\nqzGZTMTExLQEgqNf0dHRnVG32ymubqBXuL/RZYiIiBv56eTAawd2M7ia4502BPTp04c+ffowfvz4\nlmt79+5l27ZtbN26lezsbBYuXEhlZSUmk+mkSwi7suoGO9UNdk0KFBGRYyR3d+8VAme0hVFSUhJJ\nSUlcddVVLdf2799Pdna2ywrzJCU1WhkgIiLHiwywEhfi17VCwIkkJiaSmJjoqofzKNotUERETia5\nWxDLcso4VNNAVKB7vU/ozFsX0B4BIiJyMj/OC3C/7YMVAlxAuwWKiMjJDHbjeQEKAS6gngARETmZ\noz0BW91w50CFABfQboEiInIyvcIDCLSa1RPQVZWoJ0BERE7CYjYx0BbMrkM1NNgdRpdzDIUAFyip\nbsBqNhHu77LFFiIi0oUkdwui0eHk+1L3mhyoEOACxVUN2IJ8MZm0ZbCIiBzPXVcIKAS0k9PppLi6\nQUMBIiJyUu66QkAhoJ0q6+3U250KASIiclIDbUGYgGw3WyGgENBOxdojQERETiPY14deEQFkF1fh\ndDqNLqeFQkA7/bg80GpwJSIi4s6SuwVRXtdE0Q9bzbsDhYB20vJAERFpjaSIAAD2VdQaXMmPFALa\nSbsFiohIaySENYeA3MN1BlfyI4WAdtJugSIi0hoJ4f4A5FYoBHQZGg4QEZHWSPwhBOzXcEDXoeEA\nERFpjZ6h/pjQcECXUlzVQKDVTJCvxehSRETEjflazMSF+ikEnE5OTg4jRoygvr4egI0bN5Kens6U\nKVOYP39+y+3mz5/PxIkTmTJlCps3bzak1uLqBqIDtWWwiIicXkKYP/lH6t3mICG3O/GmqqqKefPm\n4ev7Y/f67Nmzef7554mPj2f69Ols27YNp9PJmjVryMzMpKCggIyMDLKysjq1VofTSWlNI8NjQzr1\neUVExDMlhgewKu8wBw7X0Tsy0Ohy3KsnwOl08sgjj3DPPfcQENC8lKKqqoqGhgYSEhIwmUykpaWx\ncuVK1q1bR1paGiaTibi4OOx2O2VlZZ1ab3ltI00ObRksIiKtkxD2w+RANxkSMKwnIDMzk0WLFh1z\nLS4ujiuvvJIBAwa0XKuqqiI4OLjlz0FBQeTl5eHn50d4ePgx1ysrK4mMjDzpc0ZEBOLj0/6xe5ut\n+ZN/UcERABKiglquSeupzVxHbek6akvXUVseb3B88/tWud3ZpvbpqLY0LASkp6eTnp5+zLVLLrmE\nrKwssrKyKCkp4ZZbbuHll1+muvrHoxerq6sJDQ3FarUedz0k5NSNVF5e0+66bbYQSkoqAfj+QAUA\nIRZTyzVpnZ+2o7SP2tJ11Jauo7Y8sfAf+t+35lVQ0rd17dPetjxVgHCr4YBly5axePFiFi9ejM1m\nY8GCBQQHB2O1WsnNzcXpdLJixQpSU1NJSUlhxYoVOBwO8vPzcTgcp+wF6AhaHigiIm3RK9y9dg10\nu4mBJzJnzhzuu+8+7HY7aWlpDBs2DIDU1FQmT56Mw+Fg1qxZnV5XcZVCgIiItF63YF/8LCa32TXQ\nbUPA559/3vL/w4cP5+233z7uNhkZGWRkZHRmWcfQCYIiItIWZpOJ+DB/cg+7x66BbjUc4Gm0ZbCI\niLRVQngAZbVNVNY3GV2KQkB76PAgERFpqx/PEDB+SEAhoB2KqxsI8bMQYNWWwSIi0jpH9wpwhyEB\nhYB2KKlu0FCAiIi0SUJY8woB9QR4sCaHg0M1jQoBIiLSJr1+GA5whxUCCgFnqKbBgROIDfEzuhQR\nEfEgCeHuMxzgtksE3V2ovw9/vWYQw2K0LaaIiLRemL+VMD8ft9gwSCGgHa4d2M3oEkRExAMlhPuz\n+1ANTqfT0KPoNRwgIiLSyRLD/altcrQsNTeKQoCIiEgnO7pCwOghAYUAERGRTpbgJisEFAJEREQ6\nWWLY0V0DjV0hoBAgIiLSyRLd5EhhhQAREZFO1jOseY8ZDQeIiIh4GX8fCzHBvhoOEBER8UYJ4f4c\nrKyn0e4wrAaFABEREQMkhgXgcMLBynrDalAIEBERMYA7LBNUCBARETGA7YdTaA/VGLdroEKAiIiI\nASIDrACU1TYZVoNCgIiIiAF+DAGNhtWgECAiImKAiIDmg3zLFQJERES8i3oCREREvFSEQoCIiIh3\nCrRaCPAxazhARETEG0UEWCnX6gARERHvExHgo+EAERERbxQZYKWqwU6DQecHKASIiIgY5OgKAaPm\nBbhlCMjJyWHEiBHU1zcfqrBs2TLGjRvH1KlTmTp1KmvWrAFg/vz5TJw4kSlTprB582YjSxYREWkz\no1cI+BjyrKdQVVXFvHnz8PX1bbm2detW7r//fi677LKWa9nZ2axZs4bMzEwKCgrIyMggKyvLiJJF\nRETOSMteATXqCcDpdPLII49wzz33EBAQ0HI9OzubrKwsrr/+ep544gmamppYt24daWlpmEwm4uLi\nsNvtlJWVGVi9iIhI2xh9foBhPQGZmZksWrTomGtxcXFceeWVDBgw4JjrF1xwAePGjaNnz57Mnj2b\nJUuWUFVVRXh4eMttgoKCqKysJDIy8qTPGRERiI+Ppd2122wh7X4MUTu6ktrSddSWrqO2PL3E7s1t\n1ORjOWV7dVRbGhYC0tPTSU9PP+baJZdcQlZWFllZWZSUlHDLLbfwxhtvMGHCBEJDQwEYO3YsH3/8\nMQMGDKC6urrlvtXV1YSEnLqRystr2l23zRZCSUllux/H26kdXUdt6TpqS9dRW7aOtam5ByC3pPKk\n7dXetjxVgHCr4YBly5axePFiFi9ejM1mY8GCBTidTsaPH09hYSEAq1atIjk5mZSUFFasWIHD4SA/\nPx+Hw3HKXgARERF3c3Ri4CGD5gS43cTA/2YymXj00UeZOXMm/v7+9OnTh0mTJmG1WklNTWXy5Mk4\nHA5mzZpldKkiIiJtEmHwEkG3DQGff/55y/+npaWRlpZ23G0yMjLIyMjozLJERERcJtL/aAgwZmKg\nWw0HiIiIeJMQPws+ZhOHtFmQiIiIdzGZTEQE+GjHQBEREW8UGWBVCBAREfFGkQFWKuqasDucnf7c\nCgEiIiIGigiw4gQq6jq/N0AhQERExEA/niTY+SsEFAJEREQMFGngSYIKASIiIgaKCGjeskchQERE\nxMtEGrhroEKAiIiIgSI0HCAiIuKdNDFQRETES2lioIiIiJfSxEAREREvFe5vxYQmBoqIiHgdi9lE\nuL+PegJERES8UUSAVSFARETEG0UEWCmvbcLp7NxDhBQCREREDBYVYKXJ4aSqwd6pz6sQICIiYrCj\nKwQOdfKQgEKAiIiIwSIM2jpYIUBERMRgUYEKASIiIl7px/MDOnfrYIUAERERg0X4/xACatQTICIi\n4oWzaE4AAAlJSURBVFWODgd09l4BCgEiIiIGO7o6oLxOIUBERMSrtJwkqOEAERER7xJh0HHCCgEi\nIiIG87WYCfa1UN7JqwN8OvXZTsPpdDJ69Gh69eoFwPDhw7n33nvZuHEjjz32GBaLhbS0NGbOnAnA\n/Pnz+fLLL/Hx8eHhhx9m6NChBlYvIiJy5iIDrJ0+J8CtQkBubi7Jycm89NJLx1yfPXs2zz//PPHx\n8UyfPp1t27bhdDpZs2YNmZmZFBQUkJGRQVZWlkGVi4iItE9EgA87S2s69TndKgRkZ2dTVFTE1KlT\n8ff353e/+x3dunWjoaGBhIQEANLS0li5ciW+vr6kpaVhMpmIi4vDbrdTVlZGZGSkwT+FiIhI20UG\nWKltclDbaCfAaumU5zQsBGRmZrJo0aJjrs2aNYvp06dzxRVX8N1333H//ffzwgsvEBwc3HKboKAg\n8vLy8PPzIzw8/JjrlZWVCgEiIuKRfnp+QJcPAenp6aSnpx9zrba2Foul+QdPTU2luLiYoKAgqqur\nW25TXV1NaGgoVqv1uOshISGnfM6IiEB8fNrfsDbbqZ9HWkft6DpqS9dRW7qO2rJtekQGAeD09z2u\n7TqqLd1qOGD+/PmEh4dz2223sWPHDmJjYwkJCcFqtZKbm0t8fDwrVqxg5syZWCwW/vSnP/GrX/2K\nwsJCHA7HaXsBysvbP9Zis4VQUlLZ7sfxdmpH11Fbuo7a0nXUlm0XgBOAnPzD9PD9cfFee9vyVAHC\nrULA9OnTuf/++/nqq6+wWCw8/vjjAMyZM4f77rsPu91OWloaw4YNA5p7CyZPnozD4WDWrFlGli4i\nItIuRhwn7FYhICwsjL/+9a/HXR8+fDhvv/32cdczMjLIyMjojNJEREQ6VOQPWwd35oZB2ixIRETE\nDSR3C/7/7d1NSFR7HMbxJ8eXLE2L6A2yULRFEfnSIrySUNJCCXSwMdMKAnEVSIsiaHYFYdBCDHKj\nIFlIuEjpBUwikIIsMxorUMpdGuSQLzEyzf8uLs29Xr1n7sLxzHC+n92cs/DhN8P8n/mf44zWJiZo\na1ryqv3NmNoJAADAqfI2r9d40x9Kcq3e53N2AgAAiBGrWQAkSgAAAI5FCQAAwKEoAQAAOBQlAAAA\nh6IEAADgUJQAAAAcihIAAIBDUQIAAHAoSgAAAA5FCQAAwKEoAQAAONQaY4yxOwQAAFh97AQAAOBQ\nlAAAAByKEgAAgENRAgAAcChKAAAADkUJAADAoSgBywiFQvJ6vfJ4PKqvr9fExMSi8wMDA3K73fJ4\nPOru7rYpZXyINMu+vj5VV1erpqZGXq9XoVDIpqSxL9Isf7ty5Ypu3LixyuniR6Q5vnv3TrW1tTp5\n8qTOnz+vQCBgU9LYF2mWDx48UGVlpdxut7q6umxKGV9GRkZUX1+/5HjU1h2DJZ48eWIuXrxojDFm\neHjYNDY2hs8tLCyYo0ePGr/fbwKBgKmqqjLfvn2zK2rMs5rlz58/zZEjR8z8/LwxxpimpibT399v\nS854YDXL3+7evWtOnDhhmpubVzte3LCaYygUMsePHzdfvnwxxhjT3d1txsfHbckZDyK9JouLi830\n9LQJBALh9038t7a2NlNRUWGqq6sXHY/musNOwDJev36tkpISSdKBAwf0/v378Lnx8XFlZWUpIyND\nycnJKiws1KtXr+yKGvOsZpmcnKx79+4pNTVVkhQMBpWSkmJLznhgNUtJevPmjUZGRuTxeOyIFzes\n5vj582dlZmaqo6NDdXV18vv9ys7OtitqzIv0mtyzZ49mZma0sLAgY4zWrFljR8y4kZWVpZaWliXH\no7nuUAKWMTs7q7S0tPBjl8ulYDAYPpeenh4+t379es3Ozq56xnhhNcuEhARt3rxZktTZ2an5+XkV\nFxfbkjMeWM1yampKra2t8nq9dsWLG1ZznJ6e1vDwsOrq6tTe3q6XL1/qxYsXdkWNeVazlKTc3Fy5\n3W6Vl5ertLRUGzZssCNm3Dh27JgSExOXHI/mukMJWEZaWprm5ubCj0OhUPiJ+fe5ubm5RU8OFrOa\n5e/H169f1+DgoFpaWvikYMFqlo8fP9b09LQaGhrU1tamvr4+9fT02BU1plnNMTMzU7t27VJOTo6S\nkpJUUlKy5NMt/mY1y48fP+rZs2d6+vSpBgYG9P37dz169MiuqHEtmusOJWAZBQUFev78uSTp7du3\nysvLC5/LycnRxMSE/H6/FhYWNDQ0pPz8fLuixjyrWUqS1+tVIBDQrVu3wpcFsDyrWZ4+fVo9PT3q\n7OxUQ0ODKioqVFVVZVfUmGY1x507d2pubi58g9vQ0JByc3NtyRkPrGaZnp6utWvXKiUlRS6XS5s2\nbdKPHz/sihrXornuLN13gMrKyjQ4OKiamhoZY3Tt2jX19vZqfn5eHo9Hly5d0rlz52SMkdvt1tat\nW+2OHLOsZrlv3z7dv39fRUVFOnPmjKS/FrOysjKbU8emSK9L/D+R5nj16lVduHBBxhjl5+ertLTU\n7sgxK9IsPR6PamtrlZSUpKysLFVWVtodOa6sxrrDrwgCAOBQXA4AAMChKAEAADgUJQAAAIeiBAAA\n4FCUAAAAHIoSAACAQ1ECAABwKL4sCEDUVFdXa8eOHcrOzlZvb68mJye1e/duXb58WYcOHbI7HuB4\nfFkQgKgIBoMqKChQQkKC9u/fr7Nnz+rXr1+6efOmJicn1d/fr40bN9odE3A0dgIARMXY2JgCgYAO\nHjyo9vZ2uVwuSVJGRobq6+s1NDTEV0QDNuOeAABRMTo6KklqamoKFwBJys7OliT5/X5bcgH4GyUA\nQFT4fD5t2bJFhYWFi45PTU1JEj+8BcQASgCAqPD5fMsu9A8fPlRqaqqKiopsSAXgn7gnAMCKC4VC\n+vTpk1JTUxUMBpWY+NdbzeTkpLq6unTq1CmtW7fO5pQA+O8AACtubGxM5eXl2r59u4qKiuR2u/X1\n61e1trYqMzNTd+7cUUpKit0xAcfjcgCAFefz+SRJbW1tmpmZUWNjo5qbm3X48GF1dHRQAIAYweUA\nACtudHRU27ZtU15enm7fvm13HAD/gZ0AACvO5/Np7969dscAEAElAMCKMsbow4cPlAAgDnBjIAAA\nDsVOAAAADkUJAADAoSgBAAA4FCUAAACHogQAAOBQlAAAAByKEgAAgENRAgAAcKg/ATUhZO4dNPzQ\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe50d9e0c18>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "p = np.linspace(0, 1, 100)[1:-1]\n",
    "ls = []\n",
    "for i in p:\n",
    "    dist = stats.binom(n=1, p=i)\n",
    "    l = np.sum(dist.logpmf(monedas_obs))\n",
    "    ls.append(l)\n",
    "\n",
    "plt.plot(p, ls)\n",
    "p_hat = p[np.argmax(ls)]\n",
    "plt.plot(p[np.argmax(ls)], np.max(ls), 'bo',\n",
    "         label='$\\hat p ={:.2f}$'.format(p_hat))\n",
    "\n",
    "plt.ylabel('$log\\mathcal{L(p;\\mathbf{x})}$', fontsize=16)\n",
    "plt.xlabel('$p$', fontsize=16)\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Intervalos de confianza\n",
    "\n",
    "$\\hat p$ es un ejemplo de estimación puntual, un número que usamos como estimación de un parámetro de una población. En este caso particular la media de la muestra como aproximación de la media de la población.\n",
    "\n",
    "Suele ser buena idea acompañar la estimación puntual con una medida de la incertidumbre asociada a esa estimación. Una forma de medir esta incerteza es calculando lo que se conoce como **intervalo de confianza**. \n",
    "\n",
    "Este intervalo se puede construir a partir de la distribución de muestreo. A medida que $n$ aumenta el teorema del límite central garantiza que dada una población con media finita $\\mu$ y varianza finita (y distinta de cero) $\\sigma^2$: \n",
    "\n",
    "$$\\hat \\mu \\sim \\mathcal{N} \\left(\\mu, \\frac{\\sigma^2}{n}\\right)$$\n",
    "\n",
    "Entonces podemos definir un intervalo de confianza en términos de la estimación puntual y el error estándar (SE) de esa media:\n",
    "\n",
    "$$IC = [\\hat p - SE, \\hat p + SE]$$\n",
    "\n",
    "Este sería un intervalo de confianza de $\\pm 1 SE$, si quisiéramos ampliar o achica el rango podría usar $X$ errores estándar, por ejemplo:\n",
    "\n",
    "$$IC = [\\hat p - X SE, \\hat p + X SE]$$\n",
    "\n",
    "En la práctica es común que la gente defina los $IC$ en términos de porcentajes y no $SE$, de hecho el porcentaje más comúnmente usado es 95%. Convertir entre porcentajes y $SE$ es simple si asumimos que la distribución es Gaussiana. En ese caso es conocido que un intervalo del 95% estará de forma proximidad dentro de 1.96 errores estándar (o lo que es lo mismo un Z-score de 1.96).\n",
    "\n",
    "Podemos usar el siguiente código para pasar de porcentajes a _cantidades_ de $SE$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def IC_to_a(ic=95):\n",
    "    a = (100 - ic) / 100\n",
    "    gaussian = stats.norm(0, 1)\n",
    "    l, u = gaussian.ppf([(a / 2), (1 - a / 2)])\n",
    "    x = np.linspace(-4, 4, 300)\n",
    "    y = gaussian.pdf(x)\n",
    "    plt.plot(x, y);\n",
    "    plt.fill_between(x[x < l], y[x < l], alpha=0.5, color='C2')\n",
    "    plt.fill_between(x[x > u], y[x > u], alpha=0.5, color='C2',\n",
    "                    label='$\\\\alpha$ = {:.3f}'.format(a))\n",
    "\n",
    "    plt.vlines([l, u], 0, gaussian.pdf([l,u]),\n",
    "              label='{:.3f} {:.3f}'.format(l, u))\n",
    "\n",
    "    plt.legend(fontsize=14)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Widget Javascript not detected.  It may not be installed or enabled properly.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "be816dd81e6f410a9bf36d2bed2ab7a1"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "interact(IC_to_a, ic=(0, 100, 5));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dada la definición de $IC$ y el hecho que usando $\\pm 1.96$ errores estándard podemos definir un intervalo de confianza del 95%, usemos Python para calcularlo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcwAAAFJCAYAAAAWph3tAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEOJJREFUeJzt3X2Q1XW9wPEPuyskD4HoQndUEExSpmkczIqGRCagvPKQ\nRD5QPFSIIASomFwaUJTxjnVvhZMYemV8SpuxtMnGwuzhKgRegnRUCMLrHbtosCKGu6Asu7/7h8Ne\nkQU/nF2WZXu9/sFzfr9zzvezC7/3+Z1z3G1XFEURAMAhlR3tBQDAsUAwASBBMAEgQTABIEEwASBB\nMAEgoeJQG6uq3mzyA5xwQsfYsWNXk+/nWGDWtsmsbZNZ257mmrOyskuj1x/xM8yKivIj/RCthlnb\nJrO2TWZte470nF6SBYAEwQSABMEEgATBBIAEwQSABMEEgATBBIAEwQSABMEEaANefHFzjBr1uXjx\nxc1HeyltlmACtAH33rssbr/9rrj33mVHeylt1iF/liwAx4aFC2/e789SrF79h1i8+N+ivr4+Roz4\nQowfP+mg+9bV1cXkyeOjsrJHfPvb34+IiLFjR0bHjh2jrKw8ysvL46677juia3nvPldf/Y2IiHj5\n5f+JBQvmNez3yitbYvLkK+Lii8eVvJ4IwQQg3gngd797S3zve7dFjx49Y/LkCTFo0HnRp0/fRvd/\n6KEHo3fvPrFrV81+199669Lo1q3bIR9r3bo/xi9/+Yv41rduKHktje0zatQF0bVrz+jV67S4++4H\nGva76KJ/jvPOG3IYX43GeUkW4Bh3/fXzYsGCf4nLL58QX/ziiPjDH1Yc9n1s2PBCnHLKqXHyyafE\ncccdF0OHDo8VK/6z0X23bdsaq1atjJEjv9DUpZe8lsb2+c1vfnPAfa1duyZOPvnk+NCH/qnJ63KG\nCdBCHlm/Nb6/6uXY9FpN9DupU8we2Csu6t+zyff74ot/iUGDBseNN/5rPPvsM/GDH3w3Pv3pQQ3b\nx40bF3//+4G/rnH69Flx7rmfjIiIqqpt0aPH/6+lsrJHrF//fKOPd+ut/x7Tps084OyyXbt2MXv2\nlVFWVhajR4+J0aPH7Lf98ssnRm1tbezevSt27twZkya98xLptGnfiE9+cmDDfpm1NLbPSy9tOmCt\nTzyxPIYO/VyjcxwuwQRoAY+s3xpX/HxDw+UNVTUNl5sSzbfffjveeGNHfPWrl0dERJ8+feLNN/eP\n4wMPPNAsv984ImLlyqeiW7fuceaZZ8W6dX/cb9uSJf8RlZU9YseO12P27OnRu/dpcfbZAxq233nn\nPRHx/i/JNpfa2tpYufLJmDp1RrPcn2ACHEHnLFkVERF/q97T6PYZv/hzLPr9f8faKwc2uv39vPTS\ni3HKKb2iQ4cOERGxceOf48MfPmO/fTJnmJWVPWLbtq0N26qqtkVlZY8DbvPcc8/GypVPxurVK2PP\nnj1RU1MdN944PxYsuKlh/xNO6B7nnXd+rF//wn7BzMqspbF9evbc/4nH6tUro1+/M6N79xMPew2N\nEUyAFlBbXxzW9VmbN2+KrVv/Fm+//XbU19fHsmVLY9q0mfvtkznDPPPM/vHXv/41XnllS1RW9ogn\nnng8rr9+0QH7TZ06o+GMbd26P8aPf3x/LFhwU+zevTuKoj46duwUu3fvjjVrno5JkyY3+lgDBnw8\nBgz4eJPW0tg+ixd/b799mvPl2AjBBDii9p05Dr5rTWyoqjlge//KTvH7r59b8v1v3vyXGDx4SEyZ\nMinq6vbG+PFfjY997OzDvp+Kioq4+upr4+qrvxH19XVx4YWjom/f0xu2z5kzM+bOnR8nnVTZ6O1f\nf317zJt3bUS888nUYcM+F5/61Kf322ffe5jv9d73MA+1lnev4737nHHGGQ1PDN6J9n/Ftdd+67C/\nFgfTriiKgz69aY7XvCsruzTba+etnVnbJrO2TS0963vfw9xn6aizmvQe5owZU+Kb35wXvXqddtB9\n/lG+r801Z2Vll0avd4YJ0AL2RXHxqpdj0/Zd0e/EjjGrGT4lu2XL/8Ypp/RqjiXyPgQToIVc1L9n\ns/xvJO/2yCOPNev9cXB+cAEAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgm\nACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYA\nJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAk\nCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQI\nJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgm\nACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYA\nJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAk\nCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQI\nJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkVBztBQD/\n75H1W+P7q16OTa/VRL+TOsXsgb3iov49j/aygBBMaDUeWb81rvj5hobLG6pqGi6LJhx97YqiKA62\nsarqzSY/QGVll4b7OWfJqibfX2tWVl4W9XX1R3sZLcKsze9v1Xuitv7Af47HlbWLD3Vuf8QfP8L3\nta1qy7OuvXJgw3+/uzdNUVnZpdHrvYcJrURjsTzU9UDLatGXZN/9TKAtaq5nN8cCsza/wXetiQ1V\nNQdc37+yU/z+6+ce8ceP8H1tq/6RZj2SnGFCKzF7YK9Gr591kOuBluVDP9BK7Ptgz+JVL8em7bui\n34kdY5ZPyUKrIZjQilzUv6dAQivlJVkASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABI\nEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQ\nTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBM\nAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwA\nSBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABI\nEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQ\nTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBM\nAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwA\nSBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABI\nEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQ\nTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASGhXFEVxtBcBAK2dM0wASBBMAEgQTABIEEwA\nSBBMAEgQTABIqGjKjevr6+OGG26IjRs3Rvv27WPRokXRu3fvhu3Lly+PO+64I9q1axcjR46MiRMn\nRm1tbcydOze2bNkSZWVlcdNNN8Xpp5/e5EGOtPebdZ/58+dH165dY86cOenbtDalzFpbWxvz5s2L\nLVu2xJ49e2LatGnx2c9+9iis/vCUMus+27dvjzFjxsSyZcta/d/hUudcunRp/Pa3v43a2tq47LLL\n4ktf+lJLL/2wlfr3ty0el+6+++546KGHonv37hERsXDhwjjttNPa5HGpsVlPPfXU5j0uFU2wfPny\n4rrrriuKoij+9Kc/FVOnTm3Ytnfv3mLYsGHFzp07i7179xbDhw8vtm/fXvz6178uZs6cWRRFUaxY\nsaKYMWNGU5bQYg416z4PPvhgcfHFFxff+c530rdpjUqZ9Sc/+UmxaNGioiiKYseOHcXgwYNbbL1N\nUcqsRVEUe/bsKa688spi+PDhxebNm1tsvaUqZc7Vq1cXV1xxRVFXV1dUV1cXt956a4uuuVSlzNpW\nj0vXXHNN8dxzzx3WbVqrUmZt7uNSk16SXbt2bXzmM5+JiIizzz47nn/++YZt5eXl8dhjj0WXLl3i\njTfeiPr6+mjfvn306dMn6urqor6+Pqqrq6OiokknuS3mULNGRKxbty6effbZuOSSS9K3aa1KmfXz\nn/98zJo1KyIiiqKI8vLylltwE5Qya0TELbfcEpdeemn06NGjxdbaFKXMuWLFiujXr19Mnz49pk6d\nGueff35LLrlkpczaVo9LL7zwQtxxxx1x2WWXxdKlS1O3aa1KmbW5j0tNCmZ1dXV07ty54XJ5eXns\n3bu34XJFRUU8/vjjMXr06PjEJz4Rxx9/fHTs2DG2bNkSF1xwQcyfPz/Gjx/flCW0mEPNum3btrjt\ntttiwYIF6du0ZqXM2qlTp+jcuXNUV1fHzJkzY/bs2S265lKVMuvDDz8c3bt3b/jHeywoZc4dO3bE\n888/H4sXL46FCxfGnDlzojgGfjBYKbO2xeNSRMSFF14YN9xwQ9xzzz2xdu3a+N3vftcmj0sRjc/a\n3MelJj2N6ty5c9TU1DRcrq+vP+CZ2fDhw2Po0KExd+7c+NnPfhabNm2KQYMGxTXXXBOvvvpqTJw4\nMR599NHo0KFDU5ZyxB1q1l/96lexY8eOmDJlSlRVVcVbb70Vffv2TX19WqNSZh0zZky8+uqrMX36\n9Bg3blyMHDnyaC3/sJQy609/+tNo165drFq1KjZs2BDXXXdd3H777VFZWXm0xnhfpczZrVu36Nu3\nb7Rv3z769u0bHTp0iNdffz1OPPHEozVGSimzbty4sc0dl4qiiIkTJ0aXLl0iImLw4MGxfv36Nnlc\nOtisQ4YMadbjUpPOMAcMGBBPPvlkREQ888wz0a9fv4Zt1dXV8ZWvfCX27NkTZWVlcfzxx0dZWVl8\n8IMfbBiqa9eusXfv3qirq2vKMlrEoWadMGFCPPzww3HffffFlClTYsSIETFmzJhD3qY1K2XW1157\nLb72ta/FtddeG2PHjj1aSz9spcz6ox/9KO6///6477774qyzzopbbrmlVccyorQ5zznnnHjqqaei\nKIrYunVr7N69O7p163a0RkgrZda2eFyqrq6OESNGRE1NTRRFEU8//XR89KMfbZPHpYPN2tzHpSY9\nrRg2bFisXLkyLr300iiKIm6++eZ49NFHY9euXXHJJZfEyJEj48tf/nJUVFTERz7ykRg1alS89dZb\nMW/evBg3blzU1tbGVVddFR07dmzyIEfa+82avc2xoJRZf/jDH8bOnTtjyZIlsWTJkoiIuPPOO+MD\nH/hASy79sJUy67GolDmHDBkSa9asibFjx0ZRFLFgwYJj4r3pUmadNGlSmzwuXXXVVTFhwoRo3759\nDBw4MAYPHhz19fVt8rjU2KyLFi1q1uOS31YCAAl+cAEAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYA\nJAgmACT8H27wZHFML8F/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe50d821550>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "p_hat = caras_obs/n\n",
    "dist_monedas = stats.binom(n=1, p=p_hat)\n",
    "std_error =  dist_monedas.std() / n**0.5\n",
    "r = std_error * 1.96  \n",
    "\n",
    "plt.errorbar(x=p_hat, y=0,\n",
    "             xerr=r,\n",
    "             fmt='o',\n",
    "            label='$\\hat p = {:.2f} \\pm {:.2f}$'.format(p_hat, r))\n",
    "plt.yticks([])\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Como ya dijimos es muy común usar un intervalo con un *nivel de confianza* del 95%. Pero este valor es arbitrario, tan arbitrario como el valor de $\\alpha$ en las pruebas de hipótesis. De hecho ambos valores están relacionados.\n",
    "\n",
    "Esencialmente un intervalo de confianza contiene a todos los valores de $\\hat p$ para los cuales $H_0$ NO sería rechazada, para un dado valor de $\\alpha$. Los intervalos de confianza se suelen expresar como $(1 - \\alpha) 100 %$, por lo que un valor de $\\alpha=0.05$ se relaciona con un intervalo de confianza de 95%.\n",
    "\n",
    "Como sucede con los _valores p_ resulta que es fácil interpretar los $IC$ de forma errónea. Por ello veamos cual es la interpretación adecuada.\n",
    "\n",
    "Si repetimos un experimento (infinitas veces) y cada vez calculamos un intervalo de confianza, entonces obtendremos que el X% de esos intervalos contendrán el valor _real_ del parámetro de interés.\n",
    "Es decir el nivel de confianza, digamos el 95%, NO es sobre un intervalo en particular es sobre todo el universo posible de intervalos que se podrían calcular con muestras similares a la nuestra. Por lo tanto NO es posible decir que tenemos una confianza del 95% que NUESTRO intervalo contenga el parámetro, en sentido estricto nuestro intervalo o contiene o no contiene al parámetro _real_. Si quisiéramos hacer enunciados como tengo un 95% de confianza que el parámetro real se encuentra en tal o cual rango, deberé hacer uso de estadística Bayesiana!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bootstrapping\n",
    "\n",
    "EL bootstrapping es una técnica computacional que permite aproximar la _distribución de muestreo_. Su principal ventaja es su simplicidad lo que permite usarla con estimadores sencillos como el que hemos usado hasta ahora y también con estimadores más complejos. Esta técnica es parte de una familia más grande de métodos que usualmente se conocen como métodos de _remuestreo_.\n",
    "\n",
    "La idea del bootstrapping consiste en generar un número _grande_ de muestras a partir de (re)muestrear los propios datos! El muestreo se hace con reemplazo:\n",
    "\n",
    "1. _saco_ un elemento de los datos\n",
    "2.  _anoto_ su valor\n",
    "3. lo _vuelvo a poner_ con el resto de los datos\n",
    "4. repito de 1 a 3 hasta obtener una (re)muestra del mismo tamaño que la muestra original\n",
    "5. a partir de esa (re)muestra calculo lo que sea que me interese (la media, la desviación estándar, etc)\n",
    "6. repito de 4-5 miles de veces.\n",
    "\n",
    "De esta forma habré logrado calcular miles de veces el estimador de interés por lo que no solo tendré una estimación puntual, si no una distribución de valores para el estimador. Usando Python se puede escribir de forma más breve incluso que usando español ;-)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAFJCAYAAABU5W56AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHSlJREFUeJzt3X1clHW+//H3MOMgwiC6h043hklJmw+XFSP8dUqM7qi8\nyTxpMiey1Uc+pMwDqyVZQFtuYh7ZHpu53qybJ0qRtG1rj92czLLUB8fIm7SwlS0ru5ESc2ZUbuT6\n/dFDinRAx4H5Mryef8Vc1zXX9/rwxXff71zzvWyWZVkCAABGigh1AwAAgH8ENQAABiOoAQAwGEEN\nAIDBCGoAAAxGUAMAYDBHqBtwMjU1nlA3wUi9evVQbe3hUDfDSNTGP2rjH7VpHfXxL9i1iY93+d3G\niLoTcTjsoW6CsaiNf9TGP2rTOurjX0fWhqAGAMBgBDUAAAYjqAEAMBhBDQCAwQhqAAAMRlADAGAw\nghoAAIMR1AAAGIygBgDAYAQ1AAAGI6gBADAYQQ0ACIrq6j0aNSpT1dV7Qt2UsGLk07MA/Ois4rdO\n+5j9+VcFvR1AW5555i/605+WacmShfrd7x4LdXPCBkENAAiK4+FMSAcXQQ0ACLmmpibNn1+sPXv+\noW7duik/v0B9+px/0n1raw9o0qRs/eEPT6lv3ws0ceJ/qEePaEnSueeep1mzivye5/3339Mrr/xd\nDz74cEDtOL59795/Sopo3r527ctau/ZlSVJ9fb327PlYf/vba3K5/D9n+lQR1EAYCmS6XGLKHIEr\nKpoly7L01Vf7dODAAU2fnq9/+7crT/n4d955S/X19Vq8+Gnt3PmBFiz4g4qLS07Yr7GxUY8//pic\nzkhJUl1dnSzL0oIFS4JyHW214/j2VatWaf36Tc3bb7pppG66aaQkaf78uRo+fFRQQloiqAEAQVBd\n/Q9deeUwPfLIHG3fvk0LFpS0COr169/QmjXlLY65++5pGjBgoCRpx45tGjLkcknSwIG/UlXVRyc9\nz4IFT2j06H9XaenTkqQ9e/6ho0ePKi/vHh07dkyTJ9+jgQN/dcJxd901QQ0NDTpy5LAOHTqkO+90\nS5Jycu5tPu+ptKOt7VVVH+qTT6o1ffrMNip26ghqAMAZqaur08GDtfrNb+6SJPXr108ej6fFPhkZ\n1yoj41q/7+Hz+RQdHdP8c0REhBobG+Vw/BhTa9e+rLi4OA0ZcnlzUHfv3l1ZWdkaOXK0Pv/8M82Y\nMU0rVqxpcZwkLV3635Lanvpuqx1tbX/mmac1ceJkv9cZCIIaAMJIoB97+HMqH4d88km1+vRJUGTk\nD9PRu3dX6aKL+rfYp60RdXR0tA4fPty8zbKsE8L2f/7nJdlsNr333v9pz56PNXt2oR59dK4yM2+U\nzWZTQkJf9ezZU999963+9V/PDuRy22xHa9s9Ho8++2yvBg9ODejc/hDUAIAzsmfPx/rmm69VV1en\npqYm/eUvi5WTM63FPm2NqH/1q19r48Z3dM0112nnzg+UmHjRCfs89dTS5v+eOnWy7rtvljZv3qjq\n6j2aMSNf335bI5/Pp1/84l/8nmfw4NRWg7StdhzfPn78mBO2b9/+vlJTL/P73oEiqAEAZ2TPnn9o\n2LAMTZ58p44da1R29m+UnDzotN4jPT1DW7ZUaMqUibIsq/nO7UOHvldx8Ww99ti8kx43YsTN+v3v\nH1ZOziTZbDY98EDhCSNx6cfPqH/u559Rn6wdP23D8e3jx49XfX1jizvMP/tsr84997zTuu5TYbMs\nywr6u56hmhpP2zt1QfHxLmrjRzjXJthTma3pand9h2O/CcXU99Spk3X//bOUkHBBUM9tsmD3nfh4\n/3eIs4QoAOCM7Nv3hfr0SQh1M8IWU98AgDPy17+uDXUTwhojagAADMaIGgDCSDDvMwjHz/A7I4Ia\nQDOe1AWYh6lvAAAMxogaABAUa9e+rL17P1VOzr2qq6vT0qV/0ocf7pTNZlNUVJTuu29WqyuGvfvu\nBi1f/mfZ7XYNHz5Ko0bd0mL7t99+q0cfLVBDQ4NiY2NVWPioevSI1uuvv6KysmcVEfHDcbfccutp\nv/ehQ98rK2uM+vW7UNIP36ceNy5L//u/r6q8fKUcDrsSEy/S9On5iojo2DEuQQ0ACLo//nG++va9\nQFOn/lmS9Pbb61VY+IAWL376pPs3NjbqySdLtHTpM4qKilJOziRdeWW6evf+RfM+zz3337rhhuG6\n8cYRWrZssV5++UXddtt/6KmnnlBpabmionro9tvH6pprrldsbOxpvffu3VW69tpM5eXd3/xaXd1R\nLV36Jz3zzCp1795dRUWztGnTO7ryymHBLlerCGoAQFA1NDTo3Xff1owZDzS/NmxYhgYNSpEk3X9/\nbov1si+4IFGjR/+7zjvv/OaATU7+tbZt26qrr/5x2dFp034ry7LU1NSk/fu/0dlnnyNJuvDC/vJ6\nvbLb7bIsSzabrUV7Pv30kzbfe/fuj7R7d5WmTp2suLheys29T71799aiRX9R9+7dJUnHjh1rfrxm\nRyKoAQBB9f3336t371+cEJg9e8ZJkh5//IkTjtm+fZtiYn58KlWPHtHy+bwt9rHZbDp27JjuvDNL\ndXX1P3la14WaNClbUVFRSk/POOE50D6fr8337tv3Al188SW67LIhev31V/TEE49r9uzHm0fdq1eX\n6ciRI7rssiGnW44zRlADAIIqLi5OXq/3hNHt66+/oquvvk6zZs046Yj68GFf82uHD7cM1+McDoee\nffZ5bdlSodmzi5Sbe582b35Xzz//kqKiovTIIwV68803tGfPx9qxY5sk6T//c0ab733ppZcpMvKH\nkXN6eob+/OdFkqSmpiYtXPhHff75Xv3+94+f8D8fHYGgBgAElcPhUFra/9Pq1as0dux4SdKbb76h\n559fqeuvv/GkI+rGxkZ98cXnOnToe0VF9dC2bVuVlZXdYp//+q9iXX31tRo8OFU9ekTLZrMpJiZG\nkZGRioyMlN1uV69eveXxHNLkyXef1nsXF8/WsGFX65prrtN77/2fLr74EknSvHmPqVu3bpozZ36H\n30R2HEENAAi6e+/9rZ58skRTpkyUzSa5XLGaPftxv/s7HA5NnZqn3/72XjU1NWn48FGKjz+rxZOr\nxo4dr3nzHtPTTy9VRESEpk/P19lnn6Obbx6ju++eJIejm847r49uumnkab/3lClTNWfOI/rrX59X\nVFSUZs4s0O7dVfr73/+mX/86RdOmTZEkjR2bpWHDMtq1dj/H07M6EVYJ8i+ca9ORT88KRGde8CSc\n+00wUB//eHoWAACQRFADAGA0ghoAAIMR1AAAGIygBgDAYAQ1AAAGO6Wg3r59u7Kzf/hy+N69e5WV\nlSW3262ioiI1NTVJksrLyzVmzBiNGzdO69evlyQdPXpU9957r9xut+666y4dOHCgnS4DAIDw1GZQ\nL126VA899JDq6uokSXPmzFFubq5WrFghy7K0bt061dTUqLS0VGVlZVq2bJlKSkpUX1+vlStXKikp\nSStWrNDo0aO1cOHCdr8gAADCSZtBnZCQoCeffLL55127diktLU2SlJ6erk2bNmnHjh1KSUmR0+mU\ny+VSQkKCqqqqVFlZqaFDhzbvu3nz5na6DAAAwlObS4hmZmbqiy++aP75p4usR0dHy+PxyOv1tnha\nSXR0tLxeb4vXj+97Knr16iGHw35aF9JVtLZ6TVdHbUKjs9e9s7e/vVEf/zqqNqe91vdPFyX3+XyK\njY1VTEyMfD5fi9ddLleL14/veypqaw+3vVMXxHJ+/lGb0OnMdafftI76+Gf0EqIDBgxQRUWFJGnD\nhg1KTU1VcnKyKisrVVdXJ4/Ho+rqaiUlJWnw4MF6++23m/e99NJLA7wEAAC6ptMeUc+cOVMFBQUq\nKSlRYmKiMjMzZbfblZ2dLbfbLcuylJeXp8jISGVlZWnmzJnKyspSt27dNH/+/Pa4BgAAwhZPz+pE\nmIbyrzPUxvSnYAWKp2eFL+rjn9FT3wAAoOMQ1AAAGIygBgDAYAQ1AAAGI6gBADAYQQ0AgMEIagAA\nDEZQAwBgMIIaAACDEdQAABiMoAYAwGAENQAABjvtp2cBwE8F+rCRzvwwD6AjMaIGAMBgBDUAAAYj\nqAEAMBhBDQCAwQhqAAAMRlADAGAwghoAAIMR1AAAGIygBgDAYAQ1AAAGI6gBADAYQQ0AgMEIagAA\nDEZQAwBgMIIaAACDEdQAABiMoAYAwGAENQAABiOoAQAwGEENAIDBCGoAAAxGUAMAYDBHqBsAdEZn\nFb8V6iYA6CIYUQMAYDCCGgAAgxHUAAAYjKAGAMBgBDUAAAYjqAEAMFhAX89qaGhQfn6+9u3bp4iI\nCD366KNyOBzKz8+XzWZT//79VVRUpIiICJWXl6usrEwOh0M5OTnKyMgI9jUAABC2Agrqt99+W42N\njSorK9PGjRv1xBNPqKGhQbm5uRoyZIgKCwu1bt06DRo0SKWlpVqzZo3q6urkdrt1xRVXyOl0Bvs6\nAAAISwFNfffr10/Hjh1TU1OTvF6vHA6Hdu3apbS0NElSenq6Nm3apB07diglJUVOp1Mul0sJCQmq\nqqoK6gUAABDOAhpR9+jRQ/v27dONN96o2tpaLVq0SFu2bJHNZpMkRUdHy+PxyOv1yuVyNR8XHR0t\nr9fb5vv36tVDDoc9kKaFvfh4V9s7dVHUpnMx5fdlSjtMRX3866jaBBTUy5cv15VXXqnp06frq6++\n0oQJE9TQ0NC83efzKTY2VjExMfL5fC1e/2lw+1NbeziQZoW9+HiXamo8oW6GkahN52PC74t+0zrq\n41+wa9Na6Ac09R0bG9scuD179lRjY6MGDBigiooKSdKGDRuUmpqq5ORkVVZWqq6uTh6PR9XV1UpK\nSgrklAAAdEkBjajvvPNOzZo1S263Ww0NDcrLy9PAgQNVUFCgkpISJSYmKjMzU3a7XdnZ2XK73bIs\nS3l5eYqMjAz2NQAAELZslmVZoW7EzzHVcnJMQ/nX0bXh6Vlnbn/+VaFuAn9TbaA+/hk/9Q0AADoG\nQQ0AgMEIagAADEZQAwBgMIIaAACDEdQAABiMoAYAwGAENQAABiOoAQAwGEENAIDBCGoAAAxGUAMA\nYDCCGgAAgxHUAAAYjKAGAMBgjlA3AEDXFMgzvU14hjXQ0RhRAwBgMIIaAACDEdQAABiMoAYAwGAE\nNQAABiOoAQAwGEENAIDBCGoAAAxGUAMAYDCCGgAAgxHUAAAYjKAGAMBgBDUAAAYjqAEAMBhBDQCA\nwQhqAAAMRlADAGAwghoAAIMR1AAAGIygBgDAYAQ1AAAGI6gBADAYQQ0AgMEIagAADEZQAwBgMIIa\nAACDEdQAABjMEeiBixcv1ptvvqmGhgZlZWUpLS1N+fn5stls6t+/v4qKihQREaHy8nKVlZXJ4XAo\nJydHGRkZwWw/AABhLaARdUVFhbZu3aqVK1eqtLRUX3/9tebMmaPc3FytWLFClmVp3bp1qqmpUWlp\nqcrKyrRs2TKVlJSovr4+2NcAAEDYCiio3333XSUlJemee+7RlClTdNVVV2nXrl1KS0uTJKWnp2vT\npk3asWOHUlJS5HQ65XK5lJCQoKqqqqBeAAAA4Sygqe/a2lp9+eWXWrRokb744gvl5OTIsizZbDZJ\nUnR0tDwej7xer1wuV/Nx0dHR8nq9bb5/r1495HDYA2la2IuPd7W9UxdFbcJfe/yO6Tetoz7+dVRt\nAgrquLg4JSYmyul0KjExUZGRkfr666+bt/t8PsXGxiomJkY+n6/F6z8Nbn9qaw8H0qywFx/vUk2N\nJ9TNMBK16RqC/Tum37SO+vgX7Nq0FvoBTX1feumleuedd2RZlr755hsdOXJEl19+uSoqKiRJGzZs\nUGpqqpKTk1VZWam6ujp5PB5VV1crKSkpsKsAAKALCmhEnZGRoS1btujWW2+VZVkqLCxUnz59VFBQ\noJKSEiUmJiozM1N2u13Z2dlyu92yLEt5eXmKjIwM9jUAABC2bJZlWaFuxM8x1XJyTEP519G1Oav4\nrQ47F360P/+qoL4ff1Otoz7+GT/1DQAAOgZBDQCAwQhqAAAMRlADAGAwghoAAIMR1AAAGIygBgDA\nYAQ1AAAGC/h51ADQ0QJdaCbYC6UAHYkRNQAABiOoAQAwGEENAIDBCGoAAAxGUAMAYDDu+kbYCOSO\nYO4GBmA6RtQAABiMoAYAwGAENQAABiOoAQAwGEENAIDBCGoAAAxGUAMAYDCCGgAAgxHUAAAYjKAG\nAMBgBDUAAAYjqAEAMBhBDQCAwQhqAAAMRlADAGAwghoAAIMR1AAAGIygBgDAYAQ1AAAGI6gBADAY\nQQ0AgMEIagAADEZQAwBgMIIaAACDEdQAABiMoAYAwGAENQAABjujoP7uu+80bNgwVVdXa+/evcrK\nypLb7VZRUZGampokSeXl5RozZozGjRun9evXB6XRAAB0FQEHdUNDgwoLC9W9e3dJ0pw5c5Sbm6sV\nK1bIsiytW7dONTU1Ki0tVVlZmZYtW6aSkhLV19cHrfEAAIQ7R6AHzp07V+PHj9eSJUskSbt27VJa\nWpokKT09XRs3blRERIRSUlLkdDrldDqVkJCgqqoqJScnB6f1CEtnFb8V6iYAgDECCuoXXnhBvXv3\n1tChQ5uD2rIs2Ww2SVJ0dLQ8Ho+8Xq9cLlfzcdHR0fJ6vW2+f69ePeRw2ANpWtiLj3e1vRNOGfXs\nGlr7PdMHWkd9/Ouo2gQU1GvWrJHNZtPmzZv10UcfaebMmTpw4EDzdp/Pp9jYWMXExMjn87V4/afB\n7U9t7eFAmhX24uNdqqnxhLoZYYV6dg3+fs/8TbWO+vgX7Nq0FvoBfUb93HPP6dlnn1VpaakuueQS\nzZ07V+np6aqoqJAkbdiwQampqUpOTlZlZaXq6urk8XhUXV2tpKSkwK4CAIAuKODPqH9u5syZKigo\nUElJiRITE5WZmSm73a7s7Gy53W5ZlqW8vDxFRkYG65QAAIS9Mw7q0tLS5v9+9tlnT9g+btw4jRs3\n7kxPAwBAl8SCJwAAGIygBgDAYAQ1AAAGI6gBADAYQQ0AgMEIagAADEZQAwBgMIIaAACDEdQAABiM\noAYAwGAENQAABgvaQzmAzuis4rdC3QQAaBUjagAADEZQAwBgMIIaAACDEdQAABiMoAYAwGAENQAA\nBiOoAQAwGEENAIDBCGoAAAzGymQA4EcgK9ftz78q6O1A18aIGgAAgxHUAAAYjKAGAMBgBDUAAAYj\nqAEAMBhBDQCAwQhqAAAMRlADAGAwghoAAIMR1AAAGIygBgDAYAQ1AAAGI6gBADAYQQ0AgMEIagAA\nDEZQAwBgMEeoG4DwdlbxW6FuAgB0aoyoAQAwGEENAIDBmPoGEPb4CAadWUBB3dDQoFmzZmnfvn2q\nr69XTk6OLrroIuXn58tms6l///4qKipSRESEysvLVVZWJofDoZycHGVkZAT7GgAACFsBBfVLL72k\nuLg4zZs3TwcPHtTo0aP1y1/+Urm5uRoyZIgKCwu1bt06DRo0SKWlpVqzZo3q6urkdrt1xRVXyOl0\nBvs6AAAISwEF9Q033KDMzExJkmVZstvt2rVrl9LS0iRJ6enp2rhxoyIiIpSSkiKn0ymn06mEhARV\nVVUpOTk5eFcAAEAYCyioo6OjJUler1fTpk1Tbm6u5s6dK5vN1rzd4/HI6/XK5XK1OM7r9bb5/r16\n9ZDDYQ+kaWEvPt7V9k4AQibc/kbD7XqCqaNqE/DNZF999ZXuueceud1ujRw5UvPmzWve5vP5FBsb\nq5iYGPl8vhav/zS4/amtPRxos8JafLxLNTWeUDcDQCvC6W+Uf3P8C3ZtWgv9gL6e9e2332rixIm6\n7777dOutt0qSBgwYoIqKCknShg0blJqaquTkZFVWVqqurk4ej0fV1dVKSkoK5JQAAHRJAY2oFy1a\npEOHDmnhwoVauHChJOnBBx/U7NmzVVJSosTERGVmZsputys7O1tut1uWZSkvL0+RkZFBvQAAAMKZ\nzbIsK9SN+DmmWk6uM05D8f1VdDX7868KdROCpjP+m9NRjJ/6BgAAHYOgBgDAYAQ1AAAGI6gBADAY\nQQ0AgMEIagAADMZjLgEgiAL9SmI4fa0LwcWIGgAAgxHUAAAYjKAGAMBgBDUAAAbjZjKcEtbsBoDQ\nYEQNAIDBCGoAAAxGUAMAYDCCGgAAgxHUAAAYjKAGAMBgBDUAAAYjqAEAMBhBDQCAwViZDAAMEMjq\nfzwas2tgRA0AgMEIagAADEZQAwBgMIIaAACDEdQAABiMu767IJ4tDQCdByNqAAAMRlADAGAwghoA\nAIMR1AAAGIybyQCgkwr0xlCWHu1cGFEDAGAwghoAAIMR1AAAGIygBgDAYAQ1AAAG467vToylQAEg\n/BHUANDFBPI/+XylK3SY+gYAwGCMqA3BNDYA4GQIagBAm1gFLXTaPaibmpr08MMPa/fu3XI6nZo9\ne7b69u3b3qcFACAstPtn1G+88Ybq6+u1atUqTZ8+XcXFxe19SgAAwka7j6grKys1dOhQSdKgQYO0\nc+fO9j5l0PC5MQCcmY78dzRcp9nbPai9Xq9iYmKaf7bb7WpsbJTD4f/U8fGu9m7WKbHmjwx1EwAA\nhuqorGr3qe+YmBj5fL7mn5uamloNaQAA8KN2D+rBgwdrw4YNkqRt27YpKSmpvU8JAEDYsFmWZbXn\nCY7f9f3xxx/Lsiw99thjuvDCC9vzlAAAhI12D2oAABA4lhAFAMBgBDUAAAbj9msDtLV622uvvaYl\nS5bIZrNp5MiRmjBhgiTplltuaf7qW58+fTRnzpyQtL+9nerqdgUFBerZs6dmzJjRZVbEC6Q2Utfo\nO23VZvny5Xr++efVu3dvSdLvfvc7XXDBBfQbnbw2iYmJ9BtJO3bsUHFxsSzLUnx8vObNm6du3bq1\nb7+xEHKvvfaaNXPmTMuyLGvr1q3WlClTmrc1NjZa1113nXXo0CGrsbHRuv76663vvvvOOnr0qHXz\nzTeHqskdqrX6HLdy5Upr3Lhx1rx58075mHAQSG26St9pqzbTp0+3Pvjgg9M6JlwEUhv6jWU1NTVZ\no0aNsj799FPLsiyrvLzcqq6ubvd+w9S3AVpbvc1ut2vt2rVyuVw6ePCgmpqa5HQ6VVVVpSNHjmji\nxIm64447tG3btlA1v921tbrd+++/r+3bt+u222475WPCRSC16Sp9p63a7Nq1S0uWLFFWVpYWL158\nSseEi0BqQ7+RPvnkE8XFxWn58uW6/fbbdfDgQSUmJrZ7v2Hq2wBtrd7mcDj0+uuv65FHHtGwYcMU\nFRWl7t27a9KkSRo7dqw+/fRT3XXXXXr11VfDcjGZ1uqzf/9+PfXUU1qwYIFeeeWVUzomnARSm67S\nd9rqA8OHD5fb7VZMTIymTp2q9evX029aqc25557b5ftNbW2ttm7dqsLCQiUkJGjKlCkaOHBgu/eb\n8KpwJ3Uqq7ddf/31uvbaa5Wfn68XX3xRI0eOVN++fWWz2dSvXz/FxcWppqZG55xzTkc3v921Vp9X\nX31VtbW1mjx5smpqanT06FElJiZ2mRXxAqnNiBEjukTfaa02lmVpwoQJcrl+WAJy2LBh+vDDD+k3\n8l+bK664osv3m7i4OPXt27d5LZChQ4dq586d7d5vmPo2QGurt3m9Xt1+++2qr69XRESEoqKiFBER\nodWrVzc/ieybb76R1+tVfHx8SNrf3lqrzx133KEXXnhBpaWlmjx5skaMGKExY8Z0mRXxAqlNV+k7\nbf1djRgxQj6fT5ZlqaKiQgMHDqTfyH9t6DfS+eefL5/Pp71790qS3nvvPfXv37/d+w0LnhjgZKu3\nffjhhzp8+LBuu+02rVq1SqtXr5bD4dDFF1+sgoICHTt2TA888IC+/PJL2Ww2zZgxQ4MHDw71pbSL\ntupz3AsvvKB//vOfLe76DvcV8QKpTX19fZfoO23V5sUXX1RpaamcTqcuv/xyTZs2jX7TSm3oNz/U\nZvPmzZo/f74sy1JKSooeeuihdu83BDUAAAZj6hsAAIMR1AAAGIygBgDAYAQ1AAAGI6gBADAYQQ0A\ngMEIagAADEZQAwBgsP8PXZVLSdZ5mVsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe50d8ba240>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "p_hats = []\n",
    "for i in range(10000):\n",
    "    muestra = np.random.choice(monedas_obs, replace=True, size=len(monedas_obs))\n",
    "    p_hats.append(np.mean(muestra))\n",
    "p_hats = np.array(p_hats)\n",
    "p_hat = np.mean(p_hats)\n",
    "sd = np.std(p_hats) * 2\n",
    "IC = p_hat - sd , p_hat + sd\n",
    "plt.hist(p_hats, bins=30, label='$\\hat p$ = {:.2f} $\\pm$ {:.2f} \\n IC={:.2f}-{:.2f}'.format(p_hat, sd,  *IC))\n",
    "#sns.kdeplot(p_hats, bins=30, label='$\\hat p$ = {:.2f} $\\pm$ {:.2f} \\n IC={:.2f}-{:.2f}'.format(p_hat, sd,  *IC))\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La idea de (re)muestrear los propios datos parece algo mágica o al menos tramposa, y no ayuda saber que el origen del término. \"Bootstraps\" son las tiras que tienen muchas botas cuyo propósito es facilitar el colocárselas. Bootstrapping hace referencia a una expresión del inglés que quiere decir algo así como _levantarse uno mismo tirando de las tiras de las propias botas_, en definitiva algo imposible o que no tiene sentido. \n",
    "\n",
    "De todas formas existen justificaciones teóricas para este procedimiento, que garantizan que bajo ciertas condiciones las aproximaciones son buenas. Además de justificaciones empíricas que han mostrado la utilidad de este método.\n",
    "\n",
    "Existen múltiples variantes de este idea general, la que hemos presentado se llama bootstap no-paramétrico. En estadística frecuentista _no-paramétrico_ quiere decir básicamente que el método NO asume una distribución para lo datos. Fijensé que no asumimos que los datos seguían una distribución binomial (ninguna otra). Algunas de las variantes de bootstrapping son:\n",
    "\n",
    "* Bootstrapping paramétrico, a partir de un distribución (por ej la distribución binomial) se obtiene un estimador puntual (por ej usando Maximun Likelihood) y se usa esa distribución para generar muestras aleatorias (en vez de muestrear de los datos).\n",
    "\n",
    "* Bootstrapping suavizado, lo que se hace es agregar un pequeña cantidad de ruido (generalmente Gaussiano). Esto puede ayudar a \"llenar\" a cada dato de esta forma se obtiene una distribución de muestreo que es más suave. En cierta forma equivale a muestrear no de los datos si no de la _Kernel density estimation_ de los datos. Además este método se puede considerar como un intermedio entre el bootstrapping no-paramétrico y el paramétrico.\n",
    "\n",
    "* Bootstrapping de a bloques, se usa con datos correlacionados por ejemplo series temporales en la que un simple bootstrapping destruiría la correlación existente en los datos.\n",
    "\n",
    "En general se enseña al bootstrapping como un _procedimiento_ o _método_ que no asume nada sobre los datos, mas o menos como hemos hecho hasta ahora. Una alternativa sería mostrar que el bootstrapping es en realidad un modelo estadístico que asume ciertas cosas, de hecho se puede mostrar que es un caso especial de [modelo Bayesiano](http://www.sumsar.net/blog/2015/04/the-non-parametric-bootstrap-as-a-bayesian-model/). De hecho el creador del Bootstrapping se refiere al mismo como un \"poor's man Bayesian method\", algo así como la versión \"trucha/berreta\" de un método Bayesiano. Este método en realidad asume que:\n",
    "\n",
    "* Los datos son discretos\n",
    "* Los datos pueden tomar valores entre $-\\infty$ e $\\infty$\n",
    "\n",
    "Aún cuando casi ningún conjunto de datos satisface estos dos criterios en la práctica el método/modelo de bootstrapping suele ser una buena aproximación siempre y cuando se cumplan condiciones como:\n",
    "\n",
    "* El tamaño de la muestra no sea demasiado pequeño\n",
    "* Los datos muestreados _cubran_ más menos el rango _real_ de la población\n",
    "* La muestra no esté sesgada.\n",
    "\n",
    "Es decir criterios que uno espera, que en general, tenga la mayoría de las muestras con las que uno querría trabajar."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Mostrar como usar bootstrtaping con la regresión lineal?\n",
    "* Mostrar minimos cuadrados como caso particular de MLE?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ejercicios\n",
    "\n",
    "1. graficar dist_monedas.cdf() en el rango de 0 a 20. Luego lo mismo pero para un n más grande  por ejmplo 1000\n",
    "2. Modificar el ejemplo de bootstraping para eliminar el loop y la necesidad de juntar los datos en una lista. Ayuda _size_ puede ser una tupla.\n",
    "3. A partir del ejemplo de bootstaping implementá una versión \"suavizada\" y una paramétrica."
   ]
  }
 ],
 "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.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}