duality.ipynb

{
 "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",
    "<div class=\"admonition\">\n",
    "<p class=\"title\">Definition</p>\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",
    "</div>\n",
    "\n",
    "\n",
    "<div class=\"admonition\">\n",
    "<p class=\"title\">Definition</p>\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",
    "</div>\n",
    "\n",
    "\n",
    "<div class=\"admonition\">\n",
    "<p class=\"title\">Proposition</p>\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",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"admonition\">\n",
    "<p class=\"title\">Proof</p>\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",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## We will use this approach&mdash;inverting tests&mdash;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": [
       "<function __main__.binomialHiLite>"
      ]
     },
     "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": [
       "<matplotlib.figure.Figure at 0x113b16450>"
      ]
     },
     "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": [
       "<function __main__.binomialHiLite>"
      ]
     },
     "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": [
       "<matplotlib.figure.Figure at 0x115ff4390>"
      ]
     },
     "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
}