From 06b98473012640930d014924b627002c99b20914 Mon Sep 17 00:00:00 2001 From: jeremyteitelbaum Date: Wed, 2 May 2018 09:44:46 -0400 Subject: [PATCH 1/4] catching up --- BDA 5.9.3.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/BDA 5.9.3.ipynb b/BDA 5.9.3.ipynb index abfbe93..f77853a 100644 --- a/BDA 5.9.3.ipynb +++ b/BDA 5.9.3.ipynb @@ -492,7 +492,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.5" + "version": "3.6.4" } }, "nbformat": 4, From 722cff49a00a72fc7f8ecd3291d3eb6b77227937 Mon Sep 17 00:00:00 2001 From: jeremyteitelbaum Date: Wed, 2 May 2018 10:26:06 -0400 Subject: [PATCH 2/4] fixed 5.9.8 --- BDA 5.9.8.ipynb | 98 +++++++++++++++++++++++++++++++++---------- Useful Formulae.ipynb | 62 +-------------------------- 2 files changed, 78 insertions(+), 82 deletions(-) diff --git a/BDA 5.9.8.ipynb b/BDA 5.9.8.ipynb index 8ad5732..9f6a2e1 100644 --- a/BDA 5.9.8.ipynb +++ b/BDA 5.9.8.ipynb @@ -7,8 +7,8 @@ "# Discrete Mixture Models\n", "\n", " This solution differs from the one published here:\n", - "http://www.stat.columbia.edu/~gelman/book/solutions3.pdf\n", - "\n", + "http://www.stat.columbia.edu/~gelman/book/solutions3.pdf and is probably wrong\n", + "\n", "\n", "Discrete mixture models: if $p_m(\\theta)$, for $m=1,\\ldots,M$ are conjugate prior densities for the sampling model $y|\\theta$, show that the class of finite mixture prior densities given by \n", "$$\n", @@ -50,29 +50,76 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "\n", + " This is the part that was wrong -- I leave it here for historical purposes, but the correct part follows. \n", + "

\n", "In the special case under consideration, $p_1$ is normal with mean $-1$ and $\\sigma=.5$, $p_2$ is normal with mean $1$ and $\\sigma=.5$ and we can set $\\lambda_1=.1$ and $\\lambda_2=.9$. The $p_m(\\{y_{i}\\})$ can be calculated from the $t$ distribution. Drawing a sample of size $10$ from $p_1$ and getting a sample mean of $-.25$ and a sample variance of $1$ gives a $t$-statistics of \n", - "$$\\frac{(\\overline{y}-\\mu)}{s/\\sqrt{N}}=\\frac{(-.25+1)}{1/\\sqrt{10}}=\\sqrt{10}(-.25+1)$$ in the first case and $\\sqrt{10}(-.25-1)$ in the second. " + "$$\\frac{(\\overline{y}-\\mu)}{s/\\sqrt{N}}=\\frac{(-.25+1)}{1/\\sqrt{10}}=\\sqrt{10}(-.25+1)$$ in the first case and $\\sqrt{10}(-.25-1)$ in the second. \n", + "

