markov.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:9bcb143abf75bc33cf703548a89e50497eebd44e8b9ff11d378329d4648611c1"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<p class=\"title\">Lower confidence bounds for the mean of a nonnegative population: Markov's Inequality and methods based on the Empirical Distribution</p>\n",
      "\n",
      "Recall the set-up. We have:\n",
      "\n",
      "1. known numbers $N > 0$,  $ \\{\\ell_j\\}_{j=1}^N$, and $\\{u_j\\}_{j=1}^N$, with $\\ell_j \\le u_j, \\; \\forall j$\n",
      "2. an unknown population of numbers $\\{x_j\\}_{j=1}^N$ that satisfy $\\ell_j \\le x_j \\le u_j, \\; \\forall j$\n",
      "\n",
      "We want to draw inferences about \n",
      "\n",
      "$$\\mu \\equiv \\frac{1}{N} \\sum_{j=1}^N x_j.$$\n",
      "\n",
      "In this section, $\\ell_j = 0$ and $u_j = \\infty$, $\\forall j$ (the population is known to be nonnegative).\n",
      "\n",
      "We seek a _lower_ confidence bound for $\\mu$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Markov's Inequality:\n",
      "\n",
      "<div class=\"theorem\">\n",
      "   If $\\mathbb P(X \\ge 0) = 1$, then $\\mathbb P (X \\ge a) \\le \\mathbb E(x)/a$.\n",
      "</div>\n",
      "\n",
      "### Confidence bounds from Markov's inequality\n",
      "\n",
      "Suppose we draw a random sample (with or without replacement) of size $n$.  Let $X_j$ be the value of the $j$th draw, and let \n",
      "\n",
      "$$\\bar{X} \\equiv \\frac{1}{n} \\sum_{j=1}^n X_j.$$\n",
      "\n",
      "Then:\n",
      "+ $\\mathbb P(\\bar{X} \\ge 0) = 1$\n",
      "+ $\\mathbb E(\\bar{X}) = \\mu \\ge 0$\n",
      "+ by Markov, $\\mathbb P (\\bar{X} \\ge a) \\le \\mu/a$\n",
      "\n",
      "Hence, for $\\alpha \\in (0, 1)$,\n",
      "$$ \n",
      "  \\mathbb P (\\alpha \\bar{X} \\ge \\mu) \\le \\alpha;\n",
      "$$\n",
      "i.e., $[\\alpha \\bar{X}, \\infty)$ is a lower 1-sided $1-\\alpha$ confidence bound for $\\mu$.\n",
      "\n",
      "A lower confidence bound is possible because the support of $\\mathbb P$ is bounded below; no upper confidence bound for the mean is possible without additional assumptions.\n",
      "\n",
      "Note that this bound does not involve $n$ at all!  \n",
      "It only uses the fact that $\\mathbb P (\\bar{X} \\ge 0) = 1$.\n",
      "Clearly it has a great deal of slack, getting worse as $n$ grows."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## The Kolmogorov statistic \n",
      "\n",
      "We now take a completely different approach, based on the empirical cumulative distribution function (ecdf).\n",
      "\n",
      "We draw a uniform sample of size $n$ with replacement from the population $\\{x_j\\}_{j=1}^N$.\n",
      "\n",
      "Let $X_i$ be the $i$th item drawn, $1 \\le i \\le n$.\n",
      "For uniform sampling with replacement, $\\{X_i \\}$ are iid.\n",
      "\n",
      "Define $1_A = \\{1, \\mbox{ if $A$}; 0, \\mbox{ otherwise}\\}$.\n",
      "The theoretical cdf of $X_i$ is\n",
      "$$\n",
      "\tF(x) \\equiv \\mathbb P \\{ X_i \\le x \\} \n",
      "             = \\frac{1}{N}\\sum_{j=1}^N 1_{x \\ge x_j}.\n",
      "$$\n",
      "\n",
      "Define the empirical cumulative distribution function\n",
      "$$\n",
      "    \\hat{F}_n (x) \\equiv \\frac{1}{n} \\sum_{i=1}^n 1_{x \\ge X_i}\n",
      "$$\n",
      "\n",
      "We shall use the ecdf to construct a confidence set for $F$; from that confidence set we can \n",
      "construct a confidence bound for $\\mathbb E X_1$.\n",
      "\n",
      "Consider the one-sided Kolmogorov-Smirnov statistics\n",
      "$$\n",
      "\tD_n^+ = \\sup_x (\\hat{F}_n(x) - F(x) )\n",
      "$$\n",
      "and\n",
      "$$\n",
      "    D_n^- = \\sup_x (F(x) - \\hat{F}_n(x)).\n",
      "$$\n",
      "Kolmogorov (1933) and Smirnov (1944) seem to have been the first to study these\n",
      "statistics, which have the same distribution&mdash;a distribution that does not depend on $F$ if $F$ is continuous.\n",
      "\n",
      "Here's a simulation: we take a random sample from a Uniform distribution and plot the ecdf and the true cdf. The ecdf is the stair-step function. In places it is above the true cdf, and in places it is below. $D^-$ measures how far $F$ ever gets above $\\hat{F}_n$. Note that as $n$ increases, the maximum distance between $F$ and $\\hat{F}_n$ tends to be smaller."
     ]
    },
    {
     "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": 125
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def ecdf(x):\n",
      "    '''\n",
      "       calculates the empirical cdf of data x\n",
      "       returns the unique values of x in ascending order and the cumulative probabity at those values\n",
      "       NOTE: This is not an efficient algorithm: it is O(n^2), where n is the length of x. \n",
      "       A better algorithm would rely on the Collections package or something similar and could work\n",
      "       in O(n log n)\n",
      "    '''\n",
      "    theVals = sorted(np.unique(x))\n",
      "    theProbs = np.array([sum(x <= v) for v in theVals])/float(len(x))\n",
      "    if (theVals[0] > 0.0):\n",
      "        theVals = np.append(0., theVals)\n",
      "        theProbs = np.append(0., theProbs)\n",
      "    return theVals, theProbs\n",
      "\n",
      "\n",
      "def plotUniformKolmogorov(n):\n",
      "    sam = np.random.uniform(size=n)\n",
      "    v, pr = ecdf(sam)\n",
      "    fig, ax = plt.subplots(nrows=1, ncols=1)\n",
      "    ax.plot([0, 1], [0, 1], linestyle='--', color='r', label='true cdf')\n",
      "    ax.step(v, pr, color='b', where='post', label='ecdf')\n",
      "    ax.legend(loc='best')\n",
      "    dLoc = np.argmax(v-pr)\n",
      "    dMinus = (v-pr)[dLoc]\n",
      "    ax.axvline(x=v[dLoc], ymin=pr[dLoc], ymax=v[dLoc], color='g', linewidth='4')\n",
      "    ax.text(0.5, 0.1, r'$D_n^-=$' + str(round(dMinus, 3)), color='g', weight='heavy')\n",
      "\n",
      "interact(plotUniformKolmogorov, n=widgets.IntSliderWidget(min=3, max=300, step=1, value=10))\n",
      "    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 130,
       "text": [
        "<function __main__.plotUniformKolmogorov>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtUVNe9B/Dv4CuSUAQHNSQ6BuUGiAGJIHpFMr4AL8Wm\nkWXBBL2aB5I2xPpIq2mvQJrEtlYxaCjk6Y2YdEWTiGYpBb0DEsNDHsFnEh9IY3QFGCRUg8Pj3D8m\nMzIC82AeZx7fz1qz5Hi25/zYC39u9/mdvSWCIAggIiKn4iZ2AEREZHlM7kRETojJnYjICTG5ExE5\nISZ3IiInxOROROSEDCb3lStXYuzYsXj44YcHbLNhwwb4+flh2rRpOHfunEUDJCIi0xlM7itWrMDh\nw4cHPF9ZWYljx47hxIkTWLduHdatW2fRAImIyHQGk/vs2bPh5eU14PmKigokJCTA29sbSUlJOHv2\nrEUDJCIi05k9515ZWYmgoCDtsY+PDy5cuGDuZYmIyAxmJ3dBEHDnCgYSicTcyxIRkRmGmnuBiIgI\nnDlzBjExMQCApqYm+Pn59Wk3efJkjuiJiEw0adIknD9/3uQ/Z/bIPSIiAvv27UNLSwv27NmDwMDA\nfttduHBBO8p39c+mTZtEj8FePuwL9gX7otenpgZCcDCEuDgIV64AEAY9KDY4ck9KSkJJSQmam5sx\nfvx4ZGRkoLOzEwCQkpKC6dOnIzIyEmFhYfD29sbu3bsHFQgRkUvbtg147TVgyxYgORkwc3rbYHL/\n4IMPDF5k8+bN2Lx5s1mBEBG5tPBwoK4O8PW1yOXMnnMn08nlcrFDsBvsi9vYF7e5ZF9ERlr0chJB\nEGyyWYdEIoGNbkVE5BTUMzODy52ij9y9vb3R2toqdhhOz8vLC0qlUuwwiFybSgW88gowahTw299a\n9VaiJ/fW1laO6G2A7x4Qiay2Fvjv/wbGjwfy8qx+O64KSURkTSoVsGkTEBMDrF0LHDhgsYem+og+\ncicicmqrVwONjRathDGG6A9U+aDVNtjPRCJpbwfuuWdQdevmPFDltIyDSk9PR3Jysvb46NGjiIqK\ngoeHB7788ksRIyMiHR4eZr+QNBhM7npMnDgRR48eFTuMft35gPTVV1/FypUr0d7ejpCQEJGiInJh\nKhXQ0iJ2FFpM7noYmsro6uqyYTQDEwQBn3/+OWbNmiV2KESuqbZW/YbpG2+IHYkWk/sAkpOT0djY\niPj4eHh4eGDLli1oaGiAm5sbPvroI0yZMgULFixASUkJxo8fr/NnJ06ciCNHjmiPDx06hEWLFuHB\nBx/Etm3b8O9//3vA+1ZXVyM1NRVjxozB5MmTUVhYCEC92uaLL76IcePGYfHixdprqFQqeHh44Nat\nWwgNDYW/v78VeoOI+nVnJcwf/iB2RFpM7gN4//33MWHCBBw8eBDt7e062wfu2bMHBQUFOHz48IAP\niTXTJgUFBfj973+PjRs3QqFQoLy8HK+99lq/92xqaoJcLkdoaCgaGhpQWlqKiRMnAgBSU1Nx9epV\n1NXVYdGiRcjJyYFEIsHw4cO1ib6+vh7ffPONhXuCiPqlGa1XV6srYZYtE2VufSD2n9zT09Udducn\nPd349gO1HaQ1a9bAz88PI0aMMNj2H//4B373u99hxowZuPfee7FhwwZ8+umn/bbdu3cv5s2bh2ef\nfRbu7u7w9fXFgw8+iK6uLhQXFyMjIwPjxo3D8uXL8cgjj1j0eyIiE/3f/9m0bt1U9l/nnp5uWnI2\ntf0gREREGN22uLgYBw8exK9//Wvt76lUKnz//fcYM2aMTluFQtHvvPnZs2fR09OjswnKI488ghY7\nenhD5HLWrBE7Ar3sf+QuoiFDhvQ77TJ06O1/E++77z4olUp0d3cDAJqbm/Htt99qz8+dOxdvvvkm\nWltbtZ8bN270SewAMGfOHJSVlfX5/YCAALi5ueks2l9dXc0lBYhoQEzuekybNg3V1dV62/j7+0Mq\nleLdd99FU1MTNm3apJN0k5OT8Ze//AVlZWXo7u5GU1MTCgoK+r1WQkICFAoF3n77bdy4cQNXrlzB\nV199hWHDhmH+/PnIyMjAtWvXsHv3btTV1Vn0eyWiAdTWqqdgHAyTux6rVq3CwYMH4e3tja1btwLo\nfwGunJwcvPPOO5g+fTqCg4Nx//33a88tXLgQmZmZ2LFjB3x8fDBz5kxUVlb2ez+pVIojR46gvLwc\nMpkMc+bMQWNjIwDgjTfewJgxYzB16lR88sknSE1N1fmzHMUTWVjvShgHnALl8gMugv1MZII7V3AU\n6YEplx8gIrKUHTtsvoIjAHh79y308/Ia/PU4cncR7GciI1VXA/fea/PRukQC9PdXdLB/d5ncXQT7\nmci+WTq5c1qGiMgJ2f9LTERElqbZy9TNTV0R44Q4cici19J7TZhnnhE7Gqthcici1yDSXqZi4bQM\nEbmGl14Czp61+V6mYuHI3crkcjnefvtt7fHOnTsREBAAXxf44SKyK5mZTj9a740jdyvrvba7SqXC\nhg0b8OWXX+KBBx4QOTIiFzNypNgR2BRH7jZUU1MDDw8PJnYia1KpgGvXxI5CdEzuerS2tiIrKwsP\nPfQQFi5ciH/+858A1HuWHjhwAI899hhGjRqFsLAw7TK/p0+fxpIlS+Dr64uXXnpJ2/7IkSOYO3cu\nrl27Bg8PD6xcuVK074vIaWkqYV5/XexIRMfkrsdTTz2FS5cu4ejRo9i4cSNWrFiB8+fP49NPP8Xa\ntWvx9NNPQ6lU4p133oG7uzsEQcDcuXMRERGB+vp6dHR04Pjx45BIJJg3bx4OHToEX19ftLe34513\n3hH72yNyHndWwrzyitgRic7u59wttZKtqW/vtre3o7y8HPn5+Rg5ciTGjh2LJUuW4JNPPkF1dTWe\ne+45/PznPwcABAcHAwAqKysxYsQIrF27FgDwpz/9CTt37uwVA1//J7K43is4ukgljDHsPrmLlQ/L\nysrQ1NSkU9XS09ODRx99FBUVFdoE3ltFRQVCQkK0xyNHjkRAQIBN4iVyWdXV6tF6crJdbVAtNrtP\n7mKZOXMmfHx80NDQgOHDh+ucS0xMRFlZGcLDw3V+PyIiAn/961+1xz/++CPOnTtnk3iJXNbTT4sd\ngV3inPsARo0ahcjISGzcuBGXL19Gd3c3Tp06hRMnTiAxMRG5ubk4dOgQurq6UF9fD6VSibCwMNy6\ndQvbtm1DU1MT/vjHP3IqhohEweSux9///nfIZDIkJCTAx8cHzz77LH744QcsWrQIf/7zn7Fjxw6M\nHj0azzzzDDo6OuDm5obi4mJ8/vnnCAkJwYgRIzBr1iyda3I7PKJBqq0FDh4UOwqHYXA999LSUqSk\npKCrqwtpaWl4/vnndc7/+OOPWLVqFerr6/Gzn/0Ma9aswS9+8Yu+N+J67qJiP5PD0qzgmJMDZGcD\nv/qV2BFZhaXXczc45/7CCy8gNzcXMpkMMTExSEpKglQq1Z7ftWsX7r77btTW1uLy5cuYO3cuFi1a\nxBEqEZmPlTCDpndapq2tDQAQFRUFmUyG6OhoVFRU6LTx9PREe3s7Ojs7oVQq4e7uzsRORObLy3OZ\nFRytQW9yr6qq0inlCwoKQnl5uU6bpKQkdHd3QyqVIjIyEvn5+daJlIhcS2SkerS+bBlLHAfB7FLI\nHTt2YOjQobh69SpOnjyJuLg4XL58GW5uff/dSE9P134tl8shl8vNvT0R2RFvb6C11VJXC7LUhRyC\nl5f6V4VCAYVCYfb19D5QbWtrg1wuR21tLQDg+eefR2xsLOLi4rRtlixZgqeeegoxMTEA1LXeu3bt\n6vPyDh+oiov9TLYw0ENBgwSBo/MBWGWDbE9PTwDqipmGhgYUFRUhIiJCp828efNw4MAB9PT04OLF\ni1AqlXwrk4iMo1kTpp83vsk8BqdlsrKykJKSgs7OTqSlpUEqlSI3NxcAkJKSgsTERJw5cwZhYWHw\n8fHB9u3bTQrAy8uLD2BtwEvzfz4ie9G7EiYvT+xonI7BOneL3YjTAkROz6hpmd5161u2cE0YA6xW\n505EZFGvvqpe7It161bFkTsRWYzRI/dhwzhaNxJH7kTkGO5YZZWsgwuHEZF1qFRAY6PYUbgsJnci\nsjzNXqZZWWJH4rKY3InIsnrvZfq3v4kdjcvinDsRWUZtLYBQVsLYCSZ3IrKMb74BEKpewZGVMKJj\nKSQRWcyg15ahAVllbRkiInJMTO5EZJraWuDDD8WOggxgcici42hWcIyJAXp6xI6GDOADVSIyjHuZ\nOhyO3IlIv/fe416mDojVMkSk38WLwF13GZXUWS1jeYPNnUzuRGQxTO6Wx1JIIjIfM7PTYHInotuV\nMM88I3YkZCFM7kSuTrOCY3U1kJnZbxNvb/WUi6EPt+q1H5xzJ3JVJuxlyrl08XAnJiIyTXY2V3B0\nYhy5E7mqri5gyBCjVnDkyF08HLkTkWmG8q+/M+MDVSJnp1L9tNY6uRImdyJnpqmE2bZN7EjIxpjc\niZyRSgX8z/+o14RZtw7YuVPsiMjGOOlG5Gxqa4HlywGZjJUwLozVMkTO5tAhoLkZePJJi+1lymoZ\n8XDhMCKyGiZ38XDhMCIi0mJyJ3JUNTXAW2+JHQXZKSZ3IkejqYSJjQVGjhQ7GrJTrJYhciQ1Neq9\nTFkJQwZw5E7kKPLz1aP19euBggImdtKL1TJEjuK779S/ipDUWS0jHqtVy5SWliIwMBD+/v7Izs7u\nt01VVRXCw8MRGBgIuVxuchBEZARfX47WyWgGR+6hoaHYvn07ZDIZYmJiUFZWBqlUqj0vCAKCg4Ox\nbds2zJ8/H83NzTrntTfiyJ3IeD09gJv9zJpy5C4eq4zc29raAABRUVGQyWSIjo5GRUWFTpsTJ04g\nODgY8+fPB4B+EzsRGUlTCbN0qdiRkIPTm9yrqqoQEBCgPQ4KCkJ5eblOm8LCQkgkEsyePRvx8fEo\nLCy0TqREzq6mBggLU/+6davY0ZCDM7sUsqOjA3V1dSguLsbNmzexYMECnDp1CiNZf0tkHBP2MiUy\nlt7kHh4ejvXr12uPT58+jdjYWJ02M2fOxK1btzBu3DgAQFhYGEpLSxETE9Pneunp6dqv5XI5H76S\nS/D2Blpb9bUYDiBD/VkO9cfOeHmJHYHrUCgUUCgUZl/H6AeqEyZMQGxsbJ8Hqi0tLVi4cCEUCgU6\nOjowY8YM1NTU4J577tG9ER+okosy+DCyp0fdiKN16ofV9lDNyspCSkoKOjs7kZaWBqlUitzcXABA\nSkoKRo8ejRUrViAsLAw+Pj7IzMzsk9iJSA87qooh58GXmIisTDty1+xl+tBDYodEDoRL/hLZs7o6\nYPp0VsGQzXDkTmRNKhUkI4ZD8BnDShgaFKvNuRPRINXXA8uWAajjCo5kc5yWIbKWtjZgzRr110zs\nZGOcliGyMq7LQubgA1UiItJicicyV10dkJUldhREOpjciQZLpQI2bQKio4HRo8WOhkgHq2WIBqOu\nTr2X6f33sxKG7BJH7kSm2rdPPVpfswY4cICJnewSq2WITNXSAty6ZXRSZ7UMmWOwuZPJncjKmNzJ\nHHxDlchCJBm6ywMIm5iZyfFwzp2IyAkxuRMROSFOyxBpaPYy5ZCHnAB/jIk0PvoIqK4WOwoii2C1\nDJHGTz+fkkzdMY+5D1RZLUPm4MJhROYaxCbV3t63/9hAHy8vK8VLpAeTO7kelQqoqbHIpVpb1aNy\nfR+l0iK3IjIJkzu5Fu5lSi6CyZ1cQ+8VHNesAd5/X+yIiKyKpZDk/E6eVG9MzRUcyYUwuZPz6+5W\nj9aTk01+YErkqJjcyflNnar+ELkQzrkTETkhJndyHnV1wMsvix0FkV1gcifH17sSRiYTOxoiu8A5\nd3Js3MuUqF8cuZPj+uwz7mVKNACO3MlxRUVxtE40ACZ3clweHuoPEfXBaRlyDJ2dYkdA5FCY3Mm+\naSph5s/nouhEJmByJ/ulWcGxuhr44AMuHUBkAiZ3sj93ruDIShgikxlM7qWlpQgMDIS/vz+ys7MH\nbFdVVYWhQ4fi448/tmiA5IIKC9Wj9bo6YNkyjtiJBsHgHqqhoaHYvn07ZDIZYmJiUFZWBqlUqtOm\nu7sbCxYsgLu7O1asWIHFixf3vRH3UCVjaX5ORErqkgzd++rbQ5X7o5K1WWUP1ba2NgBAVFQUZDIZ\noqOjUVFR0adddnY2EhIS4OPjY3IARH0MYi9TItKlN7lXVVUhICBAexwUFITy8nKdNleuXMH+/fuR\nmpoKQP2vDJFRVCrg+HGxoyBySmY/UF29ejU2b96s/a8Dp17IKJpKmG3bOK9BZAV631ANDw/H+vXr\ntcenT59GbGysTpvq6mokJiYCAJqbm3Ho0CEMGzYMixYt6nO99PR07ddyuRxyudyM0MkhqVTAK68A\nOTnAli0OsTuSvvC8vGwXB7kGhUIBhUJh9nWMfqA6YcIExMbG9vtAVWPFihWIj4/H448/3vdGfKBK\nZ84AS5eqV3DMy7Pb8kZTHqgSWdtgc6fBtWWysrKQkpKCzs5OpKWlQSqVIjc3FwCQkpJieqTkuoYP\nd4i9TDXJnJUw5MgMjtwtdiOO3MnBMLmTPbBKKSQRETkmJneyvLo6YP16DnuJRMTkTpbTe02Yhx8W\nOxoil8bNOsgyuJcpkV3hyJ3Md+QIV3AksjOsliHz3boFtLQ4XVJntQzZg8HmTiZ3ogEwuZM9YCkk\n2UZHh9gREJERmNzJOJpKmMhIDmeJHACTOxnWey/TggK7XjqAiNSY3Glg3MuUyGGxzp0G9sUXQE0N\n69aJHBCrZYgGwGoZsgesliEiIi0md1LPrR85InYURGRBTO6uTlMJs2MH0NMjdjREZCFM7q7qzkqY\njz8G3PjjQOQsWC3jis6dAxITuYIjkRNjtYwr+u479Rz7k09a5YUkb2+gtdXil7U5Ly9AqRQ7CnJ1\nXDiM7AZLCIksh6WQRESkxeTuzOrqgFWrWAVD5IKY3J1R70qY//xPLvRF5IJYLeNsuJcpEYEjd+dy\n/DhXcCQiAKyWcS7d3UBTEzBunKhhsFqGyHJYCkl2g8mdyHJYCulqbtwQOwIismNM7o5GUwkzfbp6\nGoaIqB9M7o6k916mRUXAkCFiR0REdorJ3RFwL1MiMhHr3B3ByZPqUTvr1onISKyWIYtjtQyR5bBa\nhoiItJjc7YlKBRw8KHYUROQEmNzthaYSJi8P6OoSOxoicnBGJffS0lIEBgbC398f2dnZfc7n5+cj\nJCQEISEhWLp0Kb7++muLB+q07qyE2b8fGMrn3ERkHqMeqIaGhmL79u2QyWSIiYlBWVkZpFKp9vwX\nX3yBoKAgeHp6YteuXSguLsb777+veyM+UO3r/HkgIUG9gmNentNUwvCBKpHlWO2BaltbGwAgKioK\nMpkM0dHRqKio0Gkzc+ZMeHp6AgDi4uJQUlJiciAuafRo4MUXWbdORBZnMLlXVVUhICBAexwUFITy\n8vIB2+fl5SE+Pt4y0Tk7Ly9g6VJupkFEFmfRyd3i4mLs3r0bx48f7/d8enq69mu5XA65XG7J25MB\n3t5Aa6v17+PlZf17EDkrhUIBhUJh9nUMzrm3tbVBLpejtrYWAPD8888jNjYWcXFxOu3q6+vx+OOP\n4/Dhw5g8eXLfG7nynHtdHbBlC/Duu8CwYaKFwblwIsdjtTl3zVx6aWkpGhoaUFRUhIiICJ02jY2N\nWLx4MfLz8/tN7C6rdyVMdDSrYIjIZozKNllZWUhJSUFnZyfS0tIglUqRm5sLAEhJSUFmZiaUSiVW\nrVoFABg2bBgqKyutF7Uj4F6mRCQiri1jDbW1QEyMeiomOdluHphyWobI8XCbPXsiCEBzM+DjI3Yk\nOpjciRwPkzsZxORO5Hi4KqRYfnrJi4jInjC5D5amEuaRR9RfExHZESb3waitBcLD1XuZHjsGDB8u\ndkRERDqY3E2hGa3HxADr1nFNGCKyW3yrxhQXLgCnTrFunYjsHqtlXAirZYgcD6tliIhIi8m9PyoV\n8NFHYkdBRDRoTO530lTC/O//ArduiR0N2djWL7Zi/LbxeKvmLeRU5eDFohfRdKMJAPBd+3coPF+o\n/Xzxry/Mvl/dtTrMeGsGPDd7Yva7s/FV81cDtk07lIaxW8bCLcMN8R/EG3Xuvbr34Jbh1ufT2NZo\nduxk3/hAVUOlAl55BcjJAf72N+DJJ+1mTRiynen3TUfMpBg8/cjTAICr7VfxwuEXsGfxHvh6+MLX\nw7IP0pd8tARtt9qwOmI1siuzsezTZah4uqLfthJIkDQlCa9XvA4JJEadk0+U48OEDwEAnd2deKrg\nKXiP9MZ9HvdZ9Psg+8PkDgCXLgGPPQZMmMBKGBdX8W0Fpt83XXt8r8e9ON102ir3+lfbv3BeeR5/\niPoDMuZk4HrHdeys2om2jjZ43uXZp/32hdtx+fplvF7xutHnJo6aiImjJgIA9p7ZC1W3CitDV2KI\n2xCrfE9kPzgtAwBjxgAbNwIFBUzsLq7quyqd5C4IAto6rLPExNFLRwEAD4x6QP2r1wPoEXq0v98f\nAQNXTeg7BwC51bkY4jYEz057dhDRkqPhyB0A7r4b+NWvxI6C7MDJ708ieGyw9rj6ajWCfIJMvs6Z\npjMoulDU77nlU5dj1F2jtMeapKwpd5NYYTrwgvICjlw8grj/iMMEzwkWvz7ZHyZ3op803WiC111e\ncJPc/g/th6c+xG9n/NbkawX5BBn8R2Ge3zwAwMXWiwCAS9cvwU3ihjkT5wAAOro6MEQyBMOGmL81\nY261enOd1LBUs69FjsG1knttLfDyy8CePcBdd4kdDdmZyiuVCPcN1x5/ee1L3FDdwIJJC0y+lr6R\n+7KQZfAa6YX7f3Y//Ef7I686D24SN+yu341pvtO08+3ur7jjoTEP4WTqSQDAZ19/hlPfnwIANLY1\n4u2at/HoxEcx2Xuy3nOqbhXeq3sPslEy/Jf/f5n8vZBjcvo3VL29gdZWm9/WLnl5AUql2FHYp+P/\nOo6XS1/G6JGjMWfiHHQL3ejo6kBaRBoA4NjlYzh8/jASghJwpukMngh+wiL3rblag+c+ew5nm88i\neGww3ox/EwHSAACAW4YbpoyZgvrUegDAnF1zUNJQov27JJFI8O4v3sWykGV6z3146kM88fETeHnO\ny9g4e6NF4ibb4WYdA94XEIJD1JUwubl8YEqDcvn6ZRRfLEbE/RH4puUb/DLwl2KHRC6Cyw/056uf\nXghZt46VMGSWiisVmO83H9XfVePu4XeLHQ6RQa4xcudiWUTkoDhyJyIiLedJ7i0tYkdARGQ3HD+5\na3ZHCg0Fbt4UOxoiIrvg2Mlds4JjTQ1QXg64u4sdERGRXXDM5H7nXqashCEi0uGYb6hevQqcO8cV\nHImIBsBSSCIiO8ZSSCIi0rLv5K5SAbt2cehNRGQi+03umkqYvXtZ4khEZCL7S+79VcLczbU8iIhM\nYV/VMt9+C8TFcS9TIiIz2Ve1jEoFHDwI/PKX6jIXi9yXU/ZE5Li4nvuA92VyJyLHZbVSyNLSUgQG\nBsLf3x/Z2dn9ttmwYQP8/Pwwbdo0nDt3zuQgiIjIsgwm9xdeeAG5ubkoLi7Gzp070dzcrHO+srIS\nx44dw4kTJ7Bu3TqsW7fO8F1ra4GFC4Effhh04I5MoVCIHYLdYF/cxr64jX1hPr3Jva2tDQAQFRUF\nmUyG6OhoVFRU6LSpqKhAQkICvL29kZSUhLNnzw58wd6VMEuXAh4e5n8HDog/uLexL25jX9zGvjCf\n3uReVVWFgIAA7XFQUBDKy8t12lRWViIoKEh77OPjgwsXLvR/Qc0KjnV1QHKyxR6aEhGRLrNLIQVB\n6DPZLxkoaa9dy6RORGQLgh7Xr18Xpk6dqj3+zW9+Ixw8eFCnzeuvvy5s3bpVe+zn59fvtSZNmiQA\n4Icffvjhx4TPpEmT9KXpAekduXt6egJQV8xMmDABRUVF2LRpk06biIgIrFmzBsuWLUNhYSECAwP7\nvdb58+f13YqIiCzI4LRMVlYWUlJS0NnZibS0NEilUuTm5gIAUlJSMH36dERGRiIsLAze3t7YvXu3\n1YMmIiL9bPYSExER2Y7FFw7jS0+3GeqL/Px8hISEICQkBEuXLsXXX38tQpS2YczPBaCu0Bo6dCg+\n/vhjG0ZnO8b0Q1VVFcLDwxEYGAi5XG7bAG3IUF/8+OOPWL58OUJDQ/Hoo49i//79IkRpGytXrsTY\nsWPx8MMPD9jG5Lw5qJl6PaZOnSqUlJQIDQ0NwoMPPig0NTXpnK+oqBBmzZoltLS0CHv27BHi4uIs\nHYLdMNQXx48fF65fvy4IgiC89957wpNPPilGmDZhqC8EQRC6urqEOXPmCHFxccLevXtFiNL6DPVD\nT0+PMGXKFKGoqEgQBKHffnIWhvoiJydHSE1NFQRBEBoaGgQ/Pz+hp6dHjFCtrrS0VKipqRGmTJnS\n7/nB5E2Ljtwt/tKTAzOmL2bOnKl9aB0XF4eSkhKbx2kLxvQFAGRnZyMhIQE+Pj62DtEmjOmHEydO\nIDg4GPPnzwcASKVSm8dpC8b0haenJ9rb29HZ2QmlUgl3d/eBy6wd3OzZs+Hl5TXg+cHkTYsmd4u/\n9OTAjOmL3vLy8hAfH2+L0GzOmL64cuUK9u/fj9TUVAB63pVwYMb0Q2FhISQSCWbPno34+HgUFhba\nOkybMKYvkpKS0N3dDalUisjISOTn59s6TLsxmLxp8/XcBVNeenIRxcXF2L17N44fPy52KKJZvXo1\nNm/erF0B786fEVfR0dGBuro6FBcX4+bNm1iwYAFOnTqFkSNHih2aze3YsQNDhw7F1atXcfLkScTF\nxeHy5ctwc7O/PYasbTB506K9FB4erjPRf/r0acyYMUOnTUREBM6cOaM9bmpqgp+fnyXDsAvG9AUA\n1NfXY9WqVSgoKMCoUaNsGaLNGNMX1dXVSExMxAMPPIB9+/bhueeeQ0FBga1DtSpj+mHmzJlYuHAh\nxo0bBz8/P4SFhaG0tNTWoVqdMX1RWlqKJ554Au7u7oiIiICvr69TFx3oM5i8adHk3vulp4aGBhQV\nFSEiIqLMusjXAAABIElEQVRPkPv27UNLSwv27Nkz4EtPjs6YvmhsbMTixYuRn5+PyZMnixGmTRjT\nFxcvXsSlS5dw6dIlJCQkICcnB4sWLRIjXKsxph9mzJiBkpIS3Lx5E0qlErW1tZg1a5YY4VqVMX0x\nb948HDhwAD09Pbh48SKUSqXOVI4rGUzetPi0DF96us1QX2RmZkKpVGLVqlUAgGHDhqGyslLMkK3G\nUF+4CkP9MHr0aKxYsQJhYWHw8fFBZmYm7rnnHpGjtg5DfZGYmIgzZ85o+2L79u0iR2w9SUlJKCkp\nQXNzM8aPH4+MjAx0dnYCGHze5EtMREROyPWeTBARuQAmdyIiJ8TkTkTkhJjciYicEJM7EZETYnIn\nInJCTO5ERE6IyZ2IyAn9Pyn3RMoLXbk+AAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10411ab10>"
       ]
      }
     ],
     "prompt_number": 130
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## The Massart-Dvoretzky-Kiefer-Wolfowitz (MDKW) Inequality and confidence set for $F$\n",
      "Dvoretsky, Kiefer and Wolfowitz (1956, DKW) showed that\n",
      "$$ \n",
      "\t\\mathbb P \\{D_n^- > t\\} \\le C \\exp(- 2 n t^2)\n",
      "$$\n",
      "for some constant $C$.\n",
      "\n",
      "Massart (1990) showed that $C = 1$ is the sharp constant in the DKW inequality, provided \n",
      "$\\exp(-2 n t^2) \\le \\frac{1}{2}$.\n",
      "The distribution of $D_n^-$ is stochastically larger when \n",
      "$F$ is continuous than when $F$ has jumps (Massart 1990); thus the inequality conservative for iid \n",
      "sampling from finite populations.\n",
      "\n",
      "Moreover, $D_n^-$ is stochastically larger for sampling with replacement than for sampling \n",
      "without replacement, so the inequality is conservative for sampling from a finite population\n",
      "without replacement as well.\n",
      "\n",
      "Let's call the inequality with $C=1$ the Massart-Dvoretsky-Kiefer-Wolfowitz (MDKW) inequality.\n",
      "It follows from the MDKW inequality that\n",
      "$$ \n",
      "    \\mathbb P \\left \\{\\sup_x (F(x) - \\hat{F}_n(x)) > \\sqrt{-\\frac{\\ln \\alpha}{2n}} \\right \\} \\le \\alpha.\n",
      "$$\n",
      "That is, \n",
      "$$\n",
      "    \\mathcal C(X) \\equiv \\left \\{ \\mbox{cdfs } G: \\sup_x (G(x) - \\hat{F}_n(x)) \\le  \\sqrt{-\\frac{\\ln \\alpha}{2n}} \\right \\}\n",
      "$$\n",
      "is a nonparametric confidence set for the true population cdf $F$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Using MDKW to get a lower confidence bound for the mean of a nonnegative population\n",
      "\n",
      "We now assume all $\\{x_j\\}_{j=1}^N$ are nonnegative, so $F(0^-) = 0$.\n",
      "The population mean is\n",
      "\n",
      "$$ \\frac{1}{N} \\sum_{j=1}^N x_j = \\int_0^\\infty x dF(x).$$\n",
      "\n",
      "But with probability at least $1-\\alpha$, $F \\in \\mathcal C(X)$.\n",
      "Hence \n",
      "\n",
      "$$\n",
      "   \\inf_{G \\in \\mathcal C(X)} \\int_{0}^\\infty x dG(x)\n",
      "$$\n",
      "is a lower $1-\\alpha$ confidence bound for $\\mu = \\int_0^\\infty x dF(x)$.\n",
      "\n",
      "### Among all $G \\in \\mathcal C(X)$, which has smallest mean?\n",
      "It's the distribution function that has as much mass as possible as far \"left\" as possible, \n",
      "but never exceeds $\\hat{F}_n$ by more than $\\sqrt{-\\frac{\\ln \\alpha}{2n}}$.\n",
      "\n",
      "Let $G^-$ denote this cdf.\n",
      "\n",
      "Then, \n",
      "\n",
      "$$\n",
      "   G^-(x) = \\left ( \\hat{F}_n(x) +  \\sqrt{-\\frac{\\ln \\alpha}{2n}} \\right ) \\wedge 1,\n",
      "$$ \n",
      "\n",
      "where $a \\wedge b \\equiv \\max\\{a, b\\}$.\n",
      "\n",
      "The corresponding lower confidence bound on $\\mu$ is $\\int_0^\\infty x dG^-(x)$.\n",
      "\n",
      "The integral is really a sum because $G^-$ is a step function.\n",
      "Note that $G^-$ might have a jump at $x=0$ even if $\\hat{F}_n$ does not, but the rest of its jumps are \n",
      "at (a subset of) the data. \n",
      "\n",
      "Let's implement this in code."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# find lower confidence bound on the mean of a non-negative population using the MDKW inequality\n",
      "\n",
      "def massartOneSide(n, alpha):\n",
      "    ''' \n",
      "       tolerable misfit between the cdf and ecdf for a one-sided bound \n",
      "       at significance level alpha for an iid sample of size n\n",
      "    '''\n",
      "    return np.sqrt(-np.log(alpha)/(2.0*n))\n",
      "    \n",
      "    \n",
      "def ksLowerMean(x, c):\n",
      "    '''\n",
      "       lower confidence bound for the mean of a nonnegative population\n",
      "       x is an iid sample with replacement from the population\n",
      "       c is the Massart constant for the desired coverage\n",
      "    '''\n",
      "    # find the ecdf\n",
      "    vals, probs = ecdf(x)\n",
      "    probs = np.fmin(probs+c, 1)   # This is G^-\n",
      "    gProbs = np.diff(np.append([0.0], probs))  # pre-pend a 0 so that diff does the right thing; \n",
      "                                               # gProbs is the vector of masses\n",
      "    return (vals*gProbs).sum()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 122
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now let's compare this approach to lower one-sided Student-t intervals in a simulation using a nonstandard mixture of distributions."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Nonstandard mixture: a pointmass at 1 and a uniform[0,1]\n",
      "ns = np.array([25, 50, 100, 400])  # sample sizes\n",
      "ps = np.array([0.9, 0.99, 0.999])    # mixture fraction, weight of pointmass\n",
      "alpha = 0.05  # 1- (confidence level)\n",
      "reps = 1.0e3   # just for demonstration\n",
      "\n",
      "simTable = pd.DataFrame(columns=('mass at 1', 'sample size', 'Student-t cov',\\\n",
      "                                 'MDKW cov', 'Student-t low', 'MDKW low')\n",
      "                       )\n",
      "\n",
      "for p in ps:\n",
      "    popMean = (1-p)*0.5 + p  # population is a mixture of uniform with mass (1-p) and\n",
      "                             # pointmass at 1 with mass p\n",
      "    \n",
      "    for n in ns:\n",
      "        tCrit = sp.stats.t.ppf(q=1-alpha, df=n-1)\n",
      "        mCrit = np.sqrt(-np.log(alpha)/(2.0*n))  # the 1-sided MDKW constant\n",
      "        coverT = 0\n",
      "        coverM = 0\n",
      "        lowT = 0.0\n",
      "        lowM = 0.0\n",
      "        \n",
      "        for rep in range(int(reps)):\n",
      "            sam = np.random.uniform(size=n)  # the uniform part of the sample\n",
      "            ptMass = np.random.uniform(size=n)  # for a bunch of p-coin tosses\n",
      "            sam[ptMass < p] = 1.0   # mix in pointmass at 1, with probability p\n",
      "            # t interval\n",
      "            samMean = np.mean(sam)\n",
      "            samSD = np.std(sam, ddof=1)\n",
      "            tLo = samMean - tCrit*samSD  # lower endpoint of the t interval\n",
      "            coverT += ( tLo <= popMean )\n",
      "            lowT += tLo\n",
      "            #  MDKW interval\n",
      "            mLo = ksLowerMean(sam, mCrit) # lower endpoint of the MDKW interval\n",
      "            coverM += (mLo <= popMean )\n",
      "            lowM += mLo\n",
      "\n",
      "        simTable.loc[len(simTable)] =  p, n,\\\n",
      "            str(100*float(coverT)/float(reps)) + '%',\\\n",
      "            str(100*float(coverM)/float(reps)) + '%',\\\n",
      "            str(round(lowT/float(reps),4)),\\\n",
      "            str(round(lowM/float(reps),4))\n",
      "#\n",
      "ansStr =  '<h3>Simulated coverage probability and expected lower endpoint of ' +\\\n",
      "          'one-sided Student-t and MDKW confidence intervals for ' +\\\n",
      "          'mixture of U[0,1] and pointmass at 1 population</h3>' +\\\n",
      "          '<strong>Nominal coverage probability ' + str(100*(1-alpha)) +\\\n",
      "          '%</strong>.<br /><strong>Estimated from ' + str(int(reps)) + ' replications.</strong>'\n",
      "\n",
      "display(HTML(ansStr))\n",
      "display(simTable)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<h3>Simulated coverage probability and expected lower endpoint of one-sided Student-t and MDKW confidence intervals for mixture of U[0,1] and pointmass at 1 population</h3><strong>Nominal coverage probability 95.0%</strong>.<br /><strong>Estimated from 1000 replications.</strong>"
       ],
       "metadata": {},
       "output_type": "display_data",
       "text": [
        "<IPython.core.display.HTML at 0x10dabbc10>"
       ]
      },
      {
       "html": [
        "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
        "<table border=\"1\" class=\"dataframe\">\n",
        "  <thead>\n",
        "    <tr style=\"text-align: right;\">\n",
        "      <th></th>\n",
        "      <th>mass at 1</th>\n",
        "      <th>sample size</th>\n",
        "      <th>Student-t cov</th>\n",
        "      <th>MDKW cov</th>\n",
        "      <th>Student-t low</th>\n",
        "      <th>MDKW low</th>\n",
        "    </tr>\n",
        "  </thead>\n",
        "  <tbody>\n",
        "    <tr>\n",
        "      <th>0 </th>\n",
        "      <td>   0.9</td>\n",
        "      <td>  25</td>\n",
        "      <td>  89.7%</td>\n",
        "      <td> 100.0%</td>\n",
        "      <td> 0.6846</td>\n",
        "      <td> 0.7059</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>1 </th>\n",
        "      <td>   0.9</td>\n",
        "      <td>  50</td>\n",
        "      <td>  99.1%</td>\n",
        "      <td> 100.0%</td>\n",
        "      <td> 0.6706</td>\n",
        "      <td> 0.7775</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>2 </th>\n",
        "      <td>   0.9</td>\n",
        "      <td> 100</td>\n",
        "      <td> 100.0%</td>\n",
        "      <td> 100.0%</td>\n",
        "      <td> 0.6631</td>\n",
        "      <td> 0.8276</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>3 </th>\n",
        "      <td>   0.9</td>\n",
        "      <td> 400</td>\n",
        "      <td> 100.0%</td>\n",
        "      <td> 100.0%</td>\n",
        "      <td> 0.6615</td>\n",
        "      <td> 0.8886</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>4 </th>\n",
        "      <td>  0.99</td>\n",
        "      <td>  25</td>\n",
        "      <td>  20.7%</td>\n",
        "      <td> 100.0%</td>\n",
        "      <td> 0.9617</td>\n",
        "      <td>  0.751</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>5 </th>\n",
        "      <td>  0.99</td>\n",
        "      <td>  50</td>\n",
        "      <td>  38.4%</td>\n",
        "      <td> 100.0%</td>\n",
        "      <td> 0.9425</td>\n",
        "      <td> 0.8219</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>6 </th>\n",
        "      <td>  0.99</td>\n",
        "      <td> 100</td>\n",
        "      <td>  62.9%</td>\n",
        "      <td> 100.0%</td>\n",
        "      <td> 0.9236</td>\n",
        "      <td> 0.8723</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>7 </th>\n",
        "      <td>  0.99</td>\n",
        "      <td> 400</td>\n",
        "      <td>  98.4%</td>\n",
        "      <td> 100.0%</td>\n",
        "      <td> 0.9065</td>\n",
        "      <td> 0.9338</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>8 </th>\n",
        "      <td> 0.999</td>\n",
        "      <td>  25</td>\n",
        "      <td>   2.4%</td>\n",
        "      <td> 100.0%</td>\n",
        "      <td>  0.995</td>\n",
        "      <td> 0.7547</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>9 </th>\n",
        "      <td> 0.999</td>\n",
        "      <td>  50</td>\n",
        "      <td>   4.7%</td>\n",
        "      <td> 100.0%</td>\n",
        "      <td> 0.9926</td>\n",
        "      <td> 0.8263</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>10</th>\n",
        "      <td> 0.999</td>\n",
        "      <td> 100</td>\n",
        "      <td>  11.0%</td>\n",
        "      <td> 100.0%</td>\n",
        "      <td> 0.9897</td>\n",
        "      <td>  0.877</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>11</th>\n",
        "      <td> 0.999</td>\n",
        "      <td> 400</td>\n",
        "      <td>  32.7%</td>\n",
        "      <td> 100.0%</td>\n",
        "      <td> 0.9841</td>\n",
        "      <td> 0.9383</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
        "<p>12 rows \u00d7 6 columns</p>\n",
        "</div>"
       ],
       "metadata": {},
       "output_type": "display_data",
       "text": [
        "   mass at 1 sample size Student-t cov MDKW cov Student-t low MDKW low\n",
        "0        0.9          25         89.7%   100.0%        0.6846   0.7059\n",
        "1        0.9          50         99.1%   100.0%        0.6706   0.7775\n",
        "2        0.9         100        100.0%   100.0%        0.6631   0.8276\n",
        "3        0.9         400        100.0%   100.0%        0.6615   0.8886\n",
        "4       0.99          25         20.7%   100.0%        0.9617    0.751\n",
        "5       0.99          50         38.4%   100.0%        0.9425   0.8219\n",
        "6       0.99         100         62.9%   100.0%        0.9236   0.8723\n",
        "7       0.99         400         98.4%   100.0%        0.9065   0.9338\n",
        "8      0.999          25          2.4%   100.0%         0.995   0.7547\n",
        "9      0.999          50          4.7%   100.0%        0.9926   0.8263\n",
        "10     0.999         100         11.0%   100.0%        0.9897    0.877\n",
        "11     0.999         400         32.7%   100.0%        0.9841   0.9383\n",
        "\n",
        "[12 rows x 6 columns]"
       ]
      }
     ],
     "prompt_number": 131
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note that sometimes the average lower bound for the MDKW interval is _larger_ even though it has better coverage. To some extent, that's because we have not truncated the Student-t interval below at 0, which would be legitimate.\n",
      "\n",
      "Still, the length penalty for conservative coverage is modest, especially given that the nominal 95% Student-t interval has coverage below 2% in some of these scenarios."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## What else?\n",
      "\n",
      "The MDKW inequality is not the only way to construct a confidence set for a cdf. One can also use differences of order statistics, among other things.  I have not compared the relative efficiency of such methods."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## What's next?\n",
      "Next we consider confidence bounds for bounded populations, using common probability inequalities.\n",
      "\n",
      "+ [Next: Dollar-unit sampling and taint](dus.ipynb)\n",
      "+ [Previous: Confidence bounds from the Chebychev and Hoeffding Inequalities](hoeffding.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 0x10db1ee10>"
       ]
      }
     ],
     "prompt_number": 124
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 124
    }
   ],
   "metadata": {}
  }
 ]
}