2016-02-bias-and-coverage.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pylab\n",
    "import scipy.stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "trying to reproduce Josh's presentation 2013-06-26\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Coverage = how often the true parameter is found within the quoted interval"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "consider a measurement of a parameter with mean $\\mu$ = 0 and standard deviation $\\sigma$ = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "underlying process: Gaussian with mean $\\mu$ (set $\\mu$ = 0 for simplicity, without loss of generality), standard deviation $\\sigma$ = 1\n",
    "\n",
    "(could be a number of events if $\\mu$ large and $\\sigma = \\sqrt{\\mu}$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#----------\n",
    "# PARAMETERS\n",
    "#----------\n",
    "\n",
    "# mean and sigma of the process\n",
    "muTrue = 0.0\n",
    "sigmaTrue = 1.0\n",
    "\n",
    "# bias of the measurement, relative in units of sigmaTrue\n",
    "# i.e. muMeasured = muTrue + sigmaTrue * bias\n",
    "bias = 0.14\n",
    "\n",
    "# the number of sigmas we consider for the coverage of an interval\n",
    "numSigmas = 1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculation using Monte Carlo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# number of times the experiment is repeated\n",
    "n = 10000000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# generate experimental outcome (number of events) a large number of times\n",
    "outcomes = np.random.normal(muTrue, sigmaTrue, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# apply bias to our measurements\n",
    "measurements = outcomes + sigmaTrue * bias"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# calculate expected coverage: integral of Normal distribution between -numSigmas and +numSigmas\n",
    "normalDistribution = scipy.stats.norm(loc = 0, scale = 1)\n",
    "expectedCoverage = normalDistribution.cdf(+ numSigmas) - normalDistribution.cdf(- numSigmas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# assume our measurement uncertainty reflects the \n",
    "# true width of the Gaussian above\n",
    "#\n",
    "# calculate how often the true value 'muTrue' actually\n",
    "# lies within the +/- 1 sigma interval around the value in 'measurements'\n",
    "\n",
    "actualCoverage = sum(np.abs(muTrue - measurements) <= numSigmas * sigmaTrue) / float(n)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "expected coverage: 0.6827, actual coverage: 0.6776, actual - expected = -0.50%\n"
     ]
    }
   ],
   "source": [
    "print \"expected coverage: %.4f, actual coverage: %.4f, actual - expected = %.2f%%\" % (\n",
    "    expectedCoverage,\n",
    "    actualCoverage,\n",
    "    (actualCoverage - expectedCoverage)  * 100\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analytic calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "undercoverage for 1.0 sigma intervals at bias= 0.14 is: 0.00472716395169 (0.47%)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEPCAYAAACOU4kjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm8HfP9x/HXJ5FQayy1RkQlxNJKLaXWq4hQBLUkpRJb\nqrGFIlq1pGoLih9SP4SgIfyoLWRDbmjJIpJYkkhoRFYhC7EkIvn8/vjOjeO6yT03d+bMmbnv5+Nx\nHvfMOTPnfD4dPZ/M9zPzHXN3RERE4tQo7QBERCR/VFxERCR2Ki4iIhI7FRcREYmdiouIiMROxUVE\nRGKXeHExs/ZmNsnMJptZjxreb2pm/c1sipm9bmYtotcPMbM3zGy8mY02s4MKthkWfeZYM3vTzDZJ\nOg8RESneGkl+uJk1Au4EDgZmAaPN7Bl3n1Sw2hnAfHdvbWYnAb2AjsAnwJHuPsfMdgYGA80Ltuvk\n7mOTjF9ERFZP0kcuvwCmuPs0d18K9Ac6VFunA/Bg9PwJQiHC3ce7+5zo+bvAmmbWpGA7DemJiJSp\npH+gtwKmFyzPiF6rcR13XwYsNLONClcws+OBsVGBqnJ/NCT2l/jDFhGR+ki6uFgNr1Wfb6b6Ola4\nTjQkdj3QtWCd37r7rsD+wP5mdkoMsYqISEwS7bkQjlRaFCw3J/ReCk0HtgZmmVljYH13XwBgZs2B\nfwG/c/cPqzZw99nR3y/N7BHC8Ns/q3+5mWniNBGR1eDuNR0cFC3pI5fRQCsz28bMmhIa9c9WW+c5\noHP0/ATgZQAzawYMAC5z9xFVK5tZYzPbOHreBDgSeGdlAbh7bh9XXXVV6jEoN+Wn/PL3iEOiRy7u\nvszMzgWGEApZH3efaGY9gdHuPgDoAzxsZlOAeYQCBHAOsB1whZldSRgqawd8BQw2szWAxsCLwL1J\n5lGuPvzww7RDSEyecwPll3V5zy8OSQ+L4e6DgB2qvXZVwfMlwIk1bHctcO1KPnaPOGMUEZF46XTe\nDOvSpUvaISQmz7mB8su6vOcXB4trfK0cmZnnOT8RkSSYGV7mDX1JUGVlZdohJCbPuYHyy7q85xcH\nFRcREYmdhsVEROR7NCwmIiJlScUlw/I87pvn3ED5ZV3e84uDiouIiMROPRcREfke9VxERKQsqbhk\nWJ7HffOcGyi/rMt7fnFQcRERkdip5yIiIt+jnouIiJQlFZcMy/O4b55zA+WXdXnPLw4qLiIiEjv1\nXERE5HvUcxERkbKk4pJheR73zXNuoPyyLu/5xUHFRUREYqeei4iIfI96LiIiUpZUXDIsz+O+ec4N\nlF/W5T2/OKi4iIhI7NRzERGR71HPRUREypKKS4bledw3z7mB8su6POc3Zkw8n6PiIiIiAPTpAyec\nEM9nqeciIiI8/zyceSYMHw477FD/nssacQUmIiLZNGoUnHYaPPccbL99PJ+pYbEMy/O4b55zA+WX\ndXnKb8oU6NAB7r8f9torvs9VcRERaaDmzIH27eGaa+DII+P9bPVcREQaoM8/h4oKOPZYuOKK778X\nx3UuKi4iIg3MkiXw619D69bQuzdYtTKiiygbuDyN+1aX59xA+WVdlvNbvhy6dIH11oM77/xhYYmL\nzhYTEWkg3KF7d5g5EwYPhsaNk/uuxIfFzKw9cBvhKKmPu99Y7f2mwEPA7sCnwEnu/pGZHQLcADQB\nvgEudfdh0Ta7AX2BtYAX3L37Sr5bw2IiIpHrroPHHgvXsjRrtvL1yn5YzMwaAXcChwE7A53MrE21\n1c4A5rt7a0IR6hW9/glwpLvvCnQBHi7Y5h/Ame6+PbC9mR2WXBYiItl3333hMWjQqgtLXJLuufwC\nmOLu09x9KdAf6FBtnQ7Ag9HzJ4CDAdx9vLvPiZ6/C6xpZk3MbHNgPXcfFW3zEHBMwnmUpSyP+9Ym\nz7mB8su6rOX39NPhjLDBg2GLLUrznUkXl62A6QXLM6LXalzH3ZcBC81so8IVzOx4YGxUoLaKPmdV\nnykiIkBlJXTtCgMGhLPDSiXphn5NY3bVmyDV17HCdcxsZ+B64NA6fOYKXbp0oWXLlgA0a9aMtm3b\nUlFRAXz3r4+sLle9Vi7xxLlcUVFRVvEoP+WXxfymTIHLL6/gscdg0aJKKitrXr+yspK+ffsCrPi9\nrK9EG/pmtjdwtbu3j5YvA7ywqW9mA6N1RppZY2C2u28avdcceAno7O4jotc2B4a5+47RckfgQHf/\nQw3fr4a+iDRI778PBxwQTjc+7ri6bVv2DX1gNNDKzLaJzgrrCDxbbZ3ngM7R8xOAlwHMrBkwALis\nqrAARH2Yz83sF2ZmwKnAM8mmUZ6q/uWRR3nODZRf1pV7fjNnQrt20LNn3QtLXBItLlEP5VxgCPAu\n0N/dJ5pZTzOrmsmmD7CJmU0BugOXRa+fA2wHXGFmY83sTTPbJHqvW7TdZMIJA4OSzENEJCvmz4fD\nDoPf/x7OOiu9ODT9i4hITnzxBRxySBgO69Wr9vVXRnOL1ULFRUQaiiVLwszG22wD995bv2ldstBz\nkQSV+7hvfeQ5N1B+WVdu+X37LXTqBBtsAHffndx8YXWhucVERDJs+fJwe+Ivv4Rnn4U1yuRXXcNi\nIiIZ5Q4XXgijR8OQIbDOOvF8bhzDYmVS40REpK6uvjpcgT9sWHyFJS7quWRYuY37xinPuYHyy7py\nyO/mm8MMx0OGwIYbph3ND+nIRUQkY+65B+66C155BTbdNO1oaqaei4hIhjz6KFxySRgOa9Uqme9Q\nz0VEpAF5+mm46CJ48cXkCktc1HPJsHIY901KnnMD5Zd1aeQ3eHCY0uX552HnnUv+9XWmIxcRkTI3\nfDiccgo88wzstlva0RRHPRcRkTI2YgQcfXTotRx8cGm+U9O/iIjk2JtvhsLSt2/pCktcVFwyLM/j\n2nnODZRf1pUiv7ffhiOOgP/93/A3a1RcRETKzKRJ4Z4st90Gxx6bdjSrRz0XEZEyMmUKHHQQXHst\ndO5c+/pJUM9FRCRHpk4NN/u66qr0CktcVFwyLM/j2nnODZRf1iWR30cfhab9pZeme3viuKi4iIik\nbMYM+NWv4Lzz4Jxz0o4mHuq5iIikaNYsqKiArl3h4ovTjiZQz0VEJMNmzw5HLKefXj6FJS4qLhmW\n53HtPOcGyi/r4shvzpzQYznlFLjssvrHVG5UXERESuzjj8MRS8eO8Je/pB1NMtRzEREpoY8/Dtex\ndOwIV16ZdjQ1U89FRCRDqo5YTjqpfAtLXFRcMizP49p5zg2UX9atTn5z5oQjlhNOCBdJ5p3u5yIi\nkrCqwtKpU/6PWKrUqediZuu4+5cJxhMr9VxEJG1Vpxv/9rdwxRVpR1OckvVczGwfM5sATIyWdzWz\n3vX5YhGRvJs5M1wgecop2SkscSm253IrcBgwD8DdxwMHJBWUFCfP49p5zg2UX9YVk9/06aGwnH46\nXH554iGVnaIb+u4+vdpLy2KORUQkF6ZNgwMPhD/8AXr0SDuadBTVczGzJ4C/A3cCewPnA3u4e8dk\nw6sf9VxEpNT++99w5f2FF8L556cdzeop5XUuZwPnAFsBM4C20bKIiEQmTw5DYT16ZLewxKWo4uLu\nn7r7ye6+mbtv6u6nuPu8pIOTVcvzuHaecwPll3U15TdhQjjd+Oqr4eyzSx5S2SnqOhcz+58aXv4M\neMPdn4k3JBGRbBk/Htq3h5tuCmeGSfE9l3uANsD/RS/9BpgKbAz81927r2Lb9sBthKOkPu5+Y7X3\nmwIPAbsDnwInuftHZrYR8ASwJ/CAu59fsM0wYAvga8CBdu7+aQ3frZ6LiCRq1Cg46ii44w448cS0\no4lHHD2XYq/Q/xmwr7svi774H8CrwH7A26sIsBHhJICDgVnAaDN7xt0nFax2BjDf3Vub2UlAL6Aj\nsBj4C7BL9Kiuk7uPLTJ+EZHY/fvfcNxx0KdPKDDynWIb+hsC6xYsrwNsFBWbJavY7hfAFHef5u5L\ngf5Ah2rrdAAejJ4/QShEuPtX7v7aKj6/wc+Lludx7TznBsov6yorK3npJTj2WOjXT4WlJsX+QPcC\nxpnZA2bWFxgL3GRm6wAvrmK7rYDC62NmRK/VuE5UrBZGQ2K1ud/M3jSznN4NQUTK1WuvhXnCnnwS\nDj007WjKU9Fzi5nZFoQjEQNGufusIrY5ntAP6RotnwLs6e4XFKzzTrTOrGj5/WidBdFyZ2D3aj2X\nLdx9dlTc/gU87O7/rOH71XMRkVg99hhccAE89xzsuWfa0SSjlD0XCD2Q2cBaQCsza+Xur9SyzQyg\nRcFyc0LvpdB0YGtglpk1BtavKiwr4+6zo79fmtkjhKL3g+IC0KVLF1q2bAlAs2bNaNu2LRUVFcB3\nh+5a1rKWtVzM8gsvQL9+FQwdCvPmVVJZWV7xre5yZWUlffv2BVjxe1lv7l7rAziT0LhfAAwjnKX1\nchHbNQbeB7YBmgLjgB2rrdMN6B097wj0r/Z+Z+COap+5cfS8CeEMtq4r+X7Ps2HDhqUdQmLynJu7\n8suiW291b9HC/b338plfoei3s6j6sLJHsUcuFxBOCR7h7geZWRvguiIK1zIzOxcYwnenIk80s57A\naHcfAPQBHjazKYSJMVdMKWNmU4H1gKZm1gFoB3wEDDazNaJC8yJwb5F5iIjUiTv07AmPPgqvvgot\nWsCsWpsCUux1LqPdfU8zGwfs5e5LzOxdd985+RBXn3ouIlIfy5fDH/8Iw4bB4MGw2WZpR1Qapey5\nzDCzZsDTwFAzWwBMq88Xi4iUs2+/hbPOgvfeC8Vlww3Tjihbip1b7Fh3X+juVwNXEIayjkkyMKld\nVUMuj/KcGyi/crd4cbjX/ezZMHToDwtL1vMrhVqLi5k1MrMVV9S7+3B3f9bdv0k2NBGR0lu0CI44\nAtZcE559FtZZJ+2IsqnYnsszwHnu/lHyIcVHPRcRqYtPPgmFZffd4a67oHHjtCNKRyl7LhsC75rZ\nKODLqhfd/ej6fLmISLmYNg0OOwyOPx6uuQasXj+tUuz0L1cARwJ/BW4peEiK8jzum+fcQPmVm3ff\nhf33D7cl/tvfai8sWcsvDUUdubj7cDPbBmjt7i+a2dqEa0xERDLt9dfDBJS33AInn5x2NPlRbM/l\nLKArYSbk7cysNXC3ux+cdID1oZ6LiKzKgAFw+unw0EPhZl8SxNFzKXZY7BxgX+BzAHefAmxany8W\nEUnTAw/AmWeGCShVWOJXbHFZUnjqcTT1ig4JUpbncd885wbKL03ucP31YUqXykrYa6+6f0Y551cu\nij1bbLiZ/Rn4kZkdSphs8rnkwhIRid+yZWG6/FdfDfdk2XLLtCPKr2J7Lo0ItyNuR7ify2DgvnJv\naKjnIiJVFi+GU06BefPg6adhgw3Sjqh8xdFzKba4HAu84O6ruqVx2VFxERGABQvgmGNg881D837N\nNdOOqLyVsqF/NDDZzB42s19HPRdJWZ7HffOcGyi/Upo2DfbdF/bYI0ybH0dhKaf8ylWxE1eeBrQi\n3Jjrt8AHZnZfkoGJiNTXuHGhsHTtGq5jaVTsP6el3ooaFluxslkToD1wGrC/u/84qcDioGExkYZr\nyJDQY7nrrjDDsRSvZMNiZtbezPoSbll8PHAfsEV9vlhEJCn33w+nngr/+pcKS1qKPUjsQrhR2Pbu\n3tndX3D3b5MLS4qR53HfPOcGyi8p7nDllXDttTB8OOy3XzLfk/f9F4di5xbraGabAYdamNFtlLvP\nTTQyEZE6WLIkXHE/eXKYL2xTzSGSqmJPRT4BuBmoJFznsj9wibs/kWh09aSei0jDMH9+mHxyk03g\n4Ydh7bXTjijbSnmdy3jg0KqjFTP7MfCiu+9any9PmoqLSP598AH8+tdw1FFw4406IywOpbzOpVG1\nYbB5ddhWEpLncd885wbKLy7/+U/oq1xwAdx0U+kKS973XxyKvRhykJkNBh6Nlk8CBiYTkohI7fr1\ngwsvDMNghx2WdjRSXdHXuZjZccB+hJ7LK+7+VJKBxUHDYiL54x5mNH7wwTBd/i67pB1R/pSy57It\nMNvdF0fLPwI2c/cP6/PlSVNxEcmXr78ON/eaOhWeeQY22yztiPKplD2X/wOWFywvi16TFOV53DfP\nuYHyWx1z5sBBB4X721dWpltY8r7/4lBscVmj8GZh0fOmyYQkIvJ948eHm3odcUTotay1VtoRSW2K\nHRYbCtzh7s9Gyx2A89394ITjqxcNi4lk31NPhYkn77wTTjop7WgahlL2XLYD+gFV922bAfzO3T+o\nz5cnTcVFJLvc4brr4O67Q4HZY4+0I2o4StJzie5Cubu77w3sBOzs7vuUe2FpCPI87pvn3ED51ear\nr+C3vw1N+5Ejy6+w5H3/xaHW4uLuy4FLo+dfuPuixKMSkQZr+nTYf/9wQeTw4brPfVYVOyx2A/Ap\n8BjwZdXr7j4/udDqT8NiItny2mtw/PHQvTtcckk4M0xKr5Q9l6k1vOzu/pP6fHnSVFxEsuO+++DP\nf4a+fcNZYZKekl3n4u7b1vAo68LSEOR53DfPuYHyK/TNN3DOOXDzzfDqq9koLHnff3Eo9k6Ua5vZ\nX8zsnmi5tZkdmWxoIpJ3c+fCIYfARx+Fxv0OO6QdkcSl2GGxx4AxwKnuvks0/cvr7t426QDrQ8Ni\nIuVr9Gj4zW+gc+cwV5imyi8fpZz+ZTt37wUsBXD3rwkTWNbKzNqb2SQzm2xmPWp4v6mZ9TezKWb2\nupm1iF7fyMxeNrNFZvY/1bbZzczeij7ztiJzEJEy8cADYfjr9tvhmmtUWPKo2F36TXS04rDiosol\ntW0UXSNzJ3AYsDPQyczaVFvtDGC+u7cGbgN6Ra8vBv4C/LGGj/4HcKa7bw9sb2YNcsLtPI/75jk3\naLj5VfVXbrgBXnkl3D0yi/K+/+JQbHG5ChgEbG1m/YCXiK59qcUvgCnuPs3dlwL9gQ7V1ukAPBg9\nfwI4GMDdv3L316hWxMxsc2A9dx8VvfQQcEyReYhISmbOhIoKmDEDRo2CHXdMOyJJUrFniw0FjgO6\nEG4Ytoe7Vxax6VbA9ILlGdFrNa7j7suAhWa2US2fOaOWz2wQKioq0g4hMXnODRpefsOHw557wpFH\nhqlcNtggnbjikvf9F4dV3onSzHar9tLs6G8LM2vh7m/W8vk19WWqd9irr2M1rFPXzxSRMuAOf/97\nuAXxQw9Bu3ZpRySlUtttjm+J/q4F7AGMJ/y4/wx4A/hlLdvPAFoULDcHZlVbZzqwNTDLzBoD67v7\nglo+c+taPnOFLl260LJlSwCaNWtG27ZtV/yro2rcNKvLt912W67yKVwuHNMuh3iUX93z+/JLuPFG\nWLy4gpEjYerUSioryyO+OPKrUg7xxJFP3759AVb8Xtabu9f6AP4F/LRgeRfgiSK2awy8D2xDuP/L\nOGDHaut0A3pHzzsC/au935kw3X/hayMJ/RwDXgDar+T7Pc+GDRuWdgiJyXNu7vnP7/77h/n227v/\n4Q/uixenHU388r7/ot/OourDyh7FXufyrrvvXNtrK9m2PXA7ob/Tx91vMLOewGh3H2BmawIPAz8H\n5gEdPbp9cjTtzHpRYVoItHP3SWa2O9CXcET1grtfsJLv9mLyE5H4PPggXHwx3HILnHpq2tHI6ijl\n3GKPEias/Cehv3EKsK67d6rPlydNxUWkdL7+Gs47D/79b3jiCdhll7QjktVVyosoTwPeBS4AugMT\notckRYXjvnmT59wgf/lNngy//CV8+WW48v7TTyvTDilRedt/SaitoQ+Auy8Gbo0eIiIr9O8fjlj+\n+lc4+2xNky9BscNi+wJXExrzKwqSl/nMyBoWE0nO4sVw4YXw4ovw+OPw85+nHZHEJY5hsaKOXIA+\nwIWEySuX1ecLRST7Jk2Ck06CNm1gzBhYf/20I5JyU2zP5TN3H+juc919XtUj0cikVnke981zbpDt\n/B56KNyG+JxzwpBYTYUly/kVI+/5xaHYI5dhZnYT4XqXFXN9ee1X6ItITixaBOeeG+YFe/ll+OlP\n045IylmxPZdh0dOqlY1wkc2vkgosDuq5iMRjzBjo1AkOOCBMk7/OOmlHJEkqZc+lsobX9KstknPL\nl4dicv31cMcdoc8iUoxiey5fFDy+BdoDLROKSYqU53HfPOcG2chvzhw4/HB47DEYMaJuhSUL+dVH\n3vOLQ7HXudxSuGxmNwNDEolIRFL3/PNw5pnhceWV0KRJ2hFJ1hTVc/nBRmYbEuYGaxV/SPFRz0Wk\nbr76Ci69FJ57Dh5+OPRYpOEpWc/FzN7mux5LY+DHwF/r88UiUl7GjoWTT4Zdd4Vx42DDDdOOSLKs\n2J7LkcBR0aMdsKW735lYVFKUPI/75jk3KK/8li2DXr3CjbwuvxwefbT+haWc8ktC3vOLQ7E9l2lJ\nByIipffhh2FafLMw4WRc94kSWa2eS1ao5yJSM/dw35VLLoEePcIcYY0bpx2VlItSXuciIjnx8cfQ\ntWs4annpJfjZz9KOSPKo2J6LlKE8j/vmOTdIL78nnwwN+513DtO4JFVYtP9ERy4iDcD8+XD++TBy\nJDz1VLixl0iS1HMRybkBA+D3v4fjjw/TuKy9dtoRSblTz0VEVmrBgtCof+UVeOQROPDAtCOShkQ9\nlwzL87hvnnOD5PN79lnYZRdYd114663SFxbtP9GRi0iOfPopdO8eJprU0YqkST0XkRxwD/ex794d\nOnaEa69Vb0VWn3ouIsLMmdCtG7z/fjgTbO+9045IRD2XTMvzuG+ec4N48lu+HO6+G9q2DdeuvPlm\n+RQW7T/RkYtIBk2aBGedBd9+C8OGhea9SDlRz0UkQxYvDteq3HUXXHVVGA7TnGASN/VcRBqQYcPg\n7LNhp53C/VaaN087IpGVU88lw/I87pvn3KBu+c2dC507h6nxb7wxNO3LvbBo/4mKi0iZWr4c7r03\n9FM22QQmTIBjjkk7KpHiqOciUobGjg39FAhnhO26a7rxSMMSR89FRy4iZWThQjjvPGjfHs48E/7z\nHxUWySYVlwzL87hvnnODH+a3fDn07Qs77gjffBOGwM44Axpl9P+hDW3/yQ/pbDGRlI0ZA+eeG6Zw\nee452GOPtCMSqT/1XERSMncuXH55KCjXXx/OCMvqkYrki3ouIhm0dCncemu41fB664Wr7U87TYVF\n8iXx/5zNrL2ZTTKzyWbWo4b3m5pZfzObYmavm1mLgvf+FL0+0czaFbz+oZmNN7OxZjYq6RzKVZ7H\nffOYmzu88AL89KfwyCOVvPIK/P3v0KxZ2pHFL4/7r1De84tDoj0XM2sE3AkcDMwCRpvZM+4+qWC1\nM4D57t7azE4CegEdzWwn4ERgR6A58KKZtY7GuZYDFe6+IMn4ReIyYQJcdBFMnRoKytprh+a9SF4l\n2nMxs72Bq9z98Gj5MsDd/caCdQZF64w0s8bAbHfftPq6ZjYQuDpabyqwh7vPq+X71XORVM2dG+YA\ne/LJ0F/p1g2aNEk7KpFVy0LPZStgesHyjOi1Gtdx92XAZ2a2UQ3bzizY1oHBZjbazM5KInCR+vj6\n69Ck32kn+NGPQl/lggtUWKThSPpU5JoqX/VDiZWts6pt93H3OWb2Y2ComU1093/XFECXLl1o2bIl\nAM2aNaNt27ZUVFQA342bZnX5tttuy1U+hcuFY9rlEE+xy8uWwYwZFVxxBWy7bSW33w4nn5yf/Ipd\nVn7ZWq6srKRv374AK34v683dE3sAewODCpYvA3pUW2cgsFf0vDEwt6Z1gUFV61Xb/irgopV8v+fZ\nsGHD0g4hMVnLbfly94ED3Xfd1X2ffdxfe23V62ctv7pSftkW/XbW6/c/6Z5LY+A9QkN/NjAK6OTu\nEwvW6Qbs4u7dzKwjcIy7VzX0+wF7EYbDhgKtgR8Bjdz9CzNbBxgC9HT3ITV8vyeZnwjAqFHQowfM\nng3XXQfHHgtWr9FqkXSV/f1c3H2ZmZ1LKACNgD7uPtHMegKj3X0A0Ad42MymAPOAjtG2E8zscWAC\nsBTo5u5uZpsBT5mZR/H3q6mwiCRt4kS44goYMSI07U87DdbQnBcigK7Qz7TKysoV46d5U865TZsG\nPXvCgAHwxz+GiSbXXrtun1HO+cVB+WVbFs4WE8mN2bNDIdltN9hqK5g8OQyH1bWwiDQEOnIRqcXc\nueEOkA88EIa+evSATTdNOyqR5OjIRSRBn3wCl14KbdrA4sXwzjtwyy0qLCLFUHHJsMJz7fMmzdzm\nzg1HJzvsAF98AePHw113wZZbxvcded53oPxExUVkhTlz4OKLw5HK55/DuHHQuzdsvXXakYlkj3ou\n0uB99BHcdBP06we/+x1ccgk0b552VCLpUc9FpB4mT4bTT4e2bcP8XxMmwO23q7CIxEHFJcPyPO6b\nZG5vvAEnnAD77gstWsD770OvXrD55ol95Q/ked+B8pPkJ64UKQvuMGRIGP6aPDncW+WBB2DdddOO\nTCSf1HORXPvmG3jsMbj5Zli+PPRTOnaEpk3TjkykfJX93GIiaVmwAO65B+64I5xSfOONcNhhmlBS\npFTUc8mwPI/7rm5u770H554L220XGvQDBsBLL0H79uVVWPK870D5iY5cJAeWL4ehQ8OZXmPGQNeu\n4Wr6OC96FJG6Uc9FMuuzz6Bv33D1/Nprh9sId+oEa62VdmQi2aaeizRI48bBP/4Bjz8e+igPPAD7\n7FNew14iDZ16LhmW53Hf6rl9/TU8/HAoIkcdFaa8f/dd6N8/XK+StcKS530Hyk905CJl7t13w1lf\n/frBnnuGWYqPPFJ3fBQpd+q5SNlZtChcm9KnT7jr4xlnhEfLlmlHJtIwxNFzUXGRsuAOr74aGvRP\nPQUHHghnnhlOIdZRikhpaeLKBi4P477//S/89a/QqhV06wY77QQTJ0L37pW5Hv7Kw75bFeUnOf2/\nrpSz+fPhiSdCg37SpDAdS//+sMce3zXmJ01KN0YRqR8Ni0lJfPVVuFq+Xz+orIR27cK9U9q31zxf\nIuVGPZdaqLika8kSGDw4NOeffz6c7XXyyXDccbD++mlHJyIro55LA1eO476LF8Ozz8Kpp8IWW8At\nt8B++4XOdpSwAAAJbUlEQVRp7ocOhS5diiss5ZhbnJRftuU9vzio5yL1tmgRDBwYzvIaODDc2fH4\n4+GGGzS/l0hDpWExWS1z5sBzz8Ezz8Arr4Qr5487Do4+urR3dBSR+KnnUgsVl/i4hzm9nn8+NObf\ney804zt0CH+bNUs7QhGJi3ouDVzS476ffQZPPglnnQVbbw0nngjz5sHf/gYffwyPPhpOI06isOR9\nTFv5ZVve84uDei6ywrJl8MYb4V7zQ4aEI5V994XDD4eLLw53dBQRKYaGxRow93A1/Esvwcsvw/Dh\nYbbhdu3g0EPhgAPCfVJEpGFRz6UWKi7ft3x5uPJ9+PBwIWNlZbix1sEHh8evfhVOHxaRhk09lwau\ntnHfxYvhtdfg5ptD433TTcN09SNHhqGukSPDrMP33x8ubiynwpL3MW3ll215zy8O6rnkhHsoFCNH\nhseIETB+PLRpA7/8Zbj9b+/eYdhLRCRpGhbLIHeYPh3GjoUxY0IT/o03wgzCe+313WPPPWHdddOO\nVkSyRj2XWuShuCxZEprub70VHuPGhaLSpAn8/OdhJuGqx5ZbZu92vyJSfjLRczGz9mY2ycwmm1mP\nGt5vamb9zWyKmb1uZi0K3vtT9PpEM2tX7Gdm0ZdfwptvwiOPwJVXwm9+E079bdYs9EMGDYJNNoGL\nLoJ33glXyPfoUck114R+ylZb5auw5H1MW/llW97zi0OiPRczawTcCRwMzAJGm9kz7l54t44zgPnu\n3trMTgJ6AR3NbCfgRGBHoDnwopm1BqyIzyxLn38ebo71wQfh75Qp3z3mz4fWrUNBadMmXLDYsyds\nv/3Kp6QfN24cFRUVJc2hVPKcGyi/rMt7fnFIuqH/C2CKu08DMLP+QAegsBB0AK6Knj8B3BE9Pxro\n7+7fAh+a2ZTo86yIzyy5RYtg9myYNQtmzoQZM8Ljo49Co33aNFi6FH7yk/DYdtswweMJJ4SisvXW\n0Lhx3b5z4cKFySRTBvKcGyi/rMt7fnFIurhsBUwvWJ5BKBA1ruPuy8zsMzPbKHr99YL1ZkavWRGf\nuVrc4euvQ6GoeixYAAsXhr/z58Onn4YpUD75BObODY+PPw7bbrFFeDRvHh6tWsFBB8E224THxhvn\na+hKRGRlki4uNf2UVu+wr2ydlb1eU59opV37ww8PP/zLlsG334ajh6VLQ6N88eLw96uvvns0aRLu\nN7LeeuGx4YbhscEGoThsvDFst13of2y2Wbh2ZNNNwzalLhwffvhhab+whPKcGyi/rMt7fnFI9Gwx\nM9sbuNrd20fLlwHu7jcWrDMwWmekmTUGZrv7ptXXNbNBhOEzq+0zCz4726eKiYikpL5niyV95DIa\naGVm2wCzgY5Ap2rrPAd0BkYCJwAvR68/C/Qzs1sJw2GtgFGEI5faPhOo//84IiKyehItLlEP5Vxg\nCKEo9HH3iWbWExjt7gOAPsDDUcN+HqFY4O4TzOxxYAKwFOgWXbRS42cmmYeIiNRNri+iFBGRdGRy\n4sokLswsJ6ubn5ltY2Zfmdmb0aN36aOvXRH57W9mY8xsqZkdV+29ztF275nZqaWLunj1zG9ZtO/G\nmtnTpYu6OEXkdqGZvWtm48xsqJltXfBeHvbdqvIr630HReX3ezN7K8rhFTNrU/Be3X473T1TD0JB\nfB/YBmgCjAPaVFvnD0Dv6PlJhOtlAHYCxhKGA1tGn2Np5xRjftsAb6WdQwz5tQB2AfoCxxW8viHw\nAbAB0Kzqedo5xZVf9N7naedQz9wOBNaKnp9d8N9mXvZdjfmV+76rQ37rFjw/ChgYPa/zb2cWj1xW\nXJjp7kuBqosoC3UAHoyePwH8Knq+4sJMd/8QqLows5ysTn4HF7xX7icx1Jqfu3/k7u/ww1PMDwOG\nuPtn7r6Q0HdrX4qg66A++UF5779ichvu7oujxRGEk3EgP/tuZflBee87KC6/LwoW1wWWR8/r/NuZ\nxeJS04WZ1SeS/96FmUDhhZmF286sYdu0rU5+C6P8AFpGQy7DzGy/xKOtu2LyK3bbrO6/VVnTzEaZ\n2WtmVv0fFWmra25nAANXsm0e9l1hflDe+w6KzM/MupnZ+8ANwPkr2bbW/ZfF+7kkcWFmOVmd/Cxa\nZzbQwt0XmNluwNNmtlO1f42krT77IC/7b1VauPscM9sWeNnM3nL3qTHFVl9F52ZmpwC7E4aR6rRt\niuqTH5T3voMi83P33kBvM+sIXAF0KXbbQlk8cplBGLOu0pwwgWWh6cDWANGFmRu4+4Jo261r2TZt\nq5Pf+u6+wN2/ifLE3d8kjGtvn3zIdVJMfklsWyr1itHd50R/pwKVwM/jDK6eisrNzA4B/gQcFQ2/\nFL1tyuqTX7nvO6j7PngMOKZg27r9dqbdZFqNplRjvmtKNSU0pXastk43vmt4d+SHDf2mwLaUZ0O/\nPvltAjSKnv+EUISapZ1TXfMrWPcB4DcFy4VN4arnecqvGdC0YF++R7WGa7nnRvhBfR/Yrtrrudh3\nq8ivrPddHfJrVfD8KGBU9LzOv52pJ7ya/yO1j3beFOCy6LWewJHR8zWBx6P3RwAtC7b9U/Q/zESg\nXdq5xJkfcBzwTvQfwRvAEWnnspr57UEojIuAT4C3C7btEm03GTg17VzizA/4JfBWtP/GA13SzmU1\nchtKGJ59M8rj6Zztuxrzy8K+KzK/26LfkDeBlygoPnX97dRFlCIiErss9lxERKTMqbiIiEjsVFxE\nRCR2Ki4iIhI7FRcREYmdiouIiMROxUUkZtGtD96u4fV7CqcwF8mzLM4tJpIFNc3Z1DWNQETSoCMX\nkWQ0MbN/mtkEM3vczH4UzVS9G4CZ9Y5m0H3bzK6q2sjMbii4GVWv9MIXqR8duYgkYwfgNHcfYWb3\nEeaDKzya+bO7LzSzRsBLZvYkYRrzY9y9DYCZrV/yqEVioiMXkWR85O4jouf9gOr31uloZmMIc1Ht\nFD0+B742s3vN7Fjg65JFKxIzFReRZFTvuaxYNrOWwB+Bg9x9V+AFwq1zlxHu7vckcCQwqCSRiiRA\nxUUkGduY2V7R807Aq3x3w6X1gS+ARWa2GXA4gJmtTZiGfhBwEfCz0oYsEh8VF5FkTALOMbMJhHt9\n/IPo6MXd3yLcS2Mi8E/g39E26wMDzGw88ApwYamDFomLptwXEZHY6chFRERip+IiIiKxU3EREZHY\nqbiIiEjsVFxERCR2Ki4iIhI7FRcREYmdiouIiMTu/wHNKQE1qD442AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x116ac9090>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# calculation from Josh's presentation page 3: this is the undercoverage vs. bias\n",
    "distribution = scipy.stats.norm(loc = muTrue, scale = sigmaTrue)\n",
    "\n",
    "# integrate over the true function, including our bias\n",
    "undercoverageFunc = lambda thisBias: expectedCoverage - (\n",
    "\n",
    "        # actual coverage\n",
    "        distribution.cdf(+numSigmas + thisBias * numSigmas) \n",
    "        - \n",
    "        distribution.cdf(-numSigmas + thisBias * numSigmas)\n",
    "       )\n",
    "\n",
    "xvalues = np.linspace(0, 0.3, 101)\n",
    "yvalues = [ undercoverageFunc(x) for x in xvalues ]\n",
    "\n",
    "pylab.plot(xvalues, yvalues)\n",
    "pylab.grid()\n",
    "pylab.xlabel('bias')\n",
    "pylab.ylabel('undercoverage')\n",
    "\n",
    "print \"undercoverage for\",numSigmas,\"sigma intervals at bias=\",bias,\"is:\",\\\n",
    "    undercoverageFunc(bias),\"(%.2f%%)\" % (undercoverageFunc(bias) * 100,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "given maximum 0.5% undercoverage, maximum bias is 14.4%\n",
      "given maximum 1.0% undercoverage, maximum bias is 20.4%\n"
     ]
    }
   ],
   "source": [
    "# how much bias could we allow for an undercoverage of 0.5% ?\n",
    "for undercoverage in [\n",
    "    0.005, # 0.5% undercoverage\n",
    "    0.010, # 1.0% undercoverage\n",
    "]:\n",
    "    maxBias = scipy.optimize.brentq(lambda thisBias: undercoverageFunc(thisBias) - undercoverage,\n",
    "                                    0,0.5)\n",
    "    \n",
    "    print \"given maximum %.1f%% undercoverage, maximum bias is %.1f%%\" % (\n",
    "        undercoverage * 100,\n",
    "        maxBias * 100,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## by how much would we have to inflate our sigma (uncertainty of the measurement) to have correct coverage ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class CoverageError:\n",
    "    \n",
    "    def __init__(self):\n",
    "        self.distribution = scipy.stats.norm(loc = muTrue, scale = sigmaTrue)\n",
    "\n",
    "    def __call__(self, bias, uncertScaling):\n",
    "        # inflate our measurement uncertainty\n",
    "        inflatedUncertainty = sigmaTrue * uncertScaling\n",
    "        \n",
    "        # calculate the coverage error but assume we would inflate\n",
    "        # our measurement error\n",
    "        # leave the expected coverage 'as is'\n",
    "        # \n",
    "        # note that we use 'inflatedUncertainty' in the integration boundaries\n",
    "    \n",
    "        actualCoverage = distribution.cdf(+inflatedUncertainty + bias * numSigmas) - \\\n",
    "          distribution.cdf(-inflatedUncertainty + bias * numSigmas)\n",
    "            \n",
    "        return expectedCoverage - actualCoverage\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.009815597636132711"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# now solve for uncertScaling given a bias and requiring that the coverage error is zero\n",
    "coverageError = CoverageError()\n",
    "\n",
    "func = lambda uncertScaling: coverageError(bias, uncertScaling)\n",
    "\n",
    "# get the necessary scaling for the uncertainty given out bias\n",
    "scipy.optimize.brentq(func, 0.5, 1.5) - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "maximum bias for 1.00% inflation of uncertainty is 14.1%\n",
      "maximum bias for 2.00% inflation of uncertainty is 20.0%\n"
     ]
    }
   ],
   "source": [
    "# solve for bias given a maximum allowed inflation of the uncertainty\n",
    "\n",
    "\n",
    "for maxUncertaintyScaling in [\n",
    "    0.01, # 1% inflation of uncertainty\n",
    "    0.02, # 2% inflation of uncertainty\n",
    "    ]:\n",
    "    \n",
    "    func2 = lambda bias: coverageError(bias, 1 + maxUncertaintyScaling)\n",
    "\n",
    "    maxBias = scipy.optimize.brentq(func2, 0, 0.5)\n",
    "    \n",
    "    print \"maximum bias for %.2f%% inflation of uncertainty is %.1f%%\" % (maxUncertaintyScaling * 100, maxBias * 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}