\n", + " Now we continue with what is correct \n", + "\n", + "We need to properly interpret $p_m(\\{y_{i}\\})$ and for that we should remember where it comes from. We\n", + "rewrote\n", + "$$\n", + "p(\\{y_{i}\\}|\\theta)p_{m}(\\theta)=p_m(\\{y_{i}\\},\\theta)=p_{m}(\\{y_{i}\\})p_{m}(\\theta|\\{y_{i}\\})\n", + "$$\n", + "We're dealing here with normal distributions. On the left, the quadratic form that is the log-likelihood of the\n", + "relevant bivariate normal is\n", + "$$Q=\\frac{(\\overline{y}-\\theta)^2}{\\sigma^2}+\\frac{(\\theta-\\mu)^2}{\\tau^2}\n", + "$$\n", + "where, more specifically, $\\overline{y}=-.25$, $\\mu=\\pm 1$, $\\sigma^2=1/10$, and $\\tau^2=0.25=0.5^2$. \n", + "The first term comes from $p_{m}(\\overline{y}|\\theta)$ and the second from $p_{m}(\\theta)$.\n", + "\n", + "Pure algebra (by expanding, writing $Q$ as a quadratic in $\\theta$, and completing the square) gives us\n", + "the expression\n", + "$$\n", + "Q=\\frac{(\\theta-\\mu_{1})^2}{\\tau_1^2}+\\frac{(\\overline{y}-\\mu)^2}{\\sigma^2+\\tau^2}\n", + "$$\n", + "where \n", + "$$\n", + "\\mu_{1}=\\frac{\\mu/\\tau^2+\\overline{y}/\\sigma^2}{\\tau_1^2}\n", + "$$\n", + "and \n", + "$$\n", + "\\frac{1}{\\tau_1^2}=\\frac{1}{\\sigma^2}+\\frac{1}{\\tau^2}\n", + "$$\n", + "\n", + "In the context of the problem under discussion, the two terms tell us that\n", + "$$\n", + "\\theta\\sim N(\\mu_1,\\tau_1^2)\n", + "$$\n", + "where $\\tau_1^2=1/(1/.1+1/.25)=1/14$ and $\\mu_1=(-.25/.1\\pm 1/.25)/14$ giving $6.5$ or $1.5$. \n", + "and $p_{m}(\\{y_{i}\\})=N(-.25,\\pm 1, .1+.25)$.\n", + "\n", + "\n" ] }, { "cell_type": "code", - "execution_count": 29, + "execution_count": 17, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.30191827840729224 0.07235502834102417\n" + ] + } + ], "source": [ "import numpy as np\n", "from scipy.stats import norm, t\n", "import matplotlib.pyplot as plt\n", - "t_1=np.sqrt(10)*.75\n", - "t_2=np.sqrt(10)*(-1.25)\n", - "p1y=t.pdf(t_1,df=9)\n", - "p2y=t.pdf(t_2,df=9)\n", - "lambda1,lambda2=(.1,.9)\n" + "p2y=norm.pdf(-.25,1,np.sqrt(.35))\n", + "p1y=norm.pdf(-.25,-1,np.sqrt(.35))\n", + "lambda1,lambda2=(.1,.9)\n", + "print(p1y,p2y)" ] }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 18, "metadata": {}, "outputs": [], "source": [ @@ -82,14 +129,14 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "0.600590082806086 0.3994099171939141\n" + "0.3167705292212528 0.6832294707787472\n" ] } ], @@ -99,7 +146,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 20, "metadata": {}, "outputs": [], "source": [ @@ -111,22 +158,31 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 21, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-0.4642857142857143 0.10714285714285714\n" + ] + } + ], "source": [ "post_mean1,post_var1=posterior(-1,.25,-.25,1,10)\n", - "post_mean2,post_var2=posterior(1,.25,-.25,1,10)" + "post_mean2,post_var2=posterior(1,.25,-.25,1,10)\n", + "print(post_mean1,post_mean2)" ] }, { "cell_type": "code", - "execution_count": 28, + "execution_count": 22, "metadata": {}, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VOX1+PHPSQDD6gYiAhpUUMBElriwWFFcICCIIoJVwd1atOpXf6Kt2tL67dd9qUu1VbSWgiCCKCAILkgBJW6ALMoSMIKK7IgsCc/vj5OBEBIyJPfOvXfmvF+veU1mcufeM4Gc3Dn3ec4jzjmMMcYkl7SgAzDGGOM9S+7GGJOELLkbY0wSsuRujDFJyJK7McYkIUvuxhiThCy5G2NMErLkbowxSciSuzHGJKFqFW0gIi8BPYEfnXMnlfF9AZ4EcoGtwCDn3GcV7bd+/fouMzPzgAM2xphU9umnn/7knGtQ0XYVJnfgZeBp4F/lfL870Lz4dhrwXPH9fmVmZpKXlxfH4Y0xxsSIyIp4tquwLOOcmw6s288mvYF/OTUbOEREGsUXpjHGGD94UXNvDHxb4nFB8XPGGGMC4kVylzKeK7PVpIhcLyJ5IpK3Zs0aDw5tjDGmLPHU3CtSADQt8bgJsKqsDZ1zLwAvAOTk5FivYWOS2M6dOykoKGDbtm1BhxJJGRkZNGnShOrVq1fq9V4k9/HAYBEZiV5I3eicW+3Bfo0xEVZQUEDdunXJzMxEB9WZeDnnWLt2LQUFBTRr1qxS+4hnKOQIoAtQX0QKgPuB6sUB/B2YiA6DXIIOhbyqUpEYY5LKtm3bLLFXkohw+OGHU5XydYXJ3Tk3oILvO+C3lY7AGJO0LLFXXlV/dl6UZYwJxNy5MGkStGgBvXtDms23NmY3+3UwkfTKK9CuHQwZAhddBBdcANu3Bx2Viar77ruPqVOnBh2Gpyy5m8hZtAiuuw7OOgu+/x6eegomToTbbgs6MhNFRUVFDB06lHPOOeeAXhN2ltxN5Pzud1C7NgwfDg0bws03w513wnPPwXvvBR2dCZP8/HxOPPFEBg4cSHZ2Nn379mXr1q1kZmYydOhQOnfuzOjRoxk0aBCvv/46ANOmTaNt27ZkZWVx9dVXs734I2Hp14Sd1dxNpHzxBUyZAg89BEccsef5oUNh9Gi49VbdxurvIRP7h/FSmzbwxBMVbrZ48WJefPFFOnXqxNVXX82zzz4L6DjyGTNmAPDOO+8AOsJn0KBBTJs2jRYtWnDllVfy3HPPceutt+7zmrCzXwETKX/7G9SqBddeu/fzGRnwwAMwbx6MHx9MbCacmjZtSqdOnQC4/PLLdyfnSy+9dJ9tFy9eTLNmzWjRogUAAwcOZPr06bu/X9ZrwsrO3E1kbNumZ+eXXgqHHrrv9/v1gz/8AR58UEfP2Ci8EInjDNsvpYcUxh7Xrl17n211ZHf5ynpNWNmZu4mMKVNg82ZN7mWpVg3+539g9my9GQOwcuVKZs2aBcCIESPo3LlzudueeOKJ5Ofns2TJEgBeffVVzjzzzITE6TVL7iYyxozRM/azzy5/m4EDoU4d+Oc/ExeXCbeWLVvyyiuvkJ2dzbp16/jNb35T7rYZGRkMGzaMSy65hKysLNLS0rjxxhsTGK13rCxjIsE5PXPv1g3210epTh09sx8xAh5/HOrVS1yMJpzS0tL4+9//vtdz+fn5ez1++eWXd3/dtWtXPv/88332U/o1YWdn7iYSFi7UMe1du1a87bXXwtat8Npr/sdlTFhZcjeRMG2a3seT3E87DVq3hmHD/I3JhF9mZibz588POoxAWHI3kfDee9CsGcSzproIXH45zJoFEfskbYxnLLmb0Nu1Cz78cP8XUkvr31/vR470JyZjws6Suwm9JUtg/Xro2DH+12RmQocOemHVmFRkyd2E3ief6P2ppx7Y6wYM0LbACxZ4H5MxYWfJ3YTeJ59oo7CWLQ/sdf36aY8ZO3s3VTFu3DgWVOIMYfz48fzf//2fDxHFx5K7Cb1PPoGcHEhPP7DXNWyodfoRI3ScvDGVUZnkXlhYSK9evRgyZMgBvcZLltxNqO3YAZ9/fuAlmZjLLoOlS2HOHG/jMtFQXsvf8tr6DhkyhFatWpGdnc0dd9zBzJkzGT9+PHfeeSdt2rRh6dKlLF26lG7dutG+fXvOOOMMFi1aBMCgQYO4/fbbOeuss7jrrrt4+eWXGTx4MAArVqyga9euZGdn07VrV1auXFnma7xkM1RNqM2frwk+J6dyr+/TB268Uc/eK/sHwlRdgB1/92n5+9hjj/H888/v09b3yiuvZOzYsSxatAgRYcOGDRxyyCH06tWLnj170rdvX0BnsP7973+nefPmfPzxx9x00028V7yQwNdff83UqVNJT0/fa9br4MGDufLKKxk4cCAvvfQSt9xyC+PGjdvnNV6yM3cTavPm6f3JJ1fu9YccAj166JDICCyeY3xQuuXvtGnTymzrW69ePTIyMrj22mt54403qFWr1j772rJlCzNnzuSSSy6hTZs23HDDDaxevXr39y+55JIyk/SsWbO47LLLALjiiiv26glf3muqys7cTajNm6e92o8/vvL7uOwyGDsWPvggvhmuxnsBdvzdp+VveapVq8Ynn3zCtGnTGDlyJE8//fTuM/KYXbt2ccghh/BFOR9D4m0JXDImv9oI25m7CbW5c6FVqwO/mFpSjx5Qty785z/exWWio3TL33POOafMtr5btmxh48aN5Obm8sQTT+xO4HXr1mXz5s0A1KtXj2bNmu1eZs85x5dffllhDB07dmRk8Yy64cOH77ftsFcsuZtQmzcPsrKqto+aNbX2PmYMFF83MymkdMvf2267rcy2vps3b6Znz55kZ2dz5pln8vjjjwPQv39/Hn74Ydq2bcvSpUsZPnw4L774IieffDKtW7fmzTffrDCGp556imHDhpGdnc2rr77Kk08+6ffbRipaecQvOTk5Li8vL5Bjm2j46Sdo0AAeeUQX4aiKyZO1XfDYsXDhhd7EZ/Zv4cKFtDzQyQkey8/Pp2fPnpFtHlbWz1BEPnXOVTjEwM7cTWjFLqZW9cwdtNbeoIGVZkzqsORuQiuW3LOzq76vatV0EY+33oJNm6q+PxMN1vLXmBCaOxfq19eZpl4YMEAX2R471pv9mYoFVfZNBlX92VlyN6G1aJGOlIlzJFuFOnSAY4+FEnNLjI8yMjJYu3atJfhKcM6xdu1aMjIyKr0PG+duQmvxYh3l4hURuO46uPtu+PprKJ7DYnzSpEkTCgoKWLNmTdChRFJGRgZNmjSp9OstuZtQWrdOR8t4nYAHDYJ774V//AMeftjbfZu9Va9enWbNmgUdRsqysowJpa+/1vsTTvB2v0ceCb17a2nGxrybZBZXcheRbiKyWESWiMg+PSxF5GgReV9EPheRuSKS632oJpXEkrsfpZPrr9dPBcV9m4xJShUmdxFJB54BugOtgAEi0qrUZn8ARjnn2gL9gWe9DtSklsWLteXAscd6v+9zztHFtp9+2vt9GxMW8Zy5nwoscc4tc87tAEYCvUtt44B6xV8fDKzyLkSTir7+WhN79ere7zstDW65BWbMgI8/9n7/xoRBPMm9MfBticcFxc+V9EfgchEpACYCN3sSnUlZixd7X28v6ZprtB3wI4/4dwxjghRPci9rlHHpgasDgJedc02AXOBVEdln3yJyvYjkiUieDY8y5dm1C775xt/kXreuLuLxxhu6UpMxySae5F4ANC3xuAn7ll2uAUYBOOdmARlA/dI7cs694JzLcc7lNGjQoHIRm6RXUKAzSf0eh37LLdqW4MEH/T2OMUGIJ7nPAZqLSDMRqYFeMB1fapuVQFcAEWmJJnc7NTeVsnix3vt55g7QqJGOnBk2zM7eTfKpMLk75wqBwcBkYCE6KuYrERkqIr2KN/sf4DoR+RIYAQxyNufYVJKfwyBLu+cePXsfOtT/YxmTSHHNUHXOTUQvlJZ87r4SXy8AOnkbmklVS5ZA7do64chvjRrBb38Ljz+ubQlOPNH/YxqTCDZD1YTOihWQmeldw7CK3HWXrtZ0//2JOZ4xiWDJ3YROfj4cc0zijtegAdx2G4waBZ99lrjjGuMnS+4mdPLz9cw9ke64Aw47TGvwxiQDS+4mVDZtgvXrE5/cDz5Ya+6TJ8MHHyT22Mb4wZK7CZUVK/Q+kWWZmN/+Fho31iRvY71M1FlyN6GSn6/3iT5zhz0XVWfP1rVWjYkyS+4mVGJn7kEkd4CrrtLx9ffcA0VFwcRgjBcsuZtQyc/XM+igulNUqwZ/+Qt89RX85z/BxGCMFyy5m1CJDYNM1Bj3slx8MbRrB/fdZ6s1meiy5G5CJYhhkKWlpcEDD2gsI0YEG4sxlWXJ3YTKihXBjJQp7fzzISsLHnssoiNnnNO/TgsWwC+/BB2NCYAldxMaW7bo2qZBn7mDloVuvx3mzYOpU4OO5gDs3AkPPQRNm+pagq1b66okAwbsGYpkUoIldxMaQY+UKW3AAGjYEJ54IuhI4lRQAB06aLOcVq3guedg+HBdleStt/S5114LOkqTIHF1hTQmEYKcwFSWgw6Ca6+Fv/4VvvtOJziF1ooVcNZZ+tFnzBi46KI937vsMu2vcNll0L+/TgO+7rrgYjUJYWfuJjSCnMBUnquu0mX/Xnkl6Ej2Y/NmuOACWLcOpk3bO7HHNG0K774L3bvrmfykSYmP0ySUJXcTGvn5erbcsGHQkexx3HHQpQu89FKIL6zecINeOB09Gk45pfztMjK09WV2Nvz611rGMUnLkrsJjdhImbSQ/a+85hpdhu+//w06kjK88YaO17z/fjj33Iq3r1NHE/yOHTBoUIj/YpmqCtmvkUllie7jHq8LL9ST3tGjg46klHXr4De/0RlXQ4bE/7rmzeHRR7WEYwP5k5YldxMaYZjAVJY6dbRU/frrWn8Pjb/8RS+gvvQSVK9+YK+97jot4dxxh9bsTdKx5G5C4Zdf4Mcfw5ncAfr2hVWrtGNkKCxdCk8/rVd8Tz75wF+flqavX70aHnzQ+/hM4Cy5m1AI2zDI0nr2hBo19Ow9FO69V8/Whw6t/D5OPRX69YMnn9RPACapWHI3oRDGYZAl1aun1yvHjQvBNchvvtHJSIMHw1FHVW1f998PP/+sNXiTVCy5m1AIe3IHyM2F5cs1twbqoYf0rP2226q+r1atdGLT3/6m6xuapGHJ3YTCihWarxo1CjqS8nXvrveBzv/57judUXX11XDkkd7s86679Oz9H//wZn8mFCy5m1DIz4ejjw7fGPeSmjWDE06Ad94JMIinntIhO3fe6d0+Tz4Zzj5bz9537vRuvyZQIf5VMqkkrMMgS+vWDT74IKAuutu2wYsvQu/e+pfGS7fdpjNWx4zxdr8mMJbcTSiEpY97Rbp31xw7fXoAB3/9dVi7VicueS03Vyc3Pfmk9/s2gbDkbgK3bZsOt47CmXvnzrrO6gcfBHDw557TBHz22d7vOy1N/2jMng3z53u/f5NwltxN4Fau1PsoJPfatXV4eMKT+9y5MHOmdnT068LEFVfoVe0XX/Rn/yahLLmbwIV9AlNpZ54JeXm6clTCPP+8NrgZNMi/Y9Svr/X8V1+1lcGTgCV3E7gojHEv6cwzobBQT6QTYscObfDVpw8cdpi/x7rmGq3rjx/v73GM7yy5m8Dl52sdu6qTLROlUydIT4cPP0zQASdO1AlGV1zh/7HOPVcX9rDSTOTFldxFpJuILBaRJSJSZm9REeknIgtE5CsR+Y+3YZpklp+v+aRaRBZ9rFMHcnISmNxffRWOOCK+fu1VlZ4OAwfClCl6ldtEVoXJXUTSgWeA7kArYICItCq1TXPgbqCTc641cKsPsZokFZVhkCV17qx19x07fD7Q+vXw9tu6Wnei/vpddpk20Bk1KjHHM76I58z9VGCJc26Zc24HMBLoXWqb64BnnHPrAZxzP3obpklmUZnAVNLpp+s1xy+/9PlAo0frX5BElGRiWrbUWasjRybumMZz8ST3xsC3JR4XFD9XUgughYj8V0Rmi0g3rwI0yW3HDu2THsXkDgno7/7vf8OJJ+pqS4nUv7++ueXLE3tc45l4kruU8VzppqfVgOZAF2AA8E8ROWSfHYlcLyJ5IpK3Zs2aA43VJKFvv9UKQNTKMk2aQOPGPif3lSvho490MWsp69fQR/376/1rryX2uMYz8ST3AqBpicdNgFVlbPOmc26nc245sBhN9ntxzr3gnMtxzuU0aNCgsjGbJBK1YZAlnX66z8n9jTf0/tJLfTxIOTIzoUMHW2M1wuJJ7nOA5iLSTERqAP2B0oNgxwFnAYhIfbRMs8zLQE1yinpyX7ZMlwf0xeuvQ3a2thwIQv/+OjN2wYJgjm+qpMLk7pwrBAYDk4GFwCjn3FciMlREehVvNhlYKyILgPeBO51za/0K2iSPFSt0Nn3j0ldxIiBWd//4Yx92vnq1zpK6+GIfdh6nfv20HDR6dHAxmEqLa5y7c26ic66Fc+4459wDxc/d55wbX/y1c87d7pxr5ZzLcs7ZZXYTl/x8rV9Xrx50JAeuXTsdnThrlg87HztWL0YEmdyPPBI6dtRYTOTYDFUTqCgOg4ypVUurJnPm+LDz11/XUTKtWlW8rZ/69NHxnjZqJnIsuZtARXECU0nt2sFnn3m8aPaaNTr99eKLEz9KprQ+ffR+3Lhg4zAHzJK7CczOnbr4T1TP3EGT+7p1e9oWe+LNN3Upvb59PdxpJR17rH48sdJM5FhyN4EpKNAcFuXk3r693n/2mYc7HTNGk+rJJ3u40yro0wdmzPBxWJDxgyV3E5jYMMgol2WysrTXlmfJfcsWeO89uPDC4EsyMX36aN3J2gBHiiV3E5jYIh1RPnOvWVOveXqW3KdO1Z4MPXt6tEMPZGfrgtxWmokUS+4mMPn5enLatGmFm4Za7KKqJyZMgHr1tO1kWIjo2fvUqbBpU9DRmDhZcjeByc/XBTpq1Ag6kqpp1w6+/96D9ufO6cIc550XvoH/vXvrJ4p33w06EhMnS+4mMCtWRLskExNr2Fjls/cvvtAWmT16VDkmz3XsCIccop8sTCRYcjeBifIEppLatNHKxaefVnFHb7+tO+re3ZO4PFWtGnTrpsl9166gozFxsORuAlFYqO1+kyG516kDJ5zgQXKfMAFOOQUaNvQkLs/17KnDIfPygo7ExMGSuwnEqlVQVBTtYZAltWmjDRQrbc0a+OSTcJZkYrp10y5vVpqJBEvuJhBRbvVbluxsfU8bN1ZyB5Mm6QXVMCf3ww/XHu9vvx10JCYOltxNIJIxuQPMn1/JHUyYoF0Y27b1LCZf9OypV45XlV6vx4SNJXcTiNgEpqiPcY+JJfdKlWZ27oTJkyE3V8seYRb7ZDFxYrBxmAqF/H+SSVb5+dCoEWRkBB2JN5o00ZGClUruM2dqPSfMJZmYk06Co4+20kwEWHI3gUiWYZAxInr2XqnkPmGCTlo691zP4/KciP4Revdd2LYt6GjMflhyN4HIz0+ekTIx2dkwb14leru//TaceSbUretLXJ7r2RO2btWe8ya0LLmbhCsqSp4x7iVlZcHmzXuuJ8Rl+XJYuDAaJZmYs87SjmlWmgk1S+4m4Vav1muIyZbcK3VRNTZmPErJvWZN6NpVk7unS1AZL1lyNwmXDH3cy3LSSXp/wMm9eXO9RUmPHvoPuWhR0JGYclhyNwmXDH3cy1KnDhx33AEk959/hvffj9ZZe0xurt7bbNXQsuRuEm75cr0/+uhg4/DDAY2Yee892L49msn96KP1o4qNdw8tS+4m4fLzdTJmrVpBR+K97Gz45hsdTFKhCRP0dP9Xv/I9Ll/06AEffVSFngvGT5bcTcItX558JZmY7GztiLtgQQUbOqfJ/bzzortaSW6utvecOjXoSEwZLLmbhFu+XJfkTEZZWXpfYWlm7lwoKIhmSSamY0c4+GCru4eUJXeTULE+7sma3I89VstNFSb3WEKMXZiMomrV4Pzzte5uC3iEjiV3k1DffacJPlnLMunp0Lq1zlTdrwkToH17vfgQZT16wA8/wOefBx2JKcWSu0mo2EiZZD1zhz0jZsqd37N2LcyeHe2STEy3btpvxkozoWPJ3SRUKiT3rCz46Sc9oS3TO+9oGSMZkvsRR+jSgDYkMnQsuZuEys/XE71k6eNellgbgnJLMxMmaFLMyUlYTL7q0UOXCFyzJuhITAmW3E1CLV+uvc+jOvovHrERM2Um98JCXVIvCgtzxCs3V2tQkyYFHYkpIa7/XSLSTUQWi8gSERmyn+36iogTkSQ5JTFeS+ZhkDH16+t10jJHzMyaBRs2JEdJJqZdO2jY0EozIVNhcheRdOAZoDvQChggIq3K2K4ucAvwsddBmuSRbIt0lCfW230fEyboEMIoLMwRr7Q0PXufPFk/mZhQiOfM/VRgiXNumXNuBzAS6F3Gdn8GHgJseRZTpu3bdShksp+5g5ZmFiwoI9dNmABnnKGTf5JJbq5+Ipk1K+hITLF4kntj4NsSjwuKn9tNRNoCTZ1z++3eLyLXi0ieiOStsYsvKWflSi3Npkpy37YNliwp8eTKlTB/fnKVZGLOPVc/kdiQyNCIJ7lLGc/tHsErImnA48D/VLQj59wLzrkc51xOgwYN4o/SJIXYMMhUKctAqdJMFBfmiNfBB+snEqu7h0Y8yb0AKDlwrQmwqsTjusBJwAcikg+cDoy3i6qmtNgiHalw5t6ypc5W3Su5v/02HH88nHBCYHH5KjdX3/DKlUFHYogvuc8BmotIMxGpAfQHxse+6Zzb6Jyr75zLdM5lArOBXs65PF8iNpG1fLl+cm/cuOJtoy4jQxdX2j1i5uefYdo0XVxayvownARin0js7D0UKkzuzrlCYDAwGVgIjHLOfSUiQ0Wkl98BmuSxfLmu8ZCeHnQkibHXiJnYwhw9ewYak69OPFFrbpbcQ6FaPBs55yYCE0s9d18523apelgmGeXnp0ZJJiYrC0aNgi1boM7bb0PdulqXTlYievY+bJheTc7ICDqilJYkU+RMFCxblnrJHWD+PKf19vPPT+6puaDJfetW+PDDoCNJeZbcTUJs3KitR5o3DzqSxNk9YmbCSli1KrlLMjFdukDNmjYkMgQsuZuEWLpU748/Ptg4EumYY3SJ1Lnv/qAli+7dgw7JfzVrwtlna3Ivt+exSQRL7iYhvvlG71MpuaelaWlm3sJ0OO007QSZCnJztQb39ddBR5LSLLmbhIjN1DzuuGDjSLSs47Yyb3MmrkcKlGRiYkMirTQTKEvuJiGWLIGjjoLatYOOJLGydn3JOg5n1akXBh1K4hxzjK41aEMiA2XJ3STEkiWpVZKJyf5Wz17nFe3TSDW55ebC9OmweXPQkaQsS+4mIVIyuW/fTtanLwMwb36SzkotT48esHMnTJ0adCQpy5K78d2WLfD99ymY3N9/n0O3fkfjw38pe+GOZNaxozYTs7p7YCy5G9+l4jBIAMaOhTp1yM6pUf56qsmqenU47zytu9uQyEBYcje+i42USankXlQEb74J3buTdXI6CxdqlSKl9OgBq1fDF18EHUlKsuRufJeKY9yZPRt++AH69CErC3bsSMFh39266f3b+13Dx/jEkrvx3ZIlun5y3bpBR5JAY8dqaSI3t+yFO1JBw4Zw+ukwblzQkaQkS+7Gdyk3UsY5Te5du8LBB3PiidrHPuWSO8BFF8Fnn8GKFUFHknIsuRvfpVxynzdPp9/36QNoI8gTTiD1RszA7p8BY8cGG0cKsuRufPXzz/DddymW3MeO1UZhvXvvfmqvhTtSyfHH65t/442gI0k5ltyNrxYv1vuWLYONI6HGjoVOnbTmXCwrSysTGzcGGFdQLroIZszQC8wmYSy5G18tXKj3J54YbBwJs2wZfPnlnnJEsdhF1fnzA4gpaH366HWIN98MOpKUYsnd+GrhQl0zNWUW6Xj9db0vldxjqzKlZGkmK0vbgVrdPaEsuRtfLVqkv9fJvrrcbiNHau/2UusJNm2qs/FT8qKqiJZmpk2DDRuCjiZlWHI3vlq4MIVKMl9/DZ9/Dpdeus+3RKBNGx0VmJIuukin6FqvmYSx5G58U1ios1NT5mLqa69pFu/Xr8xvt2+v5fjCwgTHFQannqoN/W3UTMJYcje+WbpUT9ZSJrmPHAlnnAGNG5f57XbtYNs2WLAgwXGFQVqaXoeYNMl6vCeIJXfjm0WL9D4lyjLz52vW7t+/3E3at9f7lC3N9O8Pv/wC48cHHUlKsORufJNSwyBHjtSz04svLneTFi2gTh349NMExhUmHTvqleURI4KOJCVYcje+WbhQy6wHHxx0JD5zTpN7165wxBHlbpaWBm3bpnByT0vTs/fJk2Ht2qCjSXqW3I1vUmakzMcf6wWGAQMq3LR9e21vnpIXVUF/RoWFe+YDGN9Ycje+2LULvvoKTjop6EgS4OWXoVYt6Nu3wk3btdOyc+x6RMpp00a7qFlpxneW3I0vli2DrVv3TLtPWr/8oiWZiy+Oq2F97KJqypZmRPTsffp07ShnfGPJ3fgiNhPz5JODjcN3b76p3cAGDYpr8xNOgNq1U3jEDGhyd07nBRjfWHI3vpg7V6+ftWoVdCQ+e+UVOPpo6NIlrs3T07UyMWeOv2GFWosWkJOjPztbPNs3cSV3EekmIotFZImIDCnj+7eLyAIRmSsi00TkGO9DNVEyd642C6tVK+hIfPTddzBlClx5pf4li9Ppp+uZ+44dPsYWdlddpf9JPv886EiSVoX/I0UkHXgG6A60AgaISOnzsc+BHOdcNvA68JDXgZpomTs3Bert//qXXjkeOPCAXnb66bB9u46aSVkDBsBBB8GwYUFHkrTiOd04FVjinFvmnNsBjAR6l9zAOfe+c25r8cPZQBNvwzRRsmWLjgxM6uReVATPPw9nn33Ay0x16KD3s2b5EFdUHHqoNhMbPlx7MhjPxZPcGwPflnhcUPxcea4BJlUlKBNtsQUpkjq5T5qkSyvddNMBv7RxY2jSJMWTO2hpZv16W8TDJ/EkdynjuTKvgojI5UAO8HA5379eRPJEJG/NmjVNjywxAAAR/UlEQVTxR2kiJVZGTeqRMs8+C40aQa9elXp5hw4we7bHMUXN2WfrxWgrzfginuReADQt8bgJsKr0RiJyDvB7oJdzbntZO3LOveCcy3HO5TRo0KAy8ZoImDMH6tfX39uktGwZvPMOXH89VK9eqV106KAn/qtXexxblKSn6xDSKVPg228r3NwcmHiS+xyguYg0E5EaQH9gr7ZuItIWeB5N7D96H6aJkrw8OOUUna+SlJ5/XkfHXHddpXdx+ul6n/Klmdj8gH/+M9AwklGFyd05VwgMBiYDC4FRzrmvRGSoiMQ+kz4M1AFGi8gXImI9PVPUzz9r24GcnKAj8cnPP2si6t273L7t8WjXTpceTPnk3qwZdO8OL7yQ4mNDvVctno2ccxOBiaWeu6/E1+d4HJeJqC++0NGBp5wSdCQ+GTYM1q2D22+v0m4OOkgXJ/rwQ4/iirKbb9YEP2ZMXM3XTHxshqrxVGzmZVKeuRcWwqOPal/yTp2qvLsuXbTHzKZNVQ8t0s47T2e8/e1vQUeSVCy5G0/NmaPVikaNgo7EB2PGQH4+/L//58nuunTRTzkzZniyu+hKS4Pf/lZrVCnbUc17ltyNp+bMSdKzdufgoYe089cFF3iyyw4dtO7+wQee7C7aBg3SjmpPPx10JEnDkrvxzA8/wDffaNUi6UyYoA1h7rzzgPrI7E+tWnDaaZbcAV2ua9AgnbFqrYA9YcndeCZWXjjjjGDj8NyuXXDvvXDccdokzEOxuvvGjZ7uNpruuEN/1o8+GnQkScGSu/HMRx9BzZp7FqRIGmPH6jCg+++v9KSl8sTq7h995OluoykzU0fLvPCCrbHqAUvuxjMffaRlhho1go7EQ0VFcN99uhjsZZd5vvsOHfQP4pQpnu86moYM0bkENnKmyiy5G09s2qQnt0lXknnpJViwAIYO1enyHqtZU1usTJhg61YA0Lq19ut56inYvDnoaCLNkrvxxKxZWl5IquS+fj3cc4++qTgWv66s3FxtV/PNN74dIlruuUd/9nb2XiWW3I0npk3TcnSsV3lS+NOfdDbqU0/52ignN1fvJ07c/3Yp47TT9Oz9wQet9l4FltyNJyZPhs6doU6doCPxyLx5Oub6hht00VMfZWbqWrMTJvh6mGh54AEtyzz4YNCRRJYld1Nlq1bpsnrdugUdiUd27tSFJA47DP7854QcMjdX+8xs2ZKQw4XfSSfBFVdoaaagIOhoIsmSu6myyZP1PmmS+0MP6eDz556Dww9PyCF79tS/KZNsDbM9/vQnvZAzZEjQkUSSJXdTZZMnay+ZrKygI/HAp59qUunXDy6+OGGH7dwZjjwSXnstYYcMv8xM7eMzfDhMnx50NJFjyd1Uyc6dOkb7/POTYHGO9et1VMyRR+oyegmUnq6HnjDBRgDu5e67dUmvwYO1K6eJmyV3UyXvv685sU+foCOpoqIiGDhQ+5qMGpWwckxJl14K27bBW28l/NDhVasWPP64XuC2oZEHxJK7qZLXX9cRMuedF3QkVeCcLr7x1lvw2GN71sBLsI4dtV3yiBGBHD68+vSBHj3g97+Hr78OOprIsORuKq2wUNuuXHABZGQEHU0VPPKIjmW/7Tb9+B+QtDQdIDJpkjVG3IuI9pvJyNDOkUVFQUcUCZbcTaVNeWcXP/0E/c7bAMuX6zTL1athw4bo1Ef/+le9aHfppZrkA3bttZq7hg0LOpKQOeooLcvMmgUPPxx0NJEgLqCGFjk5OS4vLy+QY5sDtHatNo758ksd0L50KXz7LX1XPsp0dwYFNKEGO/d+TVoaNGyodYbYLJ2TTtLbCSd41hO90oqKdIjdI49oQ7BXXoFqcS0p7Ltzz9VWBEuX+tLOJrqc0z/Cb7wB770Hv/pV0BEFQkQ+dc5VuCSOJXezr3Xr9Epp7LZgwZ7vNWoELVqwpkErGr/xFDd3/oJH+8/RDlgi8Msvelu/XmsL332nWWrp0j2dsQ4+WOvaHTvqrUMHXYUnUdasgV//Gt59V5d3e/LJUGXR0aN1JOb48Z4t+pQ8Nm3S1dc3bdLFU5JyPcf9s+RuDkxBAYwbp0X0Dz/UM9vatXUAdpcuunbeySdDgwaAzgofMkQHMZx0Uhz7/+UXWLRIz/5nz4aZM2H+fE341appE/gzz9SzsU6d4JBDvH+PzulA8ptv1uTw7LNwzTXeH6eKdu7U9aKPOgr++98kGGLqtfnztf9MmzYwdaqeWKSQeJM7zrlAbu3bt3cmYGvXOvfMM86ddppzmvqca9nSuXvuce6//3Vux44yX7Z9u3NHHeVc165VPP6GDc5NmuTc3Xc716mTc9WrawwizrVt69zvfufcmDHOrVlTteMUFTn3zjvOnX667v+UU5ybO7eKwfvrmWc01A8+CDqSkBo9Wv+f9OnjXGFh0NEkFJDn4sixltxTzc6dzk2c6Fy/fs7VqKH/BbKznfvf/3Vu4cK4dvHKK/qyiRM9jm3rVufee8+5P/7RubPPdq5mzT1/dFq1cu7GG50bMcK5ggLndu3a/7527HBuxgzn7r3XuWOP1X00bercCy9EIhls3ercEUfoH9CK3mrKeuIJ/Xf9zW9S6ocUb3K3skyqWLQIXn4ZXn1VO30dfrjWnQcNgrZt497Nzp16bbRmTa2w+Foy2LED8vK0TDR9ui7SGuusVaeOXqht0kS/rlkTtm/Xckt+PixZoiN2RLSsdPXVcMklcNBBPgbsrSefhFtv1eH3PXsGHU1I3XWX9gK68UZ45pngL9QngNXcjQ5JfO01TeqzZ+tFw9xcTeg9e1ZqPbznnoObboK339Z5JQlVWKijdmbO1GGXy5frH6qtW/V20EFQty40barL4rVrp8scHXZYggP1xs6d2q9n1y4tMyfV8oVecU5bFDz4oM4w/sc/PF/nNmys5p6qCgudmzzZuQEDnMvI0I+trVs798gjzq1eXaVdr1njXP36zv3qVyn1KThQEyfqP+H99wcdSYjt2uXc0KH6g+rSperXaEKOOMsy4RjYa6pu8WIdq/2vf+nww0MP1ZEggwbpSBQP6ie33QYbN+ogExvBkRjdu8Pll8Nf/qItlQPqjBBuInDvvdCsmc4CO/VU7eFw2mlBRxao5C9QJbPvvtNp8x06aBnioYd0eNjo0TpT9OmndQijB5n41Vfh3//WT8CtW3sQu4nb00/rpYX+/eGHH4KOJsQuv3zPMN5OneD++/W6TYqymnvUrFoFY8ZoAp8xQ2uOWVlw5ZV6gdSHSR0ffwxnnaUnQu++G5qJnCklL0+nAbRsCR98kETLGfph40a45Rb9FNu8OTz6qF5jSpKPm/HW3O3MPeyKirSfxr33anmlcWP9j7thAwwdqqNg5s6FO+7wJbF/+qmWBho1gpEjLbEHJSdHr41//jl07Qo//RR0RCF28MFaopw4UUfP9OqlF9anTNkzSzoF2Jl72BQW6hjDGTPgo490+v+6dTrSpUMHzbR9+ugpnM9GjdKSfYMG+mk3M9P3Q5oKjBsHAwbo7NXhw60GX6GdO/Ui0YMPaqmyTRuty/fvH0jPfi94OloG6AYsBpYAQ8r4/kHAa8Xf/xjIrGifNlrG6ciWhQt1Ys5ddzl3zjnO1amzZ+JOZqZzAwc699przq1bl7CwVq7UOU7gXIcOzn3/fcIObeIwc6ZzxxzjXHq6czff7NyPPwYdUQRs2+bciy/qhD3Q2dA9ezr37LPOLVsWdHQHBK8mMYlIOvA1cC5QAMwBBjjnFpTY5iYg2zl3o4j0B/o45y7d335T5sx92zb4/ntYuVIn1ixdqvdLlsDChdpzBXRsbuvW2kjrjDO0p0uTJgkLs7BQPyT8+9860CAtDf7wB+2Ga+Orw2fjRp2/889/6n+dfv20F3znzhHvrZ8IX36pIwTGjNEJbwDHHKMNyU45RedHtGihv38hnBTl2SQmEekA/NE5d37x47sBnHN/LbHN5OJtZolINeB7oIHbz85DndyLinS24/5umzbpb9iGDXofu23YoMn8++/1Y+D69Xvvu1o1rW8cf7yOcGnbVhtytWzpaxbdtQt+/lnX59y4UfuE5efrXKC8PL1ounmzzgEaOFBL+Mcc41s4xiOLFulM1uHD9d+vZk0t1WRl6blC06a6JGzDhvpvW6tWqBpgBss57a38zjvaoW3OHJ0YF5ORAccdp2u4Nmy453bEEVCvnl7VrltX7+vU0UZ71avvufn0h8HL5N4X6Oacu7b48RXAac65wSW2mV+8TUHx46XF25R72aeyyf2lq6bzyH+OAsA5Ady+XyPgYo8o8X39buz7UNY2Jfax+xVS7nN7fZ2WhpM0/e1JS8cV35OehktLh2rV9D62vxIHr+zX8WxbVKSJvax/6mrVtKtjhw7aR7x7dzvzi6Kff9ZPXu++q9ffFyzQ58pSo4b+EcjI0Pyzv5tI0gwyiU9RIWzbrkMoY7fCnfrRtrCIvTNGHESA4h+isPuHef8NP3Dp02dUKsR4k3s8Yx/K+qct/Q7j2QYRuR64HuDoo4+O49D7qn/UQZzUcM3uQ8Z+diBlfi27NyqxrUjZzyNIGruTs6QXJ+r0NP06LR2ppvfUqI7UqKG/KTVqINXS997f7ve879cVff9Av47n+3Xr6q1ePb1v3FjnfBx1lI2ASQa1a+tov1gPml274NtvdeTsDz/oB8ktW7QKGOvWsG2b/sHftWv/t9RSrfhW1voCDnbshO3bNNnvLCxO+iVubhfsKv6hlvV18RW1Q4+pl5B3UpECoGmJx02AVeVsU1BcljkYWFd6R865F4AXQM/cKxNwrwdOo9cDlXmlMakjLU3LalZa85IANYpv4RdPUWgO0FxEmolIDaA/ML7UNuOBgcVf9wXe21+93RhjjL8qPHN3zhWKyGBgMpAOvOSc+0pEhqJDcsYDLwKvisgS9Iy9v59BG2OM2b+4qq3OuYnAxFLP3Vfi623AJd6GZowxprLCN4jTGGNMlVlyN8aYJGTJ3RhjkpAld2OMSUKW3I0xJgkF1vJXRNYAKyr58vpAsnS0tvcSPsnyPsDeS1hV5b0c45xrUNFGgSX3qhCRvHh6K0SBvZfwSZb3AfZewioR78XKMsYYk4QsuRtjTBKKanJ/IegAPGTvJXyS5X2AvZew8v29RLLmbowxZv+ieuZujDFmPyKb3EXkzyIyV0S+EJEpInJU0DFVlog8LCKLit/PWBE5JOiYKkNELhGRr0Rkl4hEclSDiHQTkcUiskREhgQdT2WJyEsi8mPxKmmRJSJNReR9EVlY/H/rd0HHVFkikiEin4jIl8Xv5U++Hi+qZRkRqeec21T89S1AK+fcjQGHVSkich7aA79QRB4EcM7dFXBYB0xEWgK7gOeBO5xzIV0kt2zxLAYfFSLyK2AL8C/n3ElBx1NZItIIaOSc+0xE6gKfAhdG9N9EgNrOuS0iUh2YAfzOOTfbj+NF9sw9ltiL1eaAFzcMD+fcFOdcYfHD2ehqV5HjnFvonFscdBxVcCqwxDm3zDm3AxgJ9A44pkpxzk2njNXQosY5t9o591nx15uBhUDjYKOqHKe2FD+sXnzzLW9FNrkDiMgDIvIt8Gvgvoq2j4irgUlBB5GiGgPflnhcQEQTSTISkUygLfBxsJFUnoiki8gXwI/Au845395LqJO7iEwVkfll3HoDOOd+75xrCgwHBgcb7f5V9F6Kt/k9UIi+n1CK531EWFwLvZvEE5E6wBjg1lKf2iPFOVfknGuDfjo/VUR8K5mFet1759w5cW76H2ACcL+P4VRJRe9FRAYCPYGuYV5/9gD+TaIonsXgTYIV16fHAMOdc28EHY8XnHMbROQDoBvgy0XvUJ+574+INC/xsBewKKhYqkpEugF3Ab2cc1uDjieFxbMYvEmg4ouQLwILnXOPBR1PVYhIg9hIOBGpCZyDj3kryqNlxgAnoKMzVgA3Oue+CzaqyileWPwgYG3xU7OjOPJHRPoAfwMaABuAL5xz5wcb1YERkVzgCfYsBv9AwCFVioiMALqg3Qd/AO53zr0YaFCVICKdgY+AeejvOsA9xes6R4qIZAOvoP+30oBRzrmhvh0vqsndGGNM+SJbljHGGFM+S+7GGJOELLkbY0wSsuRujDFJyJK7McYkIUvuxhiThCy5G2NMErLkbowxSej/A2QntWHSphnaAAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VOX1+PHPIQEiCEUhbmzBCAY1LBpFBbcviKiIsoobWPdat9b2V2y/ar+0WNu6oNVaF3DBBQUj4lYEVAJIUBBUZFGUKGGRCLKJEJI8vz9OBkIIZJLcmXvnznm/XvOaZHLn3jNZTp4597nnEeccxhhjwqWe3wEYY4zxniV3Y4wJIUvuxhgTQpbcjTEmhCy5G2NMCFlyN8aYELLkbowxIWTJ3RhjQsiSuzHGhFCqXwdu0aKFy8jI8OvwxhiTkObPn/+Dcy69uu18S+4ZGRnMmzfPr8MbY0xCEpFvo9nOyjLGGBNCltyNMSaELLkbY0wI+VZzN8aE286dOyksLGT79u1+h5KQ0tLSaNWqFfXr16/V8y25G2NiorCwkCZNmpCRkYGI+B1OQnHOsX79egoLC2nXrl2t9mFlGWNMTGzfvp3mzZtbYq8FEaF58+Z1etdjyd0YEzOW2Guvrt87S+7G7ENpKUycCA88AKtW+R2NMTVjyd2YKpSVwSWXwODBcPvtcNxx8NlnfkdlYuWuu+5i2rRpfofhKUvuxlRh9GiYMAHuuQcWLYJGjeCyy6C42O/IjNdKS0sZOXIkvXr1qtFzgs6SuzGVbNgAI0fCeefBiBFw7LHwyCOa5MeP9zs6UxMFBQVkZWUxfPhwOnXqxKBBg9i2bRsZGRmMHDmSHj16MGHCBK688komTpwIwPTp0+natSvZ2dlcddVV7NixA2Cv5wSdTYU0ppKHHoLNm+HeeyFyTuuii+CYY+D+++GKK3Y/bqJ0222wcKG3++zSRd9iVWPZsmWMGTOG7t27c9VVV/Hvf/8b0Hnks2bNAuC///0voDN8rrzySqZPn06HDh0YNmwYjz32GLfddttezwk6G7kbU0FJCYwZA+eeC9nZux8X0fz02Wdg/e4SS+vWrenevTsAl19++a7kfPHFF++17bJly2jXrh0dOnQAYPjw4eTl5e36elXPCSobuRtTwZQpOjPmX//a+2sDB8KNN2ot/sQT4x9bQotihB0rlacURj5v3LjxXts65/a7r6qeE1Q2cjemgrFj4ZBDoG/fvb928MFw9tma3KvJASZAvvvuO+bMmQPASy+9RI8ePfa5bVZWFgUFBSxfvhyAcePGccYZZ8QlTq9Zcjem3LZt8M47Ov1xX+08LroICgpg2bK4hmbqoGPHjjz77LN06tSJDRs28Ktf/Wqf26alpfH0008zePBgsrOzqVevHjfccEMco/WOlWWMKTdtGvz8sybwfTn77N3bZmXFJy5TN/Xq1eM///nPHo8VFBTs8fkzzzyz6+OePXuyYMGCvfZT+TlBZyN3Y8pNmgS/+AXs7114u3Zw5JEwdWr84jKmNqpN7iIyVkTWiciifXxdRORhEVkuIp+JyPHeh2lMbJWVwZtv6tz26jqs9uwJM2boc0ywZWRksGhRlakr9KIZuT8D9NnP188F2pffrgMeq3tYxsTXp59CUZFOgazOqafCpk1WdzfBVm1yd87lARv2s8mFwHNO5QPNRORwrwI0Jh6mT9f7nj2r3/bkk/U+Pz928RhTV17U3FsCKyt8Xlj+mDEJY/p0PUF6xBHVb9uhg9bm586NfVzG1JYXyb2qC7GrnAUsIteJyDwRmVdUVOTBoY2pu+JiyMuLbtQOUK8edOtmI3cTbF4k90KgdYXPWwGrq9rQOfeEcy7HOZeTnp7uwaGNqbu5c3WOew2aAtKtG3z+OWzdGru4TDBMmjSJxYsX1/h5kydP5t57741BRNHxIrlPBoaVz5o5GdjknFvjwX6NiYv33tPeMTW5EDEnR2fLWI/38KtNci8pKaFfv36MGDGiRs/xUrUXMYnIS8CZQAsRKQTuBuoDOOf+A7wNnAcsB7YBv/Q0QmNi7MMPtUnYQQdF/5zOnfX+s8909owJpoKCAvr06UO3bt1YsGABHTp04LnnnmPOnDn87ne/o6SkhBNPPJHHHnuMhg0bMmLECCZPnkxqaiq9e/dmwIABTJ48mRkzZvDXv/6VV199FYBf//rXFBUV0ahRI5588kmysrK48sorOfjgg1mwYAHHH3882dnZzJs3j0ceeYRvv/2Wq666iqKiItLT03n66adp06bNXs+5//77PXvt1SZ359wl1XzdAb/2LCJj4qi0FObM0YU4aqJNGz2p+umnsYkrbHzs+LtXy98HHniAxx9/fK+2vsOGDeO1115j6dKliAgbN26kWbNm9OvXj759+zJo0CBAr2D9z3/+Q/v27Zk7dy433ngj7733HgBffvkl06ZNIyUlZY+rXm+66SaGDRvG8OHDGTt2LLfccguTJk3a6zlesitUTVL74gvYsgXKO8JGTQQ6dbKyTCKo3PJ3+vTpVbb1bdq0KWlpaVxzzTXk5ubSqFGjvfa1detWPvzwQwYPHkyXLl24/vrrWbNmdxV68ODBVSbpOXPmcOmllwJwxRVX7NETfl/PqSvrLWOS2uzZel+b0kqnTvDcc1p7r2fDpP3ysePvXi1/9yU1NZWPPvqI6dOnM378eB555JFdI/KIsrIymjVrxsJ9vA2JtiVwxZhi1UbYfiVNUvvwQzjsMO0ZU1OdO+uoP8H6SSWdyi1/e/XqVWVb361bt7Jp0ybOO+88Ro8evSuBN2nShC1btgDQtGlT2rVrt2uZPeccn0ZRmzv11FMZX75G4wsvvLDftsNeseRuktrs2Tpqr82yeZ066b2VZoKtcsvf3/zmN1W29d2yZQt9+/alU6dOnHHGGTz44IMADB06lH/+85907dqVr7/+mhdeeIExY8bQuXNnjj32WF5//fVqY3j44Yd5+umn6dSpE+PGjeOhhx6K9ctGqlt5JFZycnLcPFuvzPhozRq9IvW+++D222v+/C1boGlTuOceuOMO7+NLdEuWLKFjx46+xlBQUEDfvn0TtnlYVd9DEZnvnMup7rk2cjdJq/ydeo1PpkY0aaL/HJYu9S4mY7xiyd0krdmzoWFD6Nq19vvIyrLukEFmLX+NSUIffqgLXTdsWPt9ZGXpyN3WVK2aX2XfMKjr986Su0lKP/8M8+fX/erSrCzt7f79997EFSZpaWmsX7/eEnwtOOdYv349aWlptd6HzXM3SWn+fNi505vkDjp6P+ywuscVJq1ataKwsBDrAFs7aWlptGrVqtbPt+RuklKkXW9k4Y3aiiT3ZcvgzDPrtq+wqV+/Pu1qcwGB8YSVZUxSys/XC5cOPbRu+2nZEho1shkzJngsuZuklJ9f91E7aNuBo4+25G6Cx5K7STqFhbBqFZxyijf7s+RugsiSu0k6XtXbI446Cr77TpfrMyYoLLmbpJOfr3PbIwtu1FVmpnaG/PZbb/ZnjBcsuZukk58PJ5wADRp4s7/MTL3/+mtv9meMFyy5m6RSXKxz3L0qyYAldxNMltxNUvnsM9i+3dvkfvjhcMABltxNsFhyN0nF65OpoL3gjzzSkrsJFkvuJqnk52ub3jpc1V2lzExL7iZYLLmbpBK5eKk2Ky/tT2YmfPONdYc0wWHJ3SSNdet0dO1lSSYiM1M7Ta5d6/2+jakNS+4macydq/deXZlakc2YMUFjyd0kjfx8SE2F44/3ft+W3E3QWHI3SSM/X69KbdTI+323batNxCy5m6Cw5G6SQkkJfPRRbOrtoFe7tmljyd0EhyV3kxQWLoStW+G002J3jHbtYMWK2O3fmJqIKrmLSB8RWSYiy0VkRBVfbyMi74vIAhH5TETO8z5UY2pv1iy979EjdsfIyLDmYSY4qk3uIpICPAqcCxwDXCIix1Ta7H+BV5xzXYGhwL+9DtSYupg5U0fWLVvG7hht28Lq1bBjR+yOYUy0ohm5nwQsd85945wrBsYDF1baxgFNyz/+BbDauxCNqRvndOQey1E76MgdYOXK2B7HmGhEk9xbAhV/XQvLH6voz8DlIlIIvA3c7El0xnjgq6/0AqZY1ttBR+4ABQWxPY4x0YgmuVd1oXbli6wvAZ5xzrUCzgPGiche+xaR60RknojMKyoqqnm0xtTCzJl6H6+Ru9XdTRBEk9wLgdYVPm/F3mWXq4FXAJxzc4A0oEXlHTnnnnDO5TjnctLT02sXsTE1NGsWNG8OWVmxPU7LljrX3UbuJgiiSe4fA+1FpJ2INEBPmE6utM13QE8AEemIJncbmptAmDlTR+1eNwurrH597TZpI3cTBNUmd+dcCXATMAVYgs6K+UJERopIv/LNbgeuFZFPgZeAK52z/njGf6tW6YVFsa63R7RtayN3Ewyp0WzknHsbPVFa8bG7Kny8GOjubWjG1N306Xrfs2d8jpeRATNmxOdYxuyPXaFqQm36dGjRAjp1is/x2raFwkLYuTM+xzNmXyy5m9ByTpP7WWfpic54yMiAsjItBxnjJ0vuJrS+/FKTbLxKMmBz3U1wWHI3oRXvejvYXHcTHJbcTWhNn65teCMLacRD69Y65dJG7sZvltxNKJWWwvvv66g91vPbK2rYEA4/3Ebuxn+W3E0ozZ8PP/4IvXrF/9gZGTZyN/6z5G5C6a23dIbMOefE/9ht29rI3fjPkrsJpbfe0iX1mjeP/7EzMuC777Q0ZIxfLLmb0FmzRssyffv6c/y2bXXN1tV+rGrgHLzxBlx4IRxyCDRpAscfD/feq3UqkzQsuZvQebu8Ucb55/tz/Mhc9+++i/OB162Dc8+Ffv3g44/1/uqroXFjuOMObYv55ptxDsr4xZK7CZ233tLujNnZ/hw/ktzjWnf/6ivo1g3y8uDhh/U/y1NPwejR2hZz3jydxtOvHzz2WBwDM36JqnGYMYlixw6YOhUuvzy+UyAratNG7+OW3Neu1TPHW7dq17ITT9x7mxNOgDlzYMgQuPFGaNQIhg+PU4DGDzZyN6HywQea4/wqyYBWQZo3j1Ny37FDR+Pff6/1qKoSe8QBB0Burk7+v+YamD07DgEav1hyN6HyyivQtKk/89srats2TjX3O+/U+vrzz+8/sUfUrw8TJ2qAl10GmzbFPkbjC0vuJjSKi3VgeuGFkJbmbyxxmev+wQdw331www3Qv3/0z2vWTP8ZFBbCTTfFLDzjL0vuJjTefRc2boSLL/Y7kt3JPWbrkRUXw69+BUceCfffX/Pnn3wy/PGPmuTff9/7+IzvLLmb0Hj5ZTjoIDj7bL8j0ZOqP/0EGzbE6AAPPQRLl+rMmEaNarePO+6Adu109G6ri4SOJXcTCtu3w+uvw4AB0KCB39HEeK57URGMHAkXXADnnVf7/RxwgE6VXLwYHn/cu/hMIFhyN6HwzjuwZUswSjIQ47nu//gHbNum93V1wQVw+ukwahT8/HPd92cCw5K7CYWXX4b0dF1SLwhiltzXroVHH9WZLllZdd+fiL4LWLvWLm4KGUvuJuH99JO2Uxk4EFIDclle8+Za9fA8uf/tb3oy9a67vNvnGWfo3Pd779WLBEwoWHI3Ce/NN7VKEZSSDOiA2PO57j/8AE8+CcOGwVFHebhj4P/+T2v5Y8d6u1/jG0vuJuG9/LK2TTntNL8j2ZPnc90ff1zr4rff7uFOy3XvDqecoidYrVdxKFhyNwlt82a96n7wYEhJ8TuaPXma3HfsgEcegT594NhjPdppJbffDitWwKRJsdm/iStL7iahTZ6seS9IJZmINm200rFtmwc7Gz9eT3r+9rce7GwfLrqo9hdFmcCx5G4S2vjx0Lq1XnAZNJEZMytX1nFHzsGDD8Jxx8W2aU5KCtx6q3aPnD8/dscxcWHJ3SSsH3/UlgNDhuh6qUHj2XTI/Hz49FO45ZbY9zEeNkyn+dhFTQkvqj8JEekjIstEZLmIjNjHNkNEZLGIfCEiL3obpjF7e+01vWo+iCUZ8DC5P/mk9hEeOrTOMVWrWTM9zosv6gkNk7CqTe4ikgI8CpwLHANcIiLHVNqmPXAH0N05dyxwWwxiNWYPL7+sJeKcHL8jqdoRR2ilo07JffNmfaGXXKLrocbD9dfrxQMv2hgtkUUzcj8JWO6c+8Y5VwyMBy6stM21wKPOuR8BnHPrvA3TmD0VFcH06Tpq92vFpeqkpkLLlnVM7i++qGdkr73Ws7iqddJJ0LmzlmZi1tbSxFo0yb0lUPGUUGH5YxV1ADqIyGwRyReRPl4FaExVcnN1OnZQSzIRdb6Q6cknoVOn6Bbi8IqIjt4XLrQTqwksmuRe1bio8r/zVKA9cCZwCfCUiDTba0ci14nIPBGZV1RUVNNYjdllwgTo0EHzXpDVaa77J5/o7dpr4//2ZOhQaNgQxo2L73GNZ6JJ7oVA6wqftwJWV7HN6865nc65FcAyNNnvwTn3hHMuxzmXk56eXtuYTZIrKtL1JQYPDm5JJqJNG13wqKSkFk9+9lntX3zZZZ7HVa2DDtKOkS++aL3eE1Q0yf1joL2ItBORBsBQYHKlbSYBZwGISAu0TPONl4EaE5GbC2VlmtyDrm1bLR+trjwcqk5JiU7iv+ACTbR+GDZM+9n897/+HN/USbXJ3TlXAtwETAGWAK84574QkZEi0q98synAehFZDLwP/N45tz5WQZvkNmECtG8f/JIM1GHRjmnTYN06f0btEX36QIsWVppJUFE1SHXOvQ28Xemxuyp87IDflt+MiZlISWbEiOCXZGDPue49etTgiS+8oHPO67LSUl3Vr69TMJ94Qq8Y8+sdhKmVAF7XZ8y+vfZa4pRkQFsjQA1Pqv70k77QwYP1pKafrrhCm/dMmOBvHKbGLLmbhDJhgrYy79zZ70ii07ixVjZqlNxff10T/OWXxyyuqOXkwNFH2wVNCciSu0kY69cnziyZimo81/3553XIX6M6ToyI6MUEeXmwZo3f0ZgasORuEsbbb+vMk/79/Y6kZmo0172oSLuhXXZZcLqhXXyxXqk6caLfkZgaCMhvjzHVe+MNOOwwOOEEvyOpmTZtNLlHdSX/a6/pf7B4NAmL1jHHaLvhV17xOxJTA5bcTUIoLtbp1n37BmdAG622bbU9zPpoJgcHdZ7nkCEwa5ZekWUSQoL9mZhklZcHW7boNT2JJuq57j/8ENyTCpEmPlaaSRiW3E1CeOMNSEuL7UJEsRJ1X/dJk7QkM2hQzGOqsQ4doEsXbT9sEoIldxN4zmly/5//gUaN/I6m5tq00ftqk/uECZCZqUk0iIYM0VWhPFv128SSJXcTeEuWwIoViVmSAWjeXP8p7Tcnrl+vDeqDWJKJiJRm7IKmhGDJ3QTeG2/ofd++/sZRWyJRzHWPlGSCfOntkUfqVCVL7gnBkrsJvClTIDsbWrXyO5Laq3au+8SJ0K4ddO0at5hqZeBA+OgjmzWTACy5m0Dbtg1mz4bevf2OpG4ic92rtGGDdoEMckkmYsAAvZ80yd84TLUsuZtAy8vTOe5nn+13JHXTtq3OdPzppyq++Prr2r89yCWZiKOP1ouacnP9jsRUw5K7CbSpU3UxotNO8zuSuolMh1y5soovTpgAGRmJc+ntgAEwY4b+tzKBZcndBNrUqZrYE3EKZEX7nOu+caOWZAYNCn5JJmLAAO27PLnygmwmSCy5m8BaswY+/zzxSzKwn7nub76pa5QOHBj3mGqtSxd9p2GlmUCz5G4Ca9o0vQ9Dcj/iCEhJqSK55+ZCy5Zw0km+xFUrIjp6nzoVNm/2OxqzD5bcTWBNnaoLXQT1gs2aSE3VqZx7JPefftJuaP37J143tP799Uz3O+/4HYnZhwT7jTLJwjkduffqlXh5b1/2upBpyhT4+efd0wsTySmnwKGHWmkmwELyZ2PC5quvtOZ+1ll+R+Kdvea65+Zqb4JEnAqUkgIXXQRvvQXbt/sdjamCJXcTSDNn6v3pp/sbh5fatoVVq3RKO8XF2lfhwgu1ZpOIBgzQ0tLUqX5HYqpgyd0EUl4epKfrNTNh0batto9ZtQptErZ5c2KWZCLOPBOaNbPSTEBZcjeBNHOmVisSZep3NPZYtCM3F5o0gZ49fY2pTho00FadkyfrdE4TKJbcTeAUFmqL30QsRe/PrrnuK8q0N8v55+sKJIlswADtjZOX53ckphJL7iZwwlhvhwrJPe9bvXQ/kUsyEb176+XDr77qdySmEkvuJnBmztSKRefOfkfirUaNdPbgijlroWFDOPdcv0Oqu0aN9HVMmqQtCUxgWHI3gZOXB6eeqrPtwiYz0/H11w7OOQcOPNDvcLwxcKDOW83P9zsSU0FUyV1E+ojIMhFZLiIj9rPdIBFxIpLjXYgmmaxfD198Eb6STMSRv1jP1ztahaMkE3H++Xpy1UozgVJtcheRFOBR4FzgGOASETmmiu2aALcAc70O0iSP2bP1PmwnUyMyNy+kkFbs6J2gC8JWpWlTbQCUm6uXFptAiGbkfhKw3Dn3jXOuGBgPXFjFdn8B/gHY5Wqm1vLydBB44ol+RxIDzpG5fAqOeqzYdLDf0XhrwAAoKIAFC/yOxJSLJrm3BCouMVBY/tguItIVaO2ce3N/OxKR60RknojMKyoqqnGwJvxmzoRu3RJ/hmCVFi8m83t9a/L11z7H4rV+/fQkiV3QFBjRJPeqLiPZ9d5LROoBDwK3V7cj59wTzrkc51xOenp69FGapLB1K8yfH96SDLm5ZPINEMLk3qIFnHGG1d0DJJrkXgi0rvB5K2B1hc+bAMcBH4hIAXAyMNlOqpqays/Xy/PDejKV3FwOOfUoGjcOYXIHnTWzdCksWeJ3JIbokvvHQHsRaSciDYChwK71tZxzm5xzLZxzGc65DCAf6OecmxeTiE1ozZyp7X1POcXvSGLgm29g4UJk4AAyM0Oa3C+6SO9t9B4I1SZ351wJcBMwBVgCvOKc+0JERopIv1gHaJJHXp4uzNG0qd+RxMBrr+l9//5kZmquD50jjtALFKzuHghRzXN3zr3tnOvgnMt0zo0qf+wu59xeK+Q65860UbupqeJiLcuEtiTz6qv6n6tdu13JPZQXdA4YoDNmQvnfK7HYFaomEObP1zUfQnky9bvvYM4cGDwYgMxM2LEDVq+u5nmJKHJxVuSdivGNJXcTCJGmgqFM7hMn6n2F5A4hrbu3awddu1rdPQAsuZtAmDkTsrJ0gY7QeeUVTXjt2wMhT+6gs2bmzAnpW5PEYcnd+K60FGbNCumovaAA5s6FIUN2PdSmja6sF9rkbqWZQLDkbny3aBFs2hTSk6mVSjKgib1tW1i+3KeYYq1jR73ZrBlfWXI3vosszhHKkfsrr8AJJ+yuxZTr0AG+/NKnmOJhwACYMUMXJTG+sORufJeXB61b715jNDRWrICPP96jJBNx9NGa3EM5HRK07l5aaqUZH1lyN75yTkfuoSzJTJig9xVKMhFZWbBtm64XG0pduugJ5Jdf9juSpGXJ3fjq669h7doQl2ROOkmnB1Zy9NF6v2xZnGOKFxEYOhTef19/wCbuLLkbX0Xmt4du5P7113plVhUlGUiC5A6a3MvKdr+DMXFlyd34auZM7RableV3JB6LlCMGDaryy4cdpj10li6NY0zxdswx0KkTjB/vdyRJyZK78VVeHvTooe/iQ8M5GDdOX9g+zhKL6Og91CN30NH7hx/Ct9/6HUnSseRufLN6tfaXCl1J5pNPdEh++eX73SwpkvvFF+u9nViNO0vuxjeh7Sfz/PO6EOw+6u0RWVmwciX89FOc4vLDkUfqSWUrzcSdJXfjm7w8aNJEZ82FRkkJvPQS9O0LBx20300jJ1VDfTETaGlmwYIkeJsSLJbcjW9mzIDu3fVy/NCYNg2+/77akgwkyYwZ0HcwIjZ6jzNL7sYXRUWweLGuqRwq48bpiP2886rdtH17XVYw9EuOtmypP+gXXtCTzSYuLLkbX0T6yYQquW/ZopfbDxkCDRtWu3lamracWbQoDrH5bfhw+OorbQVs4sKSu/HFjBlwwAHaUys0XnsNfv4Zrrgi6qdkZ8Pnn8cwpqAYOBAaNYJnnvE7kqRhyd34YsYMXUu5QQO/I/HQmDE6FD/11Kifkp2trX+3bYthXEHQpIle0PXyy/oP0MScJXcTdz/+CJ99FrL57cuW6fSfa66p0RVZ2dlahl68OIaxBcXw4bB5M0ya5HckScGSu4m7WbM0oYWq3v7UUzrt58ora/S07Gy9T4rSzJln6jJUVpqJC0vuJu7y8rQc062b35F4pLgYnn0WLrhAm8bUQGamnlhNiuRerx4MG6bTRVet8jua0LPkbuJuxgxN7Glpfkfikddf17md115b46empGh/raRI7qClmbIynTJqYsqSu4mrzZu19UqoSjJPPqnlht69a/X07OwkmQ4JcNRR2lBt7NgQL0MVDJbcTVzNmKGrr/Xs6XckHlmxAqZOhauu0mF4LWRn63oWSbPc6PXX65z3997zO5JQs+Ru4mraNJ3ffsopfkfikUcf1aR+9dW13kXkpOpnn3kUU9ANGgTNm8Njj/kdSahFldxFpI+ILBOR5SIyooqv/1ZEFovIZyIyXUTCttSx8ci0adoFMooLOINv61adJTNoELRqVevdRBqnLVjgUVxBl5am73Ref91OrMZQtcldRFKAR4FzgWOAS0TkmEqbLQBynHOdgInAP7wO1CS+1at1PnevXn5H4pFnn4VNm+DWW+u0m0MOgdatYd48j+JKBNdfr/W5p57yO5LQimbkfhKw3Dn3jXOuGBgPXFhxA+fc+865yDV2+UDthzEmtKZP1/tQJPeyMnj4YTjxRDj55DrvLidHl1xNGpmZcM45ejK6pMTvaEIpmuTeElhZ4fPC8sf25WrgnboEZcJp2jQttXbu7HckHpgyRRux33abJ2sE5uToOcaNGz2ILVH86ldalpk82e9IQima5F7Vb26VfTtF5HIgB/jnPr5+nYjME5F5RUVF0UdpEp5zmtx79tRrWRLeAw/A4YfvcwHsmsrJ0ftPPvFkd4nh/PN1jdnRo/2OJJSi+TMrBFpX+LwVsLryRiLSC/gT0M85t6OqHTnnnnDO5Thv0bnzAAAQ3klEQVTnctLT02sTr0lQS5ZozT0UUyA/+kj/U/3mN551Pot0x0yquntqqp6vmDlTv6fGU9Ek94+B9iLSTkQaAEOBPd5HiUhX4HE0sa/zPkyT6N5+W+/79PE3Dk+MGqULctxwg2e7bN4c2rVLsuQO2mjtF7+A++/3O5LQqTa5O+dKgJuAKcAS4BXn3BciMlJE+pVv9k/gQGCCiCwUESuimT289ZbO527Txu9I6ujzz7VGfOut2sbWQzk5MHeup7sMviZNdObMxIlQUOB3NKESVfXTOfe2c66Dcy7TOTeq/LG7nHOTyz/u5Zw71DnXpfzWb/97NMlk0ybtBHn++X5H4oF77oEDD4Sbb/Z81927w3ffwcqV1W8bKjffrCdiHnrI70hCJQyntkzAvfuuznZL+OT+6ae62MRNN8HBB3u++x499H72bM93HWytWsEll+i0SJto4RlL7ibm3npLS9QeTAf31x13aH34//2/mOy+c2d9UzBrVkx2H2x//KOu0HTffX5HEhqW3E1MlZbCO+/oidTUVL+jqYMZM/SF3HGH/qeKgdRU7bmTlMk9K0tH7488AutsToYXLLmbmJo9W/9W+yXyWRjnYMQIOOKImNTaK+rRQxuIbdoU08ME0513wvbt8M8qL5MxNWTJ3cTUxInaJ6pvX78jqYNx4yA/H/7yF21pGUM9euj/kg8/jOlhgunoo+HSS7XT5vff+x1NwrPkbmKmrAxefRXOPVdryQlp40b4/e/1hEEN10etjZNP1o6Z06bF/FDBdOedumzhyJF+R5LwLLmbmJkzR69KHTzY70jq4K67dBWNRx+NS9+ERo20JfKUKTE/VDB16KAXhz3+uLYQNbVmyd3EzIQJOgpN2JJMfr4m9RtugOOPj9the/eGL76AwsK4HTJY/vxnfat3++1+R5LQLLmbmCgt1eTep4/nF3LGx7ZtMGyYzsG+5564Hvqcc/T+3XfjetjgaNFCyzP//a/eTK1YcjcxMXWqlmSGDfM7klr6wx+0B+8zz+jc9jjKztaGk0lbmgG9UCwzU1sq76iyD6GphiV3ExPPPqsXcSbkValvvKHzrX/zGzjrrLgfXkRH7+++q+cWk1LDhvozWLYs7u+cwsKSu/Hcxo3w2ms6qy3h1kr96iu4/HLtwetjUhkwQL+PkdWrklKfPnDZZfC3v+lJCFMjltyN555/Xt9JDx/udyQ1tGkT9O8P9evrHM60NN9C6d0bmjbV8xZJ7cEH9RtxzTW2HF8NWXI3nior03fTJ564e3WhhLB9O1x4oS6d9/LLukKQjxo21Kt6J02CnTt9DcVf6en6C5WfD3/9q9/RJBRL7sZT06ZpmfSWW/yOpAZ27IChQ7V/zLPPBma5qMGD4ccfk7w0A/qzGTZMrxCeOdPvaBKGJXfjqYcegkMPTaALl7Zt0xH766/Dv/6lzasC4pxz9KT02LF+RxIAjzyiS1Vddpm1BY5SIvfpM/HmHKxfDytW6Ko569fD5s1aq96xgwXr2/D227cw8qz3aThmib6lPuQQLXG0aRO8lbELC3WB648/hqeegquv9juiPTRsqOct/vUvbbVy6KF+R+SjJk20XNajBwwcqG8RPVq/NqzEOefLgXNycty8pFswMoE4p6taf/wxfPIJzJ+v7Qq3bNl725QUaNiQ/ttf5P2yMyggg2ZUamuYlgbt22tzqOOP312Ub9YsPq+nsqlTdVbMtm3aGOyii/yJoxpLl0LHjjpx5447/I4mAF58UUfv116rLQpE/I4o7kRkvnOu2jNaNnI3u333nSa96dPhvfd2d+Zr3Bi6dNG651FH6dvjjAwdlTdtCo0asfBTYVJXuPtuaPa/P+ioft063ceKFVqI//JLWLBAW0VGHHUUdOumHbNOOQU6ddLZKrHyww96gdLYsdpD/IMPNHsGVFYWnHmm5rHf/z7Be+J74dJLYdEinR7Zpg387//6HVFwOed8uZ1wwgnOBMC33zp3333OnXSSczped+6ww5y79FLnnnrKucWLnSsp2e8uysqcO+ss5w4+2LkNG6I45oYNzk2d6tw99zjXv79zhx+++9hpac716OHc73/v3KuvOrdqlTevc80a5+6807kDD3QuJcW5P/zBuW3bvNl3jE2apN+aceP8jiQgSkudGzZMvymjR/sdTdwB81wUOdbKMslo5UqdQD1hgk4xA71oZ8gQ7fLVsWON3u6++qqWrh99FG68sRbxOKf17zlzNJ78fC0DRS7PbN1aR/Unn6y3rl2rn4Me2ee778Kbb+qtpEQDHTky0KP1ysrK9I3Tzp06aE1J8TuiACgpgYsvhtxcGD0abr3V74jiJtqyjCX3ZFFYqOWQV17RJApa+x48WG+ZmbXa7aZN2gvloIM0H3tWNtixAxYu3J3s58yBb7/Vr6WmakOvVq2gZUstGzVsqNnvxx+1FLRkiZaGQLcZOhSuu05byiagiRP1xxTA877+2bFDyzS5ubpS1j33JEUN3pK7gVWrNCtMmKDr3YEOAYcM0Uxx1FF1PsSwYXqOa/ZsLZ3H1Jo1MHcuzJunib6wUF/jzz/rH3pqqv6XadFCi9XHHacF6+OOS/g/eufg9NN3n7rw6zx04JSW6tvFJ57Q3+kxYxK0DWn0ok3uVnMPm1WrnHv4Ya1bR+rYnTs7N2qUc19+6emhxo3T3d99t6e7Nfswf75z9eo598tf+h1JwJSVOff3v+s3p2NH5xYt8juimCLKmrsl9zAoKHDu/vudO/XU3Qk9O9u5v/zFuaVLY3LIWbOca9DAudNPd664OCaHMFX405/0xzthgt+RBND06c6lp+sv5qhRof3FjDa5W1kmETmnXfLefFPPZka+j1266AnDgQO1LBEjn3yija0OOkjL4c2bx+xQppKdO3UZvs8/1yvx47hAVGL4/nu4+WYtRR57LPzjH7qIb4KX5SqKtiwTsEsGzT5t3qwnjq69Vuf3ZmfrVS0i8Pe/w/LlOof8T3+KaWJ//31tcX7ggbpIjiX2+KpfX9spN2+uHXE/+cTviALm0EN10kBurjaDO/987RU0daoOipKIjdyDat06PUs5c6beFizQk0dNm0KvXjoaOeccnSYYBzt36iDorrt0wsm778bt0KYKX34JZ5+tk4PGjdP2OKaS4mI90TpqFKxdqyfWr79eJxQccojf0dWap7NlRKQP8BCQAjzlnLu30tcbAs8BJwDrgYudcwX726cl9wrWrtVpf5HbggX61ws6n7tbN+2p0bu3zveO5RWclZSVafXnj3/UStDQofr3EvIJCQmhsFDbAi9YAFdeqR1xW7b0O6oA2rFD+9KMHq3frJQUHc2ff74OlGp4XYffPEvuIpICfAmcDRQCHwOXOOcWV9jmRqCTc+4GERkK9HfOXby//SZVci8u1gS+ejV8842WUCK3r77SS+IjMjK0dn7KKZrQTzgh7ssZOadT7nJz4bnn9ON27fRv44ILEurvIPSKi+HPf4b77tOcNXSoJvru3a1VQZUWLYKXXtKa/Fdf6WOHH64DqOOP19uxx+o1FAH9BnqZ3E8B/uycO6f88zsAnHN/q7DNlPJt5ohIKrAWSHf72Xngk3tpqf7H37FDa3eRjyvetmzRq3git40bd3+8bp3Oy16zZvfFNBW1bq3zzDMz9e1ily7QuXNcJjCXlmroW7ZoKX/NGm0rU1Cgbxw+/lj/F4FeEHrzzTqFOI5vGEwNrVih7VZeegm2btV3Vqedpq16jj5axwwtWmijzoMPtp8loL/w06ZpH6X583e/WwZN7K1bw5FH6tuh9HT9Bka+ieU9lWjcWG+Rj9PS9L9sDEdAXib3QUAf59w15Z9fAXRzzt1UYZtF5dsUln/+dfk2P1S1T6h9ch/7yzzue/EIAJwTwFX4GMDhkMjDVHx1e2y/a5vy7fe1za69yj4f2+PjevVwUg/qpUBKPVy9VEipBykpuHop+oNPSa10TOL28c6d2gixKiKaCE48Ud84XHCBDmBM4vjpJ3jrLT3xPXOm5quqVnKqV0/fEEZuDRroY5GcJLLnbX+PhUZZKWzfAcU7oHgn7CzW+5ISKC2p+QnZPb5JArL747uvX8vFj5xWqzC97ApZ1Y+v8quMZhtE5DrgOoA2bdpEcei9tTiiIccdWrTrkJHvG8juXzQRRByRb6iUfz3y/dXfzgofA7Jr3pAg9UQTcr0UJDWSqFOQlHp7Pl6/PtKggf5lNKiPpKbu3l+F70iQPk5N1VFdxdvhh2vL9ZYtbUSX6Bo31vOFQ4bo5yUlWgksLNTq3w8/wIYNVb8RLSvT5+y+WGJ3PtvXY+GSAjQqv1WhtGT3N2tnecIvKd3zvrR09zeprOI3rWyPzw9qE/uTVtEk90Kg4ryIVsDqfWxTWF6W+QWwofKOnHNPAE+AjtxrE3C/Ud3oN6o2zzQm+aSm6uymBG2pEzCp5bfGfgcSlWjmuX8MtBeRdiLSABgKTK60zWQgstb9IOC9/dXbjTHGxFa1I3fnXImI3ARMQd+3jHXOfSEiI9HLYCcDY4BxIrIcHbEPjWXQxhhj9i+quT7OubeBtys9dleFj7cDibIksjHGhJ61HzDGmBCy5G6MMSFkyd0YY0LIkrsxxoSQJXdjjAkh31r+ikgR8G0tn94C2GdrgwRjryV4wvI6wF5LUNXltbR1zqVXt5Fvyb0uRGReNL0VEoG9luAJy+sAey1BFY/XYmUZY4wJIUvuxhgTQoma3J/wOwAP2WsJnrC8DrDXElQxfy0JWXM3xhizf4k6cjfGGLMfCZvcReQvIvKZiCwUkXdF5Ai/Y6otEfmniCwtfz2viUjs19qLAREZLCJfiEiZiCTkrAYR6SMiy0RkuYiM8Due2hKRsSKyrnyVtIQlIq1F5H0RWVL+u3Wr3zHVloikichHIvJp+Wv5v5geL1HLMiLS1Dm3ufzjW4BjnHM3+BxWrYhIb7QHfomI/B3AOfcHn8OqMRHpCJQBjwO/c84FeJHcvUWzGHyiEJHTga3Ac8654/yOp7ZE5HDgcOfcJyLSBJgPXJSgPxMBGjvntopIfWAWcKtzLj8Wx0vYkXsksZdrTBXL+iUK59y7zrmS8k/z0dWuEo5zbolzbpnfcdTBScBy59w3zrliYDxwoc8x1YpzLo8qVkNLNM65Nc65T8o/3gIsAVr6G1XtOLW1/NP65beY5a2ETe4AIjJKRFYClwF3Vbd9grgKeMfvIJJUS2Blhc8LSdBEEkYikgF0Beb6G0ntiUiKiCwE1gFTnXMxey2BTu4iMk1EFlVxuxDAOfcn51xr4AXgJn+j3b/qXkv5Nn8CStDXE0jRvI4EFtVC7yb+RORA4FXgtkrv2hOKc67UOdcFfXd+kojErGQW1UpMfnHO9Ypy0xeBt4C7YxhOnVT3WkRkONAX6Bnk9Wdr8DNJRNEsBm/irLw+/SrwgnMu1+94vOCc2ygiHwB9gJic9A70yH1/RKR9hU/7AUv9iqWuRKQP8Aegn3Num9/xJLFoFoM3cVR+EnIMsMQ594Df8dSFiKRHZsKJyAFAL2KYtxJ5tsyrwNHo7IxvgRucc6v8jap2yhcWbwisL38oPxFn/ohIf+BfQDqwEVjonDvH36hqRkTOA0azezH4UT6HVCsi8hJwJtp98HvgbufcGF+DqgUR6QHMBD5H/9YB/li+rnNCEZFOwLPo71Y94BXn3MiYHS9Rk7sxxph9S9iyjDHGmH2z5G6MMSFkyd0YY0LIkrsxxoSQJXdjjAkhS+7GGBNCltyNMSaELLkbY0wI/X/nkFfdAZY+dwAAAABJRU5ErkJggg==\n", "text/plain": [ "

