duality.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:4f21f4cbc93aaf2f2d233a7c028035436e86b1e7abc68e5836d2fe511e3c39d8"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<p class=\"title\"> Duality between confidence sets and hypothesis tests</p>\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=\"definition\">\n",
      "$A_\\theta \\subset \\mathcal X$ is the acceptance region for a level-$\\alpha$ test of the hypothesis $\\mu = \\theta$ iff\n",
      "$$\\mathbb P_\\theta (X \\notin A_\\theta) \\le \\alpha.$$\n",
      "</div>\n",
      "\n",
      "<p class=\"gap01\"></p>\n",
      "\n",
      "<div class=\"definition\">\n",
      "$\\mathcal I(X)$ is a $1-\\alpha$ confidence set for $\\mu$ iff \n",
      "$$ \\forall \\theta \\in \\Theta, \\;\\;\\; \\mathbb P_\\theta ( \\mathcal I(X) \\ni \\theta) \\ge 1-\\alpha.$$\n",
      "</div>\n",
      "\n",
      "<p class=\"gap01\"></p>\n",
      "\n",
      "<div class=\"proposition\">\n",
      "Suppose \n",
      "$$\\{A_\\theta: \\theta \\in \\Theta \\}$$\n",
      "is a family of level-$\\alpha$ acceptance regions. Then \n",
      "$$ \\mathcal I(X) \\equiv \\{ \\theta \\in \\Theta: X \\in A_\\theta \\}$$\n",
      "is a $1-\\alpha$ confidence set for $\\mu$.\n",
      "</div>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<div class=\"proof\">\n",
      " For any $\\theta \\in \\Theta$,\n",
      "   $$ \\mathbb P_\\theta \\left (\\{ \\eta \\in \\Theta: X \\in A_\\eta \\} \\ni \\theta \\right )  = \n",
      "        \\mathbb P_\\theta ( X \\in A_\\theta ) $$\n",
      "   $$ \\ge 1-\\alpha.$$\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",
      "$$ \\mathbb P_{n,p} (X = k) = {n \\choose k} p^k (1-p)^{n-k}$$\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",
      "$$ A_\\pi \\equiv \\ [0, a_\\pi],$$\n",
      "where \n",
      "$$ a_\\pi \\equiv \\min \\left \\{ k: \\sum_{i=k+1}^n {{n}\\choose{i}} \\pi^i (1-\\pi)^{n-i} \\le \\alpha \\right \\}.$$\n",
      "\n",
      "Let's see this graphically:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 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 IPython.html.widgets import interact, interactive, fixed\n",
      "from IPython.html import widgets\n",
      "from IPython.display import clear_output, display, HTML"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "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.IntSliderWidget(min=5, max=300, step=1, value=30),\\\n",
      "         p=widgets.FloatSliderWidget(min=0.001, max=1, step=0.001, value=.5),\\\n",
      "         alpha=widgets.FloatSliderWidget(min=0.0, max=1.0, step=0.001, value=.05)\\\n",
      "        )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "<function __main__.binomialHiLite>"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/Users/stark/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/IPython/core/formatters.py:239: FormatterWarning: Exception in image/png formatter: \n",
        "\\mathbf{P}(X \\ge 20) \\approx4.94\n",
        "             ^\n",
        "Unknown symbol: \\ge (at char 13), (line:1, col:14)\n",
        "  FormatterWarning,\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "text": [
        "<matplotlib.figure.Figure at 0x10da6e350>"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "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",
      "$$ \\min \\{ \\pi: A_\\pi \\ni X \\},$$\n",
      "that is,\n",
      "$$ \\min \\left \\{ \\pi: \\sum_{i=X}^n {{n}\\choose{i}} \\pi^i (1-\\pi)^{n-i} > \\alpha \\right \\}.$$\n",
      "\n",
      "The following code implements this, using a bisection search."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "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"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "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",
     "collapsed": false,
     "input": [
      "p_lower_bound = binoLowerCL(50, 40, cl=0.95)\n",
      "print p_lower_bound"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.68440322876\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "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",
     "collapsed": false,
     "input": [
      "interact(binomialHiLite, n=widgets.IntSliderWidget(min=5, max=300, step=1, value=50),\\\n",
      "         p=widgets.FloatSliderWidget(min=0.001, max=1, step=0.001, value=p_lower_bound),\\\n",
      "         alpha=widgets.FloatSliderWidget(min=0.0, max=1.0, step=0.001, value=.05)\\\n",
      "        )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "<function __main__.binomialHiLite>"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/Users/stark/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/IPython/core/formatters.py:239: FormatterWarning: Exception in image/png formatter: \n",
        "\\mathbf{P}(X \\ge40) \\approx5.0\n",
        "             ^\n",
        "Unknown symbol: \\ge (at char 13), (line:1, col:14)\n",
        "  FormatterWarning,\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "text": [
        "<matplotlib.figure.Figure at 0x10da639d0>"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "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",
     "collapsed": false,
     "input": [
      "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"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "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",
     "collapsed": false,
     "input": [
      "%run talkTools.py"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<style>\n",
        "\n",
        "font-family: Helvetica\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",
        ".rendered_html {\n",
        "  color: #2C5494;\n",
        "  font-family: Helvetica;\n",
        "  font-size: 140%;\n",
        "  line-height: 1.1;\n",
        "  margin: 0.5em 0;\n",
        "}\n",
        "\n",
        ".title {\n",
        "  color: #498AF3;\n",
        "  font-size: 250%;\n",
        "  font-weight:bold;\n",
        "  line-height: 1.2;\n",
        "  margin: 10px 50px 10px;\n",
        "}\n",
        "\n",
        ".subtitle {\n",
        "  color: #386BBC;\n",
        "  font-size: 180%;\n",
        "  font-weight:bold;\n",
        "  line-height: 1.2;\n",
        "  margin: 20px 50px 20px;\n",
        "}\n",
        "\n",
        ".slide-header, p.slide-header {\n",
        "  color: #498AF3;\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: #498AF3;\n",
        "  line-height: 1.2;\n",
        "  margin: 0.15em 0em 0.5em;\n",
        "  page-break-before: always;\n",
        "  text-align: center;\n",
        "}\n",
        "\n",
        "\n",
        ".rendered_html h2 {\n",
        "  color: #386BBC;\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 li {\n",
        "  line-height: 1.8;\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>"
       ],
       "metadata": {},
       "output_type": "display_data",
       "text": [
        "<IPython.core.display.HTML at 0x10da9b150>"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    }
   ],
   "metadata": {}
  }
 ]
}