markov.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lower confidence bounds for the mean of a nonnegative population\n",
    "\n",
    "### Markov's Inequality and methods based on the Empirical Distribution\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",
    "\n",
    "If $\\mathbb P(X \\ge 0) = 1$, then $\\mathbb P (X \\ge a) \\le \\mathbb E(X)/a$.\n",
    "\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",
    "$$\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^+ \\equiv \\sup_x (\\hat{F}_n(x) - F(x) )\n",
    "$$\n",
    "and\n",
    "$$\n",
    "    D_n^- \\equiv \\sup_x (F(x) - \\hat{F}_n(x)),\n",
    "$$\n",
    "and the two-sided version\n",
    "$$\n",
    "    D_n \\equiv \\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",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# This is the first cell with code: set up the Python environment\n",
    "%matplotlib inline\n",
    "from __future__ import division\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.stats\n",
    "from scipy.stats import binom\n",
    "import pandas as pd\n",
    "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__.plotUniformKolmogorov>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAEbCAYAAAB3DOvsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAMTQAADE0B0s6tTgAAIABJREFUeJzt3Xt8VNW5//HPk2CIIBMDggkSNIqxcgTLUctP0SqtNgKK\n1qqttZ4qSuVSRI7VClVO6anX2pelVkWrVKBV9Idt1VbJ6a+WBLxAewjSqtyDhEuUomZgaogk6/fH\nTuIQEphMZmbvmXzfr9e8YO9Zs+fJgswzz95rrW3OOURERPyU5XcAIiIiSkYiIuI7JSMREfGdkpGI\niPhOyUhERHynZCQiIr5TMhIREd/FlIzMbLaZVZlZo5kNPUi7i8zsXTNba2aLzOyIxIUqIiKZKtbK\n6P8CI4DN7TUws57AE8BY59xJwA5gZmcDFBGRzBdTMnLOLXPObQfsIM1GASudc+ubth8BrupkfCIi\n0gUk8prRQOC9qO3NQIGZ6bqUiIgcVDc/3tTMDOgP7Pbj/UVEJKF6AdtdJxY7TWQy2gJcELVdDOxw\nzjW20bY/sDWB7y0iIv4aAGyL98WJTEaLgV+YWYlzbh0wEVjYTtvdANXV1YRCoQSGkJlmzJjB3Xff\n7XcYgad+ip36Kjbp2k/hMBQVQXU1JPwjNhKBWbNgwQKYOZPwVVdRdOyx0MkzXTElIzObA4wBjgbK\nzGy3c67EzGYB25xzjzvn9pjZDcALZpYN/AP49sGOGwqFlIxikJOTo36Kgfopduqr2KR7P4VCCU5G\n5eUwbhz07w+rVsGJJ3qZLwFiSkbOuQnt7P+vVtt/AP6QgLhERCRoXnkFbroJpkyBrMSOTfNlAIN0\nTGlpqd8hpAX1U+zUV7FRP7Vy771JO7T5cadXMwsBtbW1tW2WwHV1ddTX16c8LmlfTk4Oubm5foch\nIjEIhyEvD2prk3DN6ID3CpOXlweQ55yL+5xd4Cqjuro6iouLqamp8TsUiVJQUEBVVZUSkkhX8MEH\n0K9fSt8ycMmovr6empoajbQLkHA4TFFREfX19UpGIpksEoHp0+HZZ6GqCnr0SNlbBy4ZNdNIOxGR\nFGoeKVdYCMuWpTQRgW4hISLStUUi3gi50aO9UXLl5d6Q7RQLbGUkIiIpMHYs7N372bwhnygZiYh0\nZU8+6S3XkJ3taxg6TdcFXXPNNYwbN65l++OPP2b06NHk5+czcOBAHyMTkZQ77jjfExEoGSXEyJEj\nmTkzfe8jOGfOHD744AN27tzJli1b/A5HRLogJaMUCfIk3o0bNzJ48GC6ddNZW5GMVFEBN94IPixy\nECslo06aOHEiS5cu5f7776dXr14tw9GfeuopioqKeOSRRyguLqZv374AFBcXM3fu3P2OkZWVxauv\nvtqyvXz5ckaOHMlRRx1FcXExM2fOpLGxrTtxeOrq6rjjjjs46aSTCIVCDBo0iAULFrQ8/5Of/IRj\njz2W3r17M378+P0S45e//GXmzZvHwoULCYVCTJo0KSH9IiIBED1SbvDgQCej9P0qvHev92hLz55t\nnwNtaPD+cVrr3t17xOHRRx9lzZo1nHPOOfzoRz9q2W9m1NTUsHr1at5++22yYlxUcO3atZx//vnM\nnTuXyy+/nK1btzJ27FgOP/xwpk+f3uZrbrjhBjZs2MBLL71ESUkJNTU1LStYPP3009x77728/PLL\nnH766cydO5fJkyfzrW99C4A///nPXHfddTQ0NDB//vy4+kBEAuqss2BACCorfR0pF4v0rYzuucdb\nfKmtx5o1bb9mzZq2299zT1JCNDNmz55Njx49Yl654JFHHmHs2LFcccUVmBlFRUXceuutB1RTzXbt\n2sXTTz/No48+SklJCeAt3fP5z38e8Cq0cePGMXz4cLKzsxk/fjxDhw5NzA8oIsHT2Ai33eb9fcIE\n3+YNdVT6VkbTp8N//mfbz/Xs2fb+z33OWzmwtTirokPp168f3Tt47PXr17NkyRJ69+7dsq+xsZH2\nFrStqqrCzFoSUWtbt27l0ksv3W9fcXFxh2ISkTSSlQUDBnh/nzgxbUqONAmzDd27f3bnqNaP9oYp\nZme33b6Tyai9U3Bt7e/VqxeRqFOF27dv3+/5goICvvnNb/Lhhx+2PD7++GNq20qiwHHHHQfAunXr\n2nx+wIABbN68eb99rbdFJMPcdJPfEXRY+iajACkoKGg3GbR2+umn88wzz1BbW0s4HGb69OmYWcvz\nkyZNYtGiRSxatIhPP/2UxsZGNm7cSFlZWZvHO+qoo7jqqquYPHlySww1NTVUVlYC8O1vf5u5c+ey\nfPlyGhoaeOKJJ3jrrbc6+ROLiCSWklEC3HLLLaxdu5bevXvvd3qtLT/+8Y8JhUIUFRVxxhlncNll\nl+33/Omnn86f/vQnfvnLX3LMMcdw1FFHccUVVxx0/s8vf/lLzj33XEaNGkWvXr04++yzeeeddwC4\n+uqrue2227jiiivo27cvK1asOOA9RSTYdu/dzRvVb+z32P1RDSxf7ndoCRO4m+s136ipvRvvSerp\n30TEX29Uv8FZc8/ab9/rrxzDmX2HwUsvHdBeN9cTEZGUiIz9BuGJ90EbH//huFOCf5SMRETS0AWP\nfg1+0P6acrm5kJOTwoA6SclIRCTo6g+c4P/738PIg0wfysnxElK6UDISEQm6nAOnn/TsmfzrQamk\n0XQiIuI7JSMRkSBpaPA7Al8oGYmIBEHzCtv/8R9+R+ILXTMSEfFbeTmMGweFhfCrX/kdjS9UGYmI\n+CX6fkNTpqTNCtvJoGQUME8++eQBq2pPmzaNo48+mlAoxNtvv+1TZCKSUP/8JwwdCitXwqpVcPPN\n7S/y3AXoNF0ARS+c+uabb/Loo4+yefNmCgoKfIxKRBKqTx948EEYM6ZLJ6FmSkYBt2HDBvr166dE\nJJJpzGDsWL+jCAydpkuAvXv3MmPGDAYNGkSfPn0477zzWLVqVcvz8+bNY9iwYRx55JEUFhYybdq0\nlufKysoYOnQooVCI888/n+rq6pbnfvjDHzJ+/Hi2b99Or169GDJkSEp/LhGRVEnLysg52L07scfs\n1cv7ohKPG2+8ke3bt7N06VKOPvpo5syZw4UXXsi6det47rnn+MEPfsDChQs599xz+eSTT1i5ciXg\n3aX1kksu4eGHH+baa69lxYoVXHLJJRxxxBGAl4yKi4u58847D3oLCREJsIoK+Ne/4MIL/Y4k0NKy\nMtq921sePZGPeJPbhx9+yPz583n44YcpLCwkKyuLSZMmEQqFeOmll5g9eza33347I0eOJCsri549\ne3LOOecA8MwzzzBkyBCuv/56srOzOfPMM/mPLjrHQCTjRCIwdao3Um7bNr+jCby0rIx69fLu05Ho\nY8Zjw4YNAAwfPrxln3OOTz/9lG3btrF582ZKSkrafO3WrVsPGDnXeltE0lBFxWfzhioru+xw7Y5I\ny2RkFpwFAgsKCjAzVq9ezYABAw54fv78+axbt44xY8Yc8NyAAQP429/+tt++qqqqpMUqIkkWicD0\n6TB3Ltx1lzd3KCstT0ClnHqpkwYOHMill17K5MmTW67r7N69m8WLF1NTU8PUqVO57777+Mtf/kJj\nYyN79uyhoqICgKuuuoq///3vzJ07l4aGBt58800WLFjg548jIp0RDsPmzV41NHWqElEHqKcS4Omn\nn+a0007jggsuIC8vj5NPPpknnngCgPHjx3PPPfdw8803k5+fT0lJCS+88ALgnZL73e9+x4MPPkh+\nfj533HEHkyZN8vNHEZHOKCyEF1/Uabk4mHMu9W9qFgJqa2trCbU639Z8P/W2nhN/6N9ExF9vVL/B\nWXPP2m/fn77+Oud/7kyfIvpM8+cDkOeci/uG56qMREQ6KhJJ/CiqLk7JSESkIyoq4NRT4eGH/Y4k\no8ScjMxskJm9ZmZrzWy5mZ3cTrvvm9nbZlZpZq+b2RmJC1dExCfR84amTIHbb/c7oozSkaHdjwFz\nnHMLzOxrwDzgC9ENzOxUYCJwsnPuEzO7GvgFMPyAo4mIpAvNG0q6mCojM+sLnAb8BsA59zxQZGbH\nt2rq8BJc8xTSI4FqRETS1Ztv6n5DKRBrZVQE7HDONUbt2wIMBDY173DOrTaznwFVZrYL2At8MVHB\nioik3PDh8O67UFTkdyQZLaErMJjZccBlwPHOuffNbDLwHHBOW+1nzJhBTk4OAKWlpZSWlrY8Fw7H\nPUJQEkz/FtKlmQUyEXU7zL/3Lisro6ysDID6+vqEHDOmeUZNp+nWA72bqyMz2wGMcM5timp3C3Ci\nc25C03YPYA+Q45zbF9Wu3XlGdXV1FBcXU1NT0+kfThKnoKCAqqoqcnNz/Q5FpMtpa57R6+Ne58yi\nzJlnFFNl5JzbaWYrgWuAeWZ2OVAdnYiabAKuNbOezrkIcDGwNjoRHUpubi5VVVUJy7aSGDk5OUpE\nkrma15S76CL4ylf8juYAp/Q7hdfHvU4kAhdcAH/6k7cvk8S8AoOZlQBPAX2AWuBa59w7ZjYL2Oac\ne7yp3V14p+rqgAgwxTlX2epY7VZGIiIpVV7ujZTr399b4DTAAxTCYe+WN7W1wVksOlGVUeCWAxIR\nSYk0XGE7k5NRWt5CQkSkU/73f+HKK71qSPOGAiHYXwNERJKhd2+46SbNGwoQVUYi0vUUF3tL+0hg\nqDISERHfKRmJSOZ67TXYs8fvKCQGSkYiknmaV9guLYW//tXvaCQGumYkIpmleYXt/v1h1SoYNMjv\niCQGqoxEJDNEIt4IueYVtpcsUSJKI6qMRCQzPP20VwmpGkpLWoFBRDJDY9MdbgK+ikJnaAUGEZGg\ny+Ak1BXoX09E0su+mG8CIGlEyUhE0kd5OQweDCtX+h2JJJiSkYgEX/NIuTFjYPJk+Pzn/Y5IEkzX\njEQk2KLvN6QVtjOWKiMRCa6f/9yrhrTCdsZTZSQiwTVqlPdQEsp4SkYiElxKQl2GTtOJiIjvVBmJ\nSIu6OqivT/GbRiLwyCPeenK5uSl+8/QSjnt9g+BTMhIRwEtE+fnen6nVE7gVfpzq901PubmQk+N3\nFImnZCQigFcR1dVBdXUK1j2LRGDWLFiwAGbOhBtv1HI+McrJycwCUslIRPYTCiU5GUXPG1pVoUEK\nAmgAg4ik2ooVmjckB1BlJCKpdeutfkcgAaTKSEREfKdkJCLJsXWr3xFIGlEyEpHEal5he8gQ+Ogj\nv6ORNKFkJCKJU14OQ4d6q2uvWOFNXBKJgZKRiHRe9P2GNFJO4qDRdCLSeVdfDbt26X5DEjclIxHp\nvDlzoF8/raIgcVMyEpHOKyjwOwJJc/oaIyIivlMyEpHYlJfDlVfCp5/6HYlkICUjETm46JFyI0bo\nupAkha4ZiUj7olfY1kg5SSJ9xRGRAzkH06Zp3pCkjCojETmQGQwerGpIUkbJSETaNn683xFIFxLz\naTozG2Rmr5nZWjNbbmYnt9OuyMxeNLM1ZvYPM5ucuHBFRCQTdeSa0WPAHOfcScD9wLx22v0OeMo5\n9znn3CnAc52MUUSSJRKBV1/1OwqR2JKRmfUFTgN+A+Ccex4oMrPjW7X7MlDnnPtt8z7n3M7EhSsi\nCdO8wva993oDFkR8FGtlVATscM41Ru3bAgxs1W4w8E8ze8bMVprZ82ZWnIhARSRBWq+wvXixN2BB\nxEeJHsDQDRgJDHfOrTGzG/FO052R4PcRkXi0mjdUV3Qi9Xu8p8Jhf0OTri3WZFQNFJpZVlR1NBCv\nOoq2Bah0zq1p2l4APGxm2c65htYHnTFjBjk5OQCUlpZSWlra4R9ARGK0bx/cfLNXDU2ZQl19Fvn5\nUFf3WZPcXGj6lRRpV1lZGWVlZQDU19cn5JjmYjxXbGavAvOcc/PM7HLgNufcF1q16QGsBr7onNtu\nZlcAM51zQ1q1CwG1tbW1hEKhhPwgIhKDhgbIzga8SigvD6qrofnXMCfHS0gisQqHw+Tl5QHkOefi\nrq87cppuAvCUmc0AaoFrAcxsFrDNOfe4c+5fZjYB+KN556BrgW/EG5yIJFhTIooWCn2WjET8EnMy\ncs6tA85qY/9/tdr+f8CwzocmInH79FM47DC/oxCJmdamE8kkzSPlLrnE70hEOkTJSCRTNM8bqqyE\n2bP9jkakQ5SMRNJd63lDWmFb0pAWShVJZ3v2wLBhUFCgFbYlrSkZiaSzI46Axx6D887THVglrSkZ\niaS7L33J7whEOk1fpURExHdKRiLpoLwcFi70OwqRpFEyEgmy6JFyH3/sdzQiSaNrRiJB1WqFbY2U\nk0ymykgkaDRvSLogVUYiQbNvH7z/vqoh6VKUjESCJi8Pnn3W7yhEUkqn6URExHdKRiJ+iURg506/\noxAJBCUjET80r7D905/6HYlIICgZiaRS65Fyd9/td0QigaABDCKponlDIu1SZSSSCm+/rXlDIgeh\nykgkFf7t32DjRjj6aL8jEQkkJSORBKurg/r6Np44/GgIpzycdoUDFIuIkpFIAtXVQX6+92c6yM2F\nnBy/oxBRMhJJnEiE+ltmUVd3P9XVEAr5HdCh5eR4CUnEb0pGIonQPFKu7wmAl4jSIRmJBIVG04l0\nRvO8odGjYcoUeOUVvyMSSUuqjETitXo1fPWrUFgIq1Z5w7U1KEAkLqqMROJVWAjTpmnekEgCqDIS\niVffvvDd7/odhUhGUGUkIiK+UzISOZSlS+GDD/yOQiSjKRmJtCcSgalTYdQo+Nvf/I5GJKPpmpFI\nWyoqvHlDhYVaYVskBVQZiURrPW9II+VEUkKVkUi0xYu9SkjVkEhKKRmJRLvsMm8ia5ZOGoikkpKR\nSDQz7yEiKaWvf9I17d3rdwQiEkXJSLqeigrvzqt/+YvfkYhIEyUj6Tqa5w2NHu0t4/PFL/odkYg0\n0TUj6Ro0b0gk0FQZSeZ78knNGxIJuJiTkZkNMrPXzGytmS03s5MP0f4pM2s0M93vUvw1apRXDU2d\nqiHbIgHVkd/Mx4A5zrmTgPuBee01NLOvAvWA61x4IgnQv7+qIZGAiykZmVlf4DTgNwDOueeBIjM7\nvo22RwPTgWmAJmyIiMghxTqAoQjY4ZxrjNq3BRgIbGrV9nHgVudcxDR5UFIlEoF77oHvfQ+OPDLm\nl9XVQX194sII67bjInFJ6Gg6M7seeM85Vx5L+xkzZpCTkwNAaWkppaWliQxHuory8s9Gyt14Y8zJ\nqK4O8vO9PxMpNxea/luLZKSysjLKysoAqE/Qtzlz7tCXdZpO060HejdXR2a2AxjhnNsU1e7XwDlA\nA94pumPxKqhLnHNvRbULAbW1tbWEQhrfIHGKRGD6dG+03F13eaPlsrNjfnk4DHl5UF0NifxvmJPj\nJSSRriAcDpOXlweQ55yL+9xATJWRc26nma0ErgHmmdnlQHV0Impq963obTNrBIY453bHG6BImyoq\n4LrrvGpo1apODVAIhRKbjESk4zoymm4CcKOZrQVuA64FMLNZZvaddl7j0CAGSYaNGzVvSCSDxHSa\nLuFvqtN0EgDNp+lqa1UZicQrUafpNANQRER8p2QkwbZhg98RiEgKKBlJMDWvsD1sGGzf7nc0IpJk\nSkYSPBUVcOqp3npylZXecj4iktGUjCQ4IhG46abPVthesgQGDfI7KhFJAd3PSIJj4kTYvNmbN6Qk\nJNKlKBlJcPzsZ95SPrrNg0iXo2QkwdG7t98RiIhP9BVURER8p2QkqVVRAaWl3mAFEZEmSkaSGs3z\nhkaP9h6HH+53RCISILpmJMlXUfHZ/YYqK7WwqYgcQJWRJNftt382b0grbItIO1QZSXKddprmDYnI\nISkZSXJdcYXfEYhIGtBpOhER8Z2SkXReJAIvvOB3FCKSxpSMpHPKy2HoUHjwQdi3z+9oRCRNKRlJ\nfJpX2B4zxvvz1Vehmy5Bikh89OnRBdXVQX19Jw6wbBlMngwFBbD0LTjhBNiTsPBSJhz2OwIRaaZk\n1MXU1UF+vvdn/M4G3oLNwL8nJCzf5OZCTo7fUYiIklEXU1/vJaLqagiF4jyIc2CW0Lj8kpPjJSQR\n8ZeSURcVCnUiGZEZiUhEgkMDGKR9nTuXJyISMyUjOVDzSLmRI71TciIiSaZkJPtrnje0ahUsWJAx\n14ZEJNiUjMTTet7QkiVa3FREUkYDGMQbYnfaadCvn1bYFhFfKBmJN7553jw44wzIUrEsIqmnZCSe\n4cP9jkBEujB9DRYREd+pMupKysvhf9cDN/gdiYjIflQZdQXRI+UaG/2ORkTkAEpGma553lBlpff4\nznf8jkhE5ABKRpnqk0/2nzdUXg4nnuh3VCIibdI1o0yVlQV79njVkJKQiAScklGm6t4d5s71OwoR\nkZjoNJ2IiPhOySjdRSLenfJERNKYklE6ax4pd999fkciItIpMScjMxtkZq+Z2VozW25mJ7fR5hQz\nKzezd8xstZk9YWbdExuyHLDC9s9/7ndEIiKd0pHK6DFgjnPuJOB+YF4bbeqAyc65wcCpwBHA9zsd\npXym9byhqVO1uKmIpL2YPsXMrC9wGvAbAOfc80CRmR0f3c45t8E594+mvzvgr8BxiQy4S3vvPRg7\nVvOGRCTjxDq0uwjY4ZyLXktmCzAQ2NTWC8ysJ94iaKqMEuXYY2HzZsjP9zsSEZGESso8IzM7DFgI\nLHbOvZiM90imujrvfnOBlJ0P4fhfHu7Ea0VEkiXWZFQNFJpZVlR1NBCvOtqPmXUDngW2OeemHeyg\nM2bMICcnB4DS0lJKS0tjDjxZ6uq8wqOuzu9Ikic317ufnohIPMrKyigrKwOgPkHf3M27tBNDQ7NX\ngXnOuXlmdjlwm3PuC63aZAPPAR8658Yf5FghoLa2tpZQKBR/9EkQDkNenjd1x5fQIhGYNQuKi2Hi\nxKS8RU6Ol5BERDorHA6Tl5cHkOeci/vcS0dO000AnjKzGUAtcC2Amc3Cq4IeB74OXAqsNrNKwAGv\nOeemxBugX0IhH5JReTmMGwf9+8Mt4yFYeVpEJGlirowS+qZpUBnV1qYwGUUiMH26t5bcXXfBlCka\nri0iacGPykiSYc0ab/Jq//5aYVtEuix9/fZbURF8//uaNyQiXZoqI7/17Km7r4pIl6fKSEREfKdk\nlCrl5bCpzcUqRES6PCWjZIteYXvVKr+jEREJJF0zSqboeUMaKSci0i5VRsnQ+n5DGiknInJQqoyS\nYcUK75ScqiERkZgoGSXDyJFw3nlg5nckIiJpQafpkkWJSEQkZkpGnfGvf4EPa/uJiGQaJaN4lZfD\nkCHwu9/5HYmISNrTNaOOar3C9qWX+h2RpInfr/k93335u1z7+WsZmDeQ9/e8z7bd23iw9EEOP+xw\nALbv3s7S95ZiZjjnMDPOKjqLAaEBnXrvNf9cw+SXJ/NG9Rv0zOnJN0/5Jj8t/Sndstr+CLj/tft5\nsvJJNny4AeccS65dwheP/WLMxzvU60VaUzLqCM0bkk64uORiJv5xIv898r+xpmuK89+azw+X/JD7\nLrgPgP69+vP1U76e0PdtaGzg4mcuZlt4G3d96S5W1qzkoRUPkZebx49G/qjN1+zdt5eLSy7m+Xef\nZ0vtlg4f72CvF2mLTtPF6plnNG9IOmXljpWc0u+UlkQE8IVjvsDijYuT+r5lG8vY+OFGxpSMYdqZ\n03jsosfoltWNh1Y81O5r7jz3Th74ygMUHFEQ1/EO9nqRtqgyitWoUaqGpFOWbVnG2UVn77evZk8N\n2Zad1Pddv2s9AANDAwHocVgPjupxFO9H3mdnZCd9e/b19XgioGQUuyOP9B4icVpWvYxJp0/ab1/5\n5nJGFI3o8LGqPqrixbUv7ldlNcvJzmHC6RMO+npHYkeBJvp40vUoGYmkyPKty1nw1QUt2845Fr27\niN9e+dsOH6s4v5ip/2dqTG1P7ONV8+/VvgdApD7Crn/tItQ91FLF1DfU45yje7fuCTmeSEcpGUWL\nROCOe4Af+x2JZJh1u9ZR2KuQHof1aNk3q3wWU4dP5cQ+J/Lm1jf547o/MqZkDIveWcQDX3ngoMfb\n9NEmXlr7UpvPde/Wfb/KqPSEUk7ofQIvr3+Zn77+UyprKmlwDUz5wpSWNiUPlbCldgtrvruGkj4l\nLH1vKWt3rWVnZCcAf1j3B9bvWs/1/359TMc72OtF2mLOh0mbZhYCamtrawmFQgc8X1cH9fUpDmrZ\nMpg8mXDfEyj662+prYU2QhPpsDe3vskvVvyCzR9vZtywcexr3MeW2i2cPfBsLhx0IQDv73mfu5fe\nzexRs/nOS9/h8YsfT2gM7+58lymvTOH16tdbhmI/8JUHOCz7MACKZxdTXVvNO5PfoaRPCde9cB3z\n35p/wHEaZjbEdLxDvV4yRzgcJi8vDyDPOReO9ziBS0Z1dZCf7/3pl9xc+Ogj70+RVNixewcL/7GQ\nCadP4Hv/8z3uPPdOjUSTtJCoZBS403T19V4iqq5OQWXyxhswYQIUFMAjj8AJJwCQk6NEJKlVWVPJ\nOceew8d1H9OnRx/q9vn4bUzEB4GrjMJhyMsjNafJXnwRqqpgyhTI0pQrEZGOytjKKKXGjvU7AhER\nQSswiIhIAHSNZPT3v+tWDyIiAZbZySgS8daSO/NM2LDB72hERKQdmZuMysth6FBvPTmtKSciEmiZ\nl4yaqyGtsC0ikjYybzTd7bfDqlWqhkRE0kjmJaO774aePTVvSEQkjWReMurVy+8IRESkg1Q+iIiI\n79IzGVVUwIgRsGuX35GIiEgCpFcyikRg6lQYPRquvNJb3ltERNJe+lwzqqiAceOgsFAj5UREMkx6\nVEYzZ3rV0JQpmjckIpKB0qMyGjECrrlGSUhEJEOlRzIqLfU7AhERSaL0OE0nIiIZLeZkZGaDzOw1\nM1trZsvN7OR22l1kZu82tVtkZkfE9AaRCPz617GGIyIiGaQjldFjwBzn3EnA/cC81g3MrCfwBDC2\nqd0OYOYhj9y8wvZjj8Enn3QgpK6hrKzM7xDSgvopduqr2KifUiemZGRmfYHTgN8AOOeeB4rM7PhW\nTUcBK51z65u2HwGuavfAzStsN4+UW7IEDj+8gz9C5tMvRGzUT7FTX8VG/ZQ6sQ5gKAJ2OOcao/Zt\nAQYCm6LjNKf+AAAEMUlEQVT2DQTei9reDBSYWVar13rOOguOOcZbZVsj5UREuixfR9Nt+8YNhCdM\nguxs2BoGYPdu77lw2MfAAqa+vp6wOuSQ1E+xU1/FRv10aInqH3POHbqRd5puPdC7ucIxsx3ACOfc\npqh2lwPXO+dGNW0PBhY75wa2Ot4xwNaE/AQiIhIEA5xz2+J9cUyVkXNup5mtBK4B5jUlneroRNRk\nMfALMytxzq0DJgIL2zjkdmAAsDvewEVEJDB64X2uxy2mygjAzEqAp4A+QC1wrXPuHTObBWxzzj3e\n1O4i4CdANvAP4NvOOSUdERFpV8zJSEREJFmSugJD0ifKZpBY+srMTjGzcjN7x8xWm9kTZtbdj3j9\nEuv/qaj2T5lZo5mFUhVjEHTgd6/IzF40szVm9g8zm5zqWP3Wgb76vpm9bWaVZva6mZ2R6lj9ZGaz\nzayq6fdp6EHaxfd57pxL2gP4M3BN09+/Bqxoo01PoAY4sWn7IeD+ZMYVxEeMfTUIOKXp74Z3PW6m\n37EHrZ+i2n4VeBxoAEJ+xx7EfgL+BlwWtd3X79iD2FfAqXhTVQ5v2r4aWO537Cnup7OB/njTeYa2\n0ybuz/NkBt4X+BjIitq3Azi+VbvLgZejtk/GGxzhe+en8B85pr5q43W3AHP9jj+I/QQcDaxo+uVo\n7ErJqAO/e18Glvkdb5r01VC8EcD9mrYnA4v8jt+nPqs6SDKK+/M8mafpDjZRNlq7E2WTGFvQxNpX\nLZqWXroB+H2SYwuSjvTT48CtzrlISiILllj7aTDwTzN7xsxWmtnzZlacsiiDIaa+cs6tBn4GVJnZ\nFuBmYErKokwfcX+ed6UP/IxhZofhnaJb7Jx70e94gsbMrgfec86V+x1LwHUDRgKznHP/DvwP8Jy/\nIQWTmR0HXIZXMQ3ES0zqqwRKZjKqBgpbZcSBeN86om0BjovaLubAbyqZLta+wsy6Ac/iDaeflqL4\ngiLWfhoJXGJmm8ysqmnfajM7NRVBBkBHfvcqnXNrmrYXAMPMLDsFMQZFrH31NWC1c+79pu1fASOa\nfh/lM3F/nictGTnndgLNE2WbV2dob6LssKZ5TND+RNmMFWtfNX1IPAvscs5NSHmgPou1n5xz33LO\nHeucO94513zaaYhz7q3URuyPDvzuvQIMMLP+TdtjgHedcw0pC9ZnHeirTXjJp2fT9sXAWufcvpQF\nmx7i/zxP8oWuEuB1YC3exeTBTftnAd+JancR8C6wDvgt0Mvvi3SpfsTSV8A38UaGVTY9VgIP+R17\n0Pqpjdd0xdF0sf7unR/1/2kJ8G9+xx7gvrqr6XOqElgGDPM79hT30xy8SrIeb5DHunb6Ka7Pc016\nFRER32kAg4iI+E7JSEREfKdkJCIivlMyEhER3ykZiYiI75SMRETEd0pGIiLiOyUjERHxnZKRiIj4\n7v8DsRLj9cWhbqgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d7b8690>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "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='3')\n",
    "    ax.text(0.5, 0.1, r'$D_n^-=$' + str(round(dMinus, 3)), color='g', weight='heavy')\n",
    "    plt.show()\n",
    "\n",
    "interact(plotUniformKolmogorov, n=widgets.IntSlider(min=3, max=300, step=1, value=10))\n",
    "    "
   ]
  },
  {
   "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",
    "A little more explanation might help: whenever $F \\in \\mathcal C(X)$ (and perhaps more often),\n",
    "\n",
    "$$\\inf_{G \\in \\mathcal C(X)} \\int_{0}^\\infty x dG(x) \\le \\int_0^\\infty x dF(x).$$\n",
    "But $\\mathbb P (C(X) \\ni F) \\ge 1-\\alpha$.\n",
    "\n",
    "Moreover, we can find confidence bounds for <em>an arbitrary collection of\n",
    "functionals, even an uncountable collection</em> and they will have simultaneous coverage\n",
    "probability $1-\\alpha$: there is no need to adjust for multiplicity, since whenever\n",
    "$F \\in C(X)$, <em>all</em> the bounds hold.\n",
    "\n",
    "This general approach&mdash;make a confidence set for an entire function, then use that\n",
    "confidence set to derive confidence bounds on properties (parameters) of that function&mdash;is\n",
    "extremely flexible, and gives conservative results.\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 \\min\\{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",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 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()"
   ]
  },
  {
   "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",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 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 = int(1.0e4) # 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)"
   ]
  },
  {
   "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: Wald's Sequential Probability Ratio Test](sprt.ipynb)\n",
    "+ [Previous: Confidence bounds from the Chebychev and Hoeffding Inequalities](hoeffding.ipynb)\n",
    "+ [Index](index.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%run talkTools.py"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}