" ] @@ -150,7 +206,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "This isn't consistent with the solutions, which I don't understand. The key issue is the meaning of the statement that \"the variance of each observation is known to be 1\". How does one compute $p_{m}(\\{y_{i}\\})$?" + "This now matches the solution published on the web!" ] }, { @@ -177,7 +233,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.5" + "version": "3.6.4" } }, "nbformat": 4, diff --git a/Useful Formulae.ipynb b/Useful Formulae.ipynb index 5a3d891..1d0ddd3 100644 --- a/Useful Formulae.ipynb +++ b/Useful Formulae.ipynb @@ -52,66 +52,6 @@ " " ] }, - { - "cell_type": "code", - "execution_count": 26, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "0.7403867575800461" - ] - }, - "execution_count": 26, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "post_sample(-.25,1,.25,-.25,1,10)+post_sample(-.25,-1,.25,-.25,1,10)" - ] - }, - { - "cell_type": "code", - "execution_count": 29, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "0.05095226579074726" - ] - }, - "execution_count": 29, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - ".1*post_sample(-.25,-1,.25,-.25,1,10)/(post_sample(-.25,1,.25,-.25,1,10)+post_sample(-.25,-1,.25,-.25,1,10))" - ] - }, - { - "cell_type": "code", - "execution_count": 30, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "0.4414296078832747" - ] - }, - "execution_count": 30, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - ".9*post_sample(-.25,1,.25,-.25,1,10)/(post_sample(-.25,1,.25,-.25,1,10)+post_sample(-.25,-1,.25,-.25,1,10))" - ] - }, { "cell_type": "code", "execution_count": null, @@ -136,7 +76,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.5" + "version": "3.6.4" } }, "nbformat": 4, From 6a7223b767323034fb5cb9e1cc10898bbf3a7f5b Mon Sep 17 00:00:00 2001 From: Jeremy Teitelbaum Date: Mon, 15 Oct 2018 07:06:03 -0400 Subject: [PATCH 3/4] more notes on likelihood ratios --- LikelihoodRatioNotes.ipynb | 377 +++++++++++++++++++++++++++++++++++++ 1 file changed, 377 insertions(+) create mode 100644 LikelihoodRatioNotes.ipynb diff --git a/LikelihoodRatioNotes.ipynb b/LikelihoodRatioNotes.ipynb new file mode 100644 index 0000000..e189011 --- /dev/null +++ b/LikelihoodRatioNotes.ipynb @@ -0,0 +1,377 @@ +{ + "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": [ + "
" + ] + }, + "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": [ + "
" + ] + }, + "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 xs])/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}" + ] + }, + "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 xs])/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": [ + "
" + ] + }, + "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(T1)=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 +} From 62ca64acb75f6ef4c6632ec3050eab32e7f425fc Mon Sep 17 00:00:00 2001 From: Jeremy Teitelbaum Date: Mon, 15 Oct 2018 11:19:39 -0400 Subject: [PATCH 4/4] changes to 5.9.8 --- BDA 5.9.8.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/BDA 5.9.8.ipynb b/BDA 5.9.8.ipynb index 9f6a2e1..c3ad5cd 100644 --- a/BDA 5.9.8.ipynb +++ b/BDA 5.9.8.ipynb @@ -88,7 +88,7 @@ "$$\n", "\\theta\\sim N(\\mu_1,\\tau_1^2)\n", "$$\n", - "where $\\tau_1^2=1/(1/.1+1/.25)=1/14$ and $\\mu_1=(-.25/.1\\pm 1/.25)/14$ giving $6.5$ or $1.5$. \n", + "where $\\tau_1^2=1/(1/.1+1/.25)=1/14$ and $\\mu_1=(-.25/.1\\pm 1/.25)/14$ giving $-.46$ and $.11$. \n", "and $p_{m}(\\{y_{i}\\})=N(-.25,\\pm 1, .1+.25)$.\n", "\n", "\n"