{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Orange Is the New Stat\n", "----------------------\n", "\n", "This notebook contains code that demonstrates the effect of the inspection paradox on observations of a prison population, which I wrote about in this article on my blog, Probably Overthinking It:\n", "\n", "http://allendowney.blogspot.com/2015/08/orange-is-new-stat.html\n", "\n", "Copyright 2015 Allen Downey\n", "\n", "MIT License: http://opensource.org/licenses/MIT" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from __future__ import print_function, division\n", "\n", "import numpy as np\n", "import thinkplot\n", "import thinkstats2\n", "\n", "from collections import Counter\n", "\n", "%matplotlib inline\n", "formats = ['png', 'pdf']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To model the distribution of sentences, I use random values from a gamma distribution, rounded to the nearest integer. All sentences are in units of months." ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Writing orange1.png\n", "Writing orange1.pdf\n" ] }, { "data": { "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sentences = np.random.gamma(shape=2, scale=60, size=1000).astype(int)\n", "cdf1 = thinkstats2.Cdf(sentences, label='actual')\n", "thinkplot.PrePlot(2)\n", "thinkplot.Cdf(cdf1)\n", "thinkplot.Config(xlabel='sentence (months)', ylabel='CDF')\n", "thinkplot.Save('orange1', formats=formats)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I chose parameters that very roughly match the histogram of sentences in this report: http://www.ussc.gov/sites/default/files/pdf/research-and-publications/quick-facts/Quick-Facts_BOP.pdf\n", "\n", "In my distribution about 28% of sentences are less than 5 years, and 40% are more than 10." ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(0.25800000000000001, 0.41599999999999998)" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(sentences < 5*12).mean(), (sentences > 10*12).mean()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we imagine a series of sentences served consecutively, we can compute the release date for each prisoner:" ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "collapsed": false }, "outputs": [], "source": [ "releases = sentences.cumsum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For example, here are the sentences of the first 10 prisoners:" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([193, 69, 148, 107, 123, 189, 61, 219, 104, 83])" ] }, "execution_count": 45, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sentences[:10]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And the release dates of the first 10 prisoners:" ] }, { "cell_type": "code", "execution_count": 46, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 193, 262, 410, 517, 640, 829, 890, 1109, 1213, 1296])" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "releases[:10]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we arrive during month 500, we can figure out which prisoner we would observe using `searchsorted`, which uses bisection search. The result is an array index, which we can think of as a prisoner ID." ] }, { "cell_type": "code", "execution_count": 47, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "releases.searchsorted(500)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following function chooses a random prisoner by choosing a random arrival time, `i`, looking up the prisoner who would be observed at `i`, and then looking up that prisoner's sentence." ] }, { "cell_type": "code", "execution_count": 48, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def random_sample(sentences, releases):\n", " i = np.random.random_integers(1, releases[-1])\n", " prisoner = releases.searchsorted(i)\n", " sentence = sentences[prisoner]\n", " return i, prisoner, sentence" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we call `random_sample` a few times, we get a sample of sentences as seen by random arrivals." ] }, { "cell_type": "code", "execution_count": 49, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(78300, 633, 54)\n", "(37764, 310, 57)\n", "(106242, 871, 110)\n", "(84510, 678, 229)\n", "(38545, 318, 95)\n", "(85459, 684, 201)\n", "(54732, 453, 100)\n", "(77412, 623, 120)\n", "(82512, 667, 212)\n", "(120811, 995, 50)\n" ] } ], "source": [ "for _ in range(10):\n", " print(random_sample(sentences, releases))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can be a little more efficient by generating 1000 random arrival times and finding the corresponding sentences. Then we compute the CDF of the resulting sample." ] }, { "cell_type": "code", "execution_count": 50, "metadata": { "collapsed": false }, "outputs": [], "source": [ "arrivals = np.random.random_integers(1, releases[-1], 10000)\n", "prisoners = releases.searchsorted(arrivals)\n", "sample = sentences[prisoners]\n", "cdf2 = thinkstats2.Cdf(sample, label='biased')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's what the biased sample looks like, compared to the actual distribution. Due to the inspection paradox, we oversample long sentences." ] }, { "cell_type": "code", "execution_count": 60, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Writing orange2.png\n", "Writing orange2.pdf\n" ] }, { "data": { "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "thinkplot.PrePlot(2)\n", "thinkplot.Cdf(cdf1)\n", "thinkplot.Cdf(cdf2)\n", "thinkplot.Config(xlabel='sentence (months)', ylabel='CDF', loc='lower right')\n", "thinkplot.Save('orange2', formats=formats)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As expected, the sample mean is substantially higher than the actual mean:" ] }, { "cell_type": "code", "execution_count": 52, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(121.14100000000001, 182.9058)" ] }, "execution_count": 52, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sentences.mean(), sample.mean()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Instead of using a random sample, we could have computed the biased distribution directly by applying an operation on the PMF. The following function takes the actual distribution of sentences, weights each sentence length, `x`, by `x`, and renormalizes.\n", "\n", "This operation is discussed in Section 3.4 of Think Stats: http://greenteapress.com/thinkstats2/html/thinkstats2004.html#toc29" ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def bias_pmf(pmf, t=0, label=None):\n", " label = label if label is not None else pmf.label\n", " new_pmf = pmf.Copy(label=label)\n", "\n", " for x, p in pmf.Items():\n", " new_pmf[x] *= (x + t)\n", " \n", " new_pmf.Normalize()\n", " return new_pmf\n", "\n", "pmf = thinkstats2.Pmf(sentences)\n", "biased = bias_pmf(pmf, t=0, label='BiasPmf').MakeCdf()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This figure compares the biased distribution computed by random sampling with the distribution computed by `BiasPmf`. They should be similar." ] }, { "cell_type": "code", "execution_count": 54, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEPCAYAAACp/QjLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8VPW57/HPTMIlwUBuEm5KuFW57Crao4AXYguIFgSR\nLV66NwoVTo+49WxPt0Xj2WkPtlp3q9vY46Va0N2DRNlA3JoqVZOKV8QSFDExhIpyMQUJVBGQmDl/\n/GYmayaTIZnMmllr5vt+vfLKrJk1a56JOM/8bs8PREREREREREREREREREREREREREREUtLvgCbg\n/SjnPAA0AFuA8YkISkREkusCzAd+R8nhUqDKf/tc4K1EBCUiIslXTMfJ4WFgnuW4DiiyOyAREYnO\nm+TXHwx8ajneBQxJUiwiIuKX7OQA4Ak79iUlChERCcpM8uvvBk6xHA/x3xdixIgRvsbGxoQFJSKS\nIhqBkbE8MdnJ4VlgCbAKmAAcxMxuCtHY2IjP594GRVlZGWVlZckOI2Zujt/NsYPi74rSilrW7dhP\nSxyvuf+lJymc8o9xvKI9MoGtS6e0u9/j8YzozjXt9BQwGSjEjC38K9DD/9gjmJlKlwLbgcPA9TbH\nIyIJZseHtrTJBGYPL7Tluna6uhPnLLE5BhEhOR/S+zfsoLDX/gS+YscCH6LL5p3Z6eeUHXuNsgjf\nyNNBsruV0kJJSUmyQ+gWN8fv5tgh9vid8m09e/gZCX29WBJANG7/99Md4TOFnMrn5jEHkURxSlII\nF+8Pbekcj8cDMX7Oq+Ug4jLdSQD6kJbOUnIQcbDKqnrKqxvYk9uLVm/XvgAqEUh3qFtJJIns6AZS\nUpAAdSuJuFBpRS2rd8Q2k0cJQOym5CCSBF1JDF6vh6K8LG66cARzTlddSkkMJQeRBDhR99FctQLE\nYZQcROJoTV0TK7fu5UjLNxw8dJTmg0dpbY0+XqbEIE6k5CASR+WvNtLUfOSECQE0biDOpuQgEqOu\nzjRSMhA3UXIQ6aLOJAVvq4/T65vp3TuTebPHMuvS0xIWn0g8aJ2DSCd1tqWgFoI4hdY5iMRZZVU9\nFes+YHd2JvtOzupwdbISgaQqtRxEwpRW1LJm+76o5SqUFMQN1HIQ6YaI3UVqKUiaU3KQtBPoMjp6\n1KSDutPyOmwlFORnMbCgD9eMG6jVyZJWlBwk7ZRXN7BnaI66jUSi0JiDpIUTzTQK1C+qXjQpoXGJ\n2EljDiIRdCYhDBuaC0BWZgbXjBuYuOBEHE7JQVJWtMSQCcwuLlC3kUgHlBwkJZVW1LZLDBpHEOk8\nJQdJSesseyVkAluXTkleMCIu5E12ACLxFt5qmD28MGmxiLiVZitJSgnfYU2tBkln3ZmtpJaDpIxI\nW2+q1SASGyUHSQmVVfWs2b4v5D7tsCYSOyUHcb3Kqnp+vb4+ZMWzEoNI92jMQVwtkBiairKD9w0s\nyNZKZxE05iBpKlJiKMjP4qYLRyQxKpHUoOQgrtRRYrh18khVTxWJAy2CE1ex7tCmxCBiH7UcxDUq\nq+p5YtUWJQaRBNCAtLhCoBvJup+z1+shL7e3EoNIB1SyW1JaR+MLuf16s/DMIUoMIjZQchDHK69u\naJcYtHWniL3sTg7TgfuBDOAx4J6wxwuB3wMD/LH8G7DC5pjERSqr6tmT2yt4rPEFkcSwc0A6A3gQ\nkyDGAFcDo8POWQJsBs4ESoBfodaM+IWvfPZ6PUoMIgliZ3I4B9gOfAwcB1YBs8LO2Qv09d/uC3wO\nHW7eJWkmvDupKC9LiUEkQez8lj4Y+NRyvAs4N+yc3wKvAHuAHOBKG+MRF4nUnaSVzyKJY2dy6Mzc\n09uBWkyX0gjgj8AZwBfhJ5aVlQVvl5SUUFJSEocQxWmsi9xa/a0GdSeJdE5NTQ01NTVxuZad6xwm\nAGWYMQeApUAroYPSVcBdwOv+45eB24BNYdfSOoc0cc2iNe0WuamQnkhsnFp4bxMwCigGegLzgGfD\nzqkDAtt0FQGnATtsjEkcrLKqPiQxeL0edSeJJImd3UotmNlIL2JmLj0OfAgs9j/+CPBzYDmwBZOo\n/gU4YGNM4kCR6iV5vR6GDc3VIjeRJFH5DEmKQEI4erSFz/N7h5TFAK1nEIkHlc8Q1ymvbmDP0JyQ\nhACqlyTiFEoOknCBaarWxBBICiqLIeIMSg6ScOXVDbTm9wbaxhayMjOUFEQcRMlBEip8cVtRXhbP\nX3V2EiMSkUiUHCQhOlrcpmmqIs6kneDEdh3t4KZaSSLOpZaD2Cq4g9tpee12cFOrQcS5lBzEVhXr\nPmCfZcqqdnATcQclB7FNoByGNTFoqqqIOyg5iG2s+zF4vR4GFvRh9dwzkxyViHSGBqTFFuFTVvNy\ne3PNuIFJjEhEukLJQWxRXt2g7T1FXEzJQeIu0kI3JQYRd1FykLgLbzVoyqqI+yg5SNyp1SDifkoO\nElelFbUh1VbVahBxJyUHiZvKqnrWbN8XPM4EtRpEXErJQeLGOtYAMHt4YRKjEZHuUHKQuLGONQws\nyGbZPC14E3ErJQeJi8qqeo01iKQQJQfptkDl1QCv16OxBhGXU3KQbrPWUAIzfVVE3E3JQbolfDV0\nQX6WupREUoCSg8Qs0J2kGkoiqUfJQWIWqTtJiUEkNSg5SMzUnSSSurTZj3RZZVU9Fes+oHVY3+B9\n6k4SSS1KDtJl5dUN7BmaEzzW1FWR1KPkIJ0SaC0cPdrCntPyQgahNXVVJPUoOUinBFoL1lXQXq+H\nMSMKtP2nSApScpCoAi2GSImhKC+L1XNVP0kkFSk5SIcC6xj2WRKD1+th2NBcsjIz1GIQSWFKDtKh\n8HUMgdbC81edncSoRCQRlBwkotKKWnbl9w4eF+RnMbCgj1oLImlCyUHaCe7oprIYImnL7hXS04E6\noAG4rYNzSoDNwFagxuZ45ARKK2pZuvmTkMHnOcUFSgwiacZz4lNilgHUA1OA3cA7wNXAh5ZzcoHX\ngYuBXUAhsD/CtXw+n8/GUCVg3C9eosVyPLAgm+pFk5IWj4jEzuPxQIyf83a2HM4BtgMfA8eBVcCs\nsHOuAf4TkxggcmKQBCmtqG2XGFQvSSQ92TnmMBj41HK8Czg37JxRQA+gGsgB/h34Dxtjkg6UVtSy\nekdbbs4EtRhE0pidyaEz/UA9gLOA7wHZwJvAW5gxCkmQ8AFogNnDC5MYkYgkm53JYTdwiuX4FNq6\njwI+xXQlHfH/vAqcQYTkUFZWFrxdUlJCSUlJXINNZ+XVDbRapq3OHV7Isnla+SziNjU1NdTU1MTl\nWnYOSGdiBqS/B+wBNtJ+QPp04EHMgHQv4G1gHrAt7FoakLbRmLv+GJydpAFokdTRnQFpO1sOLcAS\n4EXMzKXHMYlhsf/xRzDTXF8A3gNagd/SPjGIjUorakOmrWoAWkTA3pZDPKnlYINIg9Bbl05JXkAi\nEldOncoqDrduR+jMYQ1Ci0iAkkMas65p0CC0iFgpOaSpyqr6kGMlBhGxUnJIU+XVbbOFvV63DD2J\nSKIoOaSh8HLc2gNaRMIpOaSZ4GpoP6/Xo+mrItKOkkOaKa9uUDluETkhJYc0UllVz57cXsHjgQXZ\nGogWkYiUHNJEZVU9v15fH2w1qDtJRKJRckgDgcTQVJQdvK8oL0vdSSLSoWjJYYXl9nyb4xAblVc3\nhCSGgvwstRpEJKpoyeEMy+1b7A5E7BE+zlCQn8Wtk0eq1SAiUalbKcVZZyd5vR4lBhHplGglu4cA\nD2Aq+g223Aazy9s/2RuaxIO11aBxBhHprGjJ4ceYJOAB3g17TPWzXaCyql57NYhITKIlhxWJCkLs\nUV7dAP4yGV6vR60GEem0E405XAf8GfjK/7MJzVxyBdVPEpHuiNZymA/cDPwzsBnTvTQeuBfTrfSk\n7dFJTIL1k7TgTURiFK3l8D+AOUA1cBBoBl4BrgButD80iZXqJ4lId0VLDjnAXyLc/7H/MXGg8O4k\n1U8SkVhESw5HY3xMksi6L7S6k0QkVtHGHEYD73fwmD5xHKiyqj5kX2h1J4lIrKIlh28DRcCusPtP\nAfbaFpHELHzqqrqTRCRW0bqV7gcOYcYYrD+HgPvsDUtiEb4aWkQkVtGSQxGRu5XeA4bZE47EqrSi\nVquhRSRuoiWH3CiP9Y7ymCRY+L7QmaCxBhHplmjJYROwKML9N9C+1pIkUfi6htnDC5MYjYikAk+U\nxwYAa4GvaUsGZwO9gMtJ7KC0z+dTrb9ISitqWW2ZvjqwIJvqRZOSGJGIOIXH44Hon/MdijZb6TNg\nEnARMA5TMuM5zCppcQitaxARO0RLDmASwisoITiW1jWIiB20E1wK0boGEYkXJQcXK62oTXYIIpKi\nlBxcKtL0VRGReFFycClNXxUROyk5uFBlVX1IqQyV5RaReFNycJnKqnp+vb4+2GrQ9FURsYPdyWE6\nUAc0ALdFOe+/YWZlzrE5Htcrr26gqSg7eFyUl6XpqyISd3YmhwzgQUyCGANcjdkjItJ59wAvEONK\nvnQR3p1UkJ+lVoOI2MLO5HAOsB1T5vs4sAqYFeG8m4DVwL4Ij4mFdRDa6/Vw6+SRajWIiC3sTA6D\ngU8tx7v894WfMwt4yH+sAkpRhO/XoMQgInaxMzl05oP+fuAn/nM9qFupQ9qvQUQSyc61U7sxW4oG\nnEL7LUfPxnQ3ARQCl2C6oJ4Nv1hZWVnwdklJCSUlJfGL1AWsBfa0X4OIRFJTU0NNTU1crmXnN/VM\noB74HrAH2IgZlP6wg/OXA/8FrInwWFqX7K6sque2LW09dHOHF2pdg4ickF0lu7urBVgCvIiZkfQ4\nJjEs9j/+iI2vnVLKqxsg32y+5/V6lBhExHZ2l+T5g//HqqOkcL3NsbhS+PTVorysJEYjIulCK6Qd\nLnz6qgaiRSQRlBwcTtNXRSQZlBwcrLKqXtNXRSQplBwcKlBgL8Dr9ajVICIJo+TgUJEK7ImIJIqS\ngwOpwJ6IJJuSgwOpwJ6IJJuSg8NEWtegxCAiiabk4DBa1yAiTqDk4CBqNYiIUyg5OIhaDSLiFEoO\nDqJWg4g4hZKDQ2g1tIg4iZKDQ5RXNwRvazW0iCSbkoMDqCy3iDiNkkOSBWooaSBaRJxEySHJItVQ\nUpeSiCSbkkMSqYaSiDiVkkMSqYaSiDiVkkOSaDW0iDiZkkOSaDW0iDiZkkMSqNUgIk6n5JAEajWI\niNMpOSSYWg0i4gZKDgmmVoOIuIGSQwKp1SAibqHkkEAV6z5Qq0FEXEHJIUEqq+rZnZ0ZPM7L7a1W\ng4g4lpJDglhrKHm9HgYW9ElyRCIiHVNySIDwsYa83N5cM25gEiMSEYlOycFmkUpyq4aSiDidkoON\nAolBJblFxG2UHGwUvleDSnKLiFsoOdgk0l4N6k4SEbdQcrCJ9moQETdTcrCBVkKLiNslIjlMB+qA\nBuC2CI9fC2wB3gNeB76dgJhspZXQIuJ2mSc+pVsygAeBKcBu4B3gWeBDyzk7gAuBQ5hE8igwwea4\nbKOV0CJGfn4+zc3NyQ4jLeTl5XHgwIG4XtPu5HAOsB342H+8CphFaHJ403L7bWCIzTHZqmLdB+wb\nmgNoJbSkt+bmZnw+X7LDSAsejyfu17S7W2kw8KnleJf/vo4sBKpsjchmu7Mzg11KWgktIm5ld8uh\nK18bLgIWAOdFerCsrCx4u6SkhJKSku7EZYvKqnr2nZwVPB5Y0EddSiKSMDU1NdTU1MTlWvFvi4Sa\nAJRhxhIAlgKtwD1h530bWOM/b3uE6/ic3jwNXw3t9XpYNmuckoOkLY/Ho26lBOnob+3vborpc97u\nbqVNwCigGOgJzMMMSFudikkMPyByYnCFinUfhLQaNH1VRNzM7uTQAiwBXgS2ARWYwejF/h+A/w3k\nAQ8Bm4GNNscUd4EZSoGxBpXJEHGu4uJiXn755Xb3b9iwgdNPPz2hsaxYsYILLrggoa/ZWXaPOQD8\nwf9j9Yjl9g/9P64Vaa8GtRpEnMnj8USc3XPBBRdQV1eXhIicSSuk40B7NYhIqlFy6KbKqra9GgDV\nUBJxgY0bNzJ27Fjy8/NZsGABx44do6amhlNOOSV4zt13383IkSPp27cvY8eOZd26dcHHtm/fzuTJ\nk8nNzeXkk0/mqquuCj5WV1fH1KlTKSgo4PTTT+eZZ54JPvb5559z2WWX0a9fP84991waGxsT84Zj\nkIhupZRWXt0A+b0B06WkxCByYnP+8em4Xm/Nk1d2+lyfz8fKlStZv3492dnZzJw5k2XLljFlypSQ\n80aOHMlrr73GgAEDePrpp/nBD35AY2MjRUVF3HnnnUyfPp0//elPfP3112zatAmAw4cPM3XqVJYt\nW8aLL77Ie++9x9SpUxk3bhyjR4/mxhtvJDs7m88++4wdO3Zw8cUXM3z48Lj+LeJFLYduiFRgT0Sc\nzePxsGTJEgYPHkxeXh533HEHTz31VLvz5s6dy4ABAwC48sorGTVqFBs3mvkyPXv25OOPP2b37t30\n7NmTSZMmAfDcc88xbNgw5s+fj9fr5cwzz2TOnDk888wzfPPNN6xZs4af/exnZGVlMXbsWObPn+/Y\n6b5KDt0QXpZbM5RE3MHafXTqqaeyZ8+educ8+eSTjB8/nry8PPLy8ti6dSv79+8H4Je//CU+n49z\nzjmHcePGsXz5cgB27tzJ22+/HXxOXl4eK1eupKmpif3799PS0tLutZ1K3UoxUllukdh1pRvIDp98\n8knI7UGDBoU8vnPnThYtWsQrr7zCxIkT8Xg8jB8/Pvgtv6ioiEcffRSA119/nSlTpnDhhRdy6qmn\nMnnyZNavX9/uNb/55hsyMzP55JNPOO2009rF4TRqOcRIrQYRd/L5fPzmN79h9+7dHDhwgLvuuitk\nQBnM2IHH46GwsJDW1laWL1/O1q1bg48/88wz7Nq1C4Dc3Fw8Hg8ZGRnMmDGDjz76iN///vccP36c\n48eP884771BXV0dGRgZz5syhrKyMI0eOsG3bNp544glbiubFg5JDDNRqEHEvj8fDtddey7Rp0xgx\nYgSjRo2itLQUn88X/KAeM2YMt956KxMnTmTAgAFs3bqV888/P3iNTZs2MWHCBHJycpg1axYPPPAA\nxcXFnHTSSaxfv55Vq1YxePBgBg4cyNKlS/n6668BePDBB/nyyy8ZMGAACxYsYMGCBUn5G3SGM1NW\ne46qrTTlx8+xyzJDSTWURNpTbaXEcWNtpZRTWlEbTAygVoOIpCYlhy5at2N/8LbGGkQkVSk5dEFp\nRS0tluM5xQVqNYhISlJy6AJrqyETWDbvzOQFIyJiIyWHTgpvNcweXpi0WERE7Kbk0AmVVfWs2b4v\neKxWg4ikOiWHEwhs/2mtvKpWg4ikOiWHE7Bu5AMwsCBbrQYRSXlKDlGEr4TW9p8iqe1HP/oRy5Yt\nS3YYIa6//nry8/OZMGFCQl9XyaED4d1JXq9HG/mIpIDi4mKys7PJyckhPz+fGTNmBOskPfTQQ5SW\nlnbr+itWrCAjI4OcnBz69evH+PHjef7552O61oYNG3jppZfYs2cPb731Vrfi6iolhw6EdydpJbRI\navB4PDz33HN88cUX7N27l6KiIm666aa4vsZ5553HF198wcGDB1m4cCFXXnklhw4d6vJ1du7cSXFx\nMb179z7xyXGm5BBBeIkMdSeJpKZevXpxxRVXsG3bNgCuu+467rzzTgCam5uZMWMG/fv3Jz8/n5kz\nZ7J79+7gc1esWMGIESPo27cvw4cPZ+XKlcHHAnWOPB4P119/PUeOHKGxsZGamhqGDBnCvffeS//+\n/Rk0aBDr1q2jqqqKb33rWxQUFHD33XcD8Pjjj3PDDTfw5ptvkpOTw09/+tNE/VkA7ecQUXiJDHUn\nicTX91e9G9frPX/V2V06P/Dh/dVXX1FRUcHEiRMB82EeqMzq8/lYuHAhq1evpqWlhQULFrBkyRLW\nrl3L4cOHufnmm9m0aROjRo2iqamJzz//vN3rtLS08Nhjj5GTk8OoUaN49913aWpq4tixY+zdu5fl\ny5fzwx/+kIsvvpjNmzezc+dOvvOd73D11VezcOFCMjMzeeyxx9iwYUM3/0Jdp+QQRiUyRFKbz+dj\n9uzZZGZmcvjwYfr3788LL7wQ8jhAfn4+l19+efD+22+/ne9+97vBY6/Xy/vvv8+QIUMoKiqiqKjt\nc+Ktt94iLy+PzMxMRo0axdq1a8nJyQGgR48e3HHHHXg8HubNm8eiRYu45ZZb6NOnD2PGjGHMmDFs\n2bKFoUOHJrWqrbqVLEoralmtEhkiKc3j8VBZWUlzczPHjh2jvLycyZMn09TUFHLeV199xeLFiyku\nLqZfv35MnjyZQ4cO4fP56NOnDxUVFTz88MMMGjSIGTNmUF9fH3zuhAkTaG5uZt++fbzxxhshSaWg\noCDYOsnKMvvOWxNLVlYWX375pZ1/gk5Ry8EvPDGAFruJ2KWr3UB28Xg8XH755SxevJjXXnsteB/A\nr371Kz766CM2btxI//79qa2t5ayzzgpuCjRt2jSmTZvGsWPHuOOOO7jhhht49dVXk/l24kotByIn\nhrnDC9VqEElRge4an89HZWUlBw8eZMyYMfh8vuBjX375JVlZWfTr148DBw6EDAj/9a9/pbKyksOH\nD9OjRw/69OlDRkZGUt6LXdI+OSgxiKSfmTNnBtch3HnnnTzxxBOMHj06ZED6lltu4ciRIxQWFjJp\n0iQuueSS4GOtra3cd999DB48mIKCAjZs2MBDDz0EhA5qRxL+2InOTdYe02m9TagSg4h9tE1o4mib\n0DhSYhAR6VhaJgclBhGR6NIuOSgxiIicWFolByUGEZHOSZvkEL6bGygxiIh0JG2SQ3l1Q8hubkoM\nIiIdS/kV0qUVtazbsZ8WS5VV7eYmYr+8vLykzdFPN3l5eXG/pt3JYTpwP5ABPAbcE+GcB4BLgK+A\n64DN8XrxSGMMXq9H5bdFEuDAgQPJDkG6wc5upQzgQUyCGANcDYwOO+dSYCQwClgEPNTdF11T18RF\nj77BmHtebpcYMklOldWampqEvl68uTl+N8cOij/Z3B5/d9iZHM4BtgMfA8eBVcCssHMuA57w334b\nyAVi+uQurahl3C9e4va177P3869obQ1dLTh3eCFbl05JSneS2/+BuTl+N8cOij/Z3B5/d9jZrTQY\n+NRyvAs4txPnDAGaOIHgWMIJzsvEVFfVGIOISOfZmRw6W1QlfMQq4vO+v+pdGv/S3KkLelt9DDp4\njJsuGsWsS0/rZBgiIhJg51SCCUAZZswBYCnQSuig9MNADabLCaAOmEz7lsN2QKPIIiJd04gZ13WU\nTExgxUBPoJbIA9JV/tsTgLcSFZyIiCTPJUA95pv/Uv99i/0/AQ/6H98CnJXQ6EREREREJDVMx4xD\nNAC3JTmWjvwOM0byvuW+fOCPwEfAeswU3YClmPdTB0xLUIzRnAJUAx8AW4F/8t/vlvfQGzMNuhbY\nBvzCf79b4gezJmgz8F/+YzfF/jHwHib+jf773BR/LrAa+BDz7+dc3BP/aZi/e+DnEOb/X7fEH7MM\nTHdTMdCDyGMWTnABMJ7Q5PBL4F/8t28D7vbfHoN5Hz0w72s7ya9vNQAIzPM9CdMNOBp3vYds/+9M\nzLjV+bgr/n8G/h/wrP/YTbH/BfNhZOWm+J8AFvhvZwL9cFf8AV5gL+bLnhvj75KJwAuW45/4f5yo\nmNDkUEfbYr4B/mMwWdvaAnoBMxDvJOuAKbjzPWQD7wBjcU/8Q4CXgItoazm4JXYwyaEg7D63xN8P\n2BHhfrfEbzUN2OC/HZf4nZw1Ii2QG5ykWLqqiLbpuE20/YcahHkfAU57T8WYVtDbuOs9eDHfiJpo\n6yJzS/z3AT/GTPMOcEvsYNYlvQRsAm7w3+eW+IcB+4DlwJ+B3wJ9cE/8VlcBT/lvxyV+JyeHVNmZ\n3Ef09+KU93kS8J/AzcAXYY85/T20YrrGhgAXYr6FWzk1/hnAXzH9xR2tOXJq7AHnYb5QXALciOlm\ntXJy/JmYGZL/1//7MO17J5wcf0BPYCbwTITHYo7fyclhN6b/LOAUQrOekzVhmnMAAzEfAND+PQ3x\n35dsPTCJ4T8w3UrgvvcAZkDueeBs3BH/JEx9sb9gvvV9F/PfwA2xB+z1/94HrMXUVHNL/Lv8P+/4\nj1djksRnuCP+gEuAdzH/DcA9f/+YdWYRnVMU035AOtC39xPaDwj1xDRpG7F3lXpneIAnMd0bVm55\nD4W0zcbIAl4Fvod74g+YTNuYg1tizwZy/Lf7AK9j+r7dEj+Yfy/f8t8uw8TupvjBVJiYbzl2W/wx\nibSIzmmeAvYAX2PGSK7HzN54ichTyW7HvJ864OKERhrZ+ZhumVrapsRNxz3v4e8w/cW1mCmVP/bf\n75b4AybTNlvJLbEPw/zdazHToAP/j7olfoAzMC2HLcAazCC1m+LvA+ynLUmDu+IXERERERERERER\nERERERERERERERERSXVnYNbHJEt/zCptO4W/xzLg1i48/2nMGgWRDjm5fIZILMZjtp9NliXACptf\nI/w9drW+z2+B/xm/cERE7NMH8427FlOK5Er//WcDNZiqny/QVjOmBlMW4G3MKvrzMTWiPqGtmN3f\n+6/7O/95f8bUMgK4DrMi9g+YlaT3WGKZjqlVU4tZaUqU64TbRtv+EtdhalWtx9RPWgL8L//z3wTy\n/OedidmHIrBKN7CitTPv8UrgX4HHMRVpG4GbLDFH+pv2wKySFRFxvCuARy3HfTEfYm/Qtl/APMyH\nIJgPwnv9ty/B7HwFpsbMA5br/By41n87F/Mhm4354G7ElB3ohdnRbDBwMubDd6jlOdGuYzWA0Bpb\n12F23eqDqQF1CFjkf+zXmAq4YMp+BKqZ/pS2OledfY9lmLpGPTB/q/2Y2mSR/qYBf8K5tcrEAdSt\nJE7xHjAV8035fOBvmG0Qx2K+vW8G7iC0/vwa/+8/Y4ofgikkZi0mNg1TfGwz5sO2F3AqpivmZUx5\n8mOYb/zFmM1PXgV2+p9/MMp1rBUuwSSUvZZjn//cw5gP7IO0Fdd73/96fTH1fAIbtTyBKTvelffo\nA54DjgMTPU/lAAABcklEQVSfY1oV/Yn8Nw3YY7meSDuZyQ5AxK8B05f+fWAZ5oN7LWbjnkkdPOeY\n//c3RP+3PMd/fatzLc+3XiNa/32k64QLr3JpfY1Wy3ErkWPu6Pkneo9fW24Hzo30N/0/ltexbjAk\nEkItB3GKgcBRzF7K/4b5UKvHdPMEtjLsgSk7HM3fCK1Q+SJm0/WA8f7fkUoV+zB9/xfS9q06sD9y\nR9ex2knbmEhHrxH+2N+AZsw3e4B/wIw1RBP+HjsS/jc9K+yxnZGeJAJqOYhz/B2mf70V0z3y3/2/\n52L61/th/r3eh+kCChf4xl9NW/fPzzHflO/HdLF4MXsGX0bHO2Ttx4wLrPGf34QpbdzRdaw+88eY\nDXwV4TXCbweO5wMP+5/XiCn7Hkmk9/iLCNcOsP5NvwZ+5L+/B2ajl7oIzxERERuUYQbOnWwa8O/J\nDkJEJJ2cDFQlO4gTeBoNRouIiIiIiIiIiIiIiIiIiIiIiIiIiIgz/X8/2+WC/SyynwAAAABJRU5E\nrkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "thinkplot.PrePlot(2)\n", "thinkplot.Cdf(cdf2)\n", "thinkplot.Cdf(biased)\n", "thinkplot.Config(xlabel='sentence (months)', ylabel='CDF', loc='lower right')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we get to the scenario I describe in the article: observation made by a prisoner serving a sentence with duration `t`.\n", "\n", "The following function computes this distribution, given an array of (actual) sentences and the sentence of the observer.\n", "\n", "It starts by choosing a random arrival time between 1 and the maximum sentence. So the initial observation is subject to the inspection paradox.\n", "\n", "Each time through the loop it computes `first_prisoner`, which is the prisoner observed on arrival, and `last_prisoner`, the prisoner observed on departure. The prisoners between these endpoints (including both) are considered one observed sample, and added to the Counter object, which accumulates a histogram of observed sentences.\n", "\n", "Each time through the loop, the arrival time is advanced by 100, which is approximately the average sentence. This has the effect of thinning the observations.\n", "\n", "Finally, the function returns the CDF of the observed sentences.\n" ] }, { "cell_type": "code", "execution_count": 55, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def simulate_sentence(sentences, t):\n", " counter = Counter()\n", " \n", " releases = sentences.cumsum()\n", " last_release = releases[-1]\n", " arrival = np.random.random_integers(1, max(sentences))\n", " \n", " for i in range(arrival, last_release-t, 100):\n", " first_prisoner = releases.searchsorted(i)\n", " last_prisoner = releases.searchsorted(i+t)\n", " \n", " observed_sentences = sentences[first_prisoner:last_prisoner+1]\n", " counter.update(observed_sentences)\n", " \n", " print(sum(counter.values()))\n", " return thinkstats2.Cdf(counter, label='observed %d' % t)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following function takes an observed distribution of sentences and plots it along with the actual distribution and the biased distribution that would be seen by a single arrival." ] }, { "cell_type": "code", "execution_count": 56, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def plot_cdfs(cdf3):\n", " thinkplot.PrePlot(2)\n", " thinkplot.Cdf(cdf1)\n", " thinkplot.Cdf(cdf2)\n", " thinkplot.Cdf(cdf3, color='orange')\n", " thinkplot.Config(xlabel='sentence (months)', ylabel='CDF', loc='lower right')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If the sentence of the observer is only 11 months, the observed distribution is almost as biased as what would be seen by an instantaneous observer." ] }, { "cell_type": "code", "execution_count": 67, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1341\n", "Writing orange3.png\n", "Writing orange3.pdf\n" ] }, { "data": { "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "cdf13 = simulate_sentence(sentences, 13)\n", "#cdf13 = bias_pmf(pmf, 13).MakeCdf()\n", "plot_cdfs(cdf13)\n", "thinkplot.Save('orange3', formats=formats)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With a longer sentence, we get a less biased view. But even after 120 months, which is the average sentence, the observed sample is still quite biased." ] }, { "cell_type": "code", "execution_count": 68, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2397\n", "Writing orange4.png\n", "Writing orange4.pdf\n" ] }, { "data": { "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "cdf120 = simulate_sentence(sentences, 120)\n", "#cdf120 = bias_pmf(pmf, 120).MakeCdf()\n", "plot_cdfs(cdf120)\n", "thinkplot.Save('orange4', formats=formats)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After 600 months (50 years!) the observed distribution almost reaches the actual distribution." ] }, { "cell_type": "code", "execution_count": 69, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "7166\n", "Writing orange5.png\n", "Writing orange5.pdf\n" ] }, { "data": { "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "cdf600 = simulate_sentence(sentences, 600)\n", "#cdf600 = bias_pmf(pmf, 600).MakeCdf()\n", "plot_cdfs(cdf600)\n", "thinkplot.Save('orange5', formats=formats)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We conclude that during Kerman's 11 month sentence, she would have seen a biased sample of the distribution of sentences.\n", "\n", "Nevertheless, her observation that many prisoners are serving long sentences that do not fit their crimes is still valid, in my opinion." ] }, { "cell_type": "code", "execution_count": 71, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Writing orange6.png\n", "Writing orange6.pdf\n" ] }, { "data": { "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "thinkplot.PrePlot(cols=3)\n", "plot_cdfs(cdf13)\n", "thinkplot.SubPlot(2)\n", "plot_cdfs(cdf120)\n", "thinkplot.SubPlot(3)\n", "plot_cdfs(cdf600)\n", "thinkplot.Save('orange6', formats=formats)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "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.10" } }, "nbformat": 4, "nbformat_minor": 0 }