{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Duality between confidence sets and hypothesis tests\n", "\n", "### We will observe $X \\sim \\mathbb P_\\mu$, where $\\mu \\in \\Theta$. \n", "+ $\\Theta$ is known\n", "+ $\\theta \\rightarrow \\mathbb P_\\theta$ is known\n", "+ $\\mu$ is unknown\n", "+ $X$ takes values in $\\mathcal X$.\n", "\n", "(We will ignore issues of measurability here: tacitly assume that for all $\\theta \\in \\Theta$, $A_\\eta$ is $\\mathbb P_\\theta$-measurable\n", "and that $\\mathcal I(X)$ is set-valued $\\mathbb P_\\theta$-measurable function.)\n", "\n", "
\n", "

Definition

\n", " \n", "$A_\\theta \\subset \\mathcal X$ is the acceptance region for a level-$\\alpha$ test of the hypothesis $\\mu = \\theta$ iff\n", "\\begin{equation*} \\mathbb P_\\theta (X \\notin A_\\theta) \\le \\alpha.\\end{equation*}\n", " \n", "
\n", "\n", "\n", "
\n", "

Definition

\n", " \n", "$\\mathcal I(X)$ is a $1-\\alpha$ confidence set for $\\mu$ iff \n", "\\begin{equation*} \\forall \\theta \\in \\Theta, \\;\\;\\; \\mathbb P_\\theta ( \\mathcal I(X) \\ni \\theta) \\ge 1-\\alpha.\\end{equation*}\n", "\n", "
\n", "\n", "\n", "
\n", "

Proposition

\n", "\n", "Suppose \n", "\\begin{equation*} \\{A_\\theta: \\theta \\in \\Theta \\}\\end{equation*}\n", "is a family of level-$\\alpha$ acceptance regions. Then \n", "\\begin{equation*} \\mathcal I(X) := \\{ \\theta \\in \\Theta: X \\in A_\\theta \\}\\end{equation*}\n", "is a $1-\\alpha$ confidence set for $\\mu$.\n", " \n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "

Proof

\n", "\n", "For any $\\theta \\in \\Theta$,\n", "\\begin{equation*} \\mathbb P_\\theta \\left (\\{ \\eta \\in \\Theta: X \\in A_\\eta \\} \\ni \\theta \\right ) = \n", " \\mathbb P_\\theta ( X \\in A_\\theta ) \\end{equation*}\n", "\\begin{equation*} \\ge 1-\\alpha.\\end{equation*}\n", " \n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## We will use this approach—inverting tests—to construct confidence sets.\n", "\n", "+ We are interested in _confidence bounds_ for real-valued parameters, typically 1-sided bounds\n", "+ Typically, one-sided tests for a one-dimensional location parameter lead to one-sided confidence bounds for the parameter\n", "+ Inverting a family of tests that rejects $\\mu = \\theta$ when $X$ is sufficiently small typically leads to an _upper_ confidence bound for $\\mu$\n", "+ Inverting a family of tests that rejects $\\mu = \\theta$ when $X$ is sufficiently large typically leads to a _lower_ confidence bound for $\\mu$\n", "+ Inverting a family of tests that rejects $\\mu = \\theta$ when $X$ is sufficiently \"extreme\" typically leads to a confidence interval for $\\mu$\n", "\n", "(Exceptions arise from non-monotone likelihood ratios, etc.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example: one-sided binomial tests\n", "\n", "Consider $n$ independent, uniform draws (i.e., a random sample with replacement) from a $\\{0, 1\\}$ population containing a fraction $p$ of 1s.\n", "Let $X$ denote the number of 1s in the random sample. Then $X$ has a Binomial distribution with parameters $n$ and $p$.\n", "\n", "In particular,\n", "\\begin{equation*} \\mathbb P_{n,p} (X = k) = {n \\choose k} p^k (1-p)^{n-k}\\end{equation*}\n", "for $0 \\le k \\le n$.\n", "\n", "To find a one-sided lower confidence bound for $p$ in a Binomial$(n,p)$ distribution, with $n$ known, we would invert one-sided upper tests, that is, tests that reject the hypothesis $p = \\pi$ when $X$ is large.\n", "\n", "The form of the acceptance region for the test is:\n", "\\begin{equation*} A_\\pi := \\ [0, a_\\pi],\\end{equation*}\n", "where \n", "\\begin{equation*} a_\\pi := \\min \\left \\{ k: \\sum_{i=k+1}^n {{n}\\choose{i}} \\pi^i (1-\\pi)^{n-i} \\le \\alpha \\right \\}.\\end{equation*}\n", "\n", "Let's see this graphically:" ] }, { "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", "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", "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": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHARJREFUeJzt3X2UFPW95/H3Z5hgHjQaNTBHhMFIiIvxgiQimyfaGFdC\nsqJuHiC70WtiMidRMZGDoCceRj3Gh3i9Vw+jBOVyoibxKW5E4wN3o22OuAuDSII4CD4hoM4q4CIG\nEPG7f3Q76Znpmame6Zmhaz6vc/rYVfWrql9ZPZ/+8auuXykiMDOzdKnq7wqYmVn5OdzNzFLI4W5m\nlkIOdzOzFHK4m5mlkMPdzCyFEoW7pMmS1kpaJ2l2keWfkfSkpF2SLmiz7EBJd0tqkrRG0vHlqryZ\nmRVX3VUBSVXAPOBE4FWgUdJ9EbG2oNgW4Dzg1CKbuB54MCK+Laka+GjPq21mZp1J0nKfAKyPiA0R\nsQe4A5haWCAi3oyIp4D3CudL+jjw5YhYlC/3XkRsL0/VzcysI0nCfRiwsWB6U35eEkcAb0paJGml\npAWSPlJqJc3MrDS9fUG1GhgPNETEeODvwJxe3qeZ2YDXZZ87sBkYUTB9eH5eEpuAjRGxIj99D9Du\ngiyAJA9yY2ZWoohQsflJWu6NwChJtZIGA9OAxZ2Ub9lRRDQDGyWNzs86EXi2k0p2+Zo7d26icpXy\nStvxpPGY5ib8bFbKK23nJ43HlPR4OtNlyz0i9ko6F1hC7stgYUQ0SarLLY4FkoYCK4ADgPclnQ+M\niYgdwAzgt5I+BLwInNXVPs3MrGeSdMsQEQ8Dn2kz79cF75uB4R2s+1fguB7U0czMSlRxd6hmMpn+\nrkJZpe14IH3HlOnvCpRZ2s4PpO+YynE86qrfpq9Iin2lLmatSODPpu2DJBEdXFBN1C1jZn1j5MiR\nbNiwob+rYfuY2tpaXn755ZLWccvdrCt92HLPt8T6ZF9WOTr6XHTWcq+4PnczM+uaw93MLIUc7mZm\nKeRwNzNLIYe7mVkKOdzNrOxeeuml/q7CgOdwN9uHjRhRg6Ree40YUVP2Or/00kssW7as7NvtqVde\neYU777yzv6vRZ3wTk9k+bOPGZh57rPe2f8IJzYnL3njjjTQ0NNDU1MSUKVOQxHPPPUddXR0zZ85s\nKTd//nyuvvrq3qhuO7/73e947bXXWL58OaeddhrTpk0D4L777mPNmjUMGjSIww47jO9///uMGDGC\nP//5zzz77LOMGTOmLPt/+umneeSRR5gzp+PHVCxatIjNmzczePBgRo8ezamntn4aaZJtdIdb7maW\nyE9/+lO+853vIImbbrqJ+++/nzPOOINZs2bxyCOPAPC3v/2N4cNbjyHY2NjIt771LQ4++GD+9Kc/\nAfD73/+eYcOGcemll/LWW2+1lL3iiit49tkORwVv5YUXXmDLli3MnDmThoYGfvKTn/Dyyy+zfft2\nLrvsMi6++GJmz57NjTfeyJYtWwD43ve+x7x588rxv4OI4JJLLuHdd9/tsMwzzzzDokWL+MUvfsGF\nF15IQ0MDu3fvLmkb3eVwN7OSfXC35LBhuSdurl69GoD777+fE044oVXZ4447jl//+tfs2bOHY445\nBoA33niD1atXM3fuXA466KCWsnPmzGHZsmVcdNFFPP30053WYc2aNfzqV78C4NBDD2XUqFGsWLGC\nv/zlLxx99NEt5caOHctj+X/+7Lfffrz77rvs2LGj3fZuuukmvvnNbzJr1iy2b2/9qOcPvrwK/eEP\nf2h3rG09/PDDHHHEES3TQ4YMYenSpSVto7sc7mbWbevXr6e6uroloBobG4t2eRxyyCF8+9vfpqGh\ngYaGBqZPn87BBx/crtygQYM466yzuPLKK1m3bh0XXXQRTz75ZNF9T5kyhQcffLBl+rXXXmPUqFFs\n2rSp1RfGQQcdxPr161umx44d226bq1evZtCgQTzwwAPMmDGD888/v6XMihUr2Lp1a6vyW7Zsoaqq\nikMPPbTT/z/7778/e/bsaZnetWsXTU1NJW2juxzuZlaSiGD27NmceeaZPP/88yxZsoTPfe5zAOzc\nuROp6FAnnHPOOVx33XWcdNJJfPKTn+xyP9/97ne58sor2bp1Kz/84Q/ZtGlTq+XV1dV89rOfBeCB\nBx7g85//POPGjWPbtm18+MMfbik3ePDgVi31ww47rFXYA2zbto0f//jHAAwfPpxFixaxdOlSTjnl\nFO677z6mT5/eqvy9997L6aef3uUxnH766bzwwgsA7Nixg+eee4633367pG10ly+omllJJHHhhRdy\n7LHHtlu2d+/eDtdbt24dY8aM4fHHH2f06NEdlivU2NjIE088wVe/+tWWLqC2tm/fzm9+8xtuv/12\nAA444IBWLe2dO3dSU/OPXwUddNBBrFu3rtU2vvKVr3Dbbbdx1113cdRRR3HJJZcwa9YsZs2aBeS6\nZU4++WQAli9fzvHHH5+o/kOGDGHRokXcfPPN1NTUcMwxxzBkyJCSttFdDnezAjU1I2lubj3kbkC7\n1ujQobW8/vrLfVexfUxHI1dWVxePlHvvvZexY8dy0UUX8ctf/pIf/ehHnW5/6dKlLF68mAkTJnDV\nVVd1Wvaaa67h5ptvZv/992fDhg0ceeSRrFixomX5li1bGD9+fMv0zp07+djHPtZqG8888wy7d+/m\n/vvv55VXXuGCCy7gBz/4AV/4whdYsWIF27Ztaym7bNkydu7cyUMPPcTSpUvZtWsXixcv5pRTTila\nvzFjxrR0VV122WVcfvnlPPTQQyVtozsc7mYFcsHeNrjUbl5zc/Guh7TrajjioUOH8s4777QKzz/+\n8Y+MHDmSMWPGMHr0aGbOnMnjjz/OpEmT2q2fzWZ5+OGH+dKXvpTo55Tz5s3jtNNOY/fu3TQ2NrJz\n504mTZrE7NmzW8qsXLmy1RfE1q1bW7XkP5h39tlnAzBixAhuueUWrr32Wq666irGjh3L5Zdf3lL2\nvPPOa3l/6aWXIqkllF988UWOOOKIlsbAhg0bOOWUU/jrX/9KU1MTtbW1jBo1qtNtlE3CJ2xPBtYC\n64DZRZZ/BngS2AVcUGR5FbASWNzJPsKsvwGRG7z9H68oMq+3Pq9ttzt8+NB8nXrnNXz40MR1a2ho\niDFjxkRVVVV87Wtfi8cff7xdmYULF8ajjz4aERFLliyJadOmxbhx42Lbtm0REfHUU0/FMcccE8cd\nd1w88MAD7dZfunRp4vo88cQTUVVVFVVVVSEpqqqqYtOmTRERcdttt8Xll18el112Wdx+++2t1ps5\nc2ZLuZ6466674thjj43x48fH3XffHRERxx57bKxcubKlzLvvvhtz586NhoaG+PnPfx5bt27tchvF\ndPR5y88vmqldPqxDUlU+1E8EXgUagWkRsbagzKFALXAqsC0irmuzjZ8DnwM+HhFFv578sA7bF+Ra\nXK0/h4FQkdZ8b3xeK/1hHdu2bePaa6/liiuu6O+qdOjss8/mlltu6e9qlKQ7D+tI0i0zAVgfERvy\nG7sDmEquJQ9ARLwJvCnpm0V2fjgwBbgCuCDB/sysQn3iE5/gkEMOYcuWLRxyyCElr3/NNdewa9eu\nVvMiAkmceeaZ1NbW9qh+jY2NnHTSST3aRqVIEu7DgI0F05vIBX5S/wrMAg4sYR0zq1A/+9nPuPnm\nm6mrqyt53QsvvLAXapSzd+9eHn300Vb98WnWq79zl/QNoDkiVpG7KjUwr0KZDSBVVVXdCvbe9sYb\nbzBjxoz+rkafSdJy3wyMKJg+PD8viS8Cp0iaAnwEOEDSrRFxRrHC9fX1Le8zmQyZTCbhbszMOtf2\nFzKVKJvNks1mE5VNckF1EPAcuQuqrwHLgekR0VSk7FxgR0T8S5Flk4CZvqBq+zJfULV9Ua9cUI2I\nvZLOBZaQ68ZZGBFNkupyi2OBpKHACuAA4H1J5wNjIqL96DxmZtbrumy59xW33G1f4Ja77Yu603L3\nwGFmZink4QfM9iG1tbUdjqpoA1d3ft/vbhlLvWKDgXUuSbfMh4HdJDHQBxmz3tNZt4zD3VKvWD96\nJ6Xble2oz72Ubfqzbb3Bfe5mZgOMw93MLIUc7mZmKeRwNzNLIYe7mVkKOdzNzFLI4W5mlkIOdzOz\nFHK4m5mlkMPdzCyFHO5mZinkcDczSyGHu5lZCjnczcxSyOFuZpZCicJd0mRJayWtkzS7yPLPSHpS\n0i5JFxTMP1zSo5LWSFotaUY5K29mZsV1+bAOSVXAOuBE4FWgEZgWEWsLyhwK1AKnAtsi4rr8/Bqg\nJiJWSdofeAqYWrhuwTb8sA7rFX5Yh6VVTx/WMQFYHxEbImIPcAcwtbBARLwZEU8B77WZ/3pErMq/\n3wE0AcO6cQxmZlaCJOE+DNhYML2JbgS0pJHAOGBZqeuamVlpqvtiJ/kumXuA8/Mt+KLq6+tb3mcy\nGTKZTK/XzcysUmSzWbLZbKKySfrcJwL1ETE5Pz0HiIi4ukjZucDbH/S55+dVAw8AD0XE9Z3sx33u\n1ivc525p1dM+90ZglKRaSYOBacDizvbXZvrfgWc7C3YzMyuvLlvukPspJHA9uS+DhRFxlaQ6ci34\nBZKGAiuAA4D3gR3AGGAs8BdgNblmTgAXR8TDRfbhlrv1CrfcLa06a7knCve+4HC33uJwt7TqabeM\nmZlVGIe7mVkKOdzNzFLI4W5mlkIOdzOzFHK4m5mlkMPdzCyFHO5WsWpqRiKpy1f/2y9RPWtqRvZ3\nRS1FfBOTVazkNyf1/01MSevpvwErhW9iMjMbYBzuZmYp5HA3M0shh7uZWQo53M3MUsjhbmaWQg53\nM7MUcribmaWQw93MLIUc7mZmKZQo3CVNlrRW0jpJs4ss/4ykJyXtknRBKeuamVn5dTm2jKQqYB1w\nIvAq0AhMi4i1BWUOBWqBU4FtEXFd0nULtuGxZawkHlvGBrqeji0zAVgfERsiYg9wBzC1sEBEvBkR\nTwHvlbqumZmVX5JwHwZsLJjelJ+XRE/WNTOzbqru7woUqq+vb3mfyWTIZDL9Vhczs31NNpslm80m\nKpukz30iUB8Rk/PTc4CIiKuLlJ0LvF3Q517Kuu5zt5K4z90Gup72uTcCoyTVShoMTAMWd7a/Hqxr\nZmZl0GW3TETslXQusITcl8HCiGiSVJdbHAskDQVWAAcA70s6HxgTETuKrdtrR2NmZoAfs2cVzN0y\nNtD5MXtmZgOMw93MLIUc7mZmKeRwNzNLIYe7mVkKOdzNzFLI4W5mlkIOdzOzFHK4m5mlkMPdzCyF\nHO5mZinkcDczSyGHu5lZCjnczcxSyOFuZpZCDnczsxRyuJuZpZDD3cwshRzuZmYplCjcJU2WtFbS\nOkmzOyhzg6T1klZJGlcw/+eSnpH0N0m/lTS4XJU3M7Piugx3SVXAPOBk4GhguqSj2pT5OnBkRHwa\nqAPm5+cfBpwHjI+IfwKqgWllPQIzM2snSct9ArA+IjZExB7gDmBqmzJTgVsBImIZcKCkofllg4CP\nSaoGPgq8Wpaam5lZh5KE+zBgY8H0pvy8zspsBoZFxKvAvwCv5Oe9FRH/q/vVtbSrqRmJpESv9Nkv\n8bHX1Izs78raPq66Nzcu6SByrfpa4P8B90j6XkT8rlj5+vr6lveZTIZMJtOb1bN9UHPzBiASlk5b\nwO8m6bE3N6ft2C2JbDZLNptNVFYRnX+YJE0E6iNicn56DhARcXVBmfnAYxFxZ356LTAJ+DJwckT8\nKD//+8DxEXFukf1EV3Wx9Mu1yEsJ9yRle7bNQKjd+v1fT/+9mCQioug3fZJumUZglKTa/C9dpgGL\n25RZDJyR39lEct0vzeS6YyZK+rByf7UnAk3dPA4zM0uoy26ZiNgr6VxgCbkvg4UR0SSpLrc4FkTE\ng5KmSHoeeAc4K7/uckn3AE8De/L/XdBbB2NmZjlddsv0FXfLGLhbxt0yVoqedsuYmVmFcbibmaWQ\nw93MLIUc7mZmKeRwNzNLIYe7mVkKOdzNzFLI4W5mlkIOdzOzFHK4m5mlkMPdzCyFHO5mZinkcDcz\nSyGHu5lZCjnczcxSyOFuZpZCDnczsxRyuJuZpZDD3cwshRKFu6TJktZKWidpdgdlbpC0XtIqSeMK\n5h8o6W5JTZLWSDq+XJU3M7Piugx3SVXAPOBk4GhguqSj2pT5OnBkRHwaqAPmFyy+HngwIv4TMBZo\nKlPdzcysA0la7hOA9RGxISL2AHcAU9uUmQrcChARy4ADJQ2V9HHgyxGxKL/svYjYXr7qm5lZMUnC\nfRiwsWB6U35eZ2U25+cdAbwpaZGklZIWSPpITypsZmZdq+6D7Y8HzomIFZL+DZgDzC1WuL6+vuV9\nJpMhk8n0cvXMzCpHNpslm80mKquI6LyANBGoj4jJ+ek5QETE1QVl5gOPRcSd+em1wKT84v8dEZ/K\nz/8SMDsi/muR/URXdbH0kwQk/RwkLduzbQZC7dbv/3r678UkEREqtixJt0wjMEpSraTBwDRgcZsy\ni4Ez8jubCLwVEc0R0QxslDQ6X+5E4NnuHISZmSXXZbdMROyVdC6whNyXwcKIaJJUl1scCyLiQUlT\nJD0PvAOcVbCJGcBvJX0IeLHNMjMz6wVddsv0FXfLGLhbxt0yVoqedsuYmVmFcbibmaWQw93MLIUc\n7mZmKeRwNzNLIYe7mVkKOdzNzFLI4W5mlkIOd+sTNTUjkdTly5LaL9H/z5qakf1dUesnvkPV+kTy\nO0/7/87PSrlDNek2/XeVXr5D1cxsgHG4m5mlkMPdzCyFHO5mZinkcDczSyGHu5lZCjnczcxSyOFu\nZpZCDnczsxRKFO6SJktaK2mdpNkdlLlB0npJqySNa7OsStJKSYvLUWkzM+tcl+EuqQqYB5wMHA1M\nl3RUmzJfB46MiE8DdcD8Nps5H3i2LDU2M7MuJWm5TwDWR8SGiNgD3AFMbVNmKnArQEQsAw6UNBRA\n0uHAFOCWstXazMw6lSTchwEbC6Y35ed1VmZzQZl/BWaRfEQkMzProere3LikbwDNEbFKUobcUHYd\nqq+vb3mfyWTIZDK9WT0zs4qSzWbJZrOJynY55K+kiUB9REzOT88BIiKuLigzH3gsIu7MT68FJpHr\na/8fwHvAR4ADgHsj4owi+/GQvynmIX895K+VX0+H/G0ERkmqlTQYmAa0/dXLYuCM/M4mAm9FRHNE\nXBwRIyLiU/n1Hi0W7GZmVl5ddstExF5J5wJLyH0ZLIyIJkl1ucWxICIelDRF0vPAO8BZvVttMzPr\njJ/EZH3C3TLulrHy85OYzMwGGIe7mVkKOdzNzFLI4W5mlkIOdzOzFHK4m5mlkMPdzCyFHO5mZink\ncDczSyGHu5lZCjnczcxSyOFuZpZCDnfrtpqakUhK9LL+8aEPkfgcjRhR09/VtTLyqJDWbclHeoT+\nHhlxII8K+dhjybZ4wgl4BMkK41EhzcwGGIe7mVkKOdzNzFLI4W5mlkKJwl3SZElrJa2TNLuDMjdI\nWi9plaRx+XmHS3pU0hpJqyXNKGflzcysuC7DXVIVMA84GTgamC7pqDZlvg4cGRGfBuqA+flF7wEX\nRMTRwH8Gzmm7rpmZlV+SlvsEYH1EbIiIPcAdwNQ2ZaYCtwJExDLgQElDI+L1iFiVn78DaAKGla32\nZmZWVJJwHwZsLJjeRPuAbltmc9sykkYC44BlpVbSzMxK0ycXVCXtD9wDnJ9vwZuZWS+qTlBmMzCi\nYPrw/Ly2ZYYXKyOpmlyw3xYR93W2o/r6+pb3mUyGTCaToHpmZgNDNpslm80mKtvl8AOSBgHPAScC\nrwHLgekR0VRQZgpwTkR8Q9JE4N8iYmJ+2a3AmxFxQRf78fADFcbDD+xb9eyonIcfSK/Ohh/osuUe\nEXslnQssIdeNszAimiTV5RbHgoh4UNIUSc8D7wD/nN/xF4H/DqyW9DS5T+PFEfFwWY7MzMyK8sBh\n1m1uue9b9eyonFvu6eWBw8zMBhiHu5lZCjnczcxSyOFuZpZCDnczsxRyuFs7SZ+NaumS9HmrftZq\nZUhyh6oNMM3NG0j+0z1Liz17SPSzyRNOaO79yliPueVuZpZCDnczsxRyuJuZpZDD3cwshRzuZmYp\n5HA3M0shh7uZWQo53M3MUsjhbmaWQg73AcTDClg5eJiCyuDhBwYQDytg5eBhCiqDW+5mZimUKNwl\nTZa0VtI6SbM7KHODpPWSVkkaV8q6ZmZWXl2Gu6QqYB5wMnA0MF3SUW3KfB04MiI+DdQB85OuW6ps\nNtuT1fc5aTuenGx/V6Cssv1dgTJbtaq/a1B+afs7KsfxJGm5TwDWR8SGiNgD3AFMbVNmKnArQEQs\nAw6UNDThuiXxSWwt6UXSvr1Qmu3DffW+bH9XoMz6KtyTXngtx8VX50J7SS6oDgM2FkxvIhfaXZUZ\nlnBd64HkF0nBF0qtLyW98Aq++NobeuuCqlOkh/yzRRtI/PPK8lNE560+SROB+oiYnJ+eA0REXF1Q\nZj7wWETcmZ9eC0wCjuhq3YJtJG1+mplZXkQUbeUl6ZZpBEZJqgVeA6YB09uUWQycA9yZ/zJ4KyKa\nJb2ZYN1OK2hmZqXrMtwjYq+kc4El5LpxFkZEk6S63OJYEBEPSpoi6XngHeCsztbttaMxMzMgQbeM\nmZlVnoq5QzWNN0NJelnSXyU9LWl5f9enVJIWSmqW9LeCeZ+QtETSc5IekXRgf9axVB0c01xJmySt\nzL8m92cdSyHpcEmPSlojabWkGfn5FXmeihzPefn5FXmOJO0naVk+A1ZLmpuf3+PzUxEt9/zNUOuA\nE4FXyV0HmBYRa/u1Yj0k6UXgcxGxrb/r0h2SvgTsAG6NiH/Kz7sa2BIR1+S/hD8REXP6s56l6OCY\n5gJvR8R1/Vq5bpBUA9RExCpJ+wNPkbvX5Cwq8Dx1cjzfpXLP0Ucj4u+SBgFLgRnAf6OH56dSWu5l\nvxlqHyEq5xy0ExFPAG2/mKYCv8m//w1wap9Wqoc6OCao0J/3RsTrEbEq/34H0AQcToWepw6OZ1h+\ncaWeo7/n3+5H7jpoUIbzUynB0tFNUpUugP+Q1CjpR/1dmTIZEhHNkPtDBIb0c33K5Vzlxk26pVK6\nMNqSNBIYB/wfYGiln6eC41mWn1WR50hSlaSngdeB/4iIRspwfiol3NPqixExHpgCnJPvEkibfb/f\nr2s3Ap+KiHHk/gAr8Z/++wP3AOfnW7xtz0tFnacix1Ox5ygi3o+IY8n9i2qCpKMpw/mplHDfDIwo\nmD48P6+iRcRr+f++AfxP0jE0Q7Ny4wp90D/6f/u5Pj0WEW/EPy5O3Qwc15/1KZWkanJBeFtE3Jef\nXbHnqdjxVPo5AoiI7eSGMppMGc5PpYR7y41UkgaTuxlqcT/XqUckfTTf+kDSx4D/AjzTv7XqFtG6\nr3Mx8M/592cC97VdoQK0Oqb8H9cHTqfyztO/A89GxPUF8yr5PLU7nko9R5IO/aALSdJHgJPIXUfo\n8fmpiF/LQO6nkMD1/ONmqKv6uUo9IukIcq31IHcR5beVdkySfgdkgEOAZmAu8EfgbmA4sAH4TkS8\n1V91LFUHx3QCub7d94GXgboP+kP3dZK+CPwFWE3usxbAxcBy4C4q7Dx1cjzfowLPkaRjyF0wrcq/\n7oyIKyQdTA/PT8WEu5mZJVcp3TJmZlYCh7uZWQo53M3MUsjhbmaWQg53M7MUcribmaWQw93MLIUc\n7mZmKfT/AUY5OuwiqAyiAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def binomialHiLite(n, p, alpha):\n", " ''' \n", " Plots probability histogram for a binomial with parameters n and p, \n", " highlighting the upper alpha quantile in yelow.\n", " The blue region corresponds to the acceptance region for a level-alpha\n", " test about p.\n", " Plots a red vertical line at the expected value.\n", " '''\n", " fig, ax = plt.subplots(nrows=1, ncols=1)\n", " x = np.arange(n+1)\n", " val = binom.ppf(1-alpha, n, p)\n", " inx = np.searchsorted(x, val, side=\"right\")\n", " xb = x[:inx]\n", " xy = x[inx:]\n", " width = 1.0\n", " ax.bar(xb, binom.pmf(xb, n, p), width, color='b', align='center')\n", " hilit = ax.bar(xy, binom.pmf(xy, n, p), width, color='y', align='center')\n", " plt.xlim([-width,n+width])\n", " plt.axvline(x=n*p, color='r')\n", " probStr = str(round(100*(1-binom.cdf(x[inx-1],n, p)),2))\n", " label = r'$\\mathbf{P}(X \\geq ' + str(x[inx]) + r') \\approx' + probStr + '$'\n", " plt.legend([hilit[0]], [label], loc = 'best')\n", "\n", "interact(binomialHiLite, n=widgets.IntSlider(min=5, max=300, step=1, value=30),\\\n", " p=widgets.FloatSlider(min=0.001, max=1, step=0.001, value=.5),\\\n", " alpha=widgets.FloatSlider(min=0.0, max=1.0, step=0.001, value=.05)\\\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Inverting the test to find a confidence bound\n", "To turn this family of tests into a lower confidence bound, we need to find\n", "\\begin{equation*} \\min \\{ \\pi: A_\\pi \\ni X \\},\\end{equation*}\n", "that is,\n", "\\begin{equation*} \\min \\left \\{ \\pi: \\sum_{i=X}^n {{n}\\choose{i}} \\pi^i (1-\\pi)^{n-i} > \\alpha \\right \\}.\\end{equation*}\n", "\n", "The following code implements this, using a bisection search." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def bisect(lo, hi, tol, fun):\n", " mid = (lo+hi)/2.0\n", " while (hi-lo)/2 > tol:\n", " if fun(mid) == 0.0:\n", " return mid\n", " elif fun(lo)*fun(mid) < 0.0:\n", " hi = mid\n", " else:\n", " lo = mid\n", " mid = (lo+hi)/2.0\n", " return mid\n", " \n", "def binoLowerCL(n, x, cl = 0.975, inc=0.000001, p = None):\n", " \"Lower confidence level cl confidence interval for Binomial p, for x successes in n trials\"\n", " if p is None:\n", " p = float(x)/float(n)\n", " lo = 0.0\n", " if (x > 0):\n", " f = lambda q: cl - scipy.stats.binom.cdf(x-1, n, q)\n", " lo = bisect(0.0, p, inc, f)\n", " return lo" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's use the code to find a lower confidence bound for $p$ if 50 iid draws from a {0, 1} population yield 40 1s and 10 0s." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.68440322876\n" ] } ], "source": [ "p_lower_bound = binoLowerCL(50, 40, cl=0.95)\n", "print p_lower_bound" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's check that against the probability histogram. Note that reducing $p$ below $0.6844$ drops the upper tail probability below 5%; for $p > 0.6844$ the probability is at least 5%, so the confidence interval is $[0.6844, 1]$:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGPpJREFUeJzt3X+QVOWd7/H3B1BijAuKkSl+DJOoN14oF9RSWZONg0ZF\nEiGVmATiKlohoSpqqCuFg6QSJlpWJOXGjfEHGllL0Y3gjw1gLGU32hrhBkaBKxEQLJVfkbk6jLHg\noiJ+7x/djD1D93QP0zM9c/rzquqizznP6fOc0zOfeXj6PE8rIjAzs2TpU+4KmJlZ6TnczcwSyOFu\nZpZADnczswRyuJuZJZDD3cwsgYoKd0njJW2StFlSXY7tX5K0UtIHkq7LWj9M0rOSXpW0XtJPSll5\nMzPLTYXuc5fUB9gMnA/8DWgAJkfEpqwyxwMjgG8CzRHx68z6KqAqItZJ+hzwMjApe18zMyu9Ylru\nZwFbImJrROwHHgEmZReIiHcj4mXg4zbrd0XEuszzPcBGYGhJam5mZnkVE+5Dge1Zyzs4jICWVAOM\nAVZ1dF8zM+uYbvlANdMl8xgwI9OCNzOzLtSviDI7geqs5WGZdUWR1I90sC+MiCXtlPMkN2ZmHRQR\nyrW+mJZ7A3CSpBGSjgQmA0vbKd/2QP8ObIiI3xRRyUMec+fOzbm+Uh6Vfv6+Bl17/pnfvJZHtCzn\n/n1M4jXoDY9859+egi33iDgg6RpgOek/BgsiYqOk6enNca+kwcBLwDHAJ5JmACOB0cBlwHpJazM/\nNXMi4ulCxzWzcuqPlLNByODBI9i1663urY51WDHdMmTC+Ett1t2T9bwRGJ5j1xVA385U0MzK4UMO\ntuDbamzMHfrWs/T4Eaq1tbXlrkJZVfr5g69BpZ8/+BoczvkXHMTUXSRFT6mLWaVId718+nsXCBGk\nPzrL9/uogv291j0kEXk+UC2qW6acampq2Lp1a7mrYWUyYsQI3nrrrXJXw6zX6fEt98xfpjLUyHoC\nv/9dyy333q29lnuP73M3M7OOc7ibmSWQw93MLIEc7mZmCeRwNzNLIId7N3nzzTfLXQUzqyC9Ltyr\nq6uQ1GWP6uqqktf5zTffZNWqnjeN/bZt21i0aFFJXuuFF17ggw8+4MMPP+TPf/5zSV7TzA5fjx/E\n1Nb27Y0891zXvf64cY1Fl73rrru488472bhxIxMmTEASr732GtOnT2fmzJkt5ebPn8+8efO6orp5\nrV27lmeeeYbZs2cDsGTJEl599VX69u3LkCFDuPzyy6muruZPf/oTGzZsYOTIkZ063tSpU9m6dSsn\nnHAC99xzT7tlc9XFzEqr17Xce5If//jHfPe730USd999N8uWLeOKK65g1qxZPPPMMwC88sorDB/e\nek61hoYGLr30Uo477jj++Mc/AvD73/+eoUOH8otf/IL33nuvpezNN9/Mhg0bOlSviOBnP/sZH330\nEQDvv/8+N954I3PmzKGuro677rqLpqYmAL7//e9zxx13HPY1OOinP/0p27ZtY+fOnUyaNClvufbq\nYmal43AvkYMj9oYOTX8D4fr16wFYtmwZ48aNa1X2zDPP5J577mH//v2ceuqpALzzzjusX7+euXPn\nMnDgwJays2fPZtWqVdxwww2sXbu2qLo8/vjjrY75wgsvMGrUqJbl0aNH81zmvz/9+/fno48+Ys+e\nQ78g6+677+Yb3/gGs2bN4v3332+17eAfr4OOOOIIhg0bRt++7U8C2l5drGtUVdXk7Ya05Op13TI9\n3ZYtW+jXr19LuDY0NDBnzpxDyg0aNIjvfOc73HnnnVRXVzNlyhSOO+64Q8r17duXq666CoBFixax\nePFiLrnkEs4555ycx29qaqJPnz4cf/zx7N27F4AdO3a0+oMxcOBAtmzZ0rI8evRoVq5cyYUXXtiy\nbv369fTt25cnn3yS7du3M2PGDH74wx9yzjnn8NJLL7F79+5Wx129ejURQVNTEyeffDITJ07MWb9C\ndbHSa2zcSntTCVgyueVeAhFBXV0dU6dO5fXXX2f58uWcccYZAOzbty9vC+nqq6/m17/+NRdccAGf\n//znCx7ne9/7Hr/85S/ZvXs3P/jBD9ixY8chZZ544gm+9a1vtVrX3NzMZz7zmZblI488slVLfciQ\nIYcEbHNzMz/60Y8AGD58OPfffz8rVqxg4sSJLFmyhClTprQqP23aNK688kpmzpzJz3/+c/7+97/n\nPIdCdTGz0nC4l4Akrr/+eh544AEWL17cau7lAwcO5N1v8+bNjBw5kueff77oYzU0NPDiiy9y3nnn\ntXQBHbR69WrOPvvsQ/Y55phjWk30tG/fvlb/Sxg4cOAh3S5f/epXWbhwIZdccklLt8ysWbNYunQp\nN9100yHdMqNHj255fuyxx5JKpXLWv1BdzKw0HO4lkm+WvH79cvd8PfHEE4wePZobbriB3/72twVf\nf8WKFdTV1bFt2zZuueUWLrvsskP+R7Bq1Sqefvpp5s2bx+OPP86LL77I0qVLOfHEE3nnnXdayjU1\nNTFkyJCW5X379nH00Ue3eq2//vWvfPjhhyxbtoxrr72W6667jpUrVwLw0ksv0dzc3FL24YcfbnXH\ny549e/L2vReqi5mVhsO9kwpNfTp48OCWvu+D/vCHP1BTU8PIkSO59NJLaWpqytt6T6VSzJ49m+bm\nZubNm8e3v/3tvMe69tpruf7666mrq+OMM87gK1/5ChMnTuTcc89lzZo1LeXWrFnD+eef37K8e/du\nqqpa39+/e/dupk2bBkB1dTX33XcfK1eubOmWmTx5ckvZmpoapk+fDsDevXt59913Oe+88wB44403\nWl2jQnUxsxIp97d6Z32Ld+TSdv3w4YM//ar2LngMHz44Zz1yufPOO2PkyJHRp0+f+NrXvhbPP//8\nIWUWLFgQzz77bERELF++PCZPnhxjxoyJ5ubmiIh4+eWX49RTT40zzzwznnzyyUP2X7FiRdH1OWjx\n4sVx2mmnxemnnx6PPvpoREQsXLgwbrrpprjxxhvjoYcealV+5syZsWPHjg4fJ9tDDz0Ut912W8yY\nMSP+8pe/tKw/7bTTYs2aNa3KtleXtvL9XFjx0j/bkefRelu0LLe3T/+cvzuDB48o96lWnMzvR85M\n9Zd1dLHm5mZuvfVWbr755nJXJa9p06Zx3333lbsaOfX2978naPuFHG22ttpW7Jd15N7m96q79eqv\n2evtjj32WAYNGkRTUxODBg3q8P6/+tWv+OCDD1qtiwgkMXXqVEaMGNGp+jU0NHDBBRd06jXMrOdx\ny70bfPLJJ/zud79r6ZfuKQ4cOMCtt95KXV1duauSVxLe/3Jzyz252mu5O9wr2K5duxgwYABHHXVU\nuauSl9//znO4J5fD3Xotv/+d53BPLn9BtplZhXG4m5klkMPdzCyBigp3SeMlbZK0WdIht1ZI+pKk\nlZI+kHRdR/Y1M7PSK3ifu6Q+wB3A+cDfgAZJSyJiU1axJuBa4JuHsW+7RowY4XmnK1hn7+M3q1TF\nDGI6C9gSEVsBJD0CTAJaAjoi3gXelfSNju5byFtvvVVsUTMzyyimW2YosD1reUdmXTE6s6+ZmR0m\nf6BqZpZAxXTL7ASqs5aHZdYVo0P71tfXtzyvra1t9aUXZmaVLpVK5f0inLYKjlCV1Bd4jfSHom8D\nq4EpEbExR9m5wJ6I+NfD2DfnCFUz6xyPUE2uTs0KGREHJF0DLCfdjbMgIjZKmp7eHPdKGgy8BBwD\nfCJpBjAyIvbk2rdE52VmZnn0+LllzKxz3HJPLs8tY2ZWYRzuZmYJ5HA3M0sgh7uZWQI53M3MEsjh\nbmaWQA53M7MEcribmSWQw90sIaqqapB0yMMqk0eomiVE/pGoHqGaVB6hamZWYRzuZmYJ5HA3M0sg\nh7uZWQI53M3MEsjhbmaWQA53M7MEcribmSWQw93MSqR/zhGyVVU15a5YRfIIVbOE6AkjVD1ytXt5\nhKqZWYVxuJuZJZDD3cwsgRzuZmYJ5HA3M0sgh7uZWQI53M3MEsjhbmaWQA53M7MEKircJY2XtEnS\nZkl1ecrcLmmLpHWSxmSt/1+S/irpFUkPSzqyVJU3M7PcCoa7pD7AHcBFwChgiqRT2pS5GDgxIk4G\npgPzM+uHANcCp0fEPwL9gMklPQMzMztEMS33s4AtEbE1IvYDjwCT2pSZBDwIEBGrgAGSBme29QWO\nltQP+Czwt5LU3MzM8iom3IcC27OWd2TWtVdmJzA0Iv4G/CuwLbPuvYj478OvrpmZFaNfV764pIGk\nW/UjgL8Dj0n6fkT8R67y9fX1Lc9ra2upra3tyuqZmfUqqVSKVCpVVNmCU/5KGgvUR8T4zPJsICJi\nXlaZ+cBzEbEos7wJOBf4Z+CiiPhhZv3lwNkRcU2O43jKX7NO8JS/laezU/42ACdJGpG502UysLRN\nmaXAFZmDjSXd/dJIujtmrKTPKP2Tdz6w8TDPw8zMilSwWyYiDki6BlhO+o/BgojYKGl6enPcGxFP\nSZog6XVgL3BVZt/Vkh4D1gL7M//e21UnY2Zmaf4mJrOEcLdM5fE3MZmZVRiHu5lZAjnczcwSyOFu\nZpZADnczswRyuJuZJZDD3cwsgRzuZmYJ5HA3M0sgh7tZL1JVVYOknA+zbJ5+wKwXyT/FAHR8WoBD\nt3n6gd7F0w+YmVUYh7uZWQI53M3MEsjhbmaWQA53M7MEcribmSWQw93MLIEc7mZmCeRwNzNLIIe7\nmVkCOdzNzBLI4W5mXax/3snOqqpqyl25xPLEYWa9SG+dOKy91/Lv/eHzxGFmZhXG4W5mlkAOdzOz\nBHK4m5klUFHhLmm8pE2SNkuqy1PmdklbJK2TNCZr/QBJj0raKOlVSWeXqvJmZpZbwXCX1Ae4A7gI\nGAVMkXRKmzIXAydGxMnAdGB+1ubfAE9FxP8ERgMbS1R3MzPLo5iW+1nAlojYGhH7gUeASW3KTAIe\nBIiIVcAASYMl/QPwzxFxf2bbxxHxfumqb2ZmuRQT7kOB7VnLOzLr2iuzM7PuC8C7ku6XtEbSvZKO\n6kyFzcyssH7d8PqnA1dHxEuS/g2YDczNVbi+vr7leW1tLbW1tV1cPTOz3iOVSpFKpYoqW3CEqqSx\nQH1EjM8szwYiIuZllZkPPBcRizLLm4BzM5v/d0R8MbP+K0BdRFyS4zgeoWpWgEeoWrbOjlBtAE6S\nNELSkcBkYGmbMkuBKzIHGwu8FxGNEdEIbJf0PzLlzgc2HM5JmJlZ8Qp2y0TEAUnXAMtJ/zFYEBEb\nJU1Pb457I+IpSRMkvQ7sBa7KeomfAA9LOgJ4o802MzPrAp44zKwXcbeMZfPEYWZmFcbhbmaWQA53\nM7MEcribmSWQw93MLIEc7mZmCeRwNzNLIIe7mVkCOdzNzBLI4W5mlkAOd7MeqKqqBkmHPMyK5bll\nzHqg/HPIeG4Z+5TnljEzqzAOdzOzBHK4m5klkMPdzCyBHO5mZgnkcDczSyCHu5lZAjnczcwSyOFu\nZpZADnczswRyuJuZJZDD3cwsgRzuZmYJ5HA3M0sgh7uZlVH/nPPWV1XVlLtivZ7nczfrgSppPvd8\n+zgPCuv0fO6SxkvaJGmzpLo8ZW6XtEXSOklj2mzrI2mNpKUdr76ZmXVUwXCX1Ae4A7gIGAVMkXRK\nmzIXAydGxMnAdGB+m5eZAWwoSY3NzKygYlruZwFbImJrROwHHgEmtSkzCXgQICJWAQMkDQaQNAyY\nANxXslqbmVm7ign3ocD2rOUdmXXtldmZVeY2YBb5O93MzKzE+nXli0v6OtAYEesk1ZL+9CSv+vr6\nlue1tbXU1tZ2ZfXMzHqVVCpFKpUqqmzBu2UkjQXqI2J8Znk2EBExL6vMfOC5iFiUWd4EnEu6r/1f\ngI+Bo4BjgCci4oocx/HdMmYZvlvGd8sUo7N3yzQAJ0kaIelIYDLQ9q6XpcAVmYONBd6LiMaImBMR\n1RHxxcx+z+YKdjMzK62C3TIRcUDSNcBy0n8MFkTERknT05vj3oh4StIESa8De4GrurbaZmbWHg9i\nMuuB3C3jbplidHoQk5mZ9S4OdzOzBHK4m5klkMPdzCyBHO5mZgnkcDczSyCHu1kZVVXV5PyyCrPO\n8n3uZmXU8fvZfZ+7fcr3uZuZVRiHu5lZAjnczcwSyOFuZpZADnczswRyuJuZJZDD3cwsgRzuZmYJ\n5HA3M0sgh7uZWQI53M3MEsjhbmaWQA53M7MEcribWY9zxBHknApZEtXVVeWuXq/Qr9wVMDNra/9+\neO653NvGjWvs3sr0Um65m5klkMPdzCyBHO5mZgnkcDczSyCHu5lZAhUV7pLGS9okabOkujxlbpe0\nRdI6SWMy64ZJelbSq5LWS/pJKStvZma5FQx3SX2AO4CLgFHAFEmntClzMXBiRJwMTAfmZzZ9DFwX\nEaOAfwKubruvmZmVXjEt97OALRGxNSL2A48Ak9qUmQQ8CBARq4ABkgZHxK6IWJdZvwfYCAwtWe3N\nzCynYsJ9KLA9a3kHhwZ02zI725aRVAOMAVZ1tJJmvVlVVU3e0ZZmXaVbRqhK+hzwGDAj04I3qxiN\njVuByLPVAW9do5hw3wlUZy0Py6xrW2Z4rjKS+pEO9oURsaS9A9XX17c8r62tpba2tojqmZlVhlQq\nRSqVKqqsIvK1KDIFpL7Aa8D5wNvAamBKRGzMKjMBuDoivi5pLPBvETE2s+1B4N2IuK7AcaJQXcx6\no3T3S3st91zbOrq+NPsEQkSJj3N4r5V/bhlwVqRJIiJy/vevYMs9Ig5IugZYTrqPfkFEbJQ0Pb05\n7o2IpyRNkPQ6sBe4MnPgLwOXAeslrSX9Ls6JiKdLcmZmZpZTwZZ7d3HL3ZLKLXe33LtKey13j1A1\nM0sgh7uZWQI53M3MEsjhbmaWQA53M7MEcribmSWQw93MLIEc7mbWqxxxBDknYauurip31XqUbpk4\nzMysVPbvJ+cAp3HjGru/Mj2YW+5mZgnkcDczSyCHu1mJ5PtSDrNycJ+7WYnk/1IOB7x1P7fczcwS\nyOFuZpZADnczswRyuJuZJZDD3cwsgRzuZmYJ5HA3M0sgh7uZWQI53M0sEfLNFlmpM0Z6hKpZB1RV\n1WRGolpPk2+2SKjMGSMd7mYdkH+KAfA0A9aTuFvGzCyBHO5mZgnkcDczSyCHu1kOnpvderuiwl3S\neEmbJG2WVJenzO2StkhaJ2lMR/Y162k+/eC07cN6o0r8Uu2C4S6pD3AHcBEwCpgi6ZQ2ZS4GToyI\nk4HpwPxi9y0klUp1pHjiVPr5g68BpMpdgbJbt65z+x+8TbLtY/v23nGL5OH8DhTTcj8L2BIRWyNi\nP/AIMKlNmUnAgwARsQoYIGlwkfu2q9J/sSv9/KHrrkG+rpee1/2SKncFyq6z4d7bdVW4DwW2Zy3v\nyKwrpkwx+5qVRf6uF3e/VIokj2rtqkFMPa3pYwnX3sjRPn0+yyef/L9urpH1Bu2Nar3wwsac/4sb\nPnww27bt6uKadZ4i2m+lSBoL1EfE+MzybCAiYl5WmfnAcxGxKLO8CTgX+EKhfbNew80lM7MOioic\njeliWu4NwEmSRgBvA5OBKW3KLAWuBhZl/hi8FxGNkt4tYt92K2hmZh1XMNwj4oCka4DlpPvoF0TE\nRknT05vj3oh4StIESa8De4Gr2tu3y87GzMyAIrplzMys9+mxI1QrcfCTpAWSGiW9krXuWEnLJb0m\n6RlJA8pZx64kaZikZyW9Kmm9pJ9k1lfSNegvaZWktZlrMDezvmKuAaTHyEhaI2lpZrlizl/SW5L+\nT+ZnYHVmXYfPv0eGeykGP/VS95M+52yzgf+OiC8BzwI3dHutus/HwHURMQr4J+DqzPteMdcgIj4E\nxkXEacAY4GJJZ1FB1yBjBrAha7mSzv8ToDYiTouIszLrOnz+PTLcKcHgp94oIl4EmtusngQ8kHn+\nAPDNbq1UN4qIXRGxLvN8D7ARGEYFXQOAiDh432Z/0p+LBRV0DSQNAyYA92WtrpjzJ30redts7vD5\n99Rw9+CnT50QEY2QDj/ghDLXp1tIqiHdcv0LMLiSrkGmS2ItsAv4r4hooLKuwW3ALFqPJquk8w/g\nvyQ1SJqWWdfh8/c3MfU+if8EXNLngMeAGRGxJ8cYiERfg4j4BDhN0j8A/ylpFIeecyKvgaSvA40R\nsU5SbTtFE3n+GV+OiLclfR5YLuk1DuP976kt951AddbysMy6StSYmacHSVXA/y1zfbqUpH6kg31h\nRCzJrK6oa3BQRLxPemKZ8VTONfgyMFHSG8DvgfMkLQR2Vcj5ExFvZ/59B/gD6W7qDr//PTXcWwZO\nSTqS9OCnpWWuU3cRradvWApcmXk+FVjSdoeE+XdgQ0T8JmtdxVwDSccfvBNC0lHABaQ/e6iIaxAR\ncyKiOiK+SPr3/tmIuBxYRgWcv6TPZv7niqSjgQuB9RzG+99j73OXNB74DZ8OfrqlzFXqcpL+A6gF\nBgGNwFzSf7kfBYYDW4HvRsR75apjV5L0ZeAF0j/MB2fwmgOsBhZTGdfgVNIfmPXJPBZFxM2SjqNC\nrsFBks4FZkbExEo5f0lfAP6T9M9+P+DhiLjlcM6/x4a7mZkdvp7aLWNmZp3gcDczSyCHu5lZAjnc\nzcwSyOFuZpZADnczswRyuJuZJZDD3cwsgf4/FKcgAsWZif4AAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "interact(binomialHiLite, n=widgets.IntSlider(min=5, max=300, step=1, value=50),\\\n", " p=widgets.FloatSlider(min=0.001, max=1, step=0.001, value=p_lower_bound),\\\n", " alpha=widgets.FloatSlider(min=0.0, max=1.0, step=0.001, value=.05)\\\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Upper confidence bounds for $p$\n", "For completeness, here's the code to find an upper confidence bound for Binomial $p$:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def binoUpperCL(n, x, cl = 0.975, inc=0.000001, p = None):\n", " \"Upper confidence level cl confidence interval for Binomial p, for x successes in n trials\"\n", " if p is None:\n", " p = float(x)/float(n)\n", " hi = 1.0\n", " if (x < n):\n", " f = lambda q: scipy.stats.binom.cdf(x, n, q) - (1-cl)\n", " hi = bisect(p, 1.0, inc, f)\n", " return hi" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What's next?\n", "We will consider some methods for constructing confidence bounds that work where the normal approximation failed, focusing at first on (lower) one-sided confidence intervals for the mean of bounded populations\n", "\n", "- [Next: Binomial confidence intervals](binom.ipynb)\n", "- [Previous: Confidence intervals based on the normal approximation](normApprox.ipynb)\n", "- [Index](index.ipynb)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.10.4" } }, "nbformat": 4, "nbformat_minor": 4 }