LikelihoodRatioNotes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notes on the likelihood ratio test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the very simple situation in which we have a set of data which we believe is generated by a Poisson distribution.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import poisson\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "sns.set_style('darkgrid')\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\",category=FutureWarning)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our data comes from a Poisson(1) distribution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = poisson(1)\n",
    "sample=X.rvs(1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD7CAYAAABgzo9kAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAGSRJREFUeJzt3V9sVHUe9/H3mdO/004dG8ru9kGwuBpWCLthG/CiICxiiYmrbootkJpY4j4SLNsLSG1XBkxdsGmWC6n8WZK9QV3WWmO4WN2wDQ1bJK0hK4SurIlBXWjFxbbSTv8Nc+a5MPYJ2dIZpmc6099+XlfM6e/8zvebUz49c+acOVYkEokgIiLG8SS7ABERSQwFvIiIoRTwIiKGUsCLiBhKAS8iYigFvIiIoRTwIiKGUsCLiBhKAS8iYqi0ZG7ccRzC4fhupLVtK+51U40pvZjSB6iXVGVKL9PtIz3djmlcUgM+HI4wMDAc17p+vzfudVONKb2Y0geol1RlSi/T7aOgwBfTOJ2iERExlAJeRMRQCngREUPFFPDffPMNDz/8MJ999hlffPEFGzduZNOmTezevRvHcQBobm6mrKyMiooKLly4kNCiRUQkuqgBHwqFCAQCZGVlAbBv3z5qamp46623iEQitLW10d3dTVdXFy0tLezfv5+XX3454YWLiMjUogZ8Y2MjFRUVzJ07F4Du7m6WL18OwKpVq/jwww85d+4cJSUlWJZFYWEh4XCYvr6+xFYuIiJTmjLg3333XfLz81m5cuXEskgkgmVZAOTk5DA4OMjQ0BC5ubkTY75fLiIiyTPldfCtra1YlsXZs2f55JNPqK2tveXIPBgMkpeXR25uLsFg8JblPl/06zRt28Lv98ZVuG174l431ZjSiyl9gHpJVab0MlN9TBnwb7755sS/Kysr2bNnD01NTXR2drJixQpOnz7NQw89xPz582lqamLLli189dVXOI5Dfn5+1I3rRqfvmNKLKX2AeklVpvQyUzc63fGdrLW1tezatYv9+/ezcOFCSktLsW2b4uJiysvLcRyHQCBwxwXfqRujIfrHw67N5023ybRcm05EJOmsSCSStC92CIXCcf8VG7E8vH/+qmu1rFk0l7szYvt+B7fpqCT1qJfUZEov+qoCERGZFgW8iIihFPAiIoZSwIuIGEoBLyJiKAW8iIihFPAiIoZSwIuIGEoBLyJiKAW8iIihFPAiIoZSwIuIGEoBLyJiKAW8iIihFPAiIoZSwIuIGEoBLyJiKAW8iIihFPAiIoaK+tDtcDjMSy+9xOXLl7Ftm3379jE4OMjzzz/PvffeC8DGjRt57LHHaG5upr29nbS0NOrr61m6dGmi6xcRkduIGvCnTp0C4Pjx43R2drJv3z5+8Ytf8Oyzz1JVVTUxrru7m66uLlpaWujt7aW6uprW1tbEVS4iIlOKGvCPPPIIq1evBqCnp4c5c+Zw8eJFLl++TFtbGwsWLKC+vp5z585RUlKCZVkUFhYSDofp6+sjPz8/0T2IiMgkogY8QFpaGrW1tZw8eZLXXnuNa9eusWHDBpYsWcKhQ4d4/fXX8fl8+P3+iXVycnIYHBycMuBt28Lv98ZV+OjgGN7sjLjWnUxWZjr+u7Jcmw/gxmiI4Fg46rjRwTEiVvSPQ3IybfKy0t0oLSFs2xP3/kw16iU1mdLLTPURU8ADNDY2smPHDp5++mmOHz/OD37wAwDWrVtHQ0MDa9euJRgMTowPBoP4fL4p5wyHIwwMDMdVeMTyMDwyHte6kxkdCzEw4Lg2H0D/eJhTl76OOs6bnRFTL2sWzcUZDblRWkL4/d6492eqUS+pyZRepttHQcHU2fq9qIeN7733HkeOHAEgOzsby7J44YUXuHDhAgBnz55l8eLFLFu2jI6ODhzHoaenB8dxdHpGRCSJoh7BP/roo9TV1bF582Zu3rxJfX09P/rRj2hoaCA9PZ05c+bQ0NBAbm4uxcXFlJeX4zgOgUBgJuoXEZHbsCKRSCRZGw+FwnG/TRmxPLx//qprtaxZNJe7M2zX5oPEnKJxu0Y3mfL2GdRLqjKll5Q5RSMiIrOTAl5ExFAKeBERQyngRUQMpYAXETGUAl5ExFAKeBERQyngRUQMpYAXETGUAl5ExFAKeBERQyngRUQMpYAXETGUAl5ExFAKeBERQyngRUQMpYAXETGUAl5ExFBRn8kaDod56aWXuHz5MrZts2/fPiKRCC+++CKWZXH//feze/duPB4Pzc3NtLe3k5aWRn19PUuXLp2JHkREZBJRA/7UqVMAHD9+nM7OzomAr6mpYcWKFQQCAdra2igsLKSrq4uWlhZ6e3uprq6mtbU14Q2IiMjkogb8I488wurVqwHo6elhzpw5tLe3s3z5cgBWrVrFmTNnKCoqoqSkBMuyKCwsJBwO09fXR35+fkIbEBGRyUUNeIC0tDRqa2s5efIkr732GqdOncKyLABycnIYHBxkaGgIv98/sc73y6cKeNu28Pu9cRU+OjiGNzsjrnUnk5WZjv+uLNfmAxj5djSmGj0eK6ZxiajRTbbtiXt/phr1kppM6WWm+ogp4AEaGxvZsWMHTz/9NGNjYxPLg8EgeXl55ObmEgwGb1nu8/mmnDMcjjAwMBxH2RCxPAyPjMe17mRGx0IMDDiuzQcwOh6OqUZvdkZM4xJRo5v8fm/c+zPVqJfUZEov0+2joGDqbP1e1Kto3nvvPY4cOQJAdnY2lmWxZMkSOjs7ATh9+jTFxcUsW7aMjo4OHMehp6cHx3F0ekZEJImiHsE/+uij1NXVsXnzZm7evEl9fT333Xcfu3btYv/+/SxcuJDS0lJs26a4uJjy8nIcxyEQCMxE/SIichtWJBKJJGvjoVA47rcpI5aH989fda2WNYvmcneG7dp8AP3jYU5d+jrquFhP0SSiRjeZ8vYZ1EuqMqWXlDlFIyIis5MCXkTEUAp4ERFDKeBFRAylgBcRMZQCXkTEUAp4ERFDKeBFRAylgBcRMZQCXkTEUAp4ERFDKeBFRAylgBcRMZQCXkTEUAp4ERFDKeBFRAylgBcRMZQCXkTEUFM+kzUUClFfX8/Vq1cZHx9n69at/PCHP+T555/n3nvvBWDjxo089thjNDc3097eTlpaGvX19SxdunQm6hcRkduYMuBPnDiB3++nqamJ/v5+nnrqKbZt28azzz5LVVXVxLju7m66urpoaWmht7eX6upqWltbE168iIjc3pQBv379ekpLSyde27bNxYsXuXz5Mm1tbSxYsID6+nrOnTtHSUkJlmVRWFhIOBymr6+P/Pz8hDcgIiKTmzLgc3JyABgaGmL79u3U1NQwPj7Ohg0bWLJkCYcOHeL111/H5/Ph9/tvWW9wcDBqwNu2hd/vjavw0cExvNkZca07mazMdPx3Zbk2H8DIt6Mx1ejxWDGNS0SNbrJtT9z7M9Wol9RkSi8z1ceUAQ/Q29vLtm3b2LRpE48//jg3btwgLy8PgHXr1tHQ0MDatWsJBoMT6wSDQXw+X9SNh8MRBgaG4yo8YnkYHhmPa93JjI6FGBhwXJsPYHQ8HFON3uyMmMYlokY3+f3euPdnqlEvqcmUXqbbR0FB9HyFKFfRXL9+naqqKnbu3ElZWRkAW7Zs4cKFCwCcPXuWxYsXs2zZMjo6OnAch56eHhzH0ekZEZEkm/II/vDhw9y4cYODBw9y8OBBAF588UX27t1Leno6c+bMoaGhgdzcXIqLiykvL8dxHAKBwIwULyIit2dFIpFIsjYeCoXjfpsyYnl4//xV12pZs2gud2fYrs0H0D8e5tSlr6OOi/UUTSJqdJMpb59BvaQqU3pJiVM0IiIyeyngRUQMpYAXETGUAl5ExFAKeBERQyngRUQMpYAXETGUAl5ExFAKeBERQyngRUQMpYAXETGUAl5ExFAKeBERQyngRUQMpYAXETGUAl5ExFAKeBERQyngRUQMNeUzWUOhEPX19Vy9epXx8XG2bt3Kj3/8Y1588UUsy+L+++9n9+7deDwempubaW9vJy0tjfr6epYuXTpTPYiIyCSmDPgTJ07g9/tpamqiv7+fp556ikWLFlFTU8OKFSsIBAK0tbVRWFhIV1cXLS0t9Pb2Ul1dTWtr60z1ICIik5gy4NevX09paenEa9u26e7uZvny5QCsWrWKM2fOUFRURElJCZZlUVhYSDgcpq+vj/z8/MRWLyIitzVlwOfk5AAwNDTE9u3bqampobGxEcuyJn4+ODjI0NAQfr//lvUGBwejBrxtW/j93rgKHx0cw5udEde6k8nKTMd/V5Zr8wGMfDsaU40ejxXTuETU6Cbb9sS9P1ONeklNpvQyU31MGfAAvb29bNu2jU2bNvH444/T1NQ08bNgMEheXh65ubkEg8Fblvt8vqgbD4cjDAwMx1V4xPIwPDIe17qTGR0LMTDguDYfwOh4OKYavdkZMY1LRI1u8vu9ce/PVKNeUpMpvUy3j4KC6PkKUa6iuX79OlVVVezcuZOysjIAHnzwQTo7OwE4ffo0xcXFLFu2jI6ODhzHoaenB8dxdHpGRCTJpjyCP3z4MDdu3ODgwYMcPHgQgN/+9re88sor7N+/n4ULF1JaWopt2xQXF1NeXo7jOAQCgRkpXkREbs+KRCKRZG08FArH/TZlxPLw/vmrrtWyZtFc7s6wXZsPoH88zKlLX0cdF+spmkTU6CZT3j6DeklVpvSSEqdoRERk9lLAi4gYSgEvImIoBbyIiKEU8CIihlLAi4gYSgEvImIoBbyIiKGifheNSKxujIboHw+7Oqc33SbTcnVKkf8ZCnhxTXAstjt378SaRXPJTOG7d0VSmU7RiIgYSgEvImIoBbyIiKEU8CIihlLAi4gYSgEvImIoBbyIiKEU8CIihoop4M+fP09lZSUA3d3drFy5ksrKSiorK/nLX/4CQHNzM2VlZVRUVHDhwoXEVSwiIjGJeifr0aNHOXHiBNnZ2QD885//5Nlnn6WqqmpiTHd3N11dXbS0tNDb20t1dTWtra2Jq1pERKKKegQ/f/58Dhw4MPH64sWLtLe3s3nzZurr6xkaGuLcuXOUlJRgWRaFhYWEw2H6+voSWriIiEwt6hF8aWkpV65cmXi9dOlSNmzYwJIlSzh06BCvv/46Pp8Pv98/MSYnJ4fBwUHy8/OnnNu2Lfx+b1yFjw6O4c3OiGvdyWRlpuO/K8u1+QBGvh2NqUaPx4ppXCJqdJPb+wSS17Nte+L+3Uw16iX1zFQfd/xlY+vWrSMvL2/i3w0NDaxdu5ZgMDgxJhgM4vP5os4VDkcYGBi+0xIAiFgehkfG41p3MqNjIQYGHNfmAxgdD8dUozc7I6ZxiajRTW7vE0hez36/N+7fzVSjXlLPdPsoKIierxDHVTRbtmyZ+BD17NmzLF68mGXLltHR0YHjOPT09OA4TtSjdxERSaw7PoLfs2cPDQ0NpKenM2fOHBoaGsjNzaW4uJjy8nIcxyEQCCSiVhERuQMxBfy8efN4++23AVi8eDHHjx//rzHV1dVUV1e7W52IiMRNNzqJiBhKAS8iYigFvIiIoRTwIiKGUsCLiBhKAS8iYigFvIiIoRTwIiKGUsCLiBhKAS8iYigFvIiIoRTwIiKGUsCLiBhKAS8iYigFvIiIoRTwIiKGUsCLiBhKAS8iYqiYAv78+fNUVlYC8MUXX7Bx40Y2bdrE7t27cZzvnnjf3NxMWVkZFRUVEw/lFhGR5Ika8EePHuWll15ibGwMgH379lFTU8Nbb71FJBKhra2N7u5uurq6aGlpYf/+/bz88ssJL1xERKYWNeDnz5/PgQMHJl53d3ezfPlyAFatWsWHH37IuXPnKCkpwbIsCgsLCYfD9PX1Ja5qERGJKmrAl5aWkpaWNvE6EolgWRYAOTk5DA4OMjQ0RG5u7sSY75eLiEjypEUfciuP5///TQgGg+Tl5ZGbm0swGLxluc/nizqXbVv4/d47LQGA0cExvNkZca07mazMdPx3Zbk2H8DIt6Mx1ejxWDGNS0SNbnJ7n0DyerZtT9y/m6lGvaSemerjjgP+wQcfpLOzkxUrVnD69Gkeeugh5s+fT1NTE1u2bOGrr77CcRzy8/OjzhUORxgYGI6r8IjlYXhkPK51JzM6FmJgwHFtPoDR8XBMNXqzM2Ial4ga3eT2PoHk9ez3e+P+3Uw16iX1TLePgoLoB9AQR8DX1taya9cu9u/fz8KFCyktLcW2bYqLiykvL8dxHAKBwB0XLCIi7oop4OfNm8fbb78NQFFREW+88cZ/jamurqa6utrd6kREJG660UlExFAKeBERQyngRUQMpYAXETGUAl5ExFAKeBERQyngRUQMpYAXETGUAl5ExFAKeBERQyngRUQMpYAXETGUAl5ExFAKeBERQyngRUQMpYAXETGUAl5ExFAKeBERQ93xM1m/9+STT+Lzfffg13nz5lFeXs7vfvc7bNumpKSEF154wbUiRdwyFoHhUDjquJFvRxkdjz7Om26TablRmYj74gr4sbExAI4dOzax7IknnuDAgQPcc889/PrXv6a7u5vFixe7U6WIS4ZDYU5d+jrqOG92BsMj41HHrVk0l8wM243SRFwX1ymaS5cuMTIyQlVVFc888wwfffQR4+PjzJ8/H8uyKCkp4ezZs27XKiIidyCuI/isrCy2bNnChg0b+Pzzz3nuuefIy8ub+HlOTg7//ve/XStSRETuXFwBX1RUxIIFC7Asi6KiInw+HwMDAxM/DwaDtwT+7di2hd/vjacERgfH8GZnxLXuZLIy0/HfleXafPDdedxYavR4rJjGJaJGN7m9T8D9nv/X9gmAbXvi/n+WakzpZab6iCvg33nnHT799FP27NnDtWvXGBkZwev18uWXX3LPPffQ0dER04es4XCEgYHheEogYnliOkcaq9GxEAMDjmvzAYyOh2OqMdbzvYmo0U1u7xNwv+f/tX0C4Pd74/5/lmpM6WW6fRQU+GIaF1fAl5WVUVdXx8aNG7Esi7179+LxeNixYwfhcJiSkhJ++tOfxjO1iIi4JK6Az8jI4Pe///1/LX/77benXZCIiLhDNzqJiBhKAS8iYigFvIiIoRTwIiKGUsCLiBhKAS8iYigFvIiIoRTwIiKGUsCLiBhKAS8iYigFvIiIoRTwIiKGUsCLiBgq7odui8jscGM0RH8MDxCPlR40Pnso4EUMFxyL7UHjsdKDxmcPnaIRETGUAl5ExFAKeBERQyngRUQM5eqHrI7jsGfPHv71r3+RkZHBK6+8woIFC9zchIiIxMjVgP/b3/7G+Pg4f/7zn/n444959dVXOXTokJubEBEDjUVgOBT9Us6Rb0cZjeGST13K+R1XA/7cuXOsXLkSgJ/97GdcvHjRzelFxFDDodgu5fRmZzA8Mh51nNuXcsb6ByhWntGQa3NNxdWAHxoaIjc3d+K1bdvcvHmTtLTJN5OeblNQ4It7e/937QNxrzsTCoAH/o8/2WXMKO2T1KT9koJ8WQnfhKsfsubm5hIMBideO45z23AXEZHEcjXgly1bxunTpwH4+OOPeeCB1D5qEBExmRWJRCJuTfb9VTSffvopkUiEvXv3ct9997k1vYiI3AFXA15ERFKHbnQSETGUAl5ExFCzLuAdxyEQCFBeXk5lZSVffPFFskualvPnz1NZWZnsMqYlFAqxc+dONm3aRFlZGW1tbckuKW7hcJi6ujoqKirYvHkzX375ZbJLmpZvvvmGhx9+mM8++yzZpUzLk08+SWVlJZWVldTV1SW7nGk5cuQI5eXl/OpXv6KlpSWh25p11zCadLfs0aNHOXHiBNnZ2ckuZVpOnDiB3++nqamJ/v5+nnrqKdauXZvssuJy6tQpAI4fP05nZyf79u2btb9foVCIQCBAVlbir7dOpLGxMQCOHTuW5Eqmr7Ozk3/84x/86U9/YmRkhD/+8Y8J3d6sO4I36W7Z+fPnc+DAgWSXMW3r16/nN7/5zcRr2569D4N45JFHaGhoAKCnp4c5c+YkuaL4NTY2UlFRwdy5c5NdyrRcunSJkZERqqqqeOaZZ/j444+TXVLcOjo6eOCBB9i2bRvPP/88q1evTuj2Zt0R/J3eLZvKSktLuXLlSrLLmLacnBzgu32zfft2ampqklzR9KSlpVFbW8vJkyd57bXXkl1OXN59913y8/NZuXIlf/jDH5JdzrRkZWWxZcsWNmzYwOeff85zzz3HBx98MCv/z/f399PT08Phw4e5cuUKW7du5YMPPsCyEvPFObPuCF53y6am3t5ennnmGZ544gkef/zxZJczbY2Njfz1r39l165dDA8PJ7ucO9ba2sqHH35IZWUln3zyCbW1tfznP/9JdllxKSoq4pe//CWWZVFUVITf75+1vfj9fkpKSsjIyGDhwoVkZmbS19eXsO3NuoDX3bKp5/r161RVVbFz507KysqSXc60vPfeexw5cgSA7OxsLMualaec3nzzTd544w2OHTvGT37yExobGykoKEh2WXF55513ePXVVwG4du0aQ0NDs7aXn//85/z9738nEolw7do1RkZG8PsT9x08s+7Qd926dZw5c4aKioqJu2UluQ4fPsyNGzc4ePAgBw8eBL77AHk2frj36KOPUldXx+bNm7l58yb19fVkZmYmu6z/aWVlZdTV1bFx40Ysy2Lv3r2z9l37mjVr+OijjygrKyMSiRAIBBJ6AKE7WUVEDDXrTtGIiEhsFPAiIoZSwIuIGEoBLyJiKAW8iIihFPAiIoZSwIuIGEoBLyJiqP8HCKp3sHL0JZYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "j=sns.distplot(sample,kde=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'd like to consider the possibility that our data came from a Poisson distribution that has mean of $1.1$. We look at the likelihood ratio of our data with respect to the null hypothesis that the mean is $1$ and the alternative hypothesis that it is $1.1$.\n",
    "\n",
    "In the null situation, the probability $P(X=n)$ is $e^{-1}/n!$. So the log likelihood for the sample is \n",
    "$$\n",
    "\\mathcal{L}_{1}(x)=\\sum_{i=1}^{N} (-\\log(x_i!)-1).\n",
    "$$\n",
    "\n",
    "For a mean of $\\lambda$, we have $P(X=n)=\\lambda^{n}e^{-\\lambda}/n!$ and\n",
    "so the log likelhood is\n",
    "$$\n",
    "\\mathcal{L}_{\\lambda}(x)=\\sum_{i=1}^{N} x_{i}\\log(\\lambda)-\\lambda-\\log(x_{i}!).\n",
    "$$\n",
    "\n",
    "The difference\n",
    "$$\n",
    "\\mathcal{L}_{\\lambda}(x)-\\mathcal{L}_{1}(x) = +N-N\\lambda+\\log(\\lambda)\\sum x_{i}=N(-\\lambda+1+\\log(\\lambda)\\overline{x})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This difference is positive (meaning that the alternative is more likely than the null) when \n",
    "$$\n",
    "\\overline{x}>\\frac{(\\lambda -1)}{\\log(\\lambda)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see what's going on, let's draw 1000 samples from poisson(1) and poisson(1.1) distributions and see how $\\overline{x}$ is distributed in each case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [],
   "source": [
    "N=1000\n",
    "sample_means_1=[poisson(1).rvs(N).mean() for i in range(5000)]\n",
    "sample_means_2=[poisson(1.1).rvs(N).mean() for i in range(5000)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot these distributions, and include a vertical line at the critical value that gives equal likelihood either way as computed above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEFCAYAAADpIfy5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3X10U3WeP/D3vTdJm7QpsUN96EEexfFIZWZKF45OdRgFy3KmWx/otMBWl7KIroRld0WglILWAZUZdoWKD/PgmYMMyIODiHuGI0jlFFn4iQpDXUQRWJFamCmlzUOb5N77+6PT0KRpnpo0ye37dQ5H8um9N59c2ze333zvN4KqqiqIiEhzxEQ3QERE8cGAJyLSKAY8EZFGMeCJiDSKAU9EpFEMeCIijWLAU0zIsow33ngDDz30EEpKSjB9+nSsXbsWLpcr4PYvvfQSdu3aBQCoq6vDvn37etUj9fbbb2P+/PnRvQA/FRUV+NOf/hSTY/35z3/GwoULI95PVVUsWbIEv/3tb/vcpr6+HsXFxSgqKsLChQths9m8X5s0aRJKSkq8f3bv3h1V/5S6dIlugLRh1apVuHr1Kn7/+9/DbDbD4XDgqaeewvLly7F27dpe2//rv/6r9+9HjhzBLbfc0quuFXfccQfWr18f0T5nzpzBM888gxMnTuDWW28NuE1LSwuWLVuGLVu2YOTIkVi7di1++ctfYtWqVfj6669hsVjwzjvvxOIlUIpiwFO/XbhwAe+++y4aGhqQmZkJADCZTHjmmWfwySefAACWLl2K1tZWfPPNN5g8eTL++te/YuzYsUhPT8fJkyfx4osvQpIk7N+/H2PHjsXcuXNx/PhxPPfcc3A6ndDr9Xj66adx5513YseOHXjrrbfgdrtx9epVzJs3D7Nmzeqzv/LycsyZMwdFRUUA4P0H58knn8SqVatw/vx5tLa2IiMjA7/85S8xevRon9dWXFyMTz/9NODj7du3Y8uWLVAUBRaLBStWrMCYMWN8nv/IkSOora3Fnj178PHHH+P555+HoigAgPnz53v76mnz5s0oLS1Fbm5un6+roaEBd9xxB0aOHAkAmDlzJkpKSrBy5Up8+umnEEURs2bNQnt7O4qKivDEE09AkqS+/0eS5nCIhvqtsbERt9xyizfcu+Xk5PiEV0dHB9577z0sXrzYW5s9ezby8vLw9NNPY+rUqd662+3Gk08+iSeffBJ79uxBbW0tVq9eDZvNhu3bt+P111/Hrl278J//+Z8Bf0PoqbS0FG+//TaArqGk3bt3o7S0FAcPHkRWVhbeeust7N27F3l5edi8eXPYr/vo0aPYtWsXNm/ejF27duGf//mfsWDBgqD7bNiwAXPmzMHbb7+N1atX43/+538CbldTU4Pi4uKgx/ruu+9w4403eh/feOONsNlssNvtkGUZd911F37zm99g8+bNaGhowKZNm8J+baQNvIKnfhNF0XtFGsyECRPCPubp06chiiImT54MAMjLy8O7774LAHj11Vfx4Ycf4ty5czh16hQcDkfQY02fPh0vvvgiLl++jM8//xwjR470/rn55puxadMmnD9/HkePHsWPfvSjsHusr6/H+fPnUV5e7q21tbWhtbUVFosl4D5///d/j2effRYffPAB7rrrLvz7v/972M/nT1EUCILQqy6KIn7+85/71ObMmYNNmzbhn/7pn6J+Pko9vIKnfhs/fjy+/vprnzf4AKC5uRmPPfYYOjo6AHQN24RLkqRe4XX69Gl89913eOCBB/Dtt99iwoQJWLRoUchjGY1GFBUVYc+ePdi5cydKS0sBAH/4wx+wfPlypKeno7i4GD/72c/gvzSTIAg+Nbfb7f27oigoKSnBO++8g3feeQd//OMfsXPnTgwZMqTPXsrLy7F79278+Mc/RkNDA/7hH/4BnZ2dYZ0TfzfddBMuXbrkfdzc3IwhQ4bAZDJh165dOHXqlPdrqqpCp+P13GDDgKd+u+GGG1BcXIyqqipvyNtsNqxatQoWiwXp6elB95ckCR6Px6c2evRoCIKAQ4cOAegaBnr00UfxySefIDs7G//yL/+CwsJCHDhwAEDX0EswP//5z/HHP/4Rn3zyiXfYqKGhAQ8++CBKS0sxatQofPDBB72Ok5WVBbfbja+++goA8N5773m/VlhYiPfee88bslu2bMGjjz4atI/y8nL87//+Lx566CHU1taira0Nly9fDrpPXwoLC3H8+HGcO3cOALB161bcd999AIAvv/wS69evhyzL6OjowObNmzF9+vSonodSF/9Jp5hYuXIlNm7ciPLyckiSBJfLhSlTpsBqtYbc995778W6det8ro4NBgM2bNiA1atX48UXX4Rer8eGDRswbtw47N69G9OmTYMgCJg4cSKys7Nx/vz5oM+Rl5cHSZIwbdo0pKWlAQAqKytRU1ODHTt2AAB++MMf4vTp0z77mc1mLF68GPPmzUN2djamTZvm/VphYSHmzZuHyspKCIKAzMxM1NXVBRw26fbUU09h9erV+K//+i8IgoAFCxZg2LBhIc9Rtz//+c+orq7GO++8g+9973tYs2YNFi5cCLfbjeHDh+OFF14AACxYsADPPvssiouL4fF4MG3aNO9vLjR4CFwumIhImzhEQ0SkUQx4IiKNYsATEWkUA56ISKMSOotGURTIcnze45UkIW7Hjjf2PrDOX+m6UWr00IyU671bKp73buw9cnp9eEtOJHQWjdsto7U1+F2I0bJYTHE7dryx94E1/63jAIC35t+Zcr13S8Xz3o29Ry4nxxzWdhyiISLSKAY8EZFGMeCJiDSKAU9EpFEMeCIijWLAExFpFAOeiEijGPBERBrFgCci0ih+4AeRxqXDDtFt61VX9JnoQEYCOqKBEjLgZVlGdXU1zp49C0mSsGbNGrS3t+Pxxx/HyJEjAQAzZ87E9OnTUVdXh/r6euh0OlRVVWH8+PHx7p+IQhDdNiin3+9dv3UqoGfAa1nIgO/+zMutW7fiyJEjWLNmDe69917MmTMHlZWV3u0aGxtx9OhRbN++HU1NTbBardi5c2f8OicioqBCBvyUKVMwefJkAMDFixcxdOhQnDx5EmfPnsX+/fsxYsQIVFVV4dixYygsLIQgCMjNzYUsy2hpaUF2dnafx5YkARaLKWYvxvfYYtyOHW/sfWDpdF1vRaVi792C9S606QGjofcX0vQwZCX+9Wr1vCeDsMbgdTodlixZgvfffx/r169Hc3MzSktLkZeXh1deeQUvv/wyzGYzLBaLd5+MjAy0t7cHDXhZVrmaZADsfWB5PAoAQJaVlOu9W7DzbnK7oThdvepipxuOJHi9qfg9000zq0m+8MIL2Lt3L1asWIHCwkLk5eUBAKZOnYrPP/8cmZmZsNvt3u3tdjvM5vCaICKi2AsZ8Lt27cJrr70GADAajRAEAQsWLMCJEycAAIcPH8a4ceOQn5+PhoYGKIqCixcvQlGUoFfvREQUXyGHaO6//34sW7YMs2fPhsfjQVVVFW666SbU1tZCr9dj6NChqK2tRWZmJgoKClBWVgZFUVBTUzMQ/RMRUR9CBrzJZMJLL73Uq75169ZeNavVCqvVGpvOiIioX3gnKxGRRjHgiYg0iksVUEQ6VcDhln1qJr2ENCFBDRFRnxjwFBGHW8aBU5d8aj+97XqkGaQEdUREfeEQDRGRRvEKnigFBFoRkqtBUigMeKIUEGhFSK4GSaFwiIaISKMY8EREGsWAJyLSKAY8EZFGMeCJiDSKAU9EpFGcJklx4b+kAZczIBp4DHiKC/8lDbicAdHA4xANEZFG8QqegvIfapHVBDZDRBFhwFNQ/kMtd47NSWA3RBQJDtEQEWkUA56ISKM4REOUonSiCpO7GQAgtOlhcru5hDD5CBnwsiyjuroaZ8+ehSRJWLNmDVRVxdKlSyEIAsaOHYuVK1dCFEXU1dWhvr4eOp0OVVVVGD9+/EC8BqKUEOs13QW3A8qZj7oeGA1QnC4uIUw+Qgb8gQMHAABbt27FkSNHvAG/aNEiTJo0CTU1Ndi/fz9yc3Nx9OhRbN++HU1NTbBardi5c2fcXwBRquCa7jTQQgb8lClTMHnyZADAxYsXMXToUNTX12PixIkAgHvuuQeHDh3CqFGjUFhYCEEQkJubC1mW0dLSguzs7D6PLUkCLBZTbF5Jr2OLcTt2vCVT786rHTAZDd7HOkn0eQwA6Wl6WIakA7jWu/9+PbdJNjpd11tR8T7vQpse8Dt3SNPDkBX6OQPtK+gkGP9WE0Wx6+8BjhfweSN47nhLpu/3SCV772GNwet0OixZsgTvv/8+1q9fjwMHDkAQuu47z8jIQHt7O2w2GywWi3ef7nqwgJdlFa2tjn6+hMAsFlPcjh1vydR7h0uGw+nyPvbIis9jAOjodKO1VQFwrXf//Xpuk2w8nq6+ZFmJ63k3ud1Q/M6d2OmGI4znDLRvmkdG599qRqMBTqcr4PEC7RvJc8dbMn2/RypRvefkmMPaLuxZNC+88AL27t2LFStWoLOz01u32+3IyspCZmYm7Ha7T91sDq8JIiKKvZABv2vXLrz22msAAKPRCEEQkJeXhyNHjgAADh48iIKCAuTn56OhoQGKouDixYtQFCXo1TsREcVXyCGa+++/H8uWLcPs2bPh8XhQVVWFMWPGYMWKFVi3bh1Gjx6NoqIiSJKEgoIClJWVQVEU1NTUDET/lCIEQcAVl+xTM+gkuDy+Na46SRQ7IQPeZDLhpZde6lV/8803e9WsViusVmtsOiNNcXoUHP7ysk/tzrE5vWpcdZIodngnKxGRRjHgiYg0iksVECWZQHe86uBG74mORMEx4ImSTKA7XoUxdyWoG0plDHhKKoFm23BmDVF0GPCUVALNtuHMmsSK9SJpNHAY8EQUFBdJS12cRUNEpFEMeCIijWLAExFpFAOeiEijGPBERBrFWTTUbz3nrjuvdqDDJUNWE9wUETHgqf96zl03GQ1wOF24c2xOgruiUHSiCpO72afG+e3awoAnGqQEtwPKmY98apzfri0cgyci0igGPBGRRjHgiYg0imPw5NWpAg6370qOnA2TWgK9ccq15AcvBjx5OdwyDpy65FPjbJjUEuiNU64lP3gx4IkSiFfcFE9BA97tdqOqqgrffvstXC4XnnjiCdx44414/PHHMXLkSADAzJkzMX36dNTV1aG+vh46nQ5VVVUYP378QPRPlNJ4xU3xFDTgd+/eDYvFgrVr1+LKlSt48MEH8eSTT2LOnDmorKz0btfY2IijR49i+/btaGpqgtVqxc6dO+PePBER9S1owE+bNg1FRUXex5Ik4eTJkzh79iz279+PESNGoKqqCseOHUNhYSEEQUBubi5kWUZLSwuys7ODPrkkCbBYTLF5Jb2OLcbt2PGWqN6dVztgMhp8ajpJ9Kn5P/aviaIAk9EQ8X7BaulpeliGpEf/wkLQ6bomk8X7vAttesDvtQk6CcYY1ERRhNFoCHvfPutpehiyfM9BoL4DbRct/qzGT9CAz8jouqPNZrNh4cKFWLRoEVwuF0pLS5GXl4dXXnkFL7/8MsxmMywWi89+7e3tIQNellW0tjpi8DJ6s1hMcTt2vCWq9w6XDIfTd/TXIys+Nf/H/rXupQoi3S9YraPTjdZWJfoXFoLH03VsWVbiet5NbjcUv9eW5pHRGYOa0WiA0+kKe9++6mKnGw6/cxCo70DbRYs/q5HLyTGHtV3IefBNTU145JFHUFJSguLiYkydOhV5eXkAgKlTp+Lzzz9HZmYm7Ha7dx+73Q6zObwGiIgoPoIG/F/+8hdUVlZi8eLFmDFjBgBg7ty5OHHiBADg8OHDGDduHPLz89HQ0ABFUXDx4kUoihLy6p2IiOIr6BDNq6++ira2NmzcuBEbN24EACxduhSrV6+GXq/H0KFDUVtbi8zMTBQUFKCsrAyKoqCmpmZAmicior4FDfjq6mpUV1f3qm/durVXzWq1wmq1xq4ziivetUqkfbzRaZDiXasUCG+80hYGPBF58cYrbeFqkkREGsWAJyLSKA7RDAJ8Q3XgpcMO0W3zqXEsmwYaA34Q4BuqA09026Ccft+nxrFsGmgMeEp6giDgiuvabyAmvYQ0IYENEaUIBjwlPadHweEvL3sf//S265FmkBLYEVFq4JusREQaxYAnItIoBjwRkUYx4ImINIoBT0SkUQx4IiKNYsATEWkUA56ISKMY8EREGsWAJyLSKAY8EZFGcS0aon7i0sCUrBjwRP3EpYEpWQUNeLfbjaqqKnz77bdwuVx44okncMstt2Dp0qUQBAFjx47FypUrIYoi6urqUF9fD51Oh6qqKowfP36gXgMREQUQNOB3794Ni8WCtWvX4sqVK3jwwQdx2223YdGiRZg0aRJqamqwf/9+5Obm4ujRo9i+fTuamppgtVqxc+fOgXoN5Mf/E5z46U1Eg1PQgJ82bRqKioq8jyVJQmNjIyZOnAgAuOeee3Do0CGMGjUKhYWFEAQBubm5kGUZLS0tyM7ODvrkkiTAYjHF4GUEOrYYt2PHW397b7ragSPn/+p9PGHEdTAZDT7b6CQxqlqobURRgMloiHi/SGrpaXpYhqQHfvFR0Om65hpEe96FNj3g17Ogk2AcwJooijAaDWHv29/nRpoehqzY/HwN5p/VeAsa8BkZGQAAm82GhQsXYtGiRXjhhRcgCIL36+3t7bDZbLBYLD77tbe3hwx4WVbR2uro72sIyGIxxe3Y8dbf3jtcMhzOa2/xeWTF53F/aqG2MRkNcDhdEe8XSa2j043WViXwi4+Cx9N1LFlWojrvJrcbil/PaR4ZnQNYMxoNcDpdYe/b3+cWO91wxOjnazD/rEYrJ8cc1nYhp0k2NTXhkUceQUlJCYqLiyGK13ax2+3IyspCZmYm7Ha7T91sDq8BItKGdNhhcjf7/EmHPfSOFDdBA/4vf/kLKisrsXjxYsyYMQMAcPvtt+PIkSMAgIMHD6KgoAD5+floaGiAoii4ePEiFEUJefVORNrSPZuo5x//6aM0sIIO0bz66qtoa2vDxo0bsXHjRgDA8uXL8dxzz2HdunUYPXo0ioqKIEkSCgoKUFZWBkVRUFNTMyDNExFR34IGfHV1Naqrq3vV33zzzV41q9UKq9Uau86IiKhfeKMTpRxBEHDFJfvUTHoJaUL8n5t3rVIqYcBTynF6FBz+8rJP7ae3XY80gxT35+Zdq5RKuNgYEZFGMeCJiDSKAU9EpFEMeCIijWLAExFpFGfRpDj/lSMBrh5JRF0Y8CnO4ZZx4NQln9qdY3MS1A0RJRMO0RARaRQDnohIoxjwREQaxYAnItIoBjwRkUYx4ImINIrTJEkTErmE8GCkE1WY3M2+NS6bnHQY8KQJiVxCeDAS3A4oZz7yrXHZ5KTDIRoiIo3iFTxplv+wDYdsaLBhwJNm+Q/bcMiGBhsO0RARaVRYAX/8+HFUVFQAABobG3H33XejoqICFRUV+O///m8AQF1dHWbMmIHy8nKcOHEifh0TEVFYQg7R/PrXv8bu3bthNBoBAJ9//jnmzJmDyspK7zaNjY04evQotm/fjqamJlitVuzcuTN+XRMRUUghr+CHDx+ODRs2eB+fPHkS9fX1mD17NqqqqmCz2XDs2DEUFhZCEATk5uZClmW0tLTEtXEiIgou5BV8UVERLly44H08fvx4lJaWIi8vD6+88gpefvllmM1mWCwW7zYZGRlob29HdnZ20GNLkgCLxdSP9oMdW4zbseMtkt6dVztgMhp8ajpJ9Kn5P+5PLdQ2oijAZDQktIe+aulpeliGpMOfTtd1nRPOeRfa9IDf8wg6CcYE10RRhNFoCHvfgeoRaXoYsoKf08Hys5oIEc+imTp1KrKysrx/r62txX333Qe73e7dxm63w2w2hzyWLKtobXVE2kJYLBZT3I4db5H03uGS4XD63j/okRWfmv/j/tRCbWMyGuBwuhLaQ1+1jk43WlsV+PN4umqyrIQ87ya3G4rf86R5ZHQmuGY0GuB0usLed6B6FDvdcIQ4p4PlZzWWcnJC5ysQxSyauXPnet9EPXz4MMaNG4f8/Hw0NDRAURRcvHgRiqKEvHonIqL4ivgKftWqVaitrYVer8fQoUNRW1uLzMxMFBQUoKysDIqioKamJh69EhFRBMIK+GHDhmHbtm0AgHHjxmHr1q29trFarbBarbHtjoiIosYbnYiINIoBT0SkUQx4IiKNYsATEWkUA56ISKMY8EREGsWAJyLSKH7gB1Ef0mGH6Lb51PjB0pRKGPBEfRDdNiin3/ep8YOlKZUw4FNIpwo43LJPTVYT1AwRJT0GfApxuGUcOHXJp3bn2JwEdUNEyY4Bn8T8r9h5tR4/ImQIbRdgcru9NY63U6pjwCcx/yt2Xq3Hj6DKwFf7fNZ653g7pTpOkyQi0igGPBGRRjHgiYg0igFPRKRRfJOViOJGJ6owuZt9aoo+Ex3ISFBHgwsDngYNQRBwxeV7o5hJLyWom8FBcDugnPnIp2a4bQpE5doSEEKbHulIY+jHAQOeBg2nR8HhLy/71H562/UJ6mbw6hX6RgPEm38C6BnwscYxeCIijQor4I8fP46KigoAwPnz5zFz5kzMmjULK1euhKIoAIC6ujrMmDED5eXlOHHiRPw6JiKisIQM+F//+teorq5GZ2cnAGDNmjVYtGgR/vCHP0BVVezfvx+NjY04evQotm/fjnXr1uGZZ56Je+NERBRcyIAfPnw4NmzY4H3c2NiIiRMnAgDuuecefPTRRzh27BgKCwshCAJyc3MhyzJaWlri1zUREYUU8k3WoqIiXLhwwftYVVUIggAAyMjIQHt7O2w2GywWi3eb7np2dnbQY0uSAIvFFG3vIY4txu3Y8dbdu/NqB0xGg7euk0Sfx+HWot0vmmOJogCT0ZDQHiLZLz1ND51OBGQBoijC2OPrgk7yeZzMte7ew903Wfru7j0tTQ9DVur9vCZ7zkQ8i0YUr1302+12ZGVlITMzE3a73aduNptDHkuWVbS2OiJtISwWiylux4637t47XDIcPRa/8siKz+Nwa9HuF82xTEYDHE5XQnuIZL+OTjc8HgWSqkJRFDh7fD3NI6PTb/tkrRmNBjidrrD3TZa+u3vv7HTDkYI/r4nKmZyc0PkKRDGL5vbbb8eRI0cAAAcPHkRBQQHy8/PR0NAARVFw8eJFKIoS8uqdiIjiK+Ir+CVLlmDFihVYt24dRo8ejaKiIkiShIKCApSVlUFRFNTU1MSjVyIiikBYAT9s2DBs27YNADBq1Ci8+eabvbaxWq2wWq2x7Y6IiKLGG52IiDSKSxUQUcJxUbL4YMATUcIFWpRMvHUq16fpJw7REBFpFAOeBjVBEOBRVSgqcNXpgcOjwqMmuiui2GDA06Dm9ChodbjR4Zbx5aV2fPFdG1wyE560gQFPRKRRfJM1ATpVwOHu/clCaUKCGhpkbv+eigzVCQC4Qb0Ms04GP9iJtIgBnwAOt4wDpy751H562/VIMzBlBkKG6oSt8U8AAPf1Zsj2oTCYr7u2gQA4PCokteu/AGCQBOj4DzClGA7REPlxyyq++K4NVx0ufPFdG8flKWUx4ImINIpDNElCEARccXWNyzuvdqDDJYMXjf3nP97+d9kODDWpsCW4L6KBwIBPEk6PgsNfXgZwbU31O8fmJLir1Oc/3m671I6bCu5LcFdEA4NDNEREGsWAJyLSKAY8EZFGMeCJwvG3ufEeFVyvhlIGA54oDP5z4zkvnlIBA56ISKMY8EREGsWAJyLSqKhvdHrggQdgNpsBAMOGDUNZWRl+8YtfQJIkFBYWYsGCBTFrkigcPe9aBbruXO3kXauakg47RLfv/1F+dmvfogr4zs5OAMCmTZu8tZKSEmzYsAE333wzHnvsMTQ2NmLcuHGx6ZIoDD3vWgW67lzVD5+YwI4o1kS3Dcrp931r/OzWPkUV8KdOnYLT6URlZSU8Hg+sVitcLheGDx8OACgsLMThw4cZ8EQUNZ2owuRu9q3BDVeC+klFUQV8eno65s6di9LSUpw7dw7z5s1DVlaW9+sZGRn45ptvQh5HkgRYLKZoWgjj2GLcjh2Jtg437J2+H+4h6FSYjAafmk4SvTVRFGAyGnxq/ttEUot2v2iOFeveI+kBegkGw7VvaVEUIEqCtyaKXX/3rwmiAEHoWvDNYND1uV3Pmk4nwWjs+pqgk2D062ega6Iowmg0hL1vsvTd3XvA7dROSN8c8a0NnwTJ/7Wk6WHISszPerLkTF+iCvhRo0ZhxIgREAQBo0aNgtlsRmtrq/frdrvdJ/D7IssqWlsd0bQQksViituxI3HF1fvDPe4cmwOH0/c6xCMr3lr3YmM9a/7bRFKLdr9ojhXr3iPpweOW4XJ5vDVFUSHIqremKF1/V/xqqqJCVQFV/dvX+9iuZ83jkeF0KgCANI+MTr9+BrpmNBrgdLrC3jdZ+u7uXenH8cRONxwJ+llPVM7k5JjD2i6qWTQ7duzA888/DwBobm6G0+mEyWTC//3f/0FVVTQ0NKCgoCCaQxMRUYxEdQU/Y8YMLFu2DDNnzoQgCFi9ejVEUcRTTz0FWZZRWFiIH/zgB7HulYiol0Bj9ZxZ0yWqgDcYDPjVr37Vq75t27Z+N0REFAnB7YBy5iOfGmfWdOEHfhBFQ7j2gdzdH87ND+amZMM7WYmi0L34GBcgo2TGgCci0igO0cRYpwo43NfmvfOijogShQEfYw6377x3fnD2IMJxeUoyDHhKScm4sJhbVvH1pXYAwNhhLnz5XRu+f2MWdEx4ShAGPKUkLixGFBrfZCUi0igGPBGRRnGIph/8Z8wAnDVDlAy4fEEXBnw/+M+YAThrJh6ydU78Xfa1FfuS4Q1VSm5cvqALA56Snt5tS903VP2mTnpUcNokDRgGPFEc+U+dFGWV0yZpwDDgKan4D8cAgFHs/WlEKavHFT3Aq/qBFGhcXtLrIbvdPjUtjdUz4Cmp+A/HAIB497QEdRN7Pa/oAV7VD6RA4/L6MXfBreGxek6TJCLSKF7Bh4lTIoko1TDg+xBoVciDX3BKZLQCja3fbLgKj5bH2yklaWkOPQMefV+d9wx0hnn/BBpb12dP0/Tdh09tAAAHX0lEQVR4O6UmLc2hH3QBH06YAwx0GkABZtZ0qoAsc7YN9c+gC3jefRp/t39PxQ3qZZ8hGQ699C3QzBrZrXC2DfVbTANeURSsWrUKX3zxBQwGA5577jmMGDEilk8RVM+rc+fVDogqkMafh7jpXpO9Z5jr9J2w6D1wf/ERbD0CikMvlMr6GpcHTIlpKEwxDfh9+/bB5XLhrbfewmeffYbnn38er7zySiyfIqieV+cmowGTRliQZpAG7Pm1LFCYDzWqOPf/9sJ9vdkb5gaDDjnjfwI52MEoOn5DOXpce+xyeuDxqNAnqDWt63NcHr6//afDDtHtu0pSIm+mimnAHzt2DHfffTcA4Ic//CFOnjwZy8P74LTFyPX8FKTuoM7KNKLNFroWKMxvKrgvMS9kkPIfyskfq+KL79oAdP3D6nJ5UHBr7/0U+I3xSwJkWfV+rCAAfrRgFHSiCqHtAkw9wlsHN1yn6322S+TNVIKqqjGLxeXLl+P+++/HT37yEwDA5MmTsW/fPuh0g26on4go4WJ6J2tmZibsdrv3saIoDHciogSJacDn5+fj4MGDAIDPPvsMt94a4PdFIiIaEDEdoumeRXP69GmoqorVq1djzJgxsTo8ERFFIKYBT0REyYOrSRIRaRQDnohIoxjwREQalXIBrygKampqUFZWhoqKCpw/f97n66+//jpKSkowe/ZsHDhwAADQ0tKCyspKzJo1C4sWLYLT6UxE61H13traikmTJqGiogIVFRX4/e9/n4jWvY4fP46Kiope9Q8++AAPP/wwysrKsG3bNgBAR0cHrFYrZs2ahXnz5qGlpWWg2/URSe+qquLuu+/2nvdf/epXA92uj756BwCn04ny8nKcOXMGQOjvs0SIpH8AeOCBB7znftmyZQPVZkB99b5nzx6UlpaivLwcNTU1UBQl+c69mmL27t2rLlmyRFVVVf3000/Vxx9/3Pu1U6dOqcXFxWpHR4fa0dGhPvDAA6rD4VBra2vVnTt3qqqqqq+99pr6xhtvJKL1qHo/dOiQ+uyzzyakX3+vv/66+rOf/UwtLS31qbtcLnXKlClqa2ur2tnZqT700EPqpUuX1N/97nfq+vXrVVVV1T179qi1tbWJaFtV1ch7P3funDp//vwEdeurr95VVVVPnDihPvjgg+pdd92lfvXVV6qqBv8+S4RI++/o6FBLSkoGus2A+urd6XSq9913n+pwOFRVVdV/+7d/U/ft25d05z7lruCDLYdw5swZTJw4EWlpaUhLS8OIESPwxRdf+Oxzzz334KOPPgp47GTs/eTJk2hsbMQ//uM/YuHChbh06VJfh4+74cOHY8OGDb3qZ86cwfDhwzFkyBAYDAZMmDABH3/8ca/zfvjw4YFu2SvS3hsbG9Hc3IyKigrMmzcPX3/9dQK67tJX7wDgcrnw8ssvY/To0d7aQC4ZEo5I+z916hScTicqKyvxyCOP4LPPPhuoVnvpq3eDwYCtW7fCaDQCADweD9LS0pLu3KdcwNtsNmRmZnofS5IEj8cDAPj+97+Pjz/+GDabDVeuXMGnn34Kp9MJm80Gs9kMAMjIyEB7e3vAYydj76NHj8bChQvx5ptvYsqUKXjuuecS0jsAFBUVBbwzuef5BbrOsc1mS5rzDkTee05ODh577DFs2rQJ8+fPx+LFiweyXR999Q4AEyZMwE033eRTC/Z9lgiR9p+eno65c+fit7/9LZ555hk89dRTCeu/r95FUcTQoUMBAJs2bYLD4cCPf/zjpDv3KbeOQLDlEMaMGYPZs2dj3rx5GDFiBH7wgx/guuuu8+6Tnp4Ou92OrKyslOn9jjvu8F4lTJ06FevXr09I78H4vy673Q6z2exTT+R5D6av3vPy8iBJXSuRFhQUoLm5GaqqQhCSf0WuVF8yZNSoURgxYgQEQcCoUaNgsVhw+fLlXv8QJJqiKFi7di3Onj2LDRs2QBCEpDv3KXcFH2w5hJaWFly5cgVbtmzB8uXL0dTUhLFjxyI/Px8ffvghAODgwYOYMGFCyvReXV2NvXv3AgAOHz6McePGJaT3YMaMGYPz58+jtbUVLpcLH3/8MX70ox8lzXkPpq/e6+rqvG9onzp1Crm5uSkR7kDqLxmyY8cOPP/88wCA5uZm729UyaampgadnZ3YuHGj9yIs2c596vyz/jdTp07FoUOHUF5e7l0O4Y033sDw4cNx77334sKFC3j44Yeh1+vx9NNPQ5IkPPHEE1iyZAm2bduG6667LmEzIqLp/T/+4z9QVVWFLVu2wGg0JnSIxt+7774Lh8OBsrIyLF26FHPnzoWqqnj44Ydxww03YObMmViyZAlmzpwJvV6f8JkoPYXq/bHHHsPixYvx4YcfQpIkrFmzJtEte/XsPZBA32fJJFT/M2bMwLJlyzBz5kwIgoDVq1cnzW8g3b3n5eVhx44dKCgowKOPPgoAeOSRR5Lu3HOpAiIijUq5IRoiIgoPA56ISKMY8EREGsWAJyLSKAY8EZFGMeCJiDSKAU9EpFH/H4yYrfgH3tmgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax=sns.distplot(sample_means_1,kde=False)\n",
    "sns.distplot(sample_means_2,kde=False,ax=ax)\n",
    "s=.1/np.log(1.1)\n",
    "j=ax.axvline(s)\n",
    "t=ax.set_title(\"Critical value is {0:.2f}\".format(s))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a check, we can compute the \"errors\" by looking at how often a sample from poisson(1) lies to the right of the critical line, and vice versa."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fraction of draws from poisson(1.1) that lie to the left of the critical line: 0.068\n",
      "Fraction of draws from poisson(1) that lie to the right of the critical line: 0.0642\n"
     ]
    }
   ],
   "source": [
    "print('Fraction of draws from poisson(1.1) that lie to the left of the critical line:',len([x for x in sample_means_2 if x<s])/len(sample_means_2))\n",
    "print('Fraction of draws from poisson(1) that lie to the right of the critical line:',len([x for x in sample_means_1 if x>s])/len(sample_means_1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The idea behind a likelihood ratio test is that we have a parameter space $\\Omega$ and a random variable $X$ that is distributed as $P_{\\theta}$ for some $\\theta\\in\\Omega$.  We have a partition of $\\Omega$ into sets $\\Omega_{0}$ and $\\Omega_{1}$ and based on the outcome of an experiment $X$ we want to accept the hypothesis that $\\theta\\in\\Omega_{0}$ if $X\\not\\in S$.  The power function of this test is \n",
    "$$\n",
    "\\beta(\\theta)=P_{\\theta}(X\\in S)\n",
    "$$\n",
    "and the significance level is\n",
    "$$\n",
    "\\alpha = \\sup_{\\theta\\in\\Omega_{0}}P_{\\theta}(X\\in S).\n",
    "$$\n",
    "\n",
    "The power function tells us, for a given $\\theta$, how likely it is that that we will reject $H_{0}$.\n",
    "The significance level tells us the highest possible chance that we reject $H_{0}$.  \n",
    "\n",
    "In the example above the parameter space has only two values: $\\theta = 1$ and $\\theta = 1.1$.\n",
    "Since $S$ is the space where we accept $H_{1}$ and reject $H_{0}$, these two values are:\n",
    "$$\n",
    "\\beta(1)=\\textrm{ the area under the blue graph to the right of $s$}\\sim .06\n",
    "$$\n",
    "and\n",
    "$$\n",
    "\\beta(1.1)=\\textrm{ the area under the orange graph to the right of $s$}\\sim .94\n",
    "$$\n",
    "\n",
    "The significance level is $.06$.\n",
    "\n",
    "In other words: Suppose we draw $1000$ samples from our Poisson distribution and compute the mean $\\overline{x}$.  If $\\overline{x}<s$ we accept the null hypothesis and say the mean is $1$.  This choice will be incorrect $6\\%$ of the time.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One thing to notice is that the critical value we are using is chosen to make the two types of errors equally likely: that is, the chance of picking $\\theta=1$ when it's really $\\theta=1.1$ is the same as the reverse.  It doesn't control what that chance actually is (i.e. the significance level).\n",
    "\n",
    "Here's where there's a dependence on $N$.  Suppose we only drew $50$ samples.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEFCAYAAADpIfy5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3XlwFGX+P/B3d88MuZlvJIL5QTDhcNcFVmOEUsPxXcRkVUSRkBAXXUFwXQjLrijIEaIih7jsKhS6uJZlBUROAbFKCxBMBTF8UQ7JCnERsirnAiGZI3N1//4IM2buHDOZmc77VWVJPtMz83kyyXs6z3Q/LSiKooCIiFRHjHQDREQUHgx4IiKVYsATEakUA56ISKUY8EREKsWAJyJSKQY8hYTD4cC7776LsWPHYsyYMbj//vuxfPlyWK1Wn9u//vrr2LZtGwBg1apV2L17t1e9tbZu3Yqnn366bQPwMHHiRHzyyScheaxvvvkGM2bMaPX9FEXB7Nmz8c477/jdZt++fRg9ejTy8vIwY8YMGAwG121DhgzBmDFjXP/t2LGjTf1T7NJEugFSh7KyMly7dg3vvfcekpOTYTKZMGvWLMybNw/Lly/32v5Pf/qT699VVVXo27evV10tBg4ciDfeeKNV9zl16hRefPFFHDt2DP379/e5zZUrV/DCCy9g/fr1uPnmm7F8+XK89tprKCsrw/fffw+9Xo/t27eHYggUoxjw1G4//vgjPvroI1RWViIpKQkAkJCQgBdffBFff/01AGDOnDmoq6vDDz/8gBEjRuDy5cvo168f4uLicPz4cbz66quQJAl79uxBv379MHnyZBw9ehSLFi2C2WyGVqvF888/j7vuugubN2/Ghg0bYLPZcO3aNUyZMgXFxcV++ysqKsKTTz6JvLw8AHC94UybNg1lZWWora1FXV0dEhMT8dprryErK8ttbKNHj8bhw4d9fr1p0yasX78esixDr9djwYIF6NOnj9vzV1VV4eWXX8bOnTtx6NAhLF26FLIsAwCefvppV1/NrVu3DgUFBUhPT/c7rsrKSgwcOBA333wzAGDChAkYM2YMFi5ciMOHD0MURRQXF6OhoQF5eXl45plnIEmS/xeSVIdTNNRu1dXV6Nu3ryvcndLS0tzCq7GxER9//DGee+45V+2xxx7DgAED8Pzzz2PUqFGuus1mw7Rp0zBt2jTs3LkTL7/8MhYvXgyDwYBNmzZhzZo12LZtG/72t7/5/AuhuYKCAmzduhVA01TSjh07UFBQgIqKCqSkpGDDhg349NNPMWDAAKxbt67F4z548CC2bduGdevWYdu2bXjqqacwffr0gPdZuXIlnnzySWzduhWLFy/Gl19+6XO70tJSjB49OuBjnT9/Hj169HB93aNHDxgMBhiNRjgcDtx999345z//iXXr1qGyshLl5eUtHhupA/fgqd1EUXTtkQZyxx13tPgxa2pqIIoiRowYAQAYMGAAPvroIwDAW2+9hc8//xxnzpzBiRMnYDKZAj7W/fffj1dffRWXLl3Cv/71L9x8882u/3r16oXy8nLU1tbi4MGDuP3221vc4759+1BbW4uioiJXrb6+HnV1ddDr9T7v89vf/hYvvfQSPvvsM9x99934y1/+0uLn8yTLMgRB8KqLoojx48e71Z588kmUl5fj97//fZufj2IP9+Cp3QYNGoTvv//e7QM+ALhw4QKmTp2KxsZGAE3TNi0lSZJXeNXU1OD8+fN4+OGH8dNPP+GOO+7AzJkzgz5WfHw88vLysHPnTmzZsgUFBQUAgPfffx/z5s1DXFwcRo8ejQcffBCeSzMJguBWs9lsrn/LsowxY8Zg+/bt2L59Oz788ENs2bIFXbt29dtLUVERduzYgXvuuQeVlZV46KGHYLFYWvQ98XTTTTfh4sWLrq8vXLiArl27IiEhAdu2bcOJEydctymKAo2G+3OdDQOe2q179+4YPXo05s6d6wp5g8GAsrIy6PV6xMXFBby/JEmw2+1utaysLAiCgP379wNomgZ64okn8PXXXyM1NRV//OMfkZubi7179wJomnoJZPz48fjwww/x9ddfu6aNKisr8cgjj6CgoACZmZn47LPPvB4nJSUFNpsN//73vwEAH3/8seu23NxcfPzxx66QXb9+PZ544omAfRQVFeHbb7/F2LFj8fLLL6O+vh6XLl0KeB9/cnNzcfToUZw5cwYA8MEHH2DkyJEAgO+++w5vvPEGHA4HGhsbsW7dOtx///1teh6KXXxLp5BYuHAhVq9ejaKiIkiSBKvVinvvvRclJSVB7/ub3/wGK1ascNs71ul0WLlyJRYvXoxXX30VWq0WK1euxK9+9Svs2LED+fn5EAQBgwcPRmpqKmprawM+x4ABAyBJEvLz89GlSxcAwKRJk1BaWorNmzcDAG677TbU1NS43S85ORnPPfccpkyZgtTUVOTn57tuy83NxZQpUzBp0iQIgoCkpCSsWrXK57SJ06xZs7B48WL8/e9/hyAImD59Onr27Bn0e+T0zTffYP78+di+fTtuuOEGLFmyBDNmzIDNZkNGRgaWLVsGAJg+fTpeeukljB49Gna7Hfn5+a6/XKjzELhcMBGROnGKhohIpRjwREQqxYAnIlIpBjwRkUpF9CgaWZbhcETvZ7ySJER1f63RWcdSe7XpJKje/9PyY/A7klpeF7WMA4iNsWi1LVtyIqJH0dhsDtTVBT4LMZL0+oSo7q81OutYnt5wFADwj8Jfh7OlNlPL66KWcQCxMZa0tOQWbccpGiIilWLAExGpFAOeiEilGPBERCrFgCciUikGPBGRSjHgiYhUigFPRKRSDHgiIpXiBT+o04uDEaLN4FWXtUloRGIEOiIKDQY8dXqizQC5Zpd3vf8oQMuAp9jFKRoiIpViwBMRqRQDnohIpRjwREQqxYAnIlIpBjwRkUox4ImIVIoBT0SkUgx4IiKVYsATEalUi5YqePjhh5Gc3HQV7549e6KwsBCvvPIKJElCbm4upk+fDlmWUVZWhpMnT0Kn02HRokXo3bt3WJsnIiL/gga8xWIBAJSXl7tqY8aMwcqVK9GrVy9MnToV1dXV+Omnn2C1WrFhwwYcOXIES5cuxZtvvhm+zomIKKCgAX/ixAmYzWZMmjQJdrsdJSUlsFqtyMjIAADk5ubiwIEDuHTpEoYOHQoAuO2223D8+PHwdk5ERAEFDfi4uDhMnjwZBQUFOHPmDKZMmYKUlBTX7YmJifjhhx9gMBiQlJTkqkuSBLvdDo3G/1NIkgC9PqGdQwgfSRKjur/W6Kxj0WiaPmYKtL1QrwXidd43dNFClxLe75laXhe1jANQ11iCBnxmZiZ69+4NQRCQmZmJ5ORk1NXVuW43Go1ISUlBY2MjjEajqy7LcsBwBwCHQ0Fdnakd7YeXXp8Q1f21Rmcdi90uA0DA7RNsNshmq1ddtNhgCvP3TC2vi1rGAcTGWNLSklu0XdCjaDZv3oylS5cCAC5cuACz2YyEhAT85z//gaIoqKysRE5ODrKzs1FRUQEAOHLkCPr379+O9omIqL2C7sGPGzcOL7zwAiZMmABBELB48WKIoohZs2bB4XAgNzcXv/71rzFw4EDs378fRUVFUBQFixcv7oj+iYjIj6ABr9Pp8Ne//tWrvnHjRrevRVHESy+9FLrOiIioXXiiExGRSjHgiYhUigFPRKRSDHgiIpViwBMRqRQDnohIpRjwREQqxYAnIlIpBjwRkUox4ImIVIoBT0SkUgx4IiKVYsATEakUA56ISKUY8EREKsWAJyJSKQY8EZFKMeCJiFSKAU9EpFIMeCIilWLAExGpFAOeiEilGPBERCrFgCciUikGPBGRSjHgiYhUigFPRKRSDHgiIpViwBMRqRQDnohIpVoU8JcvX8bw4cNx6tQp1NbWYsKECSguLsbChQshyzIAYNWqVRg3bhyKiopw7NixsDZNRETBBQ14m82G0tJSxMXFAQCWLFmCmTNn4v3334eiKNizZw+qq6tx8OBBbNq0CStWrMCLL74Y9saJiCiwoAG/bNkyFBUV4cYbbwQAVFdXY/DgwQCAYcOG4YsvvsBXX32F3NxcCIKA9PR0OBwOXLlyJbydExFRQJpAN27duhWpqakYOnQo1qxZAwBQFAWCIAAAEhMT0dDQAIPBAL1e77qfs56amhrwySVJgF6f0N4xhI0kiVHdX2t01rFoNE37MIG2F+q1QLzO+4YuWuhSwvs9U8vropZxAOoaS8CA37JlCwRBwIEDB/Dtt99i9uzZbnvmRqMRKSkpSEpKgtFodKsnJycHfXKHQ0Fdnakd7YeXXp8Q1f21Rmcdi93e9BlRoO0TbDbIZqtXXbTYYArz90wtr4taxgHExljS0oLnKxBkimbdunVYu3YtysvL8ctf/hLLli3DsGHDUFVVBQCoqKhATk4OsrOzUVlZCVmWcfbsWciyHHTvnYiIwivgHrwvs2fPxoIFC7BixQpkZWUhLy8PkiQhJycHhYWFkGUZpaWl4eiViIhaocUBX15e7vr32rVrvW4vKSlBSUlJaLoiIqJ244lOREQqxYAnIlIpBjwRkUox4ImIVIoBT0SkUq0+TJKIwicORog2g1dd1iahEYkR6IhiGQOeKIqINgPkml3e9f6jAC0DnlqHUzRERCrFgCciUikGPBGRSjHgiYhUigFPRKRSDHgiIpViwBMRqRQDnohIpRjwREQqxTNZyY1FAUw2h1stQSuhixChhoiozRjw5MZkc2DviYtutf/9xY3oopMi1BERtRWnaIiIVIoBT0SkUgx4IiKVYsATEakUA56ISKUY8EREKsWAJyJSKQY8EZFKMeCJiFSKAU9EpFIMeCIilWLAExGpFBcbo7CKgxGizeBVl7VJaERiBDoi6jyCBrzD4cD8+fNx+vRpSJKEJUuWQFEUzJkzB4IgoF+/fli4cCFEUcSqVauwb98+aDQazJ07F4MGDeqIMVAUE20GyDW7vOv9RwFaBjxROAUN+L179wIAPvjgA1RVVbkCfubMmRgyZAhKS0uxZ88epKen4+DBg9i0aRPOnTuHkpISbNmyJewDICIi34IG/L333osRI0YAAM6ePYtu3bph3759GDx4MABg2LBh2L9/PzIzM5GbmwtBEJCeng6Hw4ErV64gNTU1rAOgzqH5hUjM1xrRaHXwQiREQbRoDl6j0WD27NnYtWsX3njjDezduxeC0PSblZiYiIaGBhgMBuj1etd9nPVAAS9JAvT6hHYOIXwkSYzq/lqjpWMxX2tEQrzOrRbXRQt917g2Pa9QrwU8Hg8A0EULXUrLv7fnrjWiqvYyAEAUBciyguH904L2pdE0HUcQaOyh6rEtPF+XSPbSHp3xdyUWtPhD1mXLlmHWrFkYP348LBaLq240GpGSkoKkpCQYjUa3enJycsDHdDgU1NWZ2tB2x9DrE6K6v9Zo6VgarQ6YzFb3msWGujq5Tc+bYLNB9ng8ABAtNpha8b1t3ldCvA4ms7VFfdntTbcHGnuoemwLz9clkr20R2f8XYmktLTA2eoU9DDJbdu24R//+AcAID4+HoIgYMCAAaiqqgIAVFRUICcnB9nZ2aisrIQsyzh79ixkWeb0DEUNiwJctTrc/rMoke6KKLyC7sHfd999eOGFF/DYY4/Bbrdj7ty56NOnDxYsWIAVK1YgKysLeXl5kCQJOTk5KCwshCzLKC0t7Yj+iVqE15qlzihowCckJOD111/3qq9du9arVlJSgpKSktB0RkRE7cIzWYmIVIoBT0SkUgx4IiKV4lo0FDHNT15yCufJS7feoCBRMbu+vkG+hDibAA1s8D4wkSj2MeApYjr6yJZExQxD9Seur+09UiBrBAh97g7L8xFFGqdoiIhUinvwpCqe0z52pelsJgdPaqJOiAFPquI57VNnsgEAbDITnjofBjzFLEEQcNXq/iGtQ3H/MHWHpun27solWBIUeF96hEi9GPDUIewKYG02T6JxKHC082gZs13Gge8uudXu6pfm9mGqw5gJALCdrIE2Y3D7npAoxjDgqUNYHQpOnq93fZ10gxGaG7pFsCMi9eNRNEREKsWAJyJSKU7RUEh5HqYoygrsdgVtu1wIEbUHA55CyvMwxTtTTTCcr0fWjS27Ag0RhQ4DnqKKr0Mfgdg4USkORog29wMxZW0SGpEYoY6os2PAU1Txdegj0HT4Y7QTbQbINbvca/1HAVoGPEUGP2QlIlIp7sET+aERFSTYLnjVOe1CsYIBT+SHYDNBPvWFV53TLhQrGPBEYcS/AiiSGPDUZr6uyBQLR7u4CIDJrkBSmv7vpJMEaEJ0VSn+FUCRxICnNvN1RaZYONrFyeZQ8P3FBvTracV3zdbJuaVHCjQBEt7fXjkv/UfRhgFPQcXysenh4G+vnJf+o2jDgKegovnYdM8LaXPdd6KfMeAppnleSNt2YzLXfSe6jic6ERGpFAOeiEilGPBERCrFgCciUqmAH7LabDbMnTsXP/30E6xWK5555hn07dsXc+bMgSAI6NevHxYuXAhRFLFq1Srs27cPGo0Gc+fOxaBBgzpqDEQxx3ksvVCvRYLN9nOdx9JTCAUM+B07dkCv12P58uW4evUqHnnkEfziF7/AzJkzMWTIEJSWlmLPnj1IT0/HwYMHsWnTJpw7dw4lJSXYsmVLR42BOoFUjRl3ppoAABqtBfZ4R4ceEmlXmi4cDsB15mt7znh1HUsfr4Ns/jnSeSw9hVLAgM/Pz0deXp7ra0mSUF1djcGDmw5DGzZsGPbv34/MzEzk5uZCEASkp6fD4XDgypUrSE1NDW/31GlobQbX4ZA6nQZWq71DD4m0OhScvH62q/PM12BnvBJFWsCAT0xsWivDYDBgxowZmDlzJpYtWwZBEFy3NzQ0wGAwQK/Xu92voaEhaMBLkgC9PqG9YwgbSRKjur/WaOlYzNcakRCvc6tpJNGr5q/uWdNoLdDpNBBFATqdplld8nn/NK0Zd3W3eD1XgiS57i8IgusxRcn9cT1rgij4rDtrOp3Gq67RSIiP10DQSIi/3p/VbHdt49zeuZ1T8+0D1ZrXRVF0u93f9uiihS4len8WO+PvSiwIeqLTuXPnMG3aNBQXF2P06NFYvny56zaj0YiUlBQkJSXBaDS61ZOTg1+D0+FQUFdnamPr4afXJ0R1f63R0rE0Wh0wmd1nge0O2avmr+5Zs8c7YLXaIcsKrFa7q66zOQAf9xetVtQd+djruYSh+a77O/fgZVmB4HB/XM+aIis+686a1WqH7FG32x0wm2V0sTtgud6f3f7zNs7tnds5Nd8+UK15PT5eB3Oz2/1tL1psMEXxz2Jn/F2JpLS0ll3jOOBRNP/9738xadIkPPfccxg3bhwA4NZbb0VVVRUAoKKiAjk5OcjOzkZlZSVkWcbZs2chyzKnZyh2XV9l0n59rt1kVyAHvxdR1Am4B//WW2+hvr4eq1evxurVqwEA8+bNw6JFi7BixQpkZWUhLy8PkiQhJycHhYWFkGUZpaWlHdI8UTj4WmUy68aW7TERRZOAAT9//nzMnz/fq7527VqvWklJCUpKSkLXGRERtQtPdCIiUikGPBGRSjHgiYhUiuvBE0WIXQGume2wN7serKQ01Xn+FIUCA54oQqwOBaevNLgdg9+vpxWiQ+EZshQSnKIhIlIpBjwRkUox4ImIVIpz8BQR3RIldFEuuZYAdooXfSy0RURtwoCniNDKjbCdrIDhYoNbXRyaH6GOiNSHAU8UA5xXgPIka5PQiMQIdESxgAFPFANcV4DyIPYfBWgZ8OQbP2QlIlIpBjwRkUox4ImIVIoBT0SkUvyQtROzKIDJ5nCrORQ/GxNRzGHAd2ImmwN7T1x0q93VLy1C3RBRqDHgidrq+sW5naTrF+nWSQKX+6WowIAnaiPnxbmdnBfpvqVHCpf7pajAD1mJiFSKAU9EpFIMeCIileIcPFEM87UIGRcgIycGPFEM87UIGRcgIydO0RARqRQDnohIpThFQxRqzU6A4slPFEncgycKMZtDwcnz9Th5vh7XTFacPF8PKxf5oQjgHnwnYFGAc9ca0WjlwmJEnQkDvhMw2Ryoqr0Mk9nqVufCYh3Ix7o1cgTboc6hRVM0R48excSJEwEAtbW1mDBhAoqLi7Fw4ULIctOP6apVqzBu3DgUFRXh2LFj4euYotKtNyi4M9WE7sol3Jlqcv3XLYF/JgDu0zbOqRuHzO8NhVfQPfi3334bO3bsQHx8PABgyZIlmDlzJoYMGYLS0lLs2bMH6enpOHjwIDZt2oRz586hpKQEW7ZsCXvzFD0SFTMM1Z/AdmMyDM0W4LopZ2QEu+qcfJ38BPAEqM4oaMBnZGRg5cqVeP755wEA1dXVGDx4MABg2LBh2L9/PzIzM5GbmwtBEJCeng6Hw4ErV64gNTU14GNLkgC9PiEEwwgPSRKjuj9f6httMFrc59oFjQJRFJAQr3OraySxzTXPukZrgU6ngSgK0Ol+/rESJcFvXfGoAYAgeNc8685/i6LgenzX43rUBFHwWXfWdDpNi+rN+w80Js96sPF7jleUBGg0EuLjPb4vGgnxPl4DX3VBsUD6ocprW/S9F3Ep4ZmWi8XfFX/UNJagAZ+Xl4cff/zR9bWiKBCEpl+axMRENDQ0wGAwQK/Xu7Zx1oMFvMOhoK7O1Nbew06vT4jq/ny5avV9EQ9ZVrzm4O0Ouc01z7o93gGr1Q5ZVmC12l3byA7Fb92zBjT9fHnWPOs6ncb1mILD43E9asr1aRB/21qtdlePgerNew00Js96sPF7jld2KLDbHTCb3Wfou9gdsPh4DXzV/W0rWmwwhennORZ/V/yJhbGkpSW3aLtWHyYpij/fxWg0IiUlBUlJSTAajW715OSWNUBEROHR6oC/9dZbUVXV9OdfRUUFcnJykJ2djcrKSsiyjLNnz0KW5aB770REFF6tPkxy9uzZWLBgAVasWIGsrCzk5eVBkiTk5OSgsLAQsiyjtLQ0HL0SEVErtCjge/bsiY0bNwIAMjMzsXbtWq9tSkpKUFJSEtruiChkQnV0TRyMEG0Gt5pQr0UcuvAonSjDE53Ir1tvUJComF3HtjsZhXj86zIXVok1vpYWBlq/vLBoM0Cu2eVejNdB7DWcyxRHGQY8+eXv2PakX+UDUMdhZERqxoCnVuuWKOFOmNz27LslKDAEuR8RdSwGPLWaVm6EoXqP2549z1iNXTzzVb0Y8ESdXKjm5in6cD14IiKVYsATEakUA56ISKUY8EREKsUPWWOYRWm6WlNzvAwfETkx4GOYyeZ7aWAiIoABTxR9PK7fCgDaCLVCsY0BTxRlbA4F3zdbGgIAcvp7hz7A4KfAGPBEMUBWmi7a7SmnfwSaoZjBgI8Bvj5MBfiBKoWXvyUMNLDB+4KAFI0Y8DHA14epAD9QpSbhmq/3t4SB0OfuED0DhRsDPoalasxu67QDQHflEkZmKKg3mN1qv06xICFN51W/M9XE9d1jmK+pG07bkBMDPoZpbQYYqj9xq9luTEZixmCcrd7jVjPVmZE2aLhX3XCxATcPfgCJStMUEJcAJlIPBjy5lv8FwCWAiVSESxUQEakU9+CJVMjzg1edJKBLhHqhyGHARxmuL0Pt5euD11t6pPDaTJ0QA74D+AptnUaC1e772PaKk1xfhmKPr+PmJa0WDpvNa1teDrBjMOA7gL9FwQ58d8lrW4Y5xSpfx81r+9wNWysuBxgHI0Sb97FbfENoGwY8UWcgAHYl+ufmRZsBcs0u7zqvD9smDHiiTsDmUHDNZMV3nJvvVBjwRJ2Zjz17nSRAwxObVYEBT9SJ+dqzv6VHCjQ+Et6uAFaHAqnZGwLfDKIbAz7EeJgjqZXV0XT4Zb+eP78h+HszCDV/K1vyw9fAGPAhFq7L6N16g4JExexWixd17X5cIi/XrygleUzdyCF8itYuRexvZUt++BpYSANelmWUlZXh5MmT0Ol0WLRoEXr37h3Kp4gIX3vlCVoJXcKw4+JvhUh9vAln/u9Tt7o4ND/0DVCn57yiVPM9dQDIujE5ZM8RqqWI/b1R+Dr+vjPu7Yc04Hfv3g2r1YoNGzbgyJEjWLp0Kd58881QPkXYNQ9z87VGNFodPk8++s0vu8OkeM+9tHc6xt8KkdqMwe17YKIIc87hA3D9ddDeOXx/bxS+jr/vjHv7IQ34r776CkOHDgUA3HbbbTh+/HgoH95Ne/eqA10lyRnmCfE6mMxWtykW51RJiu0Cfrry8552SlI86g1m/L/UBLc98JSkeCQ3W4LXqZfuGuweNYDTLhRjrk/nWM12xF8PbUkS4PAIcqBpisf5F4Hzr4NbbkqBtdlOkaQ0/W46PPaUJKXpDcLzzUCG/2vVetZThJbt7Qv1WiTYbCHb4/d18lZH/TUhKIqP3dA2mjdvHu677z4MHz4cADBixAjs3r0bGg2n+omIOlpIlwtOSkqC0Wh0fS3LMsOdiChCQhrw2dnZqKioAAAcOXIE/fvz2mFERJES0ika51E0NTU1UBQFixcvRp8+fUL18ERE1AohDXgiIooevGQfEZFKMeCJiFSKAU9EpFIMeDR9OFxaWorCwkJMnDgRtbW1brd//vnnGD9+PMaPH4+ysjJE88cWwcbyzjvvYOzYsXj00Uexa5f3hRWizdGjRzFx4kSv+meffYZHH30UhYWF2LhxYwQ6az1/Y9m5cycKCgpQVFSE0tJSyHIoV30JPX/jcFqwYAFee+21Duyo7fyN5dixYyguLsaECRMwY8YMWCyWCHQXAgopn376qTJ79mxFURTl8OHDyh/+8AfXbQ0NDcoDDzygXL58WVEURVmzZo3r39Eo0FiuXbumDB8+XLFYLEpdXZ0yYsSISLXZImvWrFEefPBBpaCgwK1utVqVe++9V6mrq1MsFosyduxY5eLFixHqsmX8jcVsNisjR45UTCaToiiK8uc//1nZvXt3JFpsEX/jcFq/fr0yfvx4Zfny5R3cWev5G4ssy8pDDz2knDlzRlEURdm4caNy6tSpSLTYbtyDR+AlFg4fPoz+/ftj2bJlKC4uRrdu3ZCamhqpVoMKNJb4+Hikp6fDbDbDbDZDEKJ7Ie+MjAysXLnSq37q1ClT9HnyAAAC7ElEQVRkZGSga9eu0Ol0uOOOO3Do0KEIdNhy/sai0+nwwQcfID4+HgBgt9vRpUs0XUTPnb9xAE2/K0ePHkVhYWEHd9U2/sZy+vRp6PV6vPfee/jd736Huro6ZGVlRaDD9mPAAzAYDEhKSnJ9LUkS7HY7AODq1auoqqrCrFmz8Pbbb+O9997D6dOnI9VqUIHGAgA33XQTHnjgATzyyCN4/PHHI9Fii+Xl5fk8E9pgMCA5+eeVDRMTE2EweF+oOZr4G4soiujWrRsAoLy8HCaTCffcc09Ht9di/sZx8eJFrFq1CqWlpRHoqm38jeXq1as4fPgwiouL8e677+LLL7/EgQMHItBh+3EdAQReYkGv12PgwIFIS2tacCwnJwfffvstMjMzI9JrMIHGUlFRgYsXL2LPnj0AgMmTJyM7OxuDBg2KSK9t5TlGo9HoFvixRpZlLF++HKdPn8bKlSuj/i8rXz755BNcvXoVU6dOxaVLl9DY2IisrCyMHTs20q21ml6vR+/evdG3b18AwNChQ3H8+HHcddddEe6s9bgHj8BLLAwYMAA1NTW4cuUK7HY7jh496nrho1GgsXTt2hVxcXHQ6XTo0qULkpOTUV9f7++holafPn1QW1uLuro6WK1WHDp0CLfffnuk22qz0tJSWCwWrF692jVVE2sef/xxbN26FeXl5Zg6dSoefPDBmAx3AOjVqxeMRqPrAIVDhw6hX79+Ee6qbbgHD2DUqFHYv38/ioqKXEssvPvuu8jIyMDIkSPx7LPP4qmnngIA5OfnR/UaO8HG8sUXX2D8+PEQRRHZ2dlRPR3g6aOPPoLJZEJhYSHmzJmDyZMnQ1EUPProo+jevXuk22sV51gGDBiAzZs3IycnB0888QSAprAcNWpUhDtsmeavSaxrPpZXXnkFzz77LBRFwe23344RI0ZEur024VIFREQqxSkaIiKVYsATEakUA56ISKUY8EREKsWAJyJSKQY8EZFKMeCJiFTq/wO+jl8GFpYqpwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N=50\n",
    "sample_means_1=[poisson(1).rvs(N).mean() for i in range(5000)]\n",
    "sample_means_2=[poisson(1.1).rvs(N).mean() for i in range(5000)]\n",
    "ax=sns.distplot(sample_means_1,kde=False)\n",
    "sns.distplot(sample_means_2,kde=False,ax=ax)\n",
    "s=.1/np.log(1.1)\n",
    "j=ax.axvline(s)\n",
    "t=ax.set_title(\"Critical value is {0:.2f}\".format(s))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fraction of draws from poisson(1.1) that lie to the left of the critical line: 0.3672\n",
      "Fraction of draws from poisson(1) that lie to the right of the critical line: 0.3476\n"
     ]
    }
   ],
   "source": [
    "print('Fraction of draws from poisson(1.1) that lie to the left of the critical line:',len([x for x in sample_means_2 if x<s])/len(sample_means_2))\n",
    "print('Fraction of draws from poisson(1) that lie to the right of the critical line:',len([x for x in sample_means_1 if x>s])/len(sample_means_1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the significance level is $36\\%$.  And if we draw the line so that the chance of falsely rejecting the null hypothesis is only $5\\%$, then the chance of correctly choosing the alternative will be low.  This is an example of an experiment with **low power** -- there's not enough information to distinguish the two cases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import uniform\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also look at the uniform distribution.  Here we look at the maximum of 100 draws from a unifrom distribution from 0 to 1, and from 0 to 1.01.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD7CAYAAABgzo9kAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAHh5JREFUeJzt3Xt0FOXBBvBnZrIbdpONayAgfDEUMFUISgsBFSNWLg2fFW9cQgKhx1BaLIZivUQCBjSAWBXPIQIWtJ6WWClKPaXKkaOo5QApCCqUVbwVUSCKSC7sLbuZeb8/OOxnSHY3O9lJNuPz+8vdd2f2SYLPznl35h1JCCFARESmI3d1ACIiMgYLnojIpFjwREQmxYInIjIpFjwRkUmx4ImITIoFT0RkUix4IiKTYsETEZlUUle+uaZpUFV9F9IqiqR7WyMxV2yYKzbMFZtEzQV0LJvForTrdVELPhgMory8HCdOnEAgEMBdd92Fyy67DA8++CAkSUJ2djaWLFkCWZbx9NNP45133kFSUhLKy8tx1VVXRdy3qgrU13vb9xNdwOm0697WSMwVG+aKDXPFJlFzAR3LlpHhaNfrohb81q1b4XQ68fjjj6Ourg633347rrjiCixYsABXX301KioqsGPHDvTr1w/79u3DSy+9hNraWpSWlmLLli26whMRUcdFLfiJEyciPz8/9FhRFLhcLowaNQoAMGbMGOzevRsDBgxAXl4eJElCv379oKoqzpw5g/T0dOPSExFRWFELPiUlBQDgdrsxf/58LFiwAI899hgkSQqNnz17Fm63G06ns8V2Z8+ejVjwiiLB6bTrCq4osu5tjcRcsWGu2DBXbBI1F9A52dr1JWttbS3mzZuHoqIiTJo0CY8//nhozOPxIC0tDampqfB4PC2edzgizxNxDr7zMFdsmCs2zBW7zpiDj3qa5OnTp1FSUoL7778fU6ZMAQAMGTIEe/fuBQDs3LkTubm5GD58OHbt2gVN03Dy5ElomsbpGSKiLhT1CP6ZZ55BY2Mj1q5di7Vr1wIAFi1ahGXLlmHVqlUYOHAg8vPzoSgKcnNzUVBQAE3TUFFRYXh4IiIKT+rKOzoFgyqnaDoJc8WGuWLDXLFLiCkaIiLqnljwREQm1aVLFRARdVSTALxBtc0x2R/s5DSJhQVPRN2aN6ji7SOn2hz732H/A1sn50kknKIhIjIpFjwRkUmx4ImITIoFT0RkUix4IiKTYsETEZkUC56IyKRY8EREJsWCJyIyKRY8EZFJseCJiEyKBU9EZFIseCIik2LBExGZFJcLJiIyQA94IAfdYcelposBWA3N0K6CP3jwIJ544gls3LgR99xzD06fPg0AOHHiBIYNG4annnoKc+fORX19PSwWC5KTk/Hss88aGpyIKJHJQTe0T94I/4IrbwKQbmiGqAW/YcMGbN26FTbbuWXzn3rqKQBAQ0MDZs2ahYULFwIAvvzyS7z22muQJMnAuERE1F5R5+CzsrJQVVXV6vmqqirMnDkTvXv3xunTp9HY2Ii5c+eisLAQb7/9tiFhiYio/aIewefn5+P48eMtnvvuu+9QU1MTOnoPBoMoKSnBrFmz0NDQgMLCQlx11VXo2bNnxH0rigSn064ruKLIurc1EnPFhrliw1yt+Rr8sNvansuWZAlOR9fkkhotQJhcACDLEpxpxmbT9SXr66+/jptvvhmKogAAevXqhenTpyMpKQk9e/bE4MGDcfTo0agFr6oC9fVePRHgdNp1b2sk5ooNc8WGuVrzB1R4fYE2x4Smv2M6yh4MQguTCwBsHciWkeFo1+t0nSZZU1ODMWPGhB7v2bMHCxYsAAB4PB58+umnGDhwoJ5dExFRnOg6gj969CguvfTS0OMbbrgBu3btwrRp0yDLMn7/+98jPd3Yb4eJiCiydhV8ZmYmNm/eHHr82muvtXrNokWL4peKiIg6jFeyEhGZFAueiMikWPBERCbFgiciMikWPBGRSbHgiYhMigVPRGRSLHgiIpNiwRMRmRQLnojIpFjwREQmxYInIjIpFjwRkUmx4ImITIoFT0RkUix4IiKTYsETEZkUC56IyKRY8EREJtWugj948CCKi4sBAC6XC9dffz2Ki4tRXFyMbdu2AQCefvppTJkyBdOnT8ehQ4eMS0xERO0S9abbGzZswNatW2Gz2QAAH374Ie68806UlJSEXuNyubBv3z689NJLqK2tRWlpKbZs2WJcaiIiiirqEXxWVhaqqqpCjw8fPox33nkHM2bMQHl5OdxuNw4cOIC8vDxIkoR+/fpBVVWcOXPG0OBERBRZ1CP4/Px8HD9+PPT4qquuwtSpUzF06FCsW7cOa9asgcPhgNPpDL0mJSUFZ8+eRXp6esR9K4oEp9OuK7iiyLq3NRJzxYa5YsNcrfka/LDbrG2OSbIEp6NrckmNFiBMLgCQZQnONGOzRS34C02YMAFpaWmh/66srMS4cePg8XhCr/F4PHA4HFH3paoC9fXeWCMAAJxOu+5tjcRcsWGu2DBXa/6ACq8v0OaY0PR3TEfZg0FoYXIBgK0D2TIyovcroOMsmtmzZ4e+RK2pqUFOTg6GDx+OXbt2QdM0nDx5EpqmRT16JyIiY8V8BL906VJUVlbCYrGgV69eqKysRGpqKnJzc1FQUABN01BRUWFEViIiioEkhBBd9ebBoMopmk7CXLFhrth0Za66gIq3j5xqc+x/h/0PbELr5ETn2IPfQPvkjbDjtitvQp2mb6bDsCkaIiLqHljwREQmxYInIjKpmL9kJSLqTE0C8AbVsONql32LmPhY8ESU0LzB8F+iAsC12RmdmKZ74RQNEZFJ8QieiEiHHvBADrrDjichiPDXsXYOFjwRkQ5y0B3xPHdp0OhOTNM2TtEQEZkUC56IyKRY8EREJsWCJyIyKRY8EZFJseCJiEyKBU9EZFI8D56IqA3d4UKmaFjwRGSoaIuF2S0KkqVODNRO3eFCpmhY8ERkqGiLhd14RW8kW5VOTPTD0a6CP3jwIJ544gls3LgRH330ESorK6EoCqxWKx577DH06tULy5Ytw3vvvYeUlBQAwNq1a+FwtO+2UkREFH9RC37Dhg3YunUrbDYbAGD58uV46KGHMHjwYGzatAkbNmzAwoUL4XK58OyzzyI9Xd89BomIKL6inkWTlZWFqqqq0ONVq1Zh8ODBAABVVZGcnAxN03Ds2DFUVFRg+vTpePnll41LTERE7RL1CD4/Px/Hjx8PPe7duzcA4L333kN1dTVeeOEFeL1ezJw5E3feeSdUVcWsWbMwdOhQXHHFFRH3rSgSnE67ruCKIuve1kjMFRvmik13zOVr8MNus4bdtkeyBc6LeoQdj7Z9kiKHHZdkCU6Hvt+X1GgBIryvlKTA1oFxWZbgTDP2b6nrS9Zt27Zh3bp1WL9+PdLT00Olfn4a55prrsGRI0eiFryqCtTXe/VEgNNp172tkZgrNswVm+6Yyx9Q4fWFP6HQ3xREfb0WfjzK9s2qFnZcaPo7xh4MQovwvsnNKpo6MG7rQLaMjPZ9vxnzhU7/+Mc/UF1djY0bN+LSSy8FAHzxxRcoKiqCqqoIBoN47733kJOTE+uuiYgojmI6gldVFcuXL0ffvn1RWloKABg5ciTmz5+PSZMmYdq0abBYLLj11luRnZ1tSGAiImqfdhV8ZmYmNm/eDADYt29fm6+ZM2cO5syZE79kRETUIVyLhojIpFjwREQmxYInIjIprkVDRD9IZlgtMhoWPBH9IJlhtchoOEVDRGRSPIInoi4lSRLqAuHXi1dFJ4YxGRY8EXUpX7OGmk+/DTt+bXZGJ6YxF07REBGZFAueiMikWPBERCbFgiciMikWPBGRSbHgiYhMigVPRGRSLHgiIpNiwRMRmRQLnojIpFjwREQm1a6CP3jwIIqLiwEAx44dQ2FhIYqKirBkyRJomgYAePrppzFlyhRMnz4dhw4dMi4xERG1S9SC37BhAxYvXoympiYAwKOPPooFCxbgr3/9K4QQ2LFjB1wuF/bt24eXXnoJq1atwsMPP2x4cCIiiixqwWdlZaGqqir02OVyYdSoUQCAMWPGYM+ePThw4ADy8vIgSRL69esHVVVx5swZ41ITEVFUUZcLzs/Px/Hjx0OPhRCQJAkAkJKSgrNnz8LtdsPpdIZec/759PT0iPtWFAlOp11XcEWRdW9rJOaKDXPFpjvm8jX4YbdZw26bpMiGjUuyBKej7VxSowWIsF8pSYHNwHFZluBMM/ZvGfN68LL8/wf9Ho8HaWlpSE1NhcfjafG8w+GIui9VFaiv98YaAQDgdNp1b2sk5ooNc8WmO+byB1R4feHvbtqsaoaNCy18x9iDQWgR9pvcrKLJwHFbhGzRZGRE71dAx1k0Q4YMwd69ewEAO3fuRG5uLoYPH45du3ZB0zScPHkSmqZFPXonIiJjxXwEX1ZWhoceegirVq3CwIEDkZ+fD0VRkJubi4KCAmiahoqKCiOyEhFRDNpV8JmZmdi8eTMAYMCAAaiurm71mtLSUpSWlsY3HRER6cYLnYiITIoFT0RkUix4IiKTYsETEZkUC56IyKRY8EREJsWCJyIyKRY8EZFJseCJiEyKBU9EZFIseCIik4p5sTEiokQwpKdAivChj/gWI9NbL7vrkWxdkCqxsOCJqFtKET64Xa8j2NsB96mzrcZTcyZ2QarEwikaIiKTYsETEZkUC56IyKRY8EREJsWCJyIyKRY8EZFJ6TpN8u9//zteeeUVAEBTUxM++ugjPPnkk/jDH/6Avn37Ajh3j9ZRo0bFLykRJaQmAdQ2+OEPqG2Oq6KTA1GIroK/4447cMcddwAAHn74YUyePBkulwv3338/8vPz4xqQiBKbN6hi77Hv4PUF2hy/NjujkxPReR2aovnPf/6Dzz77DAUFBXC5XNiyZQuKioqwcuVKNDc3xysjERHp0KErWf/4xz9i3rx5AIDrrrsO48ePR2ZmJpYsWYJNmzZh5syZEbdXFAlOp13XeyuKrHtbIzFXbJgrNomYy9fghyxLsNusbY4nKXLYsY6MJ1maYLUmQZYlWK2tqyzJokCSJTgdbf++pEYLEOF9pSQFNgPHZVmCM83Yv6Xugm9sbMR///tfXHPNNQCAyZMnIy0tDQAwbtw4bN++Peo+VFWgvr71GhLt4XTadW9rJOaKDXPFJhFz+QMqNE2EnaJpVrWwYx0Zb7apCASaoWkCgUDrGQNrUIXQwneMPRiEFuF9k5tVNBk4bouQLZqMDEe7Xqd7iubdd9/F6NGjAQBCCNxyyy34+uuvAQA1NTXIycnRu2siIooD3UfwR48eRWZmJgBAkiQsW7YMd999N3r06IFBgwZh2rRpcQtJRESx013wv/rVr1o8zsvLQ15eXocDERFRfPBCJyIik2LBExGZFAueiMikWPBERCbFgiciMikWPBGRSfGm20QUUZM4t6BYOFwtMnGx4IkoIm9QxdtHToUd52qRiYtTNEREJsWCJyIyKU7REJEp9UpRkBb4Gha0/SVBEoIIv9ajObDgiciULJof6id7oEBrc1waNLqTE3U+TtEQEZkUC56IyKQ4RUNEEc9153nu3RcLnoginuvO89y7LxY8ESWkIT0FUoQPfcS3GJne+t6lvewC7i7I1Z2w4IkoIaUIH9yu1xHs7YD71NlW431zx3VBqu5Fd8HfdtttcDjO3dk7MzMTBQUFWL58ORRFQV5eHu6+++64hSQiotjpKvimpiYAwMaNG0PP3XrrraiqqsKll16KX//613C5XMjJyYlPSiIiipmu0ySPHDkCn8+HkpISzJo1C++++y4CgQCysrIgSRLy8vJQU1MT76xERBQDXUfwPXr0wOzZszF16lR88cUXmDNnDtLS0kLjKSkp+Oqrr6LuR1EkOJ12PRGgKLLubY3EXLFhrtgYlcvX4IfdZm1zLEmRw46dH5dlqdVrfnyRCpvwoq90Gtf2aWr9npIdnzQoYfefZGmC1ZoEWZZgtbauKlmRoo5LMmBLbju7lKTAFuHnMnpcliU404z9N6ar4AcMGID+/ftDkiQMGDAADocD9fX1oXGPx9Oi8MNRVYH6+tbfjreH02nXva2RmCs2zBUbo3L5Ayq8vrZXZmlWtbBj58c1TbR6jcXmRb3rdaT3dqC+jS9JU3Mmwuuzh91/s01FINAMTRMIBJpbjWuqiDouNMAXJntys4qmCD+X0eM2TX//ZWQ42vU6XVM0L7/8MlauXAkA+Oabb+Dz+WC32/Hll19CCIFdu3YhNzdXz66JiChOdB3BT5kyBQsXLkRhYSEkScKKFSsgyzLuu+8+qKqKvLw8DBs2LN5ZiUgn3pXph0lXwVutVjz55JOtnt+8eXOHAxFR/PGuTD9MXGyMiMikWPBERCbFgiciMikWPBGRSbHgiYhMigVPRGRSLHgiIpPievBEJtAdL2TqlaJgJLy8oYeBWPBEJtAdL2SyaH64XTt4Qw8DcYqGiMikWPBERCbFgiciMikWPBGRSbHgiYhMigVPRGRSLHgiIpNiwRMRmRQLnojIpFjwREQmpWupgmAwiPLycpw4cQKBQAB33XUXLrnkEsydOxc/+tGPAACFhYW46aab4pmViIhioKvgt27dCqfTiccffxx1dXW4/fbbMW/ePNx5550oKSmJd0YiItJBV8FPnDgR+fn5oceKouDw4cM4evQoduzYgf79+6O8vBypqalxC0r0Q3Z+tUhfgx/+QOtVI7titcghPQVShA99xLcYltaEZlvLXFwNsuvpKviUlBQAgNvtxvz587FgwQIEAgFMnToVQ4cOxbp167BmzRqUlZVF3I+iSHA67XoiQFFk3dsaibliw1ztU9vgx95j30GWJWha6zYf0f9i2G3WsNsnKbLu8XBjFyWdhffDN6Gm2+Gr80GIlrmSf/ozWK1JkGUJVmvrqpEVyfBxSQZsyW3/XFKSAluE34nR47IswZlm7L8x3csF19bWYt68eSgqKsKkSZPQ2NiItLQ0AMCECRNQWVkZdR+qKlBf33od6PZwOu26tzUSc8WGudrHH1Dh9QVgt1nh9QVajTerWpvPx2M83FizTUUg0AxNExBCIBBobjGuqSI0fuFYZ40LDfCF+bmSm1U0RfidGD1u0/T3X0aGo12v03UWzenTp1FSUoL7778fU6ZMAQDMnj0bhw4dAgDU1NQgJydHz66JiChOdB3BP/PMM2hsbMTatWuxdu1aAMCDDz6IFStWwGKxoFevXu06gicyk0h3VbJbFCRLnRwoDob0FLzjUjemq+AXL16MxYsXt3p+06ZNHQ5E1F1FuqvSjVf0RrJV6eREHZcifAh+vJt3XOqmeKETEZFJseCJiEyKN90mSgCR5u+BrjnPnbo/FjxRAog0fw8A12ZndGIaMgsWPFEnkCQJdW1cgXqeUUfo37/atK0zYTySDR9+1w1P76F2YcETdQJfs4aaT78NO27UEXqK8MHteh3B3o42z4RJzZkIIHGu2KX4YsETtRPnyam7YcETtRPnyam7YcETJbAL59CTLC1XbeQcOkXCgidKYBfOoVutSS0W1uroHHqvFAUj4Y24HAF1Xyx4om6sPQUdab0Yi+aH27Uj7JewfXPHIfy3DpToWPBEXSjaaYzxKGj64WLBE3XA+YIG0GZJR5sjj3YaIwuaOoIFT9QB5wsaQJslzfPMqStxsTEiIpPiETx1G9EuNOquN9UgMgoLnrqNaBcajR3cB14R/rS+aB8A3/8A8TX44b9g7Rg9V6p29CwXoo5gwZNpRFvvJdpdlb7/AdLWza31XKnKs1yoK7HgKeH1gAdy0A1ZEzGviPj9s1x6at+iR7Dl6zRLKvxIiX9oogQQ14LXNA1Lly7Fxx9/DKvVimXLlqF///7xfAsyWEfnuY2YJ5eDbmifvIHmZgH3142txiOdqfL9s1ya+6bBfcE0S8qQ8ZDVc0fW6UDoA+T8kgBpqTY0usOfBskpFkpkcS34N998E4FAAH/729/wwQcfYOXKlVi3bl0834IMFm2eO5ZpjrZcOE9+4Vy3ng+ASPPc3y/goCrw3wumSbIzG/Dp/h0AgIHfm0Y5vyRA39xxOOk6N97WNAunWCiRxbXgDxw4gOuvvx4A8JOf/ASHDx+O5+5baPQHI95AoSvPqDD6KNiapCDQrG//HV3yVs+NK74/TZIW/AYnzvx/CfstCuzJ1tBRcs9eradLhBSEt1lAC/Oekea5WcD0QyYJEeG0gxgtWrQIP//5z3HDDTcAAH72s5/hzTffRFISp/qJiDpbXC90Sk1NhcfjCT3WNI3lTkTUReJa8MOHD8fOnTsBAB988AF+/OMfx3P3REQUg7hO0Zw/i+aTTz6BEAIrVqzAoEGD4rV7IiKKQVwLnoiIEgcXGyMiMikWPBGRSSVkwWuahoqKChQUFKC4uBjHjh1rMb5+/XrceuutmDFjBt5++20AwLfffotf/vKXKCoqwu9+9zv4fL6EyHXeu+++Gzp9NBFy1dfX4+qrr0ZxcTGKi4vx5z//OSFyeb1ePPDAAygqKsLUqVNx6NChhMi1fPny0O9q4sSJmDZtWkLkOnnyJGbOnIkZM2bgt7/9bcL8u//qq68wY8YMFBUV4b777jMk13kHDx5EcXFxq+ffeustTJ48GQUFBdi8eTMAwO/3o7S0FEVFRZgzZw7OnDmTELnOe+ONN3DvvffGL4RIQNu3bxdlZWVCCCHef/99MXfu3NDYkSNHxKRJk4Tf7xd+v1/cdtttwuv1imXLlolXXnlFCCHE6tWrxfPPP58QuYQQ4uTJk2Lu3Lli9OjRcc+kN9fu3bvFI488YkiejuRavXq1WL9+vRBCiI8++ij0N+3qXOcFAgExZcoUceTIkYTItXz5clFdXS2EEGLVqlXiL3/5S0LkKi0tFVu3bhVCCLF582axZs2auOcSQoj169eLm2++WUydOrXF84FAQIwfP17U19eLpqYmcccdd4hTp06JP/3pT2L16tVCCCFeffVVUVlZmRC5hBCisrJS5OfniwULFsQtR0IewUe6Ivbzzz/HqFGjkJycjOTkZPTv3x8ff/wxysvLccstt0DTNNTW1qJnz54JkaupqQlLlizB0qVL456nI7kOHz4Ml8uFmTNnYv78+Th1KvzyAp2Za9euXbBYLJg9ezbWrl0b2r6rc51XXV2N6667DpdffnlC5Bo8eDAaG8+tz+N2uw257kRPrs8++wxjxowBcO706QMHDsQ9FwBkZWWhqqqq1fOff/45srKycNFFF8FqtWLEiBHYv39/i59lzJgxqKmpSYhcwLnfU7x7IiEL3u12IzU1NfRYURQ0NzcDAC6//HLs378fbrcbdXV1eP/99+Hz+SBJElRVxc0334y9e/di+PDhCZHrkUceQUlJCfr06RP3PB3JNXDgQMyfPx/V1dUYP348li1blhC56urq0NjYiOeeew5jx47FY489lhC5ACAQCGDTpk2YPXt23DPpzXXJJZfghRdewC9+8Qvs3LkTEydOTIhcgwcPxltvvQUA2LFjh2FTNPn5+W1+qLndbjgcjtDjlJQUuN3uFs+npKTg7NnWSzh3RS4AuOmmmyBJ8V1fJSEvM410ReygQYMwY8YMzJkzB/3798ewYcNw8cUXAwAsFgu2bduGPXv2oKysDNXV1V2aS1EU7N+/H19++SXWrFmDhoYG3HPPPXjqqae6NNfFF1+MK6+8EjabDQAwYcIErF69Oq6Z9OZyOp0YO3YsAODGG2/E+vXrEyIXANTU1GDkyJEt/gft6lwLFy7Eo48+iuuvvx7vvPMOysrK4v4705OrrKwMlZWVePXVV3HttdeGfoed5cLMHo8HDoejxfMejwdpaWkJkcsoCXkEH+mK2DNnzqCurg4vvvgiFi1ahNraWmRnZ2Pp0qX497//DeDcp2K8Pwn15BoxYgS2b9+OjRs3YuPGjbjoooviXu56cmVnZ2Px4sXYvn07gHPFlZOTkxC5RowYgX/9618Azn0xfdlllyVELgDYs2dPaNrBCHpypaWlhQqid+/eoemars61Z88ezJs3D8899xxkWcbo0aPjniuSQYMG4dixY6ivr0cgEMD+/fvx05/+FMOHDw/9+9q5cydGjBiRELmMkpBH8BMmTMDu3bsxffr00BWxzz//PLKysjB27FgcP34ckydPhsViwQMPPABFUVBcXIylS5dizZo1kGXZkDlvPbk6g55c9957L8rLy/Hiiy/CZrMZMkWjJ9dvfvMbLF68GAUFBUhKSjJkikbv3/Ho0aO47bbb4p6nI7keeughPPLII9A0DUIIVFRUJESuAQMGoLy8HFarFdnZ2Ybkass///lPeL1eFBQU4MEHH8Ts2bMhhMDkyZPRp08fFBYWoqysDIWFhbBYLHjyyScTIpdReCUrEZFJJeQUDRERdRwLnojIpFjwREQmxYInIjIpFjwRkUmx4ImITIoFT0RkUix4IiKT+j/RpBGMZl+SEwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "j=sns.distplot([uniform(0,1).rvs(100).max() for i in range(1000)],kde=False)\n",
    "j=sns.distplot([uniform(0,1.01).rvs(100).max() for i in range(1000)],kde=False,ax=j)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So suppose we have a sample $x$ of $100$ numbers and we to use it to test whether or not   it come from the uniform distribution on $[0,1]$ (null hypothesis) or from the uniform distribution on $[0,1.01]$ (alternative hypothesis).\n",
    "\n",
    "The probability density for $[0,1]$ can be expressed in terms of the maximum $T$ and the minimum $M$ of the sample $x$. We have $P(T<r)=r^{100}$ for $0\\le r\\le 1$.  In addition,  $P(T<0)=0$ and $P(T>1)=1$.  The associated density is the derivative of this $p(r)=100r^{99}$ but is zero outside the interval.\n",
    "\n",
    "Similarly we have $p_{1.01}(r)=100r^{99}/(1.01)^{100}$ and zero outside $[0,1.01]$.\n",
    "\n",
    "So the likelihood ratio is\n",
    "$$\\frac{p_{1.01}(r)}{p(r)}=(1.01)^{-100}$$\n",
    "for $0\\le r\\le 1$ and is $\\infty$ between $1$ and $1.01$.\n",
    "\n",
    "What this means is that if we get a sample with $r>1$ then we should obviously reject the null hypothesis.  If we use $r=1$ as our criterion, then we will never falsely reject the null hypothesis, but since $(1/(1.01))^{100}=.36$,\n",
    "we will falsely accept it $36\\%$ of the time. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.3697112123291189"
      ]
     },
     "execution_count": 132,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(1/1.01)**100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}