Skip to content

Instantly share code, notes, and snippets.

@fperez
Created March 16, 2017 07:09
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save fperez/a968bacd657a504d4de1d127d65b2a80 to your computer and use it in GitHub Desktop.
Save fperez/a968bacd657a504d4de1d127d65b2a80 to your computer and use it in GitHub Desktop.
sprt test - adding an nbsp in the gap0 paragraph in 2nd cell
Display the source blob
Display the rendered blob
Raw
{
"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\">&nbsp;</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&mdash;unlike the classical Gambler's Ruin problem in which the game is fair&mdash;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