Created
March 16, 2017 07:09
-
-
Save fperez/a968bacd657a504d4de1d127d65b2a80 to your computer and use it in GitHub Desktop.
sprt test - adding an nbsp in the gap0 paragraph in 2nd cell
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"cells": [ | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"<p class=\"title\">Wald's Sequential Probability Ratio Test</p>" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## Sequential test: draw data sequentially\n", | |
"\n", | |
"Sequence of data $X_1, X_2, \\ldots$. \n", | |
"\n", | |
"Two hypotheses, $H_0$ and $H_1$, each of which completely specifies the joint distribution of the data.\n", | |
"\n", | |
"Assume that the joint distributions under $H_0$ and $H_1$ are absolutely continuous with respect to each other, relative to some underlying measure.\n", | |
"\n", | |
"Let $f_{0m}$ be the likelihood of $H_0$ for data $(X_j)_{j=1}^m$ and let $f_{1m}$ be the likelihood of $H_1$ for data $(X_j)_{j=1}^m$.\n", | |
"\n", | |
"The likelihood ratio of $H_1$ to $H_0$ is $f_{1m}/f_{0m}$; this is (loosely speaking) the probability of observing $X_1, \\ldots, X_m$ if $H_1$ is true, divided by the probability of observing $X_1, \\ldots, X_m$ if $H_0$ is true.\n", | |
"\n", | |
"The probability of observing the data actually observed will tend to be higher for whichever\n", | |
"hypothesis is in fact true, so this likelihood ratio will tend to be greater than $1$ if $H_1$ is true,\n", | |
"and will tend to be less than $1$ if $H_0$ is true.\n", | |
"The more observations we make, the more probable it is that the\n", | |
"resulting likelihood ratio will be small if $H_0$ is true.\n", | |
"Wald (1945) showed that if $H_0$ is true, then the probability is at most $\\alpha$\n", | |
"that the likelihood ratio\n", | |
"is ever greater than $1/\\alpha$, no matter how many observations are made.\n", | |
"More generally, we have:\n", | |
"\n", | |
"<p class=\"gap01\"> </p>\n", | |
"\n", | |
"<div class=\"theorem\"> Wald's SPRT\n", | |
"\n", | |
"<p>For any $\\alpha \\in (0, 1)$ and $\\beta \\in [0, 1)$, the following sequential \n", | |
" algorithm tests the hypothesis $H_0$ at level no larger than \n", | |
" $\\alpha$ and with power at least $1-\\beta$\n", | |
" against the alternative $H_1$.\n", | |
"</p>\n", | |
"\n", | |
"<p>\n", | |
" Set $m=0$.\n", | |
"</p> \n", | |
"\n", | |
" <ol>\n", | |
" <li> Increment $m$\n", | |
" <li> If $\\frac{f_{1m}}{f_{0m}} \\ge \\frac{1-\\beta}{\\alpha}$, stop and reject $H_0$.</li>\n", | |
" <li> If $\\frac{f_{1m}}{f_{0m}} \\le \\frac{\\beta}{1-\\alpha}$, stop and do not reject $H_0$.</li>\n", | |
" <li> If $\\frac{\\beta}{1-\\alpha} < \\frac{f_{1m}}{f_{0m}} < \\frac{1-\\beta}{\\alpha}$, go to step 1. </li>\n", | |
" </ol>\n", | |
" \n", | |
"</div>" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## SPRT miracle\n", | |
"\n", | |
"Don't need to know the distribution of the likelihood ratio $\\mbox{LR}=\\frac{f_{1m}}{f_{0m}}$ under the null hypothesis\n", | |
"to find the critical values for the test." | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## Connection to Gambler's ruin\n", | |
"\n", | |
"For iid data, the likelihood ratio is a product of terms. On a log scale, it's a sum. Each \"trial\" produces another term in the sum, positive or negative. But—unlike the classical Gambler's Ruin problem in which the game is fair—the terms are not equal in magnitude: the steps are not all the same size." | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## Likelihood ratio for two values of $p$ in iid Bernoulli trials\n", | |
"\n", | |
"Suppose $X_1, X_2, \\ldots,$ are iid $\\mbox{Bernoulli}(p)$ random variables\n", | |
"and let $p_1 > p_0$.\n", | |
"\n", | |
"Set $\\mbox{LR} \\leftarrow 1$ and $j \\leftarrow 0$.\n", | |
"\n", | |
" + Increment $j$\n", | |
" + If $X_j = 1$, $\\mbox{LR} \\leftarrow \\mbox{LR} \\times p_1/p_0$. \n", | |
" + If $X_j = 0$, $\\mbox{LR} \\leftarrow \\mbox{LR} \\times (1-p_1)/(1- p_0)$.\n", | |
"\n", | |
"What's $\\mbox{LR}$ at stage $m$?\n", | |
"Let $T_m \\equiv \\sum_{j=1}^m X_j$.\n", | |
"$$\n", | |
" \\frac{p_{1m}}{p_{0m}} \\equiv \n", | |
" \\frac{p_1^{T_m}(1-p_1)^{m-T_m}}{p_0^{T_m}(1-p_0)^{m-T_m}}.\n", | |
"$$\n", | |
"This is the ratio of binomial probability when $p = p_1$ to binomial probability when\n", | |
"$p = p_0$ (the binomial coefficients in the numerator and denominator cancel).\n", | |
"\n", | |
"\n", | |
"\n", | |
"\n", | |
"<div class=\"callout\">\n", | |
"\n", | |
"<p class=\"subtitle\">Wald's SPRT for $p$ in iid Bernoulli trials</p>\n", | |
"\n", | |
"Conclude $p > p_0$ if \n", | |
"$$\n", | |
" \\frac{p_{1m}}{p_{0m}} \\ge \\frac{1-\\beta}{\\alpha}.\n", | |
"$$\n", | |
"Conclude $p \\le p_0$ if\n", | |
"$$\n", | |
" \\frac{p_{1m}}{p_{0m}} \\le \\frac{\\beta}{1-\\alpha}.\n", | |
"$$\n", | |
"Otherwise, draw again.\n", | |
"\n", | |
"</div>\n", | |
"\n", | |
"The SPRT approximately minimizes the \n", | |
"expected sample size\n", | |
"when $p \\le p_0$ or $p > p_1$.\n", | |
"For values in $(p_1, p_0)$, it can have larger sample sizes than fixed-sample-size tests." | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## Simulation of SPRT for Binomial\n", | |
"\n", | |
"Let's watch the SPRT in action for iid Bernoulli trials." | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 1, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"# This is the first cell with code: set up the Python environment\n", | |
"%matplotlib inline\n", | |
"from __future__ import division\n", | |
"import matplotlib.pyplot as plt\n", | |
"import math\n", | |
"import numpy as np\n", | |
"import scipy as sp\n", | |
"import scipy.stats\n", | |
"from scipy.stats import binom\n", | |
"import pandas as pd\n", | |
"# For interactive widgets\n", | |
"from ipywidgets import interact, interactive, fixed\n", | |
"import ipywidgets as widgets\n", | |
"from IPython.display import clear_output, display, HTML" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 2, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"<function __main__.plotBinomialSPRT>" | |
] | |
}, | |
"execution_count": 2, | |
"metadata": {}, | |
"output_type": "execute_result" | |
}, | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEZCAYAAACNebLAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8FdXd+PHPl7CIlM1ooBEI+ki1uCCgVB9U4oa2lYJp\n3eCxULW1ta0KWLBgS9BalQI/2wcVoUS0bkFEhUdQUUkRBYkCBiEGKrtsQhRCFiDJ9/fHubm5SeYm\nN8ldsnzfr9d9cWfmzJlzB5jvnHNmzhFVxRhjjKmsRawLYIwxpmGyAGGMMcaTBQhjjDGeLEAYY4zx\nZAHCGGOMJwsQxhhjPFmAMCERkeEi8laE8n5GRB6sx/55ItIzfCUK6ZgniMgiEflWRNIjfKwkESkV\nEc//ryIySUT+FckymObJAoTxE5FLRORD30XvgIh8ICL9AVT1RVW9tgGUcZmI3Ba4TlXbq+q2KBfl\nZ8ApQGdVvSlwg4h09V3QTwlYN9Fj3QQRWRzi8Wp6YUl9eSaJyNYQ8wws86kiMl9EvhaRb0QkS0R+\nHpBnqYgc9n22iMj4gH1LfUH6sIjsFJFpIiK+bZ8H7FcsIoUBae+vbTlNdLWMdQFMwyAi7YFFwJ3A\nK0Br4FLgaCzL1YAlAZvU401TVd0rIpuBy4BXfasvBbIrrbsMyIhA2ery9uu/gLVAd+AYcC7QtVKe\nHVVVReQi4D0RWauq7/i2naeqW0XkdGA5sBGYo6rnlGUgIsuA51T1mTr9KhN1VoMwZb4HqKrOU+eo\nqr6rqp8DiMhIEfmgLLHvrvE3IrJJRA6JyIMicnpADeRlEWnptW/A/qdXLoSIdPI13ewXkYO+74m+\nbX/BXWhn+O5A/1E5LxHpICLP+fbfKiITA/Ie6asV/U1EckXkSxEJWisSkbN8NZZvRGS9iAzxrU8F\n/gzc7CvHLzx2/wAXAPA1DfUD/g4MClh3Me5iioj8SETW+M7ldhGZVE25eopIhi/t28DJ1aQdLyK7\nfOXMFpHLgyS9EHhWVYtUtVRVP1PVtytnB6Cqq4ANwDkB68u2bQE+BM4PVqRgZTUNjwUIU2YTUCIi\nc0XkWhHp5JGm8p3pYKAvcBEwDngaGI67Cz0XuKWafYPd5bYA0nx59AAKgCcAVPUB3IX3d6raQVXv\n9shrBtAe6AkkAz+vdAEfgLuTjwf+BszxKoQvuC0C3sI1Jd0NvCAivVQ1Ffgr8LKvHF53xMvxBQjc\nOdoIvBewrh+uBr/at3wEuFVVOwI/Bn4tIj/xPEPwIpCJCwx/AUaWbVDV7apaFiy/B/wW6K+qHYBr\ngG1B8lwJPCkiN4lI9yBpypqNBgK9gTVVEoichQvim4PkYRoRCxAGAFXNAy4BSoFZwH4ReSOwzdzD\nY6qar6rZwOfAO74LVB6wBHdhDMbzTlJVc1X1NV8NJh94hPKLarV5+e7KbwLuV9UCVd0OTANuDUi7\nXVXTfE1DzwJdRSTBI8+LgHaq+piqFqvqMuD/qBj0qvNv4BwR6YC7YH6gql8CJ/vWXQKsUtVi3+9e\nrqobfN8/B17GV9uo8ENFegAXAH9W1eOq+gEukHkpwTUVniMiLVV1h6oG65+4ARfUHgC2+GozFwQe\nGvhaRA7i/n2MV9WMgO1rROQILhAuA56q9uyYRsEChPFT1RxVvU1Ve+CaDxKBx6vZZX/A90JgX6Xl\n79S2DCLSVkSeFpFtIvIt7kLbqazTswYn4+7KdwSs2w6cGrC8t+yLqhbiLnxe5UwEdlZaVzmvoHzB\n6StccLsMV/MB+Chg3fKy9CIyQETe9zWNfYvrC/JqOvou8I2v7IHl8irDl8C9QCqwT0ReFJHvBkl7\nSFUnqOq5QBfgM+C1wCRAvKrGq+rZqvpEpSz6qup3gBuBHwDtvI5jGhcLEMaTqm4C5lLezlwf+cCJ\nZQsi0rWatPcBvYALVbUT5bWHsgBRXQfsAeA4rgO5TBLuQl1bu3HNXIF61DKvsmami3CBAWCFb91A\nAgIErtnodeBU3+9+Gu9a1h6gs4i0rVQuT6r6sqpeSvk5ebSmQqtqLjAVSBSRzgGbqgvSZX0Q84FV\nQNA+FNN4WIAwAIjImSIyRkRO9S13xzWnrAxD9p8BZ4vIeSLSBnfxCHah/w6u9nFYRE7C3f0G2gdU\n6dwGUNVSYB7wsIh8R0SSgNG4J3Rq62OgQETGiUhLEUkGrgNeqkUeHwA/B3ar6hHfuhW+dR2peG6/\ng6sZHBeRAbi+nEBlF+AdwCfAZBFpJSKXAEO8Di4i3xORy0WkNe7JpEJcE6JX2kdF5GwRiRP3RNtd\nwH9U9ZvA44foUeCXQZruTCNiAcKUycM1DXwsInm4O94s3B29l1A7nVHVzcCDuE7aTZQ3t3h5HFfb\nOOArQ+X3BP4O3OB7wqms+Svw2HfjOra34O7Qn6/hsUrPcqvqcdyF90e+sszAdSLXpvP137gO7sDf\nuw44AfhEVYsC1t8FPCQih3D9AJVfvgss53BcreQg8CdcX4qXNriL9de4GtEpwB+DpD0R16T0DfAf\nXO0psJO8uppbhW2+PpR/A3+oLp1p+CQaEwaJSDfgOVzbZikwW1X/4au+puOqv9uAG1X1UMQLZIwx\npkbRChBdga6quk5EvgN8CgwFfgEcVNUp4t7M7Kyq9nalMcY0AFFpYlLVvaq6zvf9CO459G64IFFW\nPX4WGBaN8hhjjKlZVGoQFQ7oBlXLwD0ds1NVOwdsy1XVk6JaIGOMMZ6i2knta16aD9zjq0mE3NFp\njDEmuqI2WJ9v6IL5wL9U9Q3f6n0i0kVV9/n6KfYH2dcChzHG1IGq1nn8q2jWINKAjar694B1C4FR\nvu8jgTcq71RGVe2jyqRJk2JehobysXNh58LORfWf+opKDcI3uNcIYL2IrMU1JU0AHgPmiRvffzvu\nNX1jjDENQFQChKp+CMQF2XxVNMpgjDGmduxN6kYmOTk51kVoMOxclLNzUc7ORfhE/THXuhARrTIi\njzHGmOql1q+TutEEiMZQTmOau549e7J9u+fo4yaCkpKS2LZtW5X1ImIBwhjTMPguSLEuRrMT7LzX\nN0BYH4QxxhhPFiCMMcZ4sgBhjDHGkwUIY4wxnixAGGOM8WQBwhjT7GzdujXWRWgULEAYY5qVrVu3\n8vHHH0fteDt27CA9vfIU442DBQhjTLMyc+ZMbr75Zv/yZ599xn333Rex4/Xo0YOCggI2btwYNE1e\nXh45OTkRK0NdRW0+CGOMibWsrCy6d+/uX54+fTorVqygU6dOYT3O2rVr/c1YW7Zs4e677+bee+/l\nySef9Ew/b948fvjDH4a1DOFgNQhjTLOxaNEiLr/8cv/ymDFjGDp0aFiPsW7dOg4dOkRKSgopKSks\nWbKE1q1bc+zYMY4cOeK5z65du0hMTAxrOcLBAoQxptnIzMykd+/eET3Ghg0b/CPKfvrpp5xzzjkA\n9OnTh48++qhK+pycHM4666yIlqmurInJGBM1MrnOwwJVoJPqNt5TYWEhIuEpg5edO3eSlJTE+vXr\nmTt3Lps3b+bpp58GIDExkc2bNzN48OAK+7z++uuMHj3av7xnzx7S0tI4//zzWb58OXfddRfx8fHk\n5+fTpUuXiJXdiwUIY0zU1PXCHi4lJSW1Sj9lyhSKiooqrFNVRISRI0eSlJRUYduqVatISUkhLi6O\nadOm8dRTT/HMM88wYcIEOnXqxKZNmyqkLy0tpbi4mNatWwNQUFDAsGHDWLx4MfHx8SQkJDB27FhG\njBjBddddV4dfXD/RmnJ0DnAdsE9Vz/OtmwT8EtjvSzZBVd+KRnmMMc1Ty5a1u+SNGzeuVumLioqI\niyufPDM7O5tevXoBrvbSrl27CumXLl3K1Vdf7V9OT0+nf//+xMfHA5CQkEBWVhbDhw+nVatWtSpL\nOESrD+IZ4BqP9dNVtZ/vY8HBGBNRXbp0IT8/v8r6cA1RvmLFCv/3AwcOsHLlSkaNGgVAbm4uXbt2\nrZB+1apVDBgwwL987Ngxf0AByM/PJy4ujpSUlLCUr7aiEiBUdQXwjcemyDUGGmNMJYMGDarwktyM\nGTOYM2cOGRkZTJ48mby8vDrnvW7dOoYMGcLzzz/PggULeOKJJ1iwYAHt27cH3CO2AwcO9Kc/dOgQ\nJ510UoU8brnlFg4ePMiSJUtYuHAhu3fvpm/fvsydO5fCwsI6l62uojZhkIgkAYsqNTGNAg4BnwBj\nVfVQkH1twiBjGoGGPmHQN998w9SpU3n44YfDnvdLL73ELbfcEnT7HXfcwT//+U//8uzZsxkyZEiV\nWkVdRGrCoFh2Uj8JPKiqKiJ/AaYDtwdLnJqa6v+enJxsE5MbY2qtc+fOxMfHc/DgQX87f7gE9j1U\nlpmZWaGvAdzTSuEIDoEyMjLIyMgIW34xq0GEus233WoQxjQCDb0GAe7JodmzZ3PnnXdG5XglJSVM\nnTqV8ePH+9dt2bKFrKwshg0bFpZjRKoGEc0A0RMXBM71LXdV1b2+76OBC1V1eJB9LUAY0wg0hgAR\nbXv37qVjx460bds2Ysdo1AFCRF4EkoF4YB8wCbgcOB8oBbYBd6rqviD7W4AwphGwABEbjTpA1JcF\nCGMaBwsQsRGpAGFjMRljjPFkAcIYY4wnCxDGGGM8WYAwxhjjyQKEMcYYTxYgjDHGeLIAYYwxxpMF\nCGOMqcHWrVtjXYSYsABhjDHV2Lp1a4UhwqNlx44dpKenR/24gSxAGGNMgLy8PHJycvzLM2fO5Oab\nb66Q5rPPPuO+++6LaDl69OhBQUEBGzdujOhxqmMBwhjT7Kxdu5YFCxawYMECpk6dWmHbvHnzKkzy\n07179wrbp0+fzuTJk8nNzY14OYcPH86MGTMifpxgLEAYY5qVdevWcejQIVJSUkhJSWHx4sUVtu/a\ntYvExEQAFi1axOWXX15h+5gxYxg6dGhUytqmTRuOHTvGkSNHonK8yixAGGOiRyQ8n3rYsGGDf8Kx\nNWvWcO655/q35eTkcNZZZ/mXMzMz6d27d72OV199+vTho48+ismxYzmjnDGmuYnxSK87d+4kKSmJ\n9evXM3fuXDZv3szTTz/t3/76668zevRo/3JhYSFSz4BUnT179pCWlsb555/P8uXLueuuu4iPj+fI\nkSP+2eYSExPZvHkzgwcPjlg5grEAYYxpNlatWkVKSgpxcXFMmzaNp556irS0NCZOnEhpaSnFxcW0\nbt3an76kpKTWx5gyZQpFRUUV1qkqIsLIkSNJSkoCoKCggGHDhrF48WLi4+NJSEhg7NixjBgxguuu\nu86/b6dOndi0aVMdf3H9WIAwxjQbRUVFFeaOzs7OplevXgAsXbq0yrzRLVvW/hI5bty4kNKlp6fT\nv39//9zYCQkJZGVlMXz4cFq1auVPV1hYSLt27WpdjnCISh+EiMwRkX0ikhWwrrOIvCMiOSLytoh0\njEZZjDHN14oVK/zfDxw4wMqVKxk1ahTgahcDBgyokL5Lly7k5+d75lXfiZGOHTvmD04A+fn5xMXF\nkZKSUiFdbm6uv7kp2kIOECISJyInByy3FpFfiUh2CLs/A1xTad39wLuqeibwPvDHUMtijDG1tW7d\nOoYMGcLzzz/PggULeOKJJ1iwYAHt27fn0KFDnHTSSVX2GTRoEKtXr66wbsaMGcyZM4eMjAwmT55M\nXl5encpzyy23cPDgQZYsWcLChQvZvXs3ffv2Ze7cuRQWFvrTZWVlMXDgwDodo95UtcYPcDNwCNgN\n/BsYDOwCXgP6hZhHEpAVsPwF0MX3vSvwRTX7qjGm4WvI/1dffPHFoNtmzZqle/bsqbI+NzdXJ0yY\nEMli1ej222+vMU2w8+5bH9J13usTagPbA0B/Vf2PiPQDVgI/U9VF9YhNCaq6z3f13ysiCfXIyxhj\nqhXY91DZnj17PJtxOnfuTHx8PAcPHvT3FURTZmZmlX6RaAo1QBxT1f8AqOoaEdlcz+DgpdoGvdTU\nVP/35ORk/3PMxhgTihtvvNFz/ZYtWzjvvPOC7nfvvfcye/Zs7rzzzkgVzVNJSQnvv/8+48ePD3mf\njIwMMjIywlYG0RA6WkRkFzA9YNWYwGVVnV5lp6p5JAGLVPU833I2kKyq+0SkK7BMVb8fZF8NpZzG\nmNgSkXp33hpn7969dOzYkbZt29aYNth5962v84scoXZSzwbaB3wCl78TYh7i+5RZCIzyfR8JvBFi\nPsYY0+R17do1pOAQSSHVIKrNQOReVX28hjQvAslAPLAPmAS8DrwCdAe2Azeq6rdB9rcahDGNgNUg\nYiNSNYhwBIgdqtqjXpnUfAwLEMY0AhYgYiPWTUzVidxAJcYYY2ImHAHCbheMMaYJCukxVxHJwzsQ\nCBDbXhRjTIORlJQU0dFPjbeyAQDDrd59ENFgfRDGGFN7DaEPwhhjTBNkAcIYY4wnCxDGGGM8WYAw\nxhjjyQKEMcYYTxYgjDHGeLIAYYwxxpMFCGOMMZ4sQBhjjPFkAcIYY4wnCxDGGGM8WYAwxhjjKaTR\nXCNJRLYBh4BS4LiqDohtiYwxxkADCBC4wJCsqt/EuiDGGGPKNYQmJqFhlMMYY0yAhnBhVmCpiGSK\nyC9jXRhjjDFOQ2hiGqiqe0TkFFygyFbVFbEulDHGNHcxDxCqusf359ci8howAKgSIFJTU/3fk5OT\nSU5OjlIJjTGmccjIyCAjIyNs+cV0ylERORFooapHRKQd8A4wWVXfqZTOphw1xphaqu+Uo7GuQXQB\nXhMR9ZXlhcrBwRhjTGzEtAYRKqtBGGNM7dW3BtEQnmIyxhjTAFmAMMYY48kChDHGGE8WIIwxxniy\nAGGMMcaTBQhjjDGeLEAYY4zxZAHCGGOMJwsQxhhjPFmAMMYY48kChDHGGE8WIIwxxniyABHE4cNw\n8GCsS2GMMbFjASKAKixfDiNHQo8ecMYZcOedsHmz275nD0ydChdeCLfdBuvXl+9bXAxvvQUTJ8Kn\nn8am/MYYE07NIkAUFwffVloKq1bBuHEuIPzmN9CnjwsKOTnQpQv893/DBRdA796wcSP89a8u7TXX\nwNVXw+9+B6eeCpMmwdGjkJLi9nn5ZbdsjDGNUaOZD2LfPiUhoXzdN9/AmjXuov3d71ZMX1oKWVnw\n9tvurv6jj+D734cbbnCfTp3gvfdg6VKXpmNH+OlP3YX9/PNBKo2efuQIfPghXHIJtGtXvv7oUUhP\nh5074cYboVcvt764GBYtghkzYO1aGDoUbrrJBZX33oN33oFPPoFLL4Wf/cwFmrZtI3PuSktd+Tt0\niEz+xpiGq77zQTSaANGxo9K7N5x7Lnz8MXz5pfuenQ2nnw4/+pG7yK5YAStXQkICDB7sLr6XXeYu\n1K+8Aq++CgUFkJzs7v6vvhrOPDNyZf/qK3fc9HTYsQOuvNKVq39/WLYM5s93gW7QILftqqtcMKsc\npMqUlrrazaefwtdfu5pK377Q0jc34DffwLp1rlZUdi4KC+Hss925uOYaGDAATjghcr/ZGNMwNPoA\nISLXAo/jmrvmqOpjHmm0qEjJyIAvvoCLLoJ+/aBVKzh+3NUQFi92d+4DB7pPly7exystdZ+WsZ5s\nNcDXX8O777raxXvvQX6+u4gPGOCC4K5drmlrwwb49NMMEhKSueACiI93v33bNlfz2bXL5dWnj+sn\nueQSdy5OOskFirffdrWX7Gw47zy4+GLXdHbeeS5ItmoV6zNROxkZGSQnJ8e6GA2CnYtydi7KNeoA\nISItgE3AlcBuIBO4WVW/qJSuWU05unMnZGbC6tXw+efQrZurAZx9Nrz1VipTpqRWSH/woKuFlHWs\nx8VVn39+vst/5Uq3X1aWq92ccQZ873vuz169ICkJund3n8CmtcpUXc1l3z4XfE891TXbBasFhUtq\naiqpqamRPUgjYeeinJ2LcvUNELG+jx4AbFbV7QAi8jIwFPii2r2auLKLckpK1W3Ll1ddFx/vmspC\n1a6da2ILvMkqLHS1s82b3efDD+Gll1yw2rnTNUl16eI+nTrBt99Cbq4LTgcOuDwTElxQ2L0bSkog\nMdGtS0iAU05x5ezc2dVoOnVy/SIdOkD79m7/du3gxBNdU2GLZvH4hDENW6wDxKnAzoDlXbigYaKs\nbVvXl9G3b9VtgTWEvXtdcOjUyV3w4+Ph5JOhTZuK++TluUCxf79r9tq/3wWUr7+GTZtcfnl57n2T\nw4ddrabsc/QotG7tynTCCS7vsj9btXLbWrVygevjj933li3dJy7OfUpLXZAK/BQXV/20bOn2L/vE\nxVXNKy7OnQOvPMuaKyunj4tzwbJsn2DlKft+/Hj5csuW5b+xtLS8rKrlvzXwN7ds6WqBGze6/VUr\nlinwmC1alP9GkYq/o6zcZfu1aFGePvDPFi3cvi1auGOVfaB8vUj127zSlf1bKy2tuE/l/bw+gXlk\nZJT/OwxMU7YcuL7yOq9/+8GWA8sc+D1wndc2r7IFlscrXeAxQ8kv2O+prVg3Mf0UuEZVf+Vb/h9g\ngKreXSld82lfMsaYMGrMTUxfAT0Clrv51lVQnx9ojDGmbmLd0psJnCEiSSLSGrgZWBjjMhljjCHG\nNQhVLRGR3wHvUP6Ya3Ysy2SMMcaJ+XsQxhhjGqZYNzEZY4xpoCxAGGOM8WQBwhhjjCcLEMYYYzxZ\ngDDGGOPJAoQxxhhPFiCMMcZ4sgBhjDHGkwUIY4wxnixAGGOM8WQBwhhjjKeYDtYnIm2A5UBrX1nm\nq+rkWJbJGGOME/PB+kTkRFUtEJE44EPgblVdHdNCGWOMiX0Tk6oW+L62wdUibHhZY4xpAGIeIESk\nhYisBfYCS1U1M9ZlMsYYE/spR1HVUqCviHQAXheR3qq6MTCNzUltjDF105jnpPZT1cMisgy4FthY\nJUFqtEvUQC0DLo91IRoIOxfl7FyUs3NRLrV+u8f6KaaTgeOqekhE2gJXA496pdVJVokASNVUUiel\nxroYDYKdi3J2LsrZuSgnqXWuPACxr0F8F3hWRFrg+kPSVXVxjMtkjDGGGAcIVV0P9ItlGRqb5OTk\nWBehwbBzUc7ORTk7F+ET8/cgQiEi2hjKaYwxDYmINI1OamNM5PXs2ZPt27fHuhgmzJKSkti2bVvY\n841pDUJEugHPAV2AUmC2qv7DI53VIIwJA98dZayLYcIs2N9rfWsQsQ4QXYGuqrpORL4DfAoMVdUv\nKqWzAGFMGFiAaJoiFSBi+ia1qu5V1XW+70eAbODUWJbJGGOME/OhNsqISE/gfODj2JbEGGMMNJAA\n4Wtemg/c46tJGGOMibGYP8UkIi1xweFfqvpGsHSpqan+78nJyfasszHGVJKRkUFGRkbY8ov5exAi\n8hxwQFXHVJPGOqmNCQPrpG6ammQntYgMBEYAV4jIWhFZIyLXxrJMxpjGbevWrdUum9DF+immD1U1\nTlXPV9W+qtpPVd+KZZmMMY3X1q1b+fjjj4MuR8OOHTtIT09vcHnVRYPopDbGmHCYOXMmN998c9Bl\ngM8++4z77rsvYmXo0aMHBQUFbNxYddaCMnl5eeTk5IQlr0iyAGGMaTDmz59PQkICR48erfW+WVlZ\ndO/ePegywPTp05k8eTK5ubn1LmuZtWvXsmDBAhYsWMDUqVMBGD58ODNmzAi6z7x582jfvn1I+deU\nVyTFPECIyBwR2SciWbEuizEmtgYOHEjv3r1p06ZNrfddtGgRl19+edBlgDFjxjB06NB6l7PMunXr\nOHToECkpKaSkpLB4sZutoE2bNhw7dowjR7yf2t+1axeJiYkhHaOmvCIp5o+5As8A/4sbk8kYE0My\nuX4TzJSp6wRf7777LldeeWWd9s3MzGTChAlBlyNhw4YNjBgxAoA1a9Zw7rnn+rf16dOHjz76iMGD\nB1fYJycnh7POOqtWxwmWV6TFPECo6goRSYp1OYwxsZ+58b333uPOO+8EYPfu3cyZM4d+/fqRmZnJ\nrbfeyn/9139RWlrKI488wve//3327dvHxx9/zNy5cykoKECkPMAVFhZWWA63nTt3kpSUxPr165k7\ndy6bN2/m6aef9m9PTExk8+bNVS7qr7/+OqNHj/Yv79mzh7S0NM4//3yWL1/OXXfdRXx8PEeOHKFr\n167V5hVpMQ8QxhhTZvXq1aSlpVFQUMCwYcNYsmQJ8fHxtGjRgmnTpvHkk0/ywAMPcOaZZ5KSksIL\nL7zAaaedBkBJSUmFvCovh2LKlCkUFRVVWKeqiAgjR44kKan8XnbVqlWkpKQQFxfHtGnTeOqpp0hL\nS2PixIkAdOrUiU2bNlXIq7S0lOLiYlq3bg3g/52LFy8mPj6ehIQExo4dy4gRI7juuuv8+3nlFQ0W\nIIwxDUJOTg69evWiRYsWpKenc8EFFxAfHw9AdnY2bdu2paSkhJkzZ7Jnzx7AvTl8xx13ANCqVasK\n+bVsWfvL27hx40JOW1RURFxcnH85OzubXr16+ZcLCwtp165dhX2WLl3K1Vdf7V9OT0+nf//+/t+Z\nkJBAVlYWw4cPr/B7vPKKhkYTIGyoDWOatmXLlnHFFVewaNEijh8/7r/YFhYW8uqrr/LKK6+Qn59P\nt27d/B23n3zyCbNmzQKgS5cu5Ofn+y+klZcDheNt8hUrVnDrrbcCcODAAVauXMnDDz/s356bm+tv\nIiqzatUqJk2a5F8+duxYhaCSn59PXFwcKSkpFfbzystLkxtqA/wjuS5S1XODbLehNowJg4Y81Mbb\nb7/N8uXLGTRoEBdffDGPPfYYF198MevWreP666+nd+/eAPzpT3+iT58+bNiwgS+++IKXXnoJgLS0\nNHr27MkVV1zhXz7ttNMqPMk0Y8YM5s2bx86dOxk1ahRjxowJ+XHTQOvWrWPXrl18++23nHjiiaxf\nv57bbrutwmO19913H6NHj+bUU90MBocOHeK5557j97//vT/N4cOHmTJlCgMHDuT48eOceOKJpKWl\nMXjwYG666Sbatm3rmVdlkRpqA1WN6Qd4EdgNHAV2AL/wSKPGmPpr7P+X9u7dq0VFRaqq+uijj2p6\nerp/W25urk6YMCHocji9+OKLNaa5/fbbKyzPmjVL9+zZU6fjVc6rsmB/r771db4+x7yJSVWHx7oM\nxpjGYeIfrDpQAAAaTElEQVTEifTr148OHToAcOONN/q3de7cmfj4eA4cOMDJJ5/sXz548KC/jT9c\nAvsevGRmZlboawD3tFIozUSh5BUtDaKJqSbWxGRMeDTkJqZwKC0tZfbs2f5HZSsvR0NJSQlTp05l\n/Pjx/nVbtmwhKyuLYcOG1TsvL01yTupQxSJAlD0hV8ONgjGNSlMPEA3B3r176dixo7//IBp5Ncnh\nvhuyr7+G3r0hM7N2++3fH5nyGGMah65du4YlOIQ7r7qIeYAQkWtF5AsR2SQi1dejoqhrV3jkEfjx\nj+Ef/wBVWLs2ePrt2+EnP4HERBg9GkpLo1dWY4yJhFhPGNQCmAFcA5wN3CIitRukJIJSUmDVKnju\nOejZE371Kzh+vGKaY8dgwgTo3x9+8APYtQt69YIWMQ+9xhhTP7G+jA0ANqvqdlU9DrwMhG+oxRCs\nW+cu/MHeyj/9dPjwQ5g9G1auhEova/LBB/Cf/8Dq1TBxoqt53HVXzcfdtav+ZTfGmEiKdYA4FdgZ\nsLzLty4svv46+La8PNcUdM01MGAAVDemV5s2MHgweL25f+WVMG+eCySh2L4dhg6F006DP/whtH2M\nMSYWYv4eRKhuuSWVM89032saauPgQbj/fnd3//nnFS/sqvDKKzBmjLvof/45nHJKZMsOkJ0N8+fD\n3/8O994LaWmwc2fN+9VXbi5s3AiXXBL5YxljYivcQ23E+i3qi4C3ApbvB8Z7pNMePVTHjlXdvl31\nb39TLS2t+tZgaalqWppqly6qv/2t6rffVk2zcKFqr16qy5d7vngYMc8+q3rzzapffhnd4w4ZopqY\nqHrbbar5+dE9tml4aORvUhtvwf5eqeeb1CG9ByEiFwCXAolAIfA5sFRVv6lPcBKROCAHuBLYA6wG\nblHV7Erp9MABZeRIWLrU3f0/9FDFmsGBA+6Jo9JSmDnTdRo3ZbNmwUcfwZQpkJAQPN2xY3D0KPz6\n15CV5WpPtZyrxDQh9h5E0xST9yBE5Bcisgb4I9AWdzHfD1wCvCsiz4pIj7oeXFVLgN8B7wAbgJcr\nB4cy8fGwcCHs2+ceP63cH5CTA3fc4Z46aszB4Z//hC++qDmdCJx0EvTrB//+d/B0rVtD+/bw/PNw\nzz1w6aWus92Ypmrr1q0NIo8mobrqBfBboG01288HrqxPFSaUD82oWjx7turJJ6s+/3xo6d96S7Vr\nV9Xf/151//6a03/+uWpBQf3KaBqvpv5/acuWLfrSSy/FPI+62L59u7788ss1rvMS7O+VejYxVVuD\nUNUnVLXQa5uItFPVdar6XhjjVbN3xx3w3nvw4IPu8dtCz7Nf7ppr4JNPXDPS7t0153/22RDDFzON\niaiZM2dy880312qfvLw8cnJyqs3js88+47777gtLGYPp0aMHBQUFbNy4sdp10VTjY64icqqIXCAi\nrX3LCSLyV2BzxEvXTJ13nrvoHzni3sxevrz69KeeCk8/DX36RKd8xkTK/PnzSUhI4OjRo7XeNysr\nq8J8DGXWrl3LggULWLBgAVOnTq2yfd68ef45IbzymD59OpMnTyY3N7fWZaqt4cOHM2PGjBrXRUtN\nfRD3AuuA/wVWicgdQDauP6JeLf0i8jMR+VxESkSkX33yaorat4cXXoB334X//u/IHuvwYXjggZpr\nK8ZE2sCBA+nduzdt2rSp9b6LFi2qMDkQuIl9Dh06REpKCikpKSxevLjKfrt27SIxMTFoHmPGjGHo\n0Oi8v1s2U96RI0eqXRctNdUgfgWcqaoXA8Nww2IMVtXRqrqnnsdeD1wPVNPF2ryJuA73OkytW+vj\n/Oc/cNFFEIN50Y3xe/fdd7nyyivrtG9mZqZ/1rkyGzZs8L8ztWbNGs49t+KklTk5OZwV8FifVx7R\n1qdPHz766KMa10VDTQGiSFVzAVR1B5Cjqp+G48CqmqOqm4G6T4dnwqJ9e3jpJTdEyMCB7rtppkTC\n86mj9957j6uuugqA3bt389BDD/Hmm2+SmprKl19+Cbg5Hh5++GEWLFjAU089xahRowAoKChAAo69\nc+dOkpKSWL9+PWPHjiU1NZX777+/wvFef/11rr/+ev9yYWFhhTzCbc+ePTz88MO8+eabjB8/nu3b\nt3PkyBH27t3rT5OYmMjmzRVb8L3WRUNN96bdROQfAcvfDVxW1bsjUywTbSJw551uwMEbbnCPzj75\npA062OzE+B2J1atXk5aWRkFBAcOGDWPJkiXEx8fTokULpk2bxpNPPskDDzzAmWeeSUpKCi+88AKn\nnXYa4CbXCbRq1SpSUlKIi4tj2rRpPPXUU6SlpTFx4kTABZri4mJat27t36dyHqGYMmUKRUVFFdap\nKiLCyJEjSUpKAvD/psWLFxMfH09CQgJjx45lxIgRXHfddf59O3XqxKZKVXmvddFQU4CoPFpQYO2h\nxn9JIrIU6BK4yrffRFVdFFIJfVJTU/3faxpqw9Td+efDp5+6YUEsOJhoysnJoVevXrRo0YL09HQu\nuOAC/1Sh2dnZtG3blpKSEmbOnMmePa6FOyMjgzvuuAOAVpVG0iwqKqowNWh2dja9evXyLy9durTK\nVJ4t69CeO27cuJDSpaen079/f/9vSkhIICsri+HDh1coe2FhIe3atauwr9c6L+EeaqPas6Gqzwbb\nJiJVHweoun/YJlINDBAmsjp0gNtuq/v+R4+6UW8twJjaWLZsGVdccQWLFi3i+PHj/ot5YWEhr776\nKq+88gr5+fl069bN33H7ySefMGvWLAC6dOlCfn6+/0K6YsUKbr31VgAOHDjAypUrefjhh/3HW7Vq\nFZMmTapQhsp5BNJ61q6OHTtWIUDl5+cTFxdHSkpKhXS5ublV5q72Wuel8s3z5MmT61Xm+vwXvrHm\nJCGzfohGyuvJp+nT4brr3PAnxoTqtNNOY//+/bRp04ZbbrmFgwcP8uabbzJ9+nRmz55NYmIiHTp0\nYOjQocyfP59HHnmEs846y99nMGjQIFavXg24p5eGDBnC888/z4IFC3jiiSdYsGCB/3HWQ4cOcdJJ\nJ1UpQ2AeZWbMmMGcOXPIyMhg8uTJ5OXl1en3lf2mJUuWsHDhQnbv3k3fvn2ZO3cuhQH/kbKyshg4\ncGCFfb3WRUVd37ADdtbnDT3cU1E7cWM77QGWVJO2xjcJTeRlZalu2uS+79ihev31bhDAyo4dUx03\nTrV7d9UVK6JbRlO9xv5/ae/evVpUVKSqqo8++qimp6f7t+Xm5uqECRNUVfXFF1+sNp9Zs2bpnj17\nqqwPzCNWbr/99pDWBQr290o936SutolJRKqGWN8m6nnXr6qvA6/XJw8TXevXu/Gchg9372j8/vcw\n3mOS2Fat4LHH3LhPKSkwdizcd581OZn6mzhxIv369aNDhw4A3HhjeUNG586diY+P58CBAxX6Hrzs\n2bPHs8mmLI+DBw/6+wqiKTMzs0q/iNe6qKkuegBbgS2+Pyt/ttQnMtXmQyO/62lKPv1U9Te/Ka9J\n1GT7dtWLL1adPj2y5TKhaer/l0pKSnTmzJnVpvnyyy/1tddeq1cekVBcXKyPPvpojeu8BPt7JRrD\nfceaiGhjKKfxdvy4+5x4YqxLYmy474Zr7969dOzYkbYBg6V5rfMSqeG+qw0QItJTVbdVs12AU1W1\n1jMsi8gUYAhwFPgS+IWqHg6S1gKEMWFgAaJpisl8EMDfRORVEfm5iJztG6ivh4hcISIPAR8C36/j\nsd8BzlbV83ED//2xjvkYY4yJgJqG+74B+BNwJvAE8AGwEPglbvKgK1R1aV0OrKrvqmqpb3EV0K0u\n+ZjGqbjYzQBokxcZ03A1iD4IEVmIm03uxSDbrYmpCXrjDfjlL2HcOPekUwSHwDE+1sTUNMWkDyLg\nICkeqw8B61V1fzX71TjUhohMBPqp6k+ryccCRBO1bRvcdBN06QJz57ppVE3kWIBomiIVIEIdeOR2\n4GJgmW85GTcu02ki8qCq/strJ61hqA0RGQX8CLiipgLYWExNU8+e8MEHcP/9bn7tzEw45ZRYl8qY\nxincYzGFWoN4G/i5qu7zLXcBngNuAZar6jm1PrDItcA04DJVPVhDWqtBNAOffOLmv7Cmpsjp2bMn\n27dvj3UxTJglJSWxbdu2Kuuj1cS0UVV7BywLsEFVe4vIWlXtW+sDi2wGWgNlwWGVqt4VJK0FCGOM\nqaVoNTFliMj/Aa/4ln/mW9cO+LYuB1bVXjWnMsYYEyuh1iAESAEu8a36EHg1Wrf1VoNovjZuhKVL\n4e67renJmNqK9ItygG8wD1gBvA+8h+t3sCu2ibgTT3QDA15/PXzzTaxLY0zzElKAEJEbgdW4pqUb\ngY9F5GeRLJgx4J5yWrHC/dmvH3z8caxLZEzzEWoT02fA1WXvPIjIKcC7qtqnzgcWeRAYCpQC+4BR\nqro3SFqrsBhee83Nm52aCnd5Ps4QuuJiePttSE6GEGZyjJpDh6Bjx1iXwjQVUWliAlpUeiHuYC32\nDWaKqvbxPQH1JjCpph1M83b99bBqFZx6av3yWbUKLrzQzW0xYIDr52gISkvhootg8mQoKYl1aYwJ\n/SL/loi8LSKjfC+3vQksrs+BVfVIwGI7XE3CmGqdfjoMHVr3/R94wE1idN99sHmz+/OLL8JXvvpo\n0QLefx8yMmDwYNi7F7KzwSrPJlZCHotJRH4KlE2K+oGqvlbvg4v8Bfg57lHZy4O9MCci9n/EGGNq\nSSDyL8rVOfMQxmLypRsPtFXV1CD56KRJ5S1QNtSGqWzZMujbFzp1it4xS0vdnf7EiXD55aHvd+CA\na0b685+rH1Zk/XpITITKM1/u2wc/+IGrCT36KLRuHdpxd+1y/Rvt24deVtO4VB5qY/LkyfUKEDVN\n9ZkHHPb45AGH6zOVXaXjdMcN/GdTjpo6GTdO9fTTVTMz3XJ+vuqxY3XPLycntHRvv63atavqgw+q\nFhdXn7akRPWf/1RNSFC95x7Vw4frXr6DB1WHDFH9wQ9Ut21Tzc0NnjY/X3XCBNWOHVXPPlv1+PGK\n2z/7zOX1/vt1L49pmKjnlKM1zQfRXlU7eHzaq2qHOkclQETOCFgcBmTXJz/TvD32mPv86Edu+PCz\nz4YlS+qWV3Ex3HAD3H47FBRUn3bwYPj0U3j3Xbj2Wnd372X9erj0Upg9G956Cx5/vH538ied5IZL\nv/FG9/jvhRe6aV0D7d4NDz0E55wDX37pOuMXL4aWlcZPOOUUuOQSGD7cpbcOcuNXn+hSnw8wH8gC\n1gFvAN+tJm14wqlp8jZvVv35z1WXLq1fPnl5qiNGqJ5zjrvDfvxxdycezPHjqhMnqiYmqu7dW3Hb\nV1+pdumiOnOmq0WE26ZNqgcOVF2fna16ww2q77wTWj5ffaV62WWqV1+tum9feMtoYoN61iAaxIRB\nNbH3IEwsqEJaGvzud3DVVTBnDiQkVL/P+vVw7rlV1+fnN6z3LYIpLnbvmWzaBPPmxbo0pr6iMppr\nrFmAMLFUVARt2jSvsaCOH4dWrWJdClNfFiCMMcZ4itab1BEjImNFpFREbLJJYxqwUnuVtdmJaYAQ\nkW7A1YBNcWVMA/eLX8Bf/2qBojmJdQ3i/wF/iHEZjDEhePhh9+jwD38IX38d+n5bt8LOnZErl4mc\nmAUIEfkJsFNV18eqDMaY0HXr5t5Y79fPvbW+fLl3urLgUVgIf/oT9OkDR454pzUNW6hTjtZJNUNt\nPABMwDUvBW4LKjU11f/dhtowJjZatoRHHoHLLnODJr79thsRt0xJCVxxhRuV9v33oX9/N+BgfUbg\nLSmBuLj6l705qDzURn3F5CkmETkHeBcowAWGbsBXwACtOKx4WXp7ismYBiY31z3627lzxfUHDsCD\nD7q32q+9NvT89u+v+p7J11+7t7znzHF/mtppEo+5ishWoJ+qek4qaQHCmKZF1TVRDRrkmqIeeQSe\nego+/xy6dKmYdvFiuO02uPdeN4xKi1j3nDYijf4xVx+lhiYmY0zT8fXXcMcdMGKEGysqOxvWrasa\nHMDVRDIzYdEiuO46V0NpjPLy3Mi/06c3nifBGkSAUNXTVTU31uUwxkRHQoIb5LBbN3jiCXjller7\nKbp3dxMpnXeeG/SwuDhqRQ2LN990A0ju3Amvvgo//nF4At3WrW7iq0hpEE1MNbEmJmNMmb17oWvX\nqus/+ADGj3fNVg89VHXU2lhauNDNxTFokBvG5IEHoEMHV6Ooi6IiNxfIjBmuH+gvf4Ff/arqcDBN\nog+iJhYgjDHVmT8fRo9272q8+KIbHPHll+s/f3kkqXqP73XkiHsh8YQTYMIE70B31VUu4Dz+uBuS\n/oYb3IMBw4ZVTNdoA4SITAJ+CZQ9tTRBVd8KktYChDEmqIIC9zhs+/auff+xx1yfRnp69MsS7MIf\nitdeg3vugeRk91RXfj689JJrigt04ACcfHL5cmGhG1Cycgd+Yw8Qeao6PYS0FiCMMbUS7fcn9u2D\nsWPd47i//nXd8rj/fvdocHKyC3SPPgpXX+0mhKqLxh4gjqjqtBDSWoAwxjRIJSUwaxZMmuTGq/rz\nnxvO3B/1DRCx7sb5nYjcCnwCjFXVQzEujzGmCQt3rWLNGldbaNPGvTl+zjnhy7shiGgNopqhNiYC\nq4ADqqoi8hfclKO3B8nHahDGmHq79VZISnKz5oXjKacbbnCDF44a1TBf4Gu0TUwVCiGSBCxS1fOC\nbNdJkyb5l20sJmNMXezbB//zP3DsmOv8TUyMdYnCq/JYTJMnT26cAUJEuqrqXt/30cCFqjo8SFqr\nQRhjwqKkxD1G+uST8NxzrhO4qWq0NQgReQ44HygFtgF3quq+IGktQBhjwmrZMlebeO45uPLK4OmK\niuBvf3NDg3z3u9ErXzg02k5qVf15rI5tjDGXX+7elag8Gm2gd96B3/7WDfHRHO9RG0QfRE2sBmGM\niabdu92b2ZmZ8L//68ZOaoyaymiuxhjTIBw+7CY6OuMMN/x4Yw0O4RDTGoSI/B64CygG3lTV+4Ok\nsxqEMSZqDh6E+PhYl6L+Gm0NQkSSgSHAuap6LjA1VmVpTMI5nWBjZ+einJ2LcuE4F00hOIRDLJuY\nfgM8qqrFAKraSKcBiS67EJSzc1HOzkU5OxfhE8sA8T3gMhFZJSLLROSCGJbFGGNMJRF9zLWaoTYe\n8B27s6peJCIXAvOA0yNZHmOMMaGL5Ytyi4HHVPXfvuX/AD9Q1YMeaa2H2hhj6qBRvigHvA5cAfxb\nRL4HtPIKDlC/H2iMMaZuYhkgngHSRGQ9cBSwN6uNMaYBaRRvUhtjjIm+Bv0mtYhcKyJfiMgmERkf\n6/JEk4h0E5H3RWSDiKwXkbt96zuLyDsikiMib4tIx1iXNVpEpIWIrBGRhb7lZnkuRKSjiLwiItm+\nfx8/aMbnYrSIfC4iWSLygoi0bi7nQkTmiMg+EckKWBf0t4vIH0Vks+/fzeBQjtFgA4SItABmANcA\nZwO3iMhZsS1VVBUDY1T1bOBi4Le+338/8K6qngm8D/wxhmWMtnuAjQHLzfVc/B1YrKrfB/oAX9AM\nz4WIJAK/B/r55pJpCdxC8zkXz+Cuj4E8f7uI9AZuBL4P/BB4UkRq7NttsAECGABsVtXtqnoceBkY\nGuMyRY2q7lXVdb7vR4BsoBvuHDzrS/YsMCw2JYwuEekG/Aj4Z8DqZncuRKQDcKmqPgOgqsW+qXqb\n3bnwiQPaiUhLoC3wFc3kXKjqCuCbSquD/fafAC/7/r1sAzbjrrHVasgB4lRgZ8DyLt+6ZkdEeuLm\nzlgFdCmbN8M34VJC7EoWVf8P+APuPZoyzfFcnAYcEJFnfM1ts0TkRJrhuVDV3cA0YAcuMBxS1Xdp\nhuciQEKQ3175evoVIVxPG3KAMICIfAeYD9zjq0lUfqqgyT9lICI/Bvb5alTVVYub/LnANaP0A55Q\n1X5APq5ZoTn+u+iEu2NOAhJxNYkRNMNzUY16/faGHCC+AnoELHfzrWs2fNXm+cC/VPUN3+p9ItLF\nt70rsD9W5YuigcBPRGQL8BJwhYj8C9jbDM/FLmCnqn7iW34VFzCa47+Lq4AtqpqrqiXAa8B/0zzP\nRZlgv/0roHtAupCupw05QGQCZ4hIkoi0Bm4GFsa4TNGWBmxU1b8HrFsIjPJ9Hwm8UXmnpkZVJ6hq\nD1U9Hffv4H1VvRVYRPM7F/uAnb6XSwGuBDbQDP9d4JqWLhKRE3wdrlfiHmJoTudCqFirDvbbFwI3\n+57yOg04A1hdY+YN+T0IEbkW98RGC2COqj4a4yJFjYgMBJYD63HVRAUm4P5S5+HuBrYDN6rqt7Eq\nZ7SJyCBgrKr+REROohmeCxHpg+usbwVsAX6B66xtjudiEu6m4TiwFrgDaE8zOBci8iKQDMQD+4BJ\nuBEqXsHjt4vIH4HbcefqHlV9p8ZjNOQAYYwxJnYachOTMcaYGLIAYYwxxpMFCGOMMZ4sQBhjjPFk\nAcIYY4wnCxDGGGM8WYAwJoBvKO3fVLN9RQh55IW3VMbEhgUIYyrqDNxVeaWIxAGo6iUh5GEvF5km\nIZZTjhrTED0CnC4ia3BzchThhlQ+EzhLRPJUtb2ItMMNY9AJ90bzn1S1wlAwvrFw0nFv9rYEfqOq\nH0bvpxhTP/YmtTEBRCQJWKSq5/mG9fg/4GxV3eHbflhVO/hqFG1V9YiIxAOrVLVXpTRjgDaq+ohv\nrKATVTU/Rj/NmFqzGoQx1VtdFhwqEeAREbkMKAUSRSRBVQNHDs0E5ohIK+ANVf0sCuU1JmysD8KY\n6gW74x8BnAz0VdW+uGGVTwhMoKofAJfhhlWeKyL/E8mCGhNuFiCMqSgP12cA3pMTla3rCOxX1VIR\nuRw3aU2FNCLSw5dmDm701X6RKbIxkWFNTMYEUNVcEflQRLKAQtwwyhWS+P58AVgkIp8Bn+DmDK+c\nJhn4g4gcxwWen0es4MZEgHVSG2OM8WRNTMYYYzxZgDDGGOPJAoQxxhhPFiCMMcZ4sgBhjDHGkwUI\nY4wxnixAGGOM8WQBwhhjjKf/D1NuPIndT0/NAAAAAElFTkSuQmCC\n", | |
"text/plain": [ | |
"<matplotlib.figure.Figure at 0x10d8636d0>" | |
] | |
}, | |
"metadata": {}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"def plotBinomialSPRT(n, p, p0, p1, alpha, beta):\n", | |
" '''\n", | |
" Plots the progress of the SPRT for n iid Bernoulli trials with probabiity p\n", | |
" of success, for testing the hypothesis that p=p0 against the hypothesis p=p1\n", | |
" with significance level alpha and power 1-beta\n", | |
" '''\n", | |
" fig, ax = plt.subplots(nrows=2, ncols=1, sharex=True)\n", | |
" trials = sp.stats.binom.rvs(1, p, size=n+1) # leave room to start at 1\n", | |
" terms = np.ones(n+1)\n", | |
" sfac = p1/p0\n", | |
" ffac = (1.0-p1)/(1.0-p0)\n", | |
" terms[trials == 1.0] = sfac\n", | |
" terms[trials == 0.0] = ffac\n", | |
" terms[0] = 1.0\n", | |
" logterms = np.log(terms)\n", | |
" #\n", | |
" ax[0].plot(range(n+1),np.cumprod(terms), color='b')\n", | |
" ax[0].axhline(y=(1-beta)/alpha, xmin=0, xmax=n, color='g', label=r'$(1-\\beta)/\\alpha$')\n", | |
" ax[0].axhline(y=beta/(1-alpha), xmin=0, xmax=n, color='r', label=r'$\\beta/(1-\\alpha)$')\n", | |
" ax[0].set_title('Simulation of Wald\\'s SPRT')\n", | |
" ax[0].set_ylabel('LR')\n", | |
" ax[0].legend(loc='best')\n", | |
" #\n", | |
" ax[1].plot(range(n+1),np.cumsum(logterms), color='b', linestyle='--')\n", | |
" ax[1].axhline(y=math.log((1-beta)/alpha), xmin=0, xmax=n, color='g', label=r'$log((1-\\beta)/\\alpha)$')\n", | |
" ax[1].axhline(y=math.log(beta/(1-alpha)), xmin=0, xmax=n, color='r', label=r'$log(\\beta/(1-\\alpha))$')\n", | |
" ax[1].set_ylabel('log(LR)')\n", | |
" ax[1].set_xlabel('trials')\n", | |
" ax[1].legend(loc='best')\n", | |
"\n", | |
"\n", | |
"interact(plotBinomialSPRT, n=widgets.IntSlider(min=5, max=300, step=5, value=100),\\\n", | |
" p=widgets.FloatSlider(min=0.001, max=1, step=0.01, value=.45),\\\n", | |
" p0=widgets.FloatSlider(min=0.001, max=1, step=0.01, value=.5),\\\n", | |
" p1=widgets.FloatSlider(min=0.001, max=1, step=0.01, value=.6),\\\n", | |
" alpha=widgets.FloatSlider(min=0.001, max=1, step=0.01, value=.05),\\\n", | |
" beta=widgets.FloatSlider(min=0.001, max=1, step=0.01, value=.05)\n", | |
" )" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## Behavior\n", | |
"For $p_0 < p_1$,\n", | |
"+ when $p \\ge p_1$, the likelihood ratio is likely to cross the upper (green) line eventually and unlikely to cross the lower (red) line.\n", | |
"+ when $p \\le p_0$, the likelihood ratio is likely to cross the lower (red) line eventually and unlikely to cross the upper (green) line.\n", | |
"+ the SPRT approximately minimizes the expected number of trials before rejecting one or the other hypothesis, provided $p \\notin (p_0, p_1)$.\n", | |
"\n", | |
"For $p_1 < p_0$, the directions are reversed." | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## It works for $\\beta = 0$ too: $P$-values\n", | |
"The inequalities hold when $\\beta = 0$ also. Then the likelihood ratio has chance at most $\\alpha$ of ever being greater than $1/\\alpha$, if in fact $p > p_0$.\n", | |
"\n", | |
"Hence, $1/\\mbox{LR}$ is the $P$-value of the hypothesis $p > p_0$. This can be used to construct one-sided confidence bounds for $p$ and other parameters. The next chapter does exactly that, to find a lower confidence bound for the mean of a nonnegative population." | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## Dependent observations\n", | |
"As an illustration of SPRT in for obervations that are not iid, consider the following problem.\n", | |
"\n", | |
"There is a population of $B$ items. Item $b$ has \"weight\" $N_b$ and \"value\" $a_b \\in \\{0, 1\\}$. Let $N = \\sum_{b=1}^B N_b$. (You might think of items as manufacturing lots, where $a_b = 1$ means the lot is acceptable and $a_b=0$ means the lot is defective.)\n", | |
"\n", | |
"We want to test the hypothesis $H_0$ that $\\frac{1}{N}\\sum_b N_b a_b = 1/2$ against the\n", | |
"alternative hypothesis $H_1$ that $\\frac{1}{N}\\sum_b N_b a_b = \\gamma $, for some\n", | |
"fixed $\\gamma > 1/2$.\n", | |
"\n", | |
"We will draw items sequentially, without replacement, such that the chance that item $b$ is selected in draw $i$, assuming it has not been selected already, is $N_b/\\sum_{j \\notin {\\mathcal B_{i-1}}} N_j$,\n", | |
"where ${\\mathcal B_{i-1}}$ is the sample up to and including the $i-1$st draw,\n", | |
"and ${\\mathcal B_0} \\equiv \\emptyset$. \n", | |
"\n", | |
"Let $\\mathbb B_i$ denote the item selected at random in such a manner in the $i$th draw.\n", | |
"\n", | |
"The chance that the first draw ${\\mathbb B_1}$ gives an item with value 1, i.e., \n", | |
"$\\Pr \\{a_{\\mathbb B_1} = 1\\}$, is $\\frac{1}{N}\\sum_b N_b a_b$.\n", | |
"Under $H_0$, this chance is $p_{01} = 1/2$; under $H_1$, this chance is \n", | |
"$p_{11} = \\gamma$.\n", | |
"\n", | |
"Given the values of $\\{a_{\\mathbb B_k}\\}_{k=1}^i$, the conditional\n", | |
"probability that the $i$th draw gives an item with value 1 is\n", | |
"\n", | |
"$$\n", | |
" \\Pr \\{a_{\\mathbb B_i} = 1 | {\\mathcal B_{i-1}} \\} = \\frac{ \\sum_{b \\notin {\\mathcal B_{i-1}}} N_b a_b}{\\sum_{b \\notin {\\mathcal B_{i-1}}} N_b}.\n", | |
"$$\n", | |
"\n", | |
"Under $H_0$, this chance is\n", | |
"\n", | |
"$$\n", | |
" p_{0i} = \\frac{N/2 - \\sum_{b \\in {\\mathcal B_{i-1}}} N_b a_b}{N - \\sum_{b \\in {\\mathcal B_{i-1}}} N_b}.\n", | |
"$$\n", | |
"\n", | |
"Under $H_1$, this chance is\n", | |
"\n", | |
"$$\n", | |
" p_{1i} = \\frac{N \\gamma - \\sum_{b \\in {\\mathcal B_{i-1}}} N_b a_b}{N - \\sum_{b \\in {\\mathcal B_{i-1}}} N_b}.\n", | |
"$$\n", | |
"\n", | |
"Let $X_k$ be the indicator of the event that the $k$th draw gives an item with\n", | |
"value $1$, i.e., the indicator of the event $a_{\\mathbb B_k} = 1$.\n", | |
"The likelihood ratio for a given sequence $\\{X_k\\}_{k=1}^i$ is\n", | |
"\n", | |
"$$\n", | |
" \\mbox{LR} = \\frac{\\prod_{k=1}^i p_{1k}^{X_k}(1-p_{1k})^{1-X_k}}\n", | |
" {\\prod_{k=1}^i p_{0k}^{X_k}(1-p_{0k})^{1-X_k}}.\n", | |
"$$\n", | |
"\n", | |
"This can be simplified: $p_{0k}$ and $p_{1k}$ have the same denominator,\n", | |
"$N - \\sum_{b \\in {\\mathcal B_{i-1}}} N_b$,\n", | |
"and their numerators share a term.\n", | |
"Define $N(k) \\equiv \\sum_{j = 1}^{k-1}N_b$ and\n", | |
"$A(k) \\equiv \\sum_{j = 1}^{k-1}N_b a_b$.\n", | |
"Then\n", | |
"\n", | |
"$$\n", | |
" \\mbox{LR} = \\prod_{k=1}^i \n", | |
" \\left ( \\frac{N\\gamma - A(k)}{N/2 - A(k)} \\right )^{X_k}\n", | |
" \\left ( \\frac{N-N\\gamma - (N(k) - A(k))}{N-N/2 - (N(k)-A(k))} \\right )^{1-X_k}\n", | |
"$$\n", | |
"$$\n", | |
" = \\prod_{k=1}^i\n", | |
" \\left ( {N\\gamma - A(k)}\\frac{N/2 - A(k)} \\right )^{X_k}\n", | |
" \\left ( \\frac{N(1-\\gamma) - N(k) + A(k)}{N/2 - N(k) + A(k)} \\right )^{1-X_k},\n", | |
"$$\n", | |
"where the products are defined to be infinite if the denominator vanishes anywhere.\n", | |
"\n", | |
"If $H_0$ is true, the chance that $\\mbox{LR}$ is ever greater than $1/\\alpha$\n", | |
"is at most $\\alpha$." | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"## What's next?\n", | |
"+ [Optional next: The SPRT for the population percentage, sampling without replacement](pSPRTnoReplace.ipynb)\n", | |
"+ [Next: The Kaplan-Wald Confidence Bound for a Nonnegative Mean](kaplanWald.ipynb)\n", | |
"+ [Lower confidence bounds for the mean of a nonnegative population: Markov's Inequality & methods based on the empirical distribution](markov.ipynb)\n", | |
"+ [Index](index.ipynb)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 3, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/html": [ | |
"<link href='http://fonts.googleapis.com/css?family=Lora|Open+Sans' rel='stylesheet' type='text/css'>\n", | |
"<style>\n", | |
"\n", | |
"font-family: 'Lora', serif;\n", | |
"\n", | |
".MathJax_Display {\n", | |
" padding: 10px;\n", | |
"}\n", | |
"\n", | |
"div.callout {\n", | |
" color: #000000;\n", | |
" background-color: #DDDDDD;\n", | |
" margin: 20px 20px 20px 20px;\n", | |
" border-style: solid;\n", | |
" border-width: 2px;\n", | |
" padding: 10px 10px;\n", | |
"}\n", | |
"\n", | |
".rendered_html {\n", | |
" color: #2C5494;\n", | |
" font-family: 'Lora', serif;\n", | |
" font-size: 140%;\n", | |
" line-height: 1.1;\n", | |
" margin: 0.5em 0;\n", | |
"}\n", | |
"\n", | |
"div.cell {\n", | |
" width:900px;\n", | |
" margin-left:auto;\n", | |
" margin-right:auto;\n", | |
"}\n", | |
"\n", | |
"div.text_cell_render {\n", | |
" width:800px;\n", | |
" margin-left:auto;\n", | |
" margin-right:auto;\n", | |
"}\n", | |
"\n", | |
"\n", | |
".title {\n", | |
" font-family: 'Open Sans', sans-serif;\n", | |
" color: #4773D1;\n", | |
" font-size: 250%;\n", | |
" font-weight:bold;\n", | |
" line-height: 1.2;\n", | |
" margin: 0.5em 0;\n", | |
"}\n", | |
"\n", | |
".subtitle {\n", | |
" font-family: 'Open Sans', sans-serif;\n", | |
" color: #386BBC;\n", | |
" font-size: 180%;\n", | |
" font-weight:bold;\n", | |
" line-height: 1.2;\n", | |
" margin: 0.5em 0;\n", | |
"}\n", | |
"\n", | |
".slide-header, p.slide-header {\n", | |
" color: #4773D1;\n", | |
" font-size: 200%;\n", | |
" font-weight:bold;\n", | |
" margin: 0px 20px 10px;\n", | |
" page-break-before: always;\n", | |
" text-align: center;\n", | |
"}\n", | |
"\n", | |
".rendered_html .code {\n", | |
" background-color: #999999;\n", | |
"}\n", | |
"\n", | |
".rendered_html h1 {\n", | |
" color: #7898DD;\n", | |
" font-family: 'Open Sans', sans-serif;\n", | |
" line-height: 1.2;\n", | |
" margin: 0.15em 0em 0.5em;\n", | |
" page-break-before: always;\n", | |
"}\n", | |
"\n", | |
"\n", | |
".rendered_html h2 {\n", | |
" color: #4773D1;\n", | |
" line-height: 1.2;\n", | |
" margin: 1.1em 0em 0.5em;\n", | |
"}\n", | |
"\n", | |
".rendered_html h3 {\n", | |
" font-size: 100%;\n", | |
" line-height: 1.2;\n", | |
" margin: 1.1em 0em 0.5em;\n", | |
"}\n", | |
"\n", | |
".rendered_html .definition, .proposition, .proof, .theorem {\n", | |
" padding-top: 20px;\n", | |
" color: #222299;\n", | |
" font-size: 120%;\n", | |
" font-style: italic;\n", | |
"}\n", | |
"\n", | |
".definition, .proposition, .theorem {\n", | |
" background-color: #EEEEEE;\n", | |
" border-style: solid;\n", | |
" border-width: 2px;\n", | |
" padding-left: 20px;\n", | |
" padding-right: 20px;\n", | |
"}\n", | |
"\n", | |
".rendered_html .definition::before{\n", | |
" content: \"Definition:\";\n", | |
" background-color: #DDDDDD;\n", | |
" color: #222299;\n", | |
" font-size: 110%;\n", | |
" font-weight: bold;\n", | |
" font-style: normal;\n", | |
"}\n", | |
"\n", | |
".rendered_html .proposition::before{\n", | |
" content: \"Proposition:\";\n", | |
" background-color: #DDDDDD;\n", | |
" color: #222299;\n", | |
" font-size: 110%;\n", | |
" font-weight: bold;\n", | |
" font-style: normal;\n", | |
"}\n", | |
"\n", | |
".rendered_html .proof::before{\n", | |
" content: \"Proof:\";\n", | |
" background-color: #DDDDDD;\n", | |
" color: #222299;\n", | |
" font-size: 110%;\n", | |
" font-weight: bold;\n", | |
" font-style: normal;\n", | |
"}\n", | |
"\n", | |
".rendered_html .theorem::before{\n", | |
" content: \"Theorem:\";\n", | |
" background-color: #DDDDDD;\n", | |
" color: #222299;\n", | |
" font-size: 110%;\n", | |
" font-weight: bold;\n", | |
" font-style: normal;\n", | |
"}\n", | |
"\n", | |
"\n", | |
".rendered_html ol li {\n", | |
" padding-top: 20px;\n", | |
" padding-bottom: -20px;\n", | |
" line-height: 1.5;\n", | |
"}\n", | |
"\n", | |
".rendered_html ul li {\n", | |
" padding-top: 0px;\n", | |
" padding-bottom: 0px;\n", | |
" line-height: 1.2;\n", | |
"}\n", | |
"\n", | |
"li:first-of-type {\n", | |
" padding-top: -100px;\n", | |
"}\n", | |
"\n", | |
".input_prompt, .CodeMirror-lines, .output_area {\n", | |
" font-family: Consolas;\n", | |
" font-size: 120%;\n", | |
"}\n", | |
"\n", | |
".gap-above {\n", | |
" padding-top: 200px;\n", | |
"}\n", | |
"\n", | |
".gap01 {\n", | |
" padding-top: 10px;\n", | |
"}\n", | |
"\n", | |
".gap05 {\n", | |
" padding-top: 50px;\n", | |
"}\n", | |
"\n", | |
".gap1 {\n", | |
" padding-top: 100px;\n", | |
"}\n", | |
"\n", | |
".gap2 {\n", | |
" padding-top: 200px;\n", | |
"}\n", | |
"\n", | |
".gap3 {\n", | |
" padding-top: 300px;\n", | |
"}\n", | |
"\n", | |
".emph {\n", | |
" color: #386BBC;\n", | |
"}\n", | |
"\n", | |
".warn {\n", | |
" color: red;\n", | |
"}\n", | |
"\n", | |
".center {\n", | |
" text-align: center;\n", | |
"}\n", | |
"\n", | |
".nb_link {\n", | |
" padding-bottom: 0.5em;\n", | |
"}\n", | |
"\n", | |
"</style>" | |
], | |
"text/plain": [ | |
"<IPython.core.display.HTML object>" | |
] | |
}, | |
"metadata": {}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"%run talkTools.py" | |
] | |
} | |
], | |
"metadata": { | |
"anaconda-cloud": {}, | |
"kernelspec": { | |
"display_name": "Python 2", | |
"language": "python", | |
"name": "python2" | |
}, | |
"language_info": { | |
"codemirror_mode": { | |
"name": "ipython", | |
"version": 2 | |
}, | |
"file_extension": ".py", | |
"mimetype": "text/x-python", | |
"name": "python", | |
"nbconvert_exporter": "python", | |
"pygments_lexer": "ipython2", | |
"version": "2.7.11" | |
} | |
}, | |
"nbformat": 4, | |
"nbformat_minor": 1 | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment