{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Traffic collisions in Los Angeles County\n", "\n", "This is a Jupyter Notebook for analyzing the traffic collision data available for Los Angeles County through TIMS. The analysis uses the Pandas package in Python:\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import pandas" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Analyzing data\n", "\n", "Here I will analyze data about traffic collisions in Los Angeles County. First, let's load the data set into pandas and take a look at it. This particular data set is for the collisions which occurred in Los Angeles County in January 2012." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(4146, 83)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#Load the data, show data size\n", "\n", "data = pandas.read_csv('Jan12Collisions.csv')\n", "data.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This means there are 4146 collisions recorded (instances) and 83 observations (features) for each collision." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
CASEIDPOINT_XPOINT_YYEAR_LOCATIONCHPTYPEDAYWEEKCRASHSEVVIOLCATKILLED...BICINJMCKILLMCINJURERAMP1RAMP2CITYCOUNTYSTATEX_CHPY_CHP
064138410.0000000.00000020121942023090...000--LOS ANGELESLOS ANGELESCA0.000000.00000
16348475-118.20860933.87349820121900314030...000--UNINCORPORATEDLOS ANGELESCA-118.2082133.87339
26216712-118.21056933.92539720121943144030...000--LYNWOODLOS ANGELESCA-118.2114833.92575
36071228-118.53604034.16891320121942054000...000--LOS ANGELESLOS ANGELESCA0.000000.00000
46063954-118.16833334.11479020121970044180...100--SOUTH PASADENALOS ANGELESCA0.000000.00000
\n", "

5 rows × 83 columns

\n", "
" ], "text/plain": [ " CASEID POINT_X POINT_Y YEAR_ LOCATION CHPTYPE DAYWEEK CRASHSEV \\\n", "0 6413841 0.000000 0.000000 2012 1942 0 2 3 \n", "1 6348475 -118.208609 33.873498 2012 1900 3 1 4 \n", "2 6216712 -118.210569 33.925397 2012 1943 1 4 4 \n", "3 6071228 -118.536040 34.168913 2012 1942 0 5 4 \n", "4 6063954 -118.168333 34.114790 2012 1970 0 4 4 \n", "\n", " VIOLCAT KILLED ... BICINJ MCKILL MCINJURE RAMP1 RAMP2 \\\n", "0 09 0 ... 0 0 0 - - \n", "1 03 0 ... 0 0 0 - - \n", "2 03 0 ... 0 0 0 - - \n", "3 00 0 ... 0 0 0 - - \n", "4 18 0 ... 1 0 0 - - \n", "\n", " CITY COUNTY STATE X_CHP Y_CHP \n", "0 LOS ANGELES LOS ANGELES CA 0.00000 0.00000 \n", "1 UNINCORPORATED LOS ANGELES CA -118.20821 33.87339 \n", "2 LYNWOOD LOS ANGELES CA -118.21148 33.92575 \n", "3 LOS ANGELES LOS ANGELES CA 0.00000 0.00000 \n", "4 SOUTH PASADENA LOS ANGELES CA 0.00000 0.00000 \n", "\n", "[5 rows x 83 columns]" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# What does the data look like? limit to first five rows\n", "\n", "data.iloc[:5,:]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Column labels?\n", "\n", "The many columns of this data table are labeled in ways that aren't entirely clear. One can refer to the SWITRS codebook.\n", "\n", "This allows us to interpret the meaning of each column and their codes.\n", "\n", "Next, let's look at which hours of the day have the most collisions recorded." ] }, { "cell_type": "code", "execution_count": 163, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "17 357\n", "18 333\n", "15 332\n", "16 268\n", "14 254\n", "13 242\n", "8 231\n", "12 222\n", "11 219\n", "7 212\n", "19 209\n", "10 208\n", "9 185\n", "20 141\n", "22 112\n", "21 112\n", "6 97\n", "2 78\n", "5 72\n", "1 70\n", "23 54\n", "0 51\n", "3 49\n", "4 34\n", "25 4\n", "Name: TIME_, dtype: int64\n" ] } ], "source": [ "HourCounts = (data.TIME_//100).value_counts()\n", "print HourCounts\n", "\n", "#There are apparently 4 instances where time is misrecorded as hour 25, let's remove these\n", "HourCounts = HourCounts[0:24]" ] }, { "cell_type": "code", "execution_count": 150, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [], "source": [ "#Now put the data in order by hour of day\n", "\n", "OrdHrCnts = [HourCounts.loc[i] for i in range(24)]\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Number of crashes in each hour of the day \n", "\n", "Let's now plot the number of crashes during each hour of the day for January 2012 in LA County. The rush hours of 6-9am and 3-6pm are highlighted." ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEPCAYAAAC3NDh4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucVHX9x/HXh6tykYtyE5RLKCGCQGYqmGsqoiZoBiJa\nivqLRH+hpqmZsuSv0Az7qUVWKqGpiKiB/lQQdUlQwOIaIJJCIspK4g0vCOzn98c5I8O6szvDzpkz\nM/t+Ph7z4MyZc/kwO7uf+d7N3REREalKvbgDEBGR/KUkISIiKSlJiIhISkoSIiKSkpKEiIikpCQh\nIiIp5SRJmFk9M1tsZjPD563MbLaZrTGzWWbWIunYa81srZmtNrNBuYhPRESqlquSxFhgVdLza4A5\n7t4DeA64FsDMDgGGAz2Bk4FJZmY5ilFERCqJPEmYWSfgFOCupN1DgSnh9hTg9HB7CDDV3Xe4+3pg\nLXBE1DGKiEjVclGS+A1wFZA8tLudu5cDuPsmoG24vyOwIem4jeE+ERGJQaRJwsxOBcrdfSlQXbWR\n5gYREclDDSK+/gBgiJmdAuwNNDez+4BNZtbO3cvNrD3wTnj8RuCApPM7hft2Y2ZKKiIie8DdM2rn\ntVxN8GdmxwI/dvchZvYr4F13v9nMrgZaufs1YcP1/cA3CKqZngEO8kpBmlnlXXVWaWkppaWlGZ+3\nbl3m5+RS166lGZ+zp+9FMRo7toTLLiuJO4yU9uTnu6f0udjFzDJOElGXJFK5CZhmZhcA/ybo0YS7\nrzKzaQQ9obYDY5QNRETik7Mk4e5zgbnh9hbghBTHTQAm5CouERFJTSOuC1xJSUncIeQNvRe7HHlk\nl7hDyBv6XNSOkkSB0y/ALnovdlGS2EWfi9pRkhARkZSUJEREJCUlCRERSUlJQkREUlKSEBGRlJQk\nREQkJSUJERFJSUlCRERSUpIQEZGUlCRERCQlJQkREUlJSUJERFJSkhARkZSUJEREJCUlCRERSUlJ\nQkREUlKSEBGRlCJNEmbW2MwWmtkSM1thZuPC/ePM7E0zWxw+Biedc62ZrTWz1WY2KMr4RKR27rvv\n6/z854OpqLC4Q5GINIjy4u6+zcyOc/dPzKw+MN/MngpfvtXdb00+3sx6AsOBnkAnYI6ZHeTuHmWc\nIrJnZszow6ZNzdm2rQE33vh/1KunX9ViE3l1k7t/Em42JkhKiU9RVV89hgJT3X2Hu68H1gJHRB2j\niGTuww/34tVX2zJz5h959dW2XH/9qSpRFKHIk4SZ1TOzJcAm4Bl3fzl86VIzW2pmd5lZi3BfR2BD\n0ukbw30ikmdeeqkL/fptoHXrT5g8+S9KFEUqFyWJCnfvR1B9dISZHQJMArq5e1+C5DEx6jhEJLvm\nzfsKAwe+DkCzZp8rURSpSNskkrn7h2ZWBgyu1BbxJ+DxcHsjcEDSa53CfV9SWlr6xXZJSQklJSVZ\njFZEajJvXjcmTZr2xfNEohg16lyuv/5UtVHkgbKyMsrKymp1DYuyTdjM9gO2u/sHZrY3MAu4CVjs\n7pvCYy4Hvu7uI8NSxv3ANwiqmZ4BvtRwbWZqy66ldetK4w6hWl27lsYdQkGL+uf75pstOeOMi1i4\ncOKXEsHWrY0YNepcDj74nZSJQj/feJgZ7p5RMS/q6qYOwPNmthRYCMxy9yeBX5nZ8nD/scDlAO6+\nCpgGrAKeBMYoG4jkn3nzujFgwOtVJgBVPRWXSJOEu69w9/7u3tfd+7j7L8L93w+f93X30929POmc\nCe7e3d17uvvsKOMTkT2T3B5RFSWK4qER1yKSkZ07jRdf7MrAga9Ve5wSRXFQkhCRjKxc2YH99vuY\n9u0/qvFYJYrCpyQhIhmZN69bjaWIZEoUhU1JQkQyMn/+VxgwIHV7RFW+nCgiCk6yTklCRNL26acN\nWb58f77xjfUZn5ucKC67LPuxSTSUJESKTFkZXHnl6ZFce9GizhxyyCaaNft8j85v1uxz7rrrAe66\nC7Zvz3JwEgklCZEic++98Ne/BrOzZtu8ed045pj02yOq0qLFZ3TuDKtXZykoiZSShEgR2bkTnngC\nBgx4ncceOyzr1w/GR9QuSQD07w+LF2chIImckoRIEVmwADp0gLFjy3jkkb5kc76CzZubsWnTPvTu\n/Vatr9Wvn5JEoVCSECkiM2bA0KHQr9+buMPSpZ2ydu3587tx5JHrqF+/9pmnf39YsiQLQUnklCRE\nisiMGTBkCJjBmWcuZfr0vlm7djA+IrOur6n07QtLl6KusAVASUKkSKxZAx9/DF/7WvD8jDOW89RT\nh/DZZ7VfEcA9e+0RAK1bw377wb/+lZXLSYSUJESKRHIpAqBDhw/p3fstnnnmq7W+9tq1bWjceAed\nO79X62slqF2iMChJiBSJRHtEsjPPXMojj9S+yimbpYgEtUsUBiUJkSJQXg4rV0LlBRoHDXqFZcs6\n1nrMRGL9iGxSN9jCoCQhUgSeeAIGDYLGjXffv9deOzjllJW1GjPx+ef1+fvfD+Too9fVMsrd9esX\nlCS0rFh+U5IQKQIzZ365qikhUeW0p3+MlyzpRLdu79Ky5ad7HmAVOnSAhg1hw4asXlayTElCpMB9\n8kkwX9Mpp1T9em3HTETRHpGgxuv8pyQhUuCeeQYOPxxatar69dqOmcjm+IjK1Hid/yJNEmbW2MwW\nmtkSM1thZuPC/a3MbLaZrTGzWWbWIumca81srZmtNrNBUcYnUgwSXV+rs6djJt5/f29ee60N/fpF\nUyekkkT+izRJuPs24Dh37wf0BU42syOAa4A57t4DeA64FsDMDgGGAz2Bk4FJZqZlrERSSEzol6o9\nImFPx0y8+GJXvva1N2jceGctokxNJYn8F3l1k7t/Em42BhoADgwFpoT7pwCJye+HAFPdfYe7rwfW\nAkdEHaNIoUpM6NelS83H7smYifnzM1uqNFNdugSjxN95J7JbSC1FniTMrJ6ZLQE2Ac+4+8tAO3cv\nB3D3TUDb8PCOQHK5dmO4T0SqUNUAulT2ZMxE0GgdTXsEBO0lia6wkp9qP6lLDdy9AuhnZvsAj5lZ\nL4LSxG6HZXrd0tLSL7ZLSkooqTyKSKQOmDED7r8/vWOTx0xcfPG8Go//979bsW1bAw4+ONqv+YlB\ndSedFOlt6qSysjLKyspqdY3Ik0SCu39oZmXAYKDczNq5e7mZtQcSn8KNwAFJp3UK931JcpIQqYsq\nT+iXjjPPXMpPfnI6P/zhPGpq7Zs37ysMGPB6jcfVVr9+QbKT7Kv8BXr8+PEZXyPq3k37JXoumdne\nwInAamAmcH542HlA4iMyExhhZo3MrCvQHVgUZYwiharyhH7pyGTMRND1Nbr2iARNz5Hfom6T6AA8\nb2ZLgYXALHd/ErgZONHM1gDHAzcBuPsqYBqwCngSGOOuQfsiVcmkPSIh3TETO3caCxZ0zfp8TVU5\n+GDYtAk++CDyW8keiLoL7Ap37+/ufd29j7v/Ity/xd1PcPce7j7I3d9POmeCu3d3957uPjvK+EQK\nVaoJ/dKRzpiJFSv2p127D2nbduueB5mm+vWhT59gESLJPxpxLVKAUk3ol450xkxE3aupMg2qy19K\nEiIFqLoJ/dJR05iJXLVHJGhQXf5SkhApMDVN6JeO6sZMfPxxI1au7MARR/x7z2+QIZUk8peShEiB\nqWlCv3RUt87EwoWd6d37LZo02V6LKDPTqxe8/nqQACW/KEmIFJh0JvRLR6p1JhLjI3KpcWPo0QNW\nrMjpbSUNNSYJM/uKmTUOt0vM7Edm1jL60ESksnQn9EtHqjET8+d345hjctcekaB2ifyUTkniEWCn\nmXUH/kgwIvqBSKMSkSplMqFfTaoaM7FpU3P+859m9Or1du1vkCENqstP6SSJCnffAZwB3OHuVxEM\nkhORHNuTAXTVqTxmYv78bhx11Drq18/9GFY1XuendJLEdjM7m2D6jCfCfQ2jC0lEUsl2kqg8ZiLK\npUprcthhsGoVbM9de7mkIZ0kMQo4CviFu68L51S6L9qwRKSyxIR+/ftn97rJDdjB+hG5bbROaNoU\nOncOEoXkjxqTRDif0tXA4vD5One/OerARGR3ezKhXzoSYybmzu1O06bb6NTp/ZpPiogar/NPOr2b\nTgOWAk+Hz/ua2cyoAxOR3WW7qikhMWbi+uu/HVspIkHtEvknneqmUoIlRN8HcPelQLcIYxKRSmoz\noV86zjxzKRs3toytPSJBJYn8k1bDtbtXnsS3IopgRKRqtZnQLx39+r3JBRe8xNFHr4vmBmnq2zeY\nDbZCf2HyRjpJYqWZjQTqm9lBZnYH8GLEcYlIktpO6FcTM/jZz2bRvPm26G6ShtatYb/94F//ijUM\nSZJOkvhvoBewDXgQ+BC4LMqgRGSXbEzoV0jULpFf0und9Im7X+fuX3f3w8Ptz3IRnIhkZ0K/QqKR\n1/kl9dJUITM7GLgS6JJ8vLt/K7qwRCQhWxP6FYr+/eE3v4k7CkmoMUkADwN3AncBO6MNR0SSJSb0\nu+GGuCPJnUR1k3v2x4RI5tJpk9jh7r9390Xu/o/EI52Lm1knM3vOzFaa2Qoz++9w/zgze9PMFoeP\nwUnnXGtma81stZkN2sP/l0hRyOaEfoWiQwdo1Ag2bIg7EoFqShJm1jrcfNzMxgCPETReA+DuW9K4\n/g7gCndfambNgH+Y2TPha7e6+62V7tkTGA70BDoBc8zsIPfKM96LFL+KCrj33mh7NeWrRGniwAPj\njkSqq276B+BAosB3VdJrThoD6tx9E7Ap3N5qZquBjuHLVRUkhwJTw1ln15vZWoKBfAtrupdIsaio\ngOnTYfx4aN4cHnoo7ohyLzGo7vTT445EUiYJd++azRuZWRegL8Ef/IHApWb2PeDvwI/DAXsdgZeS\nTtvIrqQiUtQqJ4eJE+Gkk+pmvXy/fvDnP8cdhUB6vZuGAU+7+0dm9jOgP3Cju6c9eD6sapoOjA1L\nFJOAn7u7m9n/ABOBizIJvLS09IvtkpISSqKar0AkYkoOX9a/P4wdG3cUha+srIyysrJaXcNqqu43\ns+Xu3sfMBgL/A9wC3ODu30jrBmYNCNaheMrdb6vi9c7A4+E9rgE8McusmT0NjHP3hZXOUTNFLa1b\nVxp3CNXq2rU07hAiVzk5lJZmLzkU+s/XPRh9/cor0K5dbmKqC8wMd8/oE5ZO76ZEt9dTgT+6+/8B\njTK4xz3AquQEYWbtk17/DvDPcHsmMMLMGoXrVnQHFmVwL5G8V1EB06ZB795w661ByeGll2Dw4Lpd\nekhmFlQ5abK/+KUzTmKjmf0BOBG42cwak15ywcwGAOcAK8xsCUGD90+BkWbWl2CiwPXAaAjWrjCz\nacAqYDswRkUGKRaqVspMovF68OCaj5XopJMkhgODgV+7+/tm1oHdezql5O7zgfpVvPR0NedMACak\nc32RQrFiBYwcGay+puSQnn794K9/jTsKSXfupkeBD8zsQIL1rV+JPDKRIvH443D88fCTn6haKRNa\nWyI/pNO7aQhB76P9gXeAAwmSRK9oQxMpbO5wyy1w221BovhGWl09JOHgg2HTJvjgA2jRIu5o6q50\n2hZuBI4EXg3HTpwALIg0KpECt20bjBoFU6cGU2soQWSufn3o0ydYhEjik+7KdO8C9cysnrs/Dxwe\ncVwiBeudd4Lqpa1b4YUX4IAD4o6ocGltifilkyTeDwfD/Q2438xuAz6ONiyRwrRiRVBqOO64oJtr\n06ZxR1TY1C4Rv3SSxFDgE+Bygl5JrwGnRRmUFIfNm5ty550DqKioG620jz8O3/oW/OIXcOONUC+t\njuJSHZUk4ldtw7WZ1QeecPfjCMY0TMlJVFLwli7tyCWXDOfzz+tz4IHvccopq+IOKTLJDdRPPKH2\nh2zq1Qteey1YwrVJk7ijqZuqTRLuvtPMKsysRTgBn0iNpk7tz8SJxzNhwkwaN97B+PEnM2jQKzRo\nUBF3aFm3bRuMHg3LlgUN1Gp/yK7GjeGrX91VjSe5l85guq0EI6afIaktwt1/FFlUUpC2bavP+PGn\n8Pe/H8hDD91Dt27v4g5t2mzlscf6MGxYcXVTeecd+M53grmF5s1T+0NUEmteK0nEI51a00eB6wka\nrv+R9BD5wttv78OIEaN4//29efTRP9Gt27tAMGjsxz9+lttvL2HbtqoG3xem5Abqhx9WgoiSGq/j\nVd3KdG2ANu4+pdL+XgSD6kQAWLiwM2PHfpfzz1/A6NHzvzSa+PDDN3DQQe/w0ENf4/vfL7z5Grdt\nC5bSXL8+eLz2Gtx1V9AGMXJk3NEVv379YPLkuKOou6qrbroDmFTF/tbAdYB+Peo4d5g8+UjuvHMg\nEyc+yjHHvJ7y2B//+DkuvPAcvvvdJTRpsj2HUdaschKo/Ni8GTp2DNaZTjxmzw7+eEn0DjsMVq2C\n7duhYcO4o6l7qksS3d39b5V3uvsLZvb7CGOSAvDppw259trTeO21Njz66F106vR+tcf36rWJww9/\ng3vvPYIf/nB+jqJMbetWuPtumDQpSASVk8CgQbu2998fGqTTeieRaNoUOncOEsVhh8UdTd1T3Ue/\neTWvKZ/XYW+80Yof/vAsevYs5+GH72avvXakdd7llz/PWWeNYuTIf7DPPp9FHGXVysvhjjvgD3+A\nkhKYMgUOP1xJIN8l2iWUJHKvuobrf5nZKZV3mtnJQOp6BSlqc+d258wzL+Sssxbz618/lnaCAPjK\nV/7Dcce9yt13HxVhhFVbswZ+8IOgO+WWLcFsrA8/DEceqQRRCDSoLj7V/XpcBvyfmQ1nV2+mw4Gj\ngG9HHZjkl4oK4/e/H8hf/vJ1fve7aRxxxBt7dJ2xY+cyZMgP+P73F7Lvvp9kOcove/HFYKDb/Plw\n8cXw6qvQpk3kt5Us699fa0vEJWVJwt3XAr2BuUCX8DEX6OPur+YiOMkfDz74NWbN6sljj/1pjxME\nQKdO7zNkyAruvHNgFqPbXUVF8AdlwAA491w44YSg3WH8eCWIQtW3bzBgsaL4xmPmvZpGXG8D1PlM\nWLiwM+edt5D27T+q9bUuueQFTjppDBdcsIAOHT7MQnSBzz6D++4LVn5r3hyuuioY7KbqpMLXujXs\ntx+sXQs9esQdTd2iKcgkLcuXd+SwwzZm5Vpt2mzlrLMW89vffjMr14OgS2rXrvDYY3DnnbBoEQwf\nrgRRTAYMgDlz4o6i7ok0SZhZJzN7zsxWmtkKM/tRuL+Vmc02szVmNsvMWiSdc62ZrTWz1WY2KMr4\nJD3vvbc3W7Y0+WIUdTaMHj2fp58+hPXrW9f6WjNmBNVKU6fCk08GvZa0PGjxOftseOCBuKOoe1Im\nCTN7Nvz35lpcfwdwhbv3ImjwvsTMvgpcA8xx9x7Ac8C14b0OAYYDPYGTgUlm+nWP2/LlHend+23q\n1fOsXbNly085//wF3HZbSa2uM3Vq0GvpySfh2GOzE5vkp0GDgo4H69bFHUndUl1JooOZHQ0MMbN+\nZtY/+ZHOxd19k7svDbe3AquBTgRrVCSm+5gCnB5uDwGmuvsOd18PrAWOyPh/JVm1fPn+9O6dnaqm\nZKNGLWDevG6sWdN2j86fPBmuuCKogjhcayUWvYYNYdgwePDBuCOpW6pLEjcQTOzXCbgVmJj0+HWm\nNzKzLkBfgvWx27l7OQSJBEj8legIbEg6bWO4T2K0fHlH+vR5K+vXbdbsc0aPns9vfnNcxuf+7ncw\nbhw8/zz07p310CRPjRwJ998fTAkjuZGyWc/dpwPTzex6d7+xNjcJlz+dDox1961mVvlHnPGPvLS0\n9IvtkpISSkpKahOipOAOy5Z1pLT0yUiuf+65LzN58pEsW5Z+w/gtt8Dvfw9z5waN1VJ3HH10MKXK\nihXQp0/c0eS/srIyysrKanWNGvt+uPuNZjYESHRFKXP3J9K9gZk1IEgQ97n7jHB3uZm1c/dyM2vP\nrlllNwLJy7Z0Cvd9SXKSkOi8/fY+uMP++0ez5tRee+3gkkv+xsSJ3+Lee++r9lh3+PnPg+qGv/0N\nOnWKJCTJY/Xq7SpNKEnUrPIX6PHjx2d8jRp7N5nZBGAssCp8jDWzX2Zwj3uAVe5+W9K+mcD54fZ5\nwIyk/SPMrJGZdQW6A4U3t3QRSXzDj7L7wLBhS9iwoRULFnRJeYw7XH01PPJIUIJQgqi7Ro4Mviho\nYF1upNMF9lTgRHe/x93vAQaT5rQcZjYAOAf4lpktMbPFZjYYuBk40czWAMcDNwG4+ypgGkEyehIY\n467axzhF1R6RrGHDCsaOLWPixG9VWddcUQGXXhq0Pzz/fLASnNRdvXtDixbBVCsSvXTHSbRM2m6R\n8qhK3H2+u9d3977u3s/d+7v70+6+xd1PcPce7j7I3d9POmeCu3d3957uPjvde0k0sjmIrjqnnbaC\njz7ai+efP2i3/Tt3woUXBlMyzJkD++4beShSAM45J6hykuilkyQmAEvM7M9mNoVgsr9fRBuW5IOK\nCmPFig707h1tSQKgfn3n8sufY+LE46moCOq2tm8P/hhs2ACzZgXfHkUARoyA6dPh88/jjqT41Zgk\n3P1B4EiCta4fAY5y94eiDkzit27dvrRq9SmtW0c/WyvAoEGv0KDBTp566hC2bQv6xG/dCk88oTWk\nZXddugTTvs9WXUPk0qpucve33X1m+NgUdVCSHzLplpoNZnDllc9y663HMWRIMHjq0Udhr71yFoIU\nEFU55YYm+JOUli/fnz59cpckAAYOfJ399/+A9u2DHiyNGuX09lJAhg2Dp54KSpsSHSUJSWnZso45\nTxJm8Oc//4UpUzSDq1Rvv/1g4MBggkeJTrVJwszqm9kruQpG8sfnn9fn1Vfbcuihb+f83vXrq9ez\npCcxsE6iU22ScPedwBozOzBH8UieWLOmLQcc8B5NmmyPOxSRlIYMCZao3bw57kiKVzrVTa2AlWb2\nrJnNTDyiDkziletGa5E90awZnHoqPPxw3JEUr3Rqfa+PPArJO7kaRCdSWyNHwi9/CWPGxB1JcUpn\nnMRcYD3QMNx+GVgccVwSMyUJKRRajCha6Uzw918Es7j+IdzVEfhrlEFJvLZubcSGDS3p0aM87lBE\napRYjGjq1LgjKU7ptElcAgwAPgRw97XsWiRIitDKlR3o0eMdGjbUNJtSGLQYUXTSSRLb3P2LGVLC\n9SH0oyhiqmqSQpO8GJFkVzpJYq6Z/RTY28xOBB4GHo82LIlTHCOtRWqjXj04+2x44IG4Iyk+6SSJ\na4DNwApgNME6Dz+LMiiJl7q/SiEaOTJIElqMKLvSWb60IpwifCFBNdMaLQRUvN59twkffLA3Xbps\niTsUkYz07g0tWwaLER1zTNzRFI90ejedCrwG3A78FviXmZ0cdWASj+XLO9K791vUq6fvAVJ4EqUJ\nyZ50qpsmAse5e4m7HwscB/wm2rAkLmqPkEJ29tlajCjb0kkSH7n7v5Kevw58FFE8ErNcrGktEpXO\nnbUYUbalTBJm9h0z+w7wdzN70szON7PzCHo2vZzOxc3sbjMrN7PlSfvGmdmbZrY4fAxOeu1aM1tr\nZqvNbFAt/l+yB9zV/VUKn6qcsqu6huvTkrbLgWPD7c3A3mlefzJwB3Bvpf23uvutyTvMrCcwHOgJ\ndALmmNlBaiTPnbfeakG9ehW0b/9h3KGI7LFhw+Daa4NxE82axR1N4UuZJNx9VG0v7u7zzKxzFS9Z\nFfuGAlPdfQew3szWAkcQ9KqSHEh0fbWqfjoiBSJ5MaJzzok7msKXTu+mrmZ2q5k9msWpwi81s6Vm\ndpeZtQj3dQQ2JB2zMdwnOaL2CCkWqnLKnnSmCv8rcDdBW0Q2hqlMAn7u7m5m/0PQe+qiTC9SWlr6\nxXZJSQklJSVZCK1uW7asIxdf/ELcYYjU2tChwdThmzdDmzZxRxOfsrIyysrKanWNdJLEZ+5+e63u\nksTdk9eQ+hO7pvjYCByQ9FqncF+VkpOE1N7OncbKlR1UkpCi0LTprsWI6vI6E5W/QI8fPz7ja6TT\nBfa2sEfSUWbWP/HI4B5GUhuEmbVPeu07wD/D7ZnACDNrZGZdge7AogzuI7Xw+uv7se++H9Oy5adx\nhyKSFapyyo50ShK9ge8B32JXdZOHz6tlZg8AJcC+ZvYGMA44zsz6htdaTzAfFO6+ysymAauA7cAY\n9WzKnWXLOmoQnRSVQYPg/PNh/Xro0iXmYApYOkliGNAtebrwdLn7yCp2T67m+AnAhEzvI7WnkdZS\nbBKLEd1/P1x3XdzRFK50qpv+CbSMOhCJV1CSUHuEFJcxY+C222D58pqPlaqlU5JoCbxiZi8D2xI7\n3X1IZFFJTm3bVp+1a9vSq9fbcYciklWHHhokiaFDYdGiut3TaU+lkyTGRR6FxOqVV9rRpcu7NGmy\nPe5QRLLu7LNh5Uo480yYMwcaNYo7osJSY3WTu8+t6pGL4CQ3tMiQFLuf/xz23TeoflJ3mMykM+L6\nIzP7MHx8ZmY7zUyT+xSRFSvUHiHFrV49uO8+ePlluD1ro77qhnRWpmue2DYzI5hj6cgog5LcWras\nI6NGLYg7DJFINWsGM2fCkUcG04mfdFLcERWGdHo3fcEDfwX09haJrVsb8dZbLTjooHfiDkUkcp07\nB6Owv/c9eOWVuKMpDDWWJMI1JRLqAYcDn0UWkeTUP/+5P1/9ajkNG2r1eKkbBg6Em26CIUNgwQJo\n3TruiPJbOr2bkteV2EEwSnpoJNFIzmmktdRFF1wQ9HgaPhyeeioYeCdVS6dNotbrSkj+Wr58f048\nUeVuqXt+9Sv49rfhiivgjjvijiZ/pUwSZnZDNee5u98YQTySY8uWdeTKK5+NOwyRnKtfH6ZODRqy\n//AHGD067ojyU3UliY+r2NcUuBDYF1CSKHCbNzdl69bGdOmyJe5QRGLRogU8/jgMGAA9eoCWpfmy\n6pYvnZjYNrPmwFhgFDCVYKEgKXDLl3ekd++3tFyp1GnduwdTio8YAS++CN26xR1Rfqm2C6yZtQ5X\nj1tOkFD6u/vV7q7+kkVg+fL9NdJaBDj+eLjhBjjtNPhQQ4V3kzJJmNktwMvAR0Bvdy919/dyFplE\nTmtai+xLZXzvAAAPNElEQVQyZgwce2ywWNHOnXFHkz+qa5P4McGsrz8DrrNddRJG0HC9T8SxSYTc\ngyRx000z4w5FJG/cdlswEvuKK+C88zI//9BDi28CweraJDIajS2F5c03W9Ko0Q7atfso7lBE8kbD\nhsGI7JEj4aKLMjt3yxY4/XT43/+NJra4pDOYToqQZn4Vqdq++8KsWZmf98Yb0K8f/OIX0LRp9uOK\ni0oLdVSiZ5OIZMeBB8IxxwQ9pYpJpEnCzO42s3IzW560r5WZzTazNWY2y8xaJL12rZmtNbPVZjYo\nytjqOq1pLZJ9Y8bApEnFtWZF1CWJyXx5xthrgDnu3gN4DrgWwMwOAYYDPYGTgUlm6sEfhZ07YeXK\nDurZJJJlJ5wAW7cGEwcWi0iThLvPAyp3mx0KTAm3pwCnh9tDgKnuvsPd1wNrgSOijK+uWr0a2rTZ\nSosWmsxXJJvq1YOLLw5KE8UijjaJtu5eDuDum4C24f6OwIak4zaG+yTLFi1CjdYiETn/fHjiCdi8\nOe5IsiMfejftUe1daWnpF9slJSWUaNKVtL38MmqPEIlI69Zwxhlwzz1w9dXxxlJWVkZZWVmtrhFH\nkig3s3buXm5m7YHEFB8bgQOSjusU7qtScpKQ9H3wAcyeDbfcoiQhEpUxY2DYMLjyymC22bhU/gI9\nfvz4jK+Ri+omCx8JM4Hzw+3zgBlJ+0eYWSMz6wp0BxblIL46Y+tWOOUUOPlk6Nv3zbjDESlahx8O\nbdrA00/HHUntRd0F9gHgReBgM3vDzEYBNwEnmtka4PjwOe6+CpgGrAKeBMa4F1NHsnh9+ikMHRos\nAH/77WjmV5GIJbrDFrpIq5vcfWSKl05IcfwEYEJ0EdVN27bBmWdCu3bwxz8GPTBEJFpnnQVXXQWv\nv17Y04/rz0WR27EDzj4b9toLpkyJt35UpC7Ze+9gksA//CHuSGpHSaKI7dwZfEg/+wwefFCLvYvk\n2g9/CJMnB7+DhUpJokhVVAQf0LffhkcegcaN445IpO7p3j2Y9O/hh+OOZM8pSRQhd7jsMli5EmbO\nDIq9IhKPQm/AVpIoMu5wzTXBWr1PPQXNmsUdkUjdduqpsHEjLF4cdyR7RkmiyNx4Izz5ZDAffosW\nNR8vItFq0ABGj4bf/z7uSPaMkkQ1tmwJRigXil//OpjLfs6cYOEUEckPF14YtEu8/37ckWROSaKS\nd9+Fu++GwYOhS5eg0ek//4k7qppNmhQ85swJxkOISP5o3z6Y6WDKlJqPzTdKEuyeGLp1C+ryR42C\nt96C4cPhu9+Fzz+PO8rUJk+Gm26CZ5+FTp3ijkZEqlKoCxLV2SSRKjFs3AjTpwejJZs1C9ar3Wcf\nuPTS/PzhTp0KP/tZUILo2jXuaEQklYEDoVEjeO65uCPJTJ1KEukmhmT168P998NLL8FvfxtP3KnM\nmxd0dZ01Cw4+OO5oRKQ6ZoXZHTYf1pOI3JYtwcCyWbPgxBODxDB9evrdQ5s3D8YbHHUU9OgBg/Jg\n9e1t2+C//iv4wB16aNzRiEg6zj0Xrrsu+GLasUCWVCv6ksTSpcG0vQccUH2JoSZdu8K0acEPec2a\naGLNxM03BwnrjDPijkRE0tW8eTCX2p/+FHck6SvqJHH//UHJ4Ze/hIkTaz+w7JvfDK41ZAi8V3nl\n7hx65ZVguu877tCU3yKF5uKLgySxfXvckaSnKJPE9u1BXf24cUEj0YgR2bv2RRcFXdlGjAhmWM21\nigr4wQ+C/9sBB9R8vIjkl0MPDeZ0mjGj5mPzQdElifJyOOEEWLs2WMu5d+/s3+PXvw56Ol15Zfav\nXZN77gnaI8aMyf29RSQ7CqkBu6iSxIIFQftDSQk8/ji0ahXNfRo0gIceCnpH5bJusbwcfvrT4J5a\nF0KkcJ1xBqxeHTzyXdEkiT/+MWgr+N3vYPz46Fdfa9UqSETXXQdz50Z7r4TLLoMLLoA+fXJzPxGJ\nRqNGQdV1IcznVPBJ4rPPgq6gt90GL7wQJIpcOfjgoHH8rLNg3bpo7/XUU7BoEdxwQ7T3EZHc+MEP\n4C9/ga1b446kerElCTNbb2bLzGyJmS0K97Uys9lmtsbMZplZtfOYbtgAxx4bTJq1YEHQJTTXTjwx\nKE0MGQIffRTNPT7+OKjDvPNOaNIkmnuISG4dcEDw9+uBB+KOpHpxliQqgBJ37+fuR4T7rgHmuHsP\n4Dng2lQnl5XBN74BZ54ZjF9o3jz6gFO59FI4+mg455xgydBsGzcuGNJ/4onZv7aIxGfMmPyvcooz\nSVgV9x8KJOZJnAKcnurkESOCGRV/8pP4xwqYBWMWPvwwmEcpmxYvhvvug1tvze51RSR+xx+f/0ub\nxpkkHHjGzF42s4vCfe3cvRzA3TcBbVOdvGBBfn2zbtQoGM09bVpQz5gNO3YE9ZY33wxt2mTnmiKS\nP+rVC8ZM5LM4524a4O5vm1kbYLaZrSFIHMlSzrv65z+XfrFdUlJCSUlJFDFmZL/9gjmejjsOXn01\n6K661157fr077ghmoD3vvOzFKCJ1R1lZGWVlZbW6hnkezH9tZuOArcBFBO0U5WbWHnje3XtWcbzn\nQ9ypbNwIY8fC8uVBfePxx2d+jfXrgzEfL70EBx2U9RBZt640+xfNoq5dS+MOoaDp5ytVMTPcPaMK\n+liqm8ysiZk1C7ebAoOAFcBM4PzwsPOAAhm4vruOHYOqp4kTg3EN3/8+bN6c/vnucMklcPnl0SQI\nEZF0xdUm0Q6YZ2ZLgAXA4+4+G7gZODGsejoeuCmm+LLitNNg5Upo2zaYr+Wee9JbuOjhh+Hf/4ar\nroo+RhGR6sTSJuHu64C+VezfApyQ+4ii06xZMNfTOefA6NFBj6w774SeX6pEC7z3XjCy+pFHgsZw\nEZE4FfyI60LRr1/QvjBsWDDl+A03BKPFK7v6ajj99GCBIxGRuClJ5FD9+sHAu6VLYdWqYA6mZ5/d\n9foLL8CTT8KECfHFKCKSTEkiBlU1bL/5ZjAm4vbboUW1k5GIiOSOkkSMkhu2DzoomDBQy5GKSD6J\nczCdsKth+8ILoV27+KcYERFJpiSRJ1L1dhIRiZOqm0REJCUlCRERSUlJQkREUlKSEBGRlJQkREQk\nJSUJERFJSUlCRERSUpIQEZGUlCRERCQlJQkREUlJSUJERFJSkhARkZTyMkmY2WAze8XMXjWzq+OO\nR0Skrsq7JGFm9YDfAicBvYCzzeyr8UaVv8rKyuIOIW/ovdhlwYL1cYeQN/S5qJ28SxLAEcBad/+3\nu28HpgJDY44pb+kXYBe9F7soSeyiz0Xt5GOS6AhsSHr+ZrhPRERyLB+ThIiI5Alz97hj2I2ZHQmU\nuvvg8Pk1gLv7zUnH5FfQIiIFwt0zWiQ5H5NEfWANcDzwNrAIONvdV8camIhIHZR3a1y7+04zuxSY\nTVAddrcShIhIPPKuJCEiIvmj4BquNdBuFzNbb2bLzGyJmS2KO55cMrO7zazczJYn7WtlZrPNbI2Z\nzTKzFnHGmCsp3otxZvammS0OH4PjjDFXzKyTmT1nZivNbIWZ/SjcX+c+G1W8F/8d7s/os1FQJYlw\noN2rBO0VbwEvAyPc/ZVYA4uJmb0OfM3d34s7llwzs4HAVuBed+8T7rsZeNfdfxV+gWjl7tfEGWcu\npHgvxgEfufutsQaXY2bWHmjv7kvNrBnwD4JxVqOoY5+Nat6Ls8jgs1FoJQkNtNudUXg/w6xw93lA\n5eQ4FJgSbk8BTs9pUDFJ8V5A8PmoU9x9k7svDbe3AquBTtTBz0aK9yIx5iztz0ah/YHRQLvdOfCM\nmb1sZv8VdzB5oK27l0PwCwK0jTmeuF1qZkvN7K66UL1SmZl1AfoCC4B2dfmzkfReLAx3pf3ZKLQk\nIbsb4O79gVOAS8JqB9mlcOpSs28S0M3d+wKbgLpW7dQMmA6MDb9FV/4s1JnPRhXvRUafjUJLEhuB\nA5Oedwr31Unu/nb472bgMYLquLqs3MzawRf1se/EHE9s3H2z72pw/BPw9TjjySUza0DwR/E+d58R\n7q6Tn42q3otMPxuFliReBrqbWWczawSMAGbGHFMszKxJ+A0BM2sKDAL+GW9UOWfsXrc6Ezg/3D4P\nmFH5hCK223sR/iFM+A5167NxD7DK3W9L2ldXPxtfei8y/WwUVO8mCLrAArexa6DdTTGHFAsz60pQ\nenCCQZH316X3wsweAEqAfYFyYBzwV+Bh4ADg38Bwd38/rhhzJcV7cRxBHXQFsB4YnaiTL2ZmNgD4\nG7CC4HfDgZ8SzNwwjTr02ajmvRhJBp+NgksSIiKSO4VW3SQiIjmkJCEiIikpSYiISEpKEiIikpKS\nhIiIpKQkISIiKSlJSFEzs48qPT/PzO7IcQzfNbNVZvZspf2dzewTM/tH+PoCMzsvl7GJ1CTvVqYT\nybKqBgLVenCQmdV3951pHn4hcJG7v1jFa/9y96+F1+wCPGZmuPuUKo4VyTmVJKTOCr/JPxvOhvmM\nmXUK9082s+8kHfdR+O+xZvY3M5sBrKziemeb2fLwMSHcdz0wELg7XO8iJXdfD1wBjA3P/bqZvRiW\nNOaZ2UHh/rlm1ifpvi+YWe/avRsiVVNJQopdEzNbHG4b0Ipd833dAUx297+Y2ajw+RlVXCO55NEP\n6OXubyQfYGYdgJvC198nmMJ9iLvfaGbfAq5w9yVpxLsY6BFurwYGunuFmR0PTAC+C9xFsIjO5WHi\naOzuK9K4tkjGVJKQYveJu/cPH/0I5jVKOAp4MNy+DxiQxvUWVU4Qoa8Dz7v7FnevAO4Hvpn0erqL\nvCQf1xKYbmYrgN8Ah4T7pwOnmll94ALgz2leWyRjShJSl6Vqm9hB+LthZgY0Snrt42qul42V4PoT\nlCAAbgSec/fewGnAXgDu/inwDMHqasMIEpJIJJQkpNhV94f7ReDscPtc4IVwez1weLg9FGiYxn0W\nAd80s9bhN/yzgbJM4gsbrm8Bbg93tWDXeimjKp13d3jcInf/II37iOwRJQkpdtX1ZPoRMMrMlgLn\nEDYYEyzEcqyZLQGOpPrSQ3CTYEnMawgSwxLgZXd/Io0YuiW6wBKs2f6/7n5v+NqvgJvM7B9U+l11\n98XAh8DkmmITqQ1NFS5SgMxsf4KqqK/GHYsUN5UkRAqMmX0PeIlgARmRSKkkISIiKakkISIiKSlJ\niIhISkoSIiKSkpKEiIikpCQhIiIpKUmIiEhK/w/wqEAgnNNrpQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "\n", "import matplotlib.pyplot as plt\n", "\n", "plt.plot(OrdHrCnts)\n", "plt.axvspan(6, 9, color='y', alpha=0.5, lw=0)\n", "plt.axvspan(15, 18, color='y', alpha=0.5, lw=0)\n", "plt.xlabel('Hour of Day')\n", "plt.ylabel('Number of Crashes')\n", "\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's now do a similar thing, but accross different days rather than hours. Weekends are highlighted." ] }, { "cell_type": "code", "execution_count": 119, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import re\n", "days=pandas.Series([re.search('-(.+)-(.+)',i).group(2) for i in data.DATE_])\n", "dayCounts = days.value_counts()\n", "ordDayCounts = [dayCounts.loc[\"{0:0=2d}\".format(i)] for i in range(1,dayCounts.shape[0]+1)]" ] }, { "cell_type": "code", "execution_count": 96, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEPCAYAAAC3NDh4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmYHGW1h98zSWayk20mCQmZhEQkCWHJgkBYWkQUlYBe\nhYhXBe9FBAVcUIzITQJiBEUuqAgoxCggFwIKIgoCjgoYCFlIIHuYrGTfZibLZDJz7h9fdaYz00tV\nd9f0dPd5n6efdH9dXXWm0l2/Oud853yiqhiGYRhGPEpybYBhGIbRfjGRMAzDMBJiImEYhmEkxETC\nMAzDSIiJhGEYhpEQEwnDMAwjIaGKhIgMFpGXReQdEVksIte1eP9bItIkIn1ixqaIyEoRWSoi54dp\nn2EYhpGcjiHv/xDwTVVdKCLdgXki8oKqLhORwcCHgbXRjUVkJHAJMBIYDLwoIu9TK+YwDMPICaF6\nEqq6WVUXes/rgKXAIO/tu4Bvt/jIRcBjqnpIVdcAK4FTw7TRMAzDSEyb5SREZChwMvC6iEwC1qvq\n4habDQLWx7zeSLOoGIZhGG1M2OEmALxQ02zgeqAR+B4u1GQYhmG0Y0IXCRHpiBOI36nq0yJyAjAU\neEtEBJd7mC8ip+I8hyExHx/sjbXcp+UoDMMw0kBVJcj2EnZOWER+C2xX1W8meL8aGKuqu0RkFPAI\n8AFcmOlvQKvEtYi0ymVXV0/zZc+wYf62C5Np06YxbVp6dvj9O/2Q7rnIxP5skcl5+N//reLrX48c\nMdYevhd+aXn+s/mdgHDPRdjfnbDPRXv47meCiAQWiVA9CRGZCHwOWCwiCwAFvqeqf43ZTAEBUNUl\nIvI4sARoAK6xmU2GYRi5I1SRUNVXgQ4ptjm2xesZwIww7TIMwzD8YRXXOSASieTahIzId/tPO21o\nrk3IiHw+//lsO+S//elgIpED8v2Llu/2m0jkjny2HfLf/nQwkTAMwzASYiJhGIZhJMREwjAMw0iI\niYRhGIaREBMJwzAMIyEmEoZhGEZCTCQMwzCMhJhIGIZhGAkxkTAMwzASYiJhGIZhJMREwjAMw0iI\niYRhGIaREBMJwzAMIyEmEoZhGEZCTCQMwzCMhJhIGIaRkL174Te/ybUVRi4xkTAMIyFLl8L//E+u\nrTBySagiISKDReRlEXlHRBaLyLXe+B0islREForIkyLSM+YzU0Rkpff++WHaZxhGcnbvhu3bc22F\nkUvC9iQOAd9U1dHA6cDXROR44AVgtKqeDKwEpgCIyCjgEmAkcAFwr4hIyDYahpGAPXtg/37Yty/X\nlhi5IlSRUNXNqrrQe14HLAUGqeqLqtrkbTYHGOw9nwQ8pqqHVHUNTkBODdNGwzASs3u3+3fbttza\nYeSONstJiMhQ4GTg9RZvfQl4zns+CFgf895Gb8wwjBywZ4/710JOxUubiISIdAdmA9d7HkV0/Cag\nQVV/3xZ2GIYRjKgnYSJRvHQM+wAi0hEnEL9T1adjxi8HPgacG7P5RuCYmNeDvbFWTJs27fDzSCRC\nZWXWTDYMw8M8ifymqqqKqqqqjPYRukgADwFLVPXu6ICIfBT4NnC2qtbHbPsM8IiI3IULM40A3oi3\n01iRAKiursqq0YZhOE+iRw8TiXwlEokQiUQOv54+fXrgfYQqEiIyEfgcsFhEFgAK3ATcA5QCf/Mm\nL81R1WtUdYmIPA4sARqAa1RVw7TRMIzE7NkDI0ZY4rqYCVUkVPVVoEOct96X5DMzgBmhGWUYhm92\n74bhw82TKGas4towjITs3u08CROJ4sVEwjCMhETDTSYSxYuJhGEYCbFwk2EiYRhGXJqaoKbGRKLY\nMZEwDCMudXXQpQsMGAA7djjRMIoPEwnDMOKyZw/06gWdOkG3bs2FdUZxYSJhGEZcdu+Go45yz/v1\ns5BTsWIiYRhGXKKeBJhIFDMmEoZhxMU8CQNMJAzDSIB5EgaYSBiGkYBYT6K83Po3FSsmEoZhxMU8\nCQNMJAzDSIDlJAwwkTAMIwG7d5snYZhIGIaRAAs3GWAiYRhGAlqGmyxxXZyYSBiGEZdYT6K83DyJ\nYsVEwjCMuMR6Er16QW0tNDTk1iaj7TGRMAwjLrGeREkJ9OkDO3fm1iaj7QlVJERksIi8LCLviMhi\nEbnOG+8tIi+IyHIReV5Ejor5zBQRWSkiS0Xk/DDtMwwjMbGeBFjyulgJ25M4BHxTVUcDpwNfFZHj\nge8CL6rq+4GXgSkAIjIKuAQYCVwA3CsiErKNhmG0oL4eDh2Crl2bx0wkipNQRUJVN6vqQu95HbAU\nGAxcBMzyNpsFXOw9nwQ8pqqHVHUNsBI4NUwbDcNozZ49zouIvUWzGU7FSZvlJERkKHAyMAfor6pb\nwAkJUOFtNghYH/Oxjd6YYRhtSGw+IorNcCpOOrbFQUSkOzAbuF5V60REW2zS8nVKpk2bdvh5JBKh\nsjIjEw3DiKFlPgIs3JSPVFVVUVVVldE+QhcJEemIE4jfqerT3vAWEemvqltEZACw1RvfCBwT8/HB\n3lgrYkUCoLq6KotWG0ZxE8+T6NcP1q7NjT1GekQiESKRyOHX06dPD7yPtgg3PQQsUdW7Y8aeAS73\nnn8ReDpmfLKIlIrIMGAE8EYb2GgYRgzmSRhRQvUkRGQi8DlgsYgswIWVvgfcDjwuIl8C1uJmNKGq\nS0TkcWAJ0ABco6qBQ1GGYWRGIk/CRKL4CFUkVPVVoEOCt89L8JkZwIzQjDIMIyWxHWCj2MJDxYlV\nXBuG0QoLNxlRTCQMw2iFhZuMKCYShmG0Ip4n0a0bNDbCvn25scnIDSYShmG0Ip4nIeK8iR07cmOT\nkRtMJAzDaEU8TwIs5FSMpBQJERkuImXe84iIXCcivVJ9zjCM/CWeJwE2w6kY8eNJPAk0isgI4AFc\nRfSjoVplGEZOMU/CiOJHJJpU9RDwSeBnqvptYGC4ZhmGkUsSeRImEsWHH5FoEJHP4tpnPOuNdQrP\nJMMwcklTE9TUQM+erd8zkSg+/IjEFbgFg25T1Wqvp9LvwjXLMIxcUVfnFhvqGKcfg4lE8ZGyLYfX\nT+lGYIj3uhrXe8kwjAIkUT4CbOGhYsTP7KYLgYXAX73XJ4vIM2EbZhhGbkiUjwBbeKgY8RNumoZb\nQnQ3gLcc6bEh2mQYRg5J5UmYSBQXvhLXqrqnxVhTGMYYhpF74nWAjWIiUXz4EYl3ROQyoIOIvE9E\nfga8FrJdhmHkiGThpr59XVsOW+WlePAjEtcCo4F64PdADfD1MI0yjFh27cq1BcVFsnBTWRl07uym\nyBrFQUqRUNV9qnqTqk5Q1fHe8wNtYZxhAHzgA7B4ca6tKB6SeRJgM5yKDT+zm44TkQdE5AUReTn6\naAvjDGPXLli5Etaty7UlxUMyTwJshlOx4Wf50ieA+4BfA43hmmMYR7Jggft306bc2lFM7NkDI0Yk\nft+S18WFn5zEIVX9paq+oarzog8/OxeRB0Vki4gsihk7SUT+LSILROQNERkf894UEVkpIktF5Pw0\n/h6jwJjnfdNMJNqOVJ6EiURxkVAkRKSPiPQB/iQi14jIwOiYN+6HmcBHWozdAUxV1VOAqcCPveON\nAi4BRgIXAPeKiAT8e4wCY/58GD8eNm/OtSXFg5+chIlE8ZAs3DQPUCB6of52zHuKj4I6VX1FRCpb\nDDcB0fuUXsBG7/kk4DGv4+waEVmJK+J7PdVxjMJl/nz47Gdh0aLU2xrZwTwJI5aEIqGqw0I65jeA\n50XkTpwAneGNDwL+HbPdRm/MKFJqamDDBvjgB+H553NtTf7Q2CjU1namV6/9aX3ejyexalWaxhl5\nR8rEtYh8BvirqtaKyPeBscCtqrogzWNeDVyvqn8UkU8DDwEfDrqTadOmHX4eiUSobOmvGHnPW2/B\nmDFwzDEWbgrCiy++n8ceG8fMmY+k9Xmb3VQ4VFVVUVVVldE+/MxuullVnxCRM4HzcDmE+4APpHnM\nL6rq9QCqOltEfu2Nb8StehdlMM2hqFbEigRAdXVVmuYY7ZX582HsWBgwwImEKliWKjWrVpWzaVOS\nq3wKLCdROEQiESKRyOHX06dPD7wPP7ObotNePw48oKp/BkoDHENozmsAbBSRcwBE5EPASm/8GWCy\niJR6a1aMAN4IcByjwIiKRNeuUFrq7nCN1FRX92Xr1u5pfba+vgOHDkGXLom3MZEoLvx4EhtF5H5c\nSOh2ESnDn7ggIo8CEaCviKzDzWa6ErhHRDoAB4Avw+F1Kx4HlgANwDWq1iGmmJk3D66/3j0fONBN\ng+3dO7c25QPV1X3Zvbsr9fUdKCsLVtrkchnJPTYTieLCj0hcAnwU+Imq7haRgRw50ykhqnpZgrfG\nxxtU1RnADD/7Ngqbfftg9WoYPdq9HjjQhZxGjcqtXe0dVXj33X5061bPtm09GDw4mPtVU9M5aagJ\nnFDv2QOHDsVfvc4oLPz2bnoK2CMiQ3DrWy8L3TKjqFm0CEaOdA3loNmTMJKza1dXVGHEiG1phZxq\nazsnTVoDdOjghGLnzjSNNPIKP72bJnk1C9XAP7x//xK2YUZxM38+jBvX/HrAABMJP1RX92XYsB30\n71/Lli09An++pqYspScBFnIqJvzkFm4FTgNWeLUT5wFzQrXKKHqiSeso5kn4IyoS5eV1bNsWXCT8\neBJgIlFM+F2ZbgdQIiIlqvp3EuQUDCNbxBMJq5VITXV1X449Nn1PIpq4ToWJRPHgRyR2i0h34J/A\nIyJyN7A3XLOMYqa+HpYtgxNPbB6zcJM/mj2J2rRyEjU15kkYR+JHJC4C9uHaafwVWA1cGKZRRnHz\n9tswfPiRc/Ut3OSPd99tzkmkF27yn5OwhYeKg6QT2LxahmdV9YO4xnyz2sQqo6hpGWoCCzf5oalJ\nWLu2D5WVbtpReonrzknXkohSXg7r1wfevZGHJPUkVLURaBKR9Gv8DSMgLWc2gZtyuX+/exjx2bSp\nJ71776dbt4OUl9eybVt6U2AtJ2HE4qcUpg5YLCJ/IyYXoarXhWaVUdTMnw+f+9yRYyLQv7/zJoaF\n1Z84z4mGmgD69t1HbW1nDh7sQGmp/6pry0kYLfEjEk95D8MInYYGl5M46aTW70XzEiYS8XFJa3fl\nLilR+vbdy7Zt3Rk0aI/vfQTJSZhIFAcJRUJEyoFyVZ3VYnw0sDVsw4ziZNky1xq8R5xwuuUlkrNm\nTbMnAVBRUcvWrT0CiYR5EkZLkuUkfgb0izPeB7g7HHOMYmfevNZJ6yg2DTY5777bL45IBMtLBMlJ\n2Oym4iCZSIxQ1X+2HFTVfwEnxtneMDIm3symKDYNNjnRGokoFRXBC+r8NPgD5+k1NNhEgmIgmUgk\n+3Z1yrYhhgGpRcLCTfGpr+/Ali1Hdn2tqAjWmqOpSdi7t5SePVNvK+K8iR07Um9r5DfJRGKViHys\n5aCIXAC8G55JRrHS2OiWLDVPIjjr1vVh0KDddOrUdHisf/9g4aa6ulK6dm2gQwd/21teojhINrvp\n68CfReQSYJ43Nh44HfhE2IYZxcfKlVBRkXjpTMtJJCbasykW15rDvydRW9uZHj0OAGW+tjeRKA4S\nehKquhIYg2sPPtR7/AM4UVVXtIVxRnGRLNQE5kkko2U+Aly4KYhI1NRERcIflrwuDpLWSahqPTCz\njWwxipxUItG/v7tzbQy2ImdRUF3dl5NO2njEmAs3BfUk6n1vb55EceBrrWrDaAuSTX8Ft1Rmnz52\n9xqP2GrrKH367GXPns40NPj7mdfWltGzp39PorzcRKIYCFUkRORBEdkiIotajF8rIktFZLGI/Chm\nfIqIrPTeOz9M24z2RVMTLFiQXCTA8hKJiBdu6tChueraD+mEm0wkCp+EIiEiL3n/3p7B/mcCH2mx\n3wiu1fgYVR0D/MQbHwlcAowELgDuFRHJ4NhGHlFdDT17urvTZFheojU1NZ05cKATFRW1rd4LskJd\nbW3nQJ6EiURxkMyTGCgiZwCTROQUERkb+/Czc1V9BdjVYvhq4EeqesjbJvo1uwh4TFUPqeoaYCVw\napA/xshfUuUjolitRGuqq/swbNgO4t1SBZkG6zwJy0kYR5Iscf0/wM3AYOCnLd5T4Nw0j3kccLaI\n/BDYD9ygqvOAQcC/Y7bb6I0ZRYBfkbBwU2uqq/sydGj8qrYgVde1tWX06bPP93FtdlNxkFAkVHU2\nMFtEblbVW7N8zN6qepqITACeAI4NupNp06Ydfh6JRKiszJp9Rg6YPx+uvTb1dgMHwgqbgH0E8fIR\nUYKEm2pqOh9esMgP5km0f6qqqqiqqspoHylbhavqrSIyCTg7elxVfTaDY67Haz2uqnNFpFFE+uI8\nhyEx2w32xuISKxIA1dVVGZhk5BLVYOGmf7bqKFbcVFf349xzl8d9r3//Wt56y59Dnm5OQpW4oS4j\n90QiESKRyOHX06dPD7yPlLObRGQGcD2wxHtc74WK/CLeI8of8UJVInIcUKqqO4BngEtFpFREhgEj\ngDcCHMfIU9avhw4dnACkwhLXrUnmSQTpBBs0J9G5M5SVQW3rfLlRQPhZdOjjwMmq2gQgIrOABcD3\nUn1QRB4FIkBfEVkHTAUeAmaKyGKgHvgCgKouEZHHcULUAFyjqhr4LzLyjqgX4edu1HISR6LqRyT8\n5ySCeBLQ7E34aQpo5Cd+RAKgFxANVvpe71pVL0vw1ucTbD8DmOF3/0ZhEG9N60REPQkLcTg2bYKu\nXQ/Ss2d8DyBIa46ammDhJmgWiWMDZxWNfMFPMd0MYIGI/MbzIuYBt4VrllFM+M1HAHTrBp06ubte\nwyXxE3kRAH377mXPni4cOpT6px60mA5shlMxkPKbo6q/B07DJZufBE5X1f8L2zCjeAgiEuBCTkHW\nSShkUolEx45N9Oq1j+3bu6XcV9DENdgMp2LAV1sOVd2kqs94DytlMrLGpk1QXw9DhqTeNsrAgQRe\nlrNQWb48uUiAv0Z/9fUdURXKyg4FOr71byp8rMGfkVOCJK2jOJEwTwJSexLgaiVSna+aGpe0Dprn\nMU+i8DGRMHJK0FATOJHw27Su0PEjEn6mwdbWdqZ792ChJjCRKAaSioSIdBCRZW1lTHvl4EEXEjGy\nT5CZTVFcTsJEoqEB1q6FIUOSV0n7ac2RTj4CTCSKgaQioaqNwHIRCRAxLjzuuANuvjnXVhQm6XoS\nFm6CNWtg0CAoK0u+ClNFRerWHEEL6aLY7KbCx0+dRG/gHRF5A9gbHVTVSaFZ1c54+23rPBoG27fD\n7t3B59hbuMmxYgUcd1zq7fr3r+Xll5NvmE4hHYSTuH7+edi/Hy6+OLv7NdLDj0gU/T30ihWwapVb\nGKfEsjhZY8ECOOWU4Oe0vYSb1q93q+X5aScSBn5Fory8NuX5SqdGAsIJN91yi6vgNpFoH/ipk/gH\nsAbo5D2fC8wP2a52g6r7MZaVwerVubamcFCFe+6Bj30s+GfbS7hp6lQXKps3LzfHX77cvyeR6nyl\nm5Po3dt5g9lad3z1aliyBObMcTdlRu7x0+DvSmA2cL83NAjXpK8o2LTJVfmefbaLnxvZ4Q9/cBeE\n668P/tm+fWH//k7U1/vtKhMOy5bBf/4nXHABPPdc2x/fryfRr99edu7sSmNj4vmt6eYkOnaEo46C\nXS2XFkuTRx5x57RfPycWRu7x8yv7Km6FuNcBVHWliFSEalU7Inq3NnasE4lLL821ReHQ0ODujPfv\nT73tZz4DZ5yR/rFqapw4PPKI89CCIuLaTWzf3o1Bg/akb0iGLF/uxO4//gM++UmYNg2uuqrtjh8V\niVR33K7qej87dnSjoqIu7ja1tWUMHep/LYlYosnrfv3S+vhhVOHhh92jthZeew1OOCGzfRqZ4yca\nXK+qB6MvRKQjbmW6omDFCnj/+5tFolBZsQJmzXKVz8keffrApz7l4vHpcvPNcP75zjtLFz8FYmGy\nfbsLsVRUwGmnwSuvwJ13wpQp6YVJ9u+HvXtTbxelrg527oRjjvG3fappsOk094uSrbzEG97CABMm\nuJuQV1/NfJ9G5vjxJP4hIt8DuojIh4FrgD+Fa1b7IXq3FhWJQu0+umqVSyJ/4xupty0rg0sugX/8\nA0pLgx3nzTfh//4P3nknPTujBFknIQyWL3c3D9HvwvDh7s73ootcuGTmTH9e0oYN8ItfwAMPwPjx\nbmaPH1atghEj/Cf9m0U1fp/1dBPXbt/ZEYmHH3bnTgQmToQf/zjzfRqZ4+cr9l1gG7AYuAp4Dvh+\nmEa1J6Ii0b+/W2Rl7dpcWxQO0YuOH264wV0Ybrwx2DEOHXLhmDvucHmFTMi1JxEViVj69YMXX3TF\nlx/5SPI4/euvw2c/CyeeCPv2OU9k1Sr4+9/9Hd9vPiJKquR1uolryI4n0dDgbh7+8z/d65Ej3T63\nbMlsv0bm+Jnd1ATMAm4FpgOzimkxoNgfYyGHnIKIREmJC0398Y8we7b/Y/ziFy7J+fm4q4kEo6Ki\nlu3bc+9JtKRLF3j8cVdFPnGiK3iLEr0Qnn46TJ4Mp54K1dVw993uonjLLS5c5efX5XdmU5RUnle6\niWvIjkg8/zy8733NNTMlJe48/fvfme3XyBw/s5s+DqwG7gF+DqwSkQvCNqw9EG17MHy4e13IIrF6\ntX+RADf18Ykn4OqrYeXK1NuvXw+33gq//GV2wnXOk2h/IgHuAnfnnfCVrziheOkluP12dwG89174\nznecKH/jG040o3z2s86reOaZ1McP6kmkWqEu3WI6yI5IPPxw65sHy0u0D/yEm+4EPqiqEVU9B/gg\ncFe4ZrUPqqvh6KObY8vjxrWNSKhCJAKLF4d/rChBPIko48fD9Onw6U+nnhV13XVw7bWJL6xBCbLi\nWhgkE4ko113nvKfLLnPTOZ9+2uVxPvlJt6Z3S0pK4Lbb4KabUtcdBBeJ5Ocrk5xEpq059uyBv/zF\nzZqLZeJEE4n2gB+RqFXVVTGv3wWKYunz6MymKNHCqbCDba+/7i4mb70V7nGiHDwIGzdCZWXwz159\nNYwa5QQgEU8/7S6S3/1u+ja2xE8VcVgcOuRuIPyI6sUXu7j6rFn+elR94hOu2vjRRxNvEy3wDCK4\nLicR/3w1NQn79pXSvXtuwk1PPQUf/GDrPNWpp7rfwIH0tMvIEglFQkQ+JSKfAt4UkedE5HIR+SJu\nZtNcPzsXkQdFZIuILIrz3rdEpElE+sSMTRGRlSKyVETOT+PvySot79YGDXI/0E3xJ4hkjfvucz+Y\ntqrwXrMGBg92y4IGRcTNzHnlFXchbEldnROQ++5LryYiEeXldTkTiepqV/XdpUv29y0CM2a4mpWD\nB+Nvs3278zqCJP/LyxOHm+rqyujatYEOHdK7+8l0dlO8UBO4ItaRI3NX0W44knkSF3qPzsAW4Bwg\ngpvp5PfnMRP4SMtBERkMfBhYGzM2ErgEGAlcANwrktvJpi1FQiT8Ngy7drmE8He+03YikU6oKZYe\nPVwC+4YbWofIpk6Fc891d4rZxFURd0taRRwWfkJNmXDOOe5796tfxX8/aKgJnKgmOl+1tWVph5og\nM09iwwZYuBA+/vH4759xhptabOSOhCKhqlcke/jZuaq+AsSbCHgX8O0WYxcBj6nqIVVdA6zEVXrn\njHg/xrCT17/7nWvzcNpp+SMS4Cpj77zTxZVrvWDkggXuLvEnP8ncxpaUljbSo8cBdu7smv2dpyBs\nkQD44Q9dfiJegV3QmU0AnTo10bPnAXbubL3WdSb5CMhMJB591FWsd+4c/33LS+QeP7ObhonIT0Xk\nKRF5JvpI94AiMglYr6ot07KDgNg63o3eWM5oa5FQhfvvd7UEw4e3nUgEndmUiC98Ac46C6680iVe\nr7oKfvSjzNs1JKK8vC4n02DbQiTGjoUzz3RNEFuSjicB0ZBT6/OVSY0EuBzKgQPBF+ZSdTdF0dqI\neEQ9ieKZdN/+8FNx/UfgQVwuIqO+jCLSBfgeLtSUEdOmTTv8PBKJpJV0TUaitgfjxvmrSk6HV15x\nF9dzznE/ipoaZ0f3kK+Dq1bBeedlZ1/33ON+2Oef72L2l1+enf3GI1pQN3Jk21ZcLV/eNj28br3V\nCcVXvuKmHEdZscJNlw1KtKBu9OgjF0fJpEYCor20nDcxKMBt3aJF7jt+5pmJtznmGOdlrFrl6iiM\nYFRVVVFVVZXRPvyIxAFVjXM/kxbDgaHAW16+YTAwX0ROxXkOsSvgDfbG4hIrEgDV1VVZMtGxcmX8\ntgdDh7oL99atrm9PNrn/fvjyl92PTgSGDYN333VVuWGyalVzLUimdOni6icuuMDN9w8zq5Sr1hxt\n4UmAO8ZFF7kK9RkzmseDzmyKkqhKPZMaiSjRkFMQkYi24UjVWmTiROdNmEgEJxKJEIlEDr+ePn16\n4H34mQJ7t4hMFZHTRWRs9BHgGOI9UNW3VXWAqh6rqsOADcApqroVeAa4VERKRWQYMAJ4I+DfkzUS\nufTR5PWCBdk93o4d8Oyz8MUvNo8NH+4u4GFy6JArGAy6OlwyRoxwIjtyZPb2GY9czHDas8fdJAS5\nGGbC1Knu5iE6o66xMf3wYKJpsJk094sSdIZTY6PLRyQLNUWxorrc4kckxgBXAj/CFdbdCfhKRYrI\no8BrwHEisk5EWia8lWYBWQI8DizB9Ye6JpftP5LFfcPIS8yaBRdeeOS0xrbIS6xf7zyiRInD9owT\nibYtqIsmjdtq3t0xx7iQ3a23utfr17sLctc08vXl5fE7wdbWZpa4huDJ67//3U0j9nMjEfUkjNzg\nJ9z0GeDY2HbhflHVy1K8f2yL1zOAGQk2b1NWrIAPfSj+e2PHwpNPZu9Y0YT1gw8eOT5iRPhV19lK\nWueCiopa5s3z2Ss7S7RVqCmWKVPg+OPhW99ynmU6SWtwnsSrr7aOK9bUdKZfv/jrTPhl4kTXAn7E\nCJe3S0WqhHUsJ54I69a56eGxuRmjbfDjSbwN9ArbkPZGW3oSVVWukG3ixCPH28KTyMb011yRi9Yc\nuRCJ8nLX4mPq1PRnNoE7X/E8iUwT1+AWkbr1VpeLuuuu5LOR9u51VfiTJ/vbd8eObo2JOXMyMjFj\n1qxJ3lnJ7dveAAAb+klEQVSgUPEjEr2AZSLyfDamwOYD0bYHiX6M73ufc62ztWRjdNpryxBGW4lE\ntpLWbU15edt3gs2FSAB885vwt7+5Fhbpi0T8VibZSFyDm/H1+uuu0+0nPuEmd8TjmWdch9cBA/zv\nuz3kJV56yXUXKLY2IX5EYirwSeCHNOck7gzTqFyzbZubcZFofn9JCZx8cnaS11u3wl//Gr8tQWWl\n66nU0JD5cRKRz55EtBNsW2auciUSPXq43ldVVekfP1pX0tR05N1IpsV0sQwbBv/6lwsRnXKKu7C2\nJEioKUp7yEvMnetapcz11ZSocPCznsQ/4j3awrhc4celz1bIaeZM1xW0V5yAXmmp60Ib5kJH+SwS\n3bsfRMT1HmoLmpoyywlkytVXuzvqdKdEl5Y20r17fasq9UyL6VrSqZObsjtrliuw/N73mm90tmxx\nF/uLLw62z9NOcxfnMG+YUjF3rvOAXnkldzbkAj8V17UiUuM9DohIo4jUtIVxuaKtRKKpybmvX/lK\n4m3CDDk1Nbk6jHwNN0HbdoNdt86t8R12cWMiOnd2IZdMpt/GCzllIycRj/POc972woVuPfM1a+Cx\nx2DSJNe8Lwi9erkapUWtWoW2DQcOwNKlLifxr3/lxoZc4ceT6KGqPVW1J66x338A94ZuWQ7xKxKZ\nNvp76SUXRjg1SYeqMEVi0ya36E2uLnrZoC1rJXIVasomFRWtp8FmKycR/3iu/uczn3Hf85/+NHio\nKUo6eYmGhhJ27cq8Xe/ChW6G2XnnOU8o1XofhYTPZdQd6vgjcTq7FhJ+ROL4410Hy5oMfKpECetY\nwiyoy+dQU5S2nOFUKCLR8nxlo04iGSUlLvH+3HOuG3CiqeWpSCcvcdttH+H66z+d3gFjmDvXzbAq\nL3cJ97ffzniXeUPKOglvTYkoJcB4oKDz+35EomNHGDPGLYpy1lnBj7Fpk/MkHnoo+XbDh4c3qyOf\nZzZFacvWHIUhEkeKan29uwSUlR0K/djjx8NvfpP+5884w63a55e33hrEn/88mgMHOtHQUEKnTum3\nnps714XMwP3e//UvOOmktHeXV/jxJC6MeXwEtyrdRWEalUuibQ/89InJZDnThx5yLnjPnsm3CzPc\nVAiehIWbgtGyNUc0H5HblVv8MXy4m120fn3qbQ8dKuGmmz7BTTc9z6BBu1myJMB82zjMndscFj7z\nzOJKXvvJScSuI3Glqt7m9VoqSNat89/2IN3kdWOjW1DmqqtSbzt8uEsuhzHNs3BEwsJNfmm5Ql2m\nCw61JSL+8xK/+c0H6N17HxddtJgJE9Yxd276baL37HHCNGqUex31JIqlfXnCcJOI/E+Sz6mq3hqC\nPTknSEXr2LHwv/8b/BgvvOBqMPy0L+jRwz02bXLTYbNJPrfkiNJW4aa9e10B5ZAhqbdtz7QMN2Wj\nuV9bEl2EKFm19saNR/HLX57Fk0/+GhEYP34df/nLKP77v/+d1jHnzXN1UR29q+WwYU4g1qxxzwud\nZJ7E3jgPgP8CbgzZrpwRRCRGj3YX2n37gh3j/vuTT3ttSRghJ9XCyEm0lSexYoUT1A4dQj9UqLQM\nN4WdtM42qZYzVYVp0z7GFVfMYejQnQBMmLCWefOGpH3nH01aRxFp9iaKgWTLl94ZfQAP4Ka/XgE8\nBmSxsXT7IohIlJa6LpZB5m5v2AD//Kf/vjUQjkjs2NGNTp3yv2FaW+UkCiHUBM1V19ELZr55EuPG\nwbJlrl17PF544XjWrOnDlVc2x6SOPrqGsrIGqqv7xv9QClqKBBRXXiJpTkJE+ojID4BFuNDUWFW9\nsZBzEkEbqAXNS0yZAl/6UrDahDBEYu3aPnkfagLo02cfdXVlHDwY7i1+oYhEWdkhunY9yK5dLukW\nViFdWJSVudDPG3FWmqmthVtuuYDbbnuWsrIjCxkyyUvEJq2jmCcBiMiPgbm42UxjVHWaqmappV37\nJejFIIhIPPmka4AWdHEoE4nElJQoffvuDd2bKBSRgCO7wYZZSBcW0bxES26+GSZOfJdTT23dx2b8\n+HW8+WbwhNKWLU58WoZlx4xxecJt2wLvMu9I5kl8Czga+D7wXkxrjtpCbcuxfz9s3kyg9bL9ToPd\nvBm++lXXzyZoS4IwCuoKRSSgbVpzFJZINJ+vbDb3ayviFdXNm+dafkyZ8kLcz0yYsJa5c4OLxNy5\nrr6j5RThDh1cH6dcd6ZtC5LlJEpUtUtsWw7v0cNr0VFwrF7tZit09LMUk8eYMS5GWp/EY1d1a1f/\n13+5L1ZQwvEkeheQSISbl4i2ji8kkWj2JPIrJwHuNzRnjus9Bm4J3i9/GW6/HXr33h/3MyNGbKem\npnPc9TSSES8fEeXMM4sj5BSoLUehk86CLl26uFkv77yTeJuZM139xdSp6dlVUeFEaPfu9D4fj0Ly\nJMJuzfHee65uJl6n3nwktjVHvuUkwP0eysthyRL3+he/cEWpX/hC4s+UlCjjxq0PHHJKJhJnnVUc\nyWsTiRjSXfUrWbO/NWvgxhvht791s6HSQST73sTatX3yfvprFFcgFp4nUUihJnCiGvW88jEnAc1F\ndRs2uBXx7rsv9brj48evCxRyUo2ftI4yYYLr4bR3b/z3C4VQRUJEHhSRLSKyKGbsDhFZKiILReRJ\nEekZ894UEVnpvX9+mLbFIxORiJeXaGpyi9jfcEP6awBEGTEieyKxe3cXDh3qQHl5dvaXayoq6kJd\noa7wRKI53JSPOQlozktcdx187Wv+/n8mTFgbyJPYsKHX4TVd4tGli5tp9frrvneZl4TtScykdcfY\nF4DRqnoysBKYAiAio4BLgJHABcC9Im3bUSbdi0EikbjnHrdIyg03ZG5bNj2JtWt7U1m5My/69fjB\nrVAXXripEEUiWoCYjzkJcJ7EE0+4MO93v+vvMyecsIk1a/pSU+NvkapFiwYlDDVFKYapsKGKhKq+\nAuxqMfaiqkbbMc4BBnvPJwGPqeohVV2DE5AkKy1kn3Q9iZNPdm5n7KpZS5bAD37gZjNlo0o3myKx\nbl0fKit3Zmdn7YCwW3MUnkjUtfAk8isnAa6ItV8/+OUv3WJMfigtbWTMmPdYsOAYX9svWnR0SpEo\nhqK6XOckvgQ85z0fBMT2d9zojbUJO3e65HD//sE/27276+mzbJl73dDgkmg/+EH2eiNl15MoLJEI\nuzVH4YmEmwKrmr85iZIS93s499xgnwsyFdaPJzFxogs3HQq/03rOCDDZM7uIyE1Ag6r+Pp3PT5s2\n7fDzSCQSqLYhHitXOi8i3RBMNOQ0Zgz88IfuLsdPl1e/ZFskJkwIceHsNqZfvzp27OhGU5NQUpLd\n1pwHDrjZTYXUyK1z50N07tzAzp1d2bevlG7d8s+TALeWdlDGj1/HvfemXgCmsVF4552BjB+ffLve\nvV1d1cKFpNw2F1RVVVFVVZXRPnIiEiJyOfAxIPY+YCMQ6wcO9sbiEisSANXVVRnZlG6oKUpUJEaP\ndlPyFixIX3Diccwxrvqzvr5jxgvErFnTh09/ekGWLMs9ZWWNdOtWz65dXejbN2C3xRSsWuXWVk7n\ngtSe6d+/lnff7UfXrgfp0KFIel4Dp5yygbffPpr6+g6tWnfEsnp1P/r2raNPn9T5i2heoj2KRCQS\nIRKJHH49PWi7B9om3CTew70Q+SjwbWCSqsbewjwDTBaRUhEZBowA4nRoCYdsiMRrr7kw0913Z7ZY\nfTw6dnQhrfXrM5+sX2jhJgivVqLQQk1RysvrWL26X16GmjKhR496hg3bwdtvJ++7v2jRIE488T1f\n+yz0orqwp8A+CrwGHCci60TkCuBnQHfgbyIyX0TuBVDVJcDjwBJcnuIa1bZb1iPTi8Epp8Cbb8IJ\nJwTr8BqE4cPdBT4T6upKqasro6IiQRvNPCWs1hyFKhL9+9eyalV5XiatM8VPXmLRoqM58cSEgYwj\niBbVFeoiRKGGm1T1sjjDM5NsPwOYEZ5FicnUk+jVC265Ba6+OrthpliGD4d16zLr7b1+fW+GDNmZ\n9dh9rgmrNcfy5c1rGxcS5eW1LF/ev+g8CXB5iaeeOglI3Hhp0aJBXHjh2772d8wxrmaikFq3xJLr\n2U3tgqYml7j2s651Mm6+2SWsw8KJRGaexJo1fQ4vxlJIWLgpGBUVdZ4nUZwiMW/eEJqa4t/N1dd3\nYOXKckaP3uR7n4XcosNEAjd7pWdP92jPjBjhCuEyYe3aPgwZUngd38MIN6kWrkj071/Lxo29itKT\nqKioo1ev/axcGb/lwLJl/ams3EnXrg1x349HIeclTCTIPNTUVmTDkyjEpDWEE27ats2FDsP0DnNF\neXktQFHmJADGj0+clwiStI5inkSBky+xxGOPdf1kGhvTT3oUqkgcd9xW5swZxptv+qum9UPUiyiU\n9iWx9O/vRKIYPQlIvghRkKR1lJEjYdcutxBRoWEigbsY5IMn0aUL9O69j82b04+LrV1bmDmJ44/f\nyp13PsXVV0/muedGZWWfhRpqAg7PbivGnAQkX8508eLgnkRJiau+LkRvwkSC/Ak3AQwZsivtabAH\nDnRkx45uDBy4J8tWtQ/OPns1s2b9jh/84KP8+tenZzwlsZBFokuXBnr0OFC0nsSwYTs4eLADGzce\ndcR4XV0p69f34v3v3xJ4n4Xa7M9EgnwTiZ1pT4Ndv743gwfvLugK21GjNvPEEw8ye/Yp3HLLBRmF\n5gpZJMD1cCrWnIRI/PUl3nlnIMcfv4VOnZoSfDIxhdrsr+hF4uBBWL/exfvzgcrK9D0JN7Op8EJN\nLRk0aA+PP/4QK1ZUcM01l7J/f3o9NQpdJIYN28GAAQW5XL0vJkxonZdIJ2kdZdw4d8NZU2CntOhF\noroaBg9Of9W4tiYTT6JQayTi0bPnAWbOfJju3eu57LIvsn17t0Cfb2iAtWuz18W3PfLAA48xbtz6\n1BsWKC55fWReIp2kdZTSUte/6bXXsmFd+6HoRSJfZjZFqazclfY02EKd2ZSI0tJGfvKTP3D22av5\n9Kf/i3ff7ev7s+++6/pvlflbn8bIQ0aN2sx77x3Frl1dDo9l4klAYU6FNZHIo3wEQGXlTtau7Z1W\nUnbdut5FJRLgYs/f+MbfufrqfzF58hW+p8gWeqjJgI4dmzjppA3Mn+++Ezt3dmX37i4MG7Yj7X0W\nYlFd0YtEvkx/jdKr135KSpRdu7oG/myxeRKxXHrpAn7ykz/wla9M5q67PphyJTsTieIgdirsokVH\nc8IJ72XU1+z002HePLeAWaFQ9CKRb54EuGmwQfMSBw92YPPmngwaVJjTX/1w9tmrmT37QXbu7Mr5\n53+Vb33rk7z99sC425pIFAexRXWZhprAtfY57jgnFIWCiUQeikRl5U7WrAmWl9i48Sj696+ltDTx\nQivFwNChO7n11j9TVXUPxx23lauumszkyZfzhz9AY8ypMZEoDk4+eQNLl/bnwIGOGSWtY7nyStc0\ntFAoWpFQhV/9yk2BzfYCQWHjPIlgIlHMoaZ49Oq1n6uuepWqqrv5/Ofn8uMfu5lMd90Fe/aYSBQL\nXbs2cNxxW1m4cHBWPAlwywWceWYWjGsnFKVI7N4Nl14KP/+5SzKV5NlZSGcarIlEfDp1auLjH3+H\n116Dxx6DuXPdetYHDsDA+JEoo8CYMGEdf/rTCQAcfXTxhmMTkWeXx8x57TW3ityAAfD6664xV76R\nTkGdiURqPvABePRRWLQIZs8uzMZ+RmuiixCdeOJG+z+PQ6gr07UnGhuF++47k4cfhgcegEmTcm1R\n+qTrSZxxRnVIFhUWgwe7h1EcjBu3jvr6TlkJNRUiRSESmzf34Fvf+hRNTcKbb+b/BWDAgFpqajqz\nb18n3wujFEtLDsMISt+++3jf+7Zy8skbcm1KuyTUcJOIPCgiW0RkUcxYbxF5QUSWi8jzInJUzHtT\nRGSliCwVkfOzYcNLLx3HpElXcdpp1Tz88Ky8FwiAkhJl8ODdvr2JxkZhw4ZeBbkinWFkg4cfnsVZ\nZ63OtRntkrBzEjOBj7QY+y7woqq+H3gZmAIgIqOAS4CRwAXAvSLpRwjr6ztwyy0fZerUj/GLXzzO\ntdf+s6C6n1ZW7vQ9w2nTpqPo23cvnTsfCtkqw8hPysv3Wj4iAaGKhKq+ArS8fb0ImOU9nwVc7D2f\nBDymqodUdQ2wEjg13WP/6Efn8957R/Hss/czYcK6dHfTbglSULd2bfG14zAMIzvkIidRoapbAFR1\ns4hUeOODgH/HbLfRG0uLb3/7Rbp0aSjYu4PKyp0JF3Jvic1sMgwjXdpD4jqtGNC0adMOP49EIlS2\nWInQb0I3XxkyZBcvveSv2sslrS0fYRjFRlVVFVVVVRntIxcisUVE+qvqFhEZAGz1xjcCsS06B3tj\ncYkVCYDq6qrsWtnOcTmJ1OGmZ58dzVNPncx99z3WBlYZhtGeiEQiRCKRw6+nT58eeB9tUUwn3iPK\nM8Dl3vMvAk/HjE8WkVIRGQaMAN5oA/vykkGDdrN5c08aGuL/F+7b14kbb5zET396LjNnPlzUi8sY\nhpE+YU+BfRR4DThORNaJyBXAj4APi8hy4EPea1R1CfA4sAR4DrhGNdOl7AuXsrJG+vWr4733jmr1\n3pIlA5g06Sqamkp45pn7OeGETTmw0DCMQiDUcJOqXpbgrfMSbD8DmBGeRYVFdJW6ykqXb1CFWbM+\nwM9/fjbf//5fufjixTm20DCMfKc9JK6NNBkyxK1Sd9ZZblWtG2+8iK1buzN79oNFs5a1YRjhUnQN\n/gqJaMvwOXOGcuGFV3Hssdt54omHTCAMw8ga5knkMZWVO3noodN4+ukx3HHH05xzzqpcm2QYRoFh\nIpHHnHLKBiKRVXznOy9SXl6Xa3MMwyhATCTymIEDa/jxj/+YazMMwyhgLCdhGIZhJMREwjAMw0iI\niYRhGIaREBMJwzAMIyEmEoZhGEZCTCQMwzCMhJhIGIZhGAkxkTAMwzASYiJhGIZhJMREwjAMw0iI\niYRhGIaREBMJwzAMIyEmEoZhGEZCciYSIvINEXlbRBaJyCMiUioivUXkBRFZLiLPi0jrBZwNwzCM\nNiMnIiEiRwPXAmNV9URcy/LPAt8FXlTV9wMvA1NyYV/YVFVV5dqEjMh3++fMWZNrEzIin89/PtsO\n+W9/OuQy3NQB6CYiHYEuwEbgImCW9/4s4OIc2RYq+f5Fy3f7TSRyRz7bDvlvfzrkRCRU9T3gTmAd\nThz2qOqLQH9V3eJtsxmoyIV9hmEYhiNX4aZeOK+hEjga51F8DtAWm7Z8bRiGYbQhotr212ER+TTw\nEVW90nv9eeA04FwgoqpbRGQA8HdVHRnn8yYehmEYaaCqEmT7XK1xvQ44TUQ6A/XAh4C5QB1wOXA7\n8EXg6XgfDvpHGoZhGOmRE08CQESmApOBBmAB8N9AD+Bx4BhgLXCJqu7OiYGGYRhG7kTCMAzDaP/k\nXcW1iHxURJaJyAoRuTHX9gRFRNaIyFsiskBE3si1PakQkQdFZIuILIoZy5uixwT2TxWRDSIy33t8\nNJc2JkJEBovIyyLyjogsFpHrvPG8OP9x7L/WG8+X818mIq97v9XFXvQjL85/EtsDn/u88iREpARY\ngcthvIfLY0xW1WU5NSwAIvIuME5Vd+XaFj+IyJm4XNFvvcJHROR2YIeq3uEJdW9V/W4u7UxEAvun\nArWq+tOcGpcCb/LGAFVdKCLdgXm4WYFXkAfnP4n9l5IH5x9ARLqq6j4R6QC8ClwH/Af5cf7j2X4B\nAc99vnkSpwIrVXWtqjYAj+G+dPmEkEfnXVVfAVoKWt4UPSawH9z/Q7tGVTer6kLveR2wFBhMnpz/\nBPYP8t5u9+cfQFX3eU/LcBN9lPw5//Fsh4DnPm8uVh6DgPUxrzfQ/KXLFxT4m4jMFZErc21MmlQU\nQNHj10RkoYj8uj2GC1oiIkOBk4E55GHRaYz9r3tDeXH+RaRERBYAm4G/qepc8uT8J7AdAp77fBOJ\nQmCiqo4FPgZ81QuH5Dv5E7N03Ascq6on435A7Trs4YVqZgPXe3fkeVV0Gsf+vDn/qtqkqqfgPLhT\nRWQ0eXL+49g+ijTOfb6JxEZgSMzrwd5Y3qCqm7x/twF/wIXQ8o0tItIfDsedt+bYnkCo6jZtTsb9\nCpiQS3uS4fU2mw38TlWjdUN5c/7j2Z9P5z+KqtYAVcBHyaPzD0fans65zzeRmAuMEJFKESnF1Vk8\nk2ObfCMiXb27KkSkG3A+8HZurfKFcGQc8xlc0SMkKXpsRxxhv/fDjvIp2vf/wUPAElW9O2Ysn85/\nK/vz5fyLSL9oOEZEugAfxuVV2v35T2D7snTOfV7NbgI3BRa4GydwD6rqj3Jskm9EZBjOe1BcIumR\n9m6/iDwKRIC+wBZgKvBH4AnyoOgxgf0fxMXHm4A1wFXRGHN7QkQmAv8EFuO+Mwp8D3iDPCg6TWL/\nZeTH+R+DS0yXeI//U9XbRKQP7fz8J7H9twQ893knEoZhGEbbkW/hJsMwDKMNMZEwDMMwEmIiYRiG\nYSTERMIwDMNIiImEYRiGkRATCcMwDCMhJhJGQSEijV4L5Le9NsnfFJFQm8mJyI+9dsy3txj/oog0\nici5MWMXe2OfSvNY54jI6TGvZ6a7L8PwQ66WLzWMsNjr9cZCRPoBvwd6AtNCPOaVuHbR8YqOFuE6\nA7zsvZ4MLMzgWBFc6/N/Z7APw/CNeRJGwaKq24EvA18D8Nq5/FNE3vQep3njs0RkUvRzIvKwiFzY\ncn8xHsNbIvIZb+xpoDswLzrWgldwzdU6eK1YRhAjEiLyIc/zecvrytnJG68WkWkiMs977zgRqQS+\nAnzd+8xEbzfniMirIrLKvAoj25hIGAWNqlYDJSJSjmvLcZ6qjsfd0f/M2+xB3EI+iEhP4HTgz7H7\n8S6+J6rqGFwfnJ+ISH9VvQjYp6pjVfWJeCYAL+Iaw11ETJ8fESkDZgKfUdWTgE7A1TGf3aqq44D7\ngBtUda33/C7veK962w1Q1YnAhcARIS/DyBQTCaMYiOYkSoFfi1vK9AlgJICq/hPXOLIv8FngSVVt\narGPM3GhK1R1K66rZrSDZrKch+IWx5qMW5Ht9zHbvx94V1VXe69nAWfHfPYP3r/zgKFJjvFHz66l\ntNO1DYz8xXISRkEjIscCh1R1m7ds6WZVPdFb0nF/zKa/BT6Pu5hf7mfXMc+TNkBT1Te9hmt1qrqq\nRR49mcDUe/82kvy3Wh/zPC9WfDPyB/MkjEIjtiV4OfBLmsNKRwGbvOdfADrEfG4W8HVAE6yZ/i/g\nUm+1r3LgLJpXWfNzYb4RuKnF2HKg0hMycCJVlWI/tbhEfCJMJIysYp6EUWh0FpH5uNBSA/BbVb3L\ne+9e4EkR+QLwV2Bv9EOqulVEltIc4jkCVf2Dl+h+C9dm+dvewlHgY2UyVX0+9qU3Vi8iVwCzPc9m\nLnB/in3+ydt+EnBtnO2srbORVaxVuGHgFoTCCcBYVa3NtT2G0V6wcJNR9IjIh4AlwD0mEIZxJOZJ\nGIZhGAkxT8IwDMNIiImEYRiGkRATCcMwDCMhJhKGYRhGQkwkDMMwjISYSBiGYRgJ+X9PUdoY3trv\nxQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "\n", "plt.plot(range(1,32),ordDayCounts)\n", "plt.axvspan(0.5, 1.5, color='y', alpha=0.5, lw=0)\n", "plt.axvspan(6.5, 8.5, color='y', alpha=0.5, lw=0)\n", "plt.axvspan(13.5, 15.5, color='y', alpha=0.5, lw=0)\n", "plt.axvspan(20.5, 22.5, color='y', alpha=0.5, lw=0)\n", "plt.axvspan(27.5, 29.5, color='y', alpha=0.5, lw=0)\n", "plt.xlabel('Day of Month')\n", "plt.ylabel('Number of Crashes')\n", "\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It looks like weekend crashes may be less common... Let's look at the data another way." ] }, { "cell_type": "code", "execution_count": 120, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[103, 132, 111, 180, 108, 114, 99, 94, 101]\n" ] } ], "source": [ "weekendCounts = [];\n", "weekendDays = [1,7,8,14,15,21,22,28,29];\n", "for i in weekendDays[::-1]: weekendCounts.append(ordDayCounts.pop(i-1))\n", "\n", "print weekendCounts" ] }, { "cell_type": "code", "execution_count": 149, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD7CAYAAACIYvgKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH4JJREFUeJzt3Xt023eZ5/H3IyVOJNuRL4mTuG7TpC1pg2lS2oSkLYxb\nBiizWcptuAywLJ3O4ewwUA4w0M5lGzgzu9O5wMyBXWZgGIZ2KSy7y6WUU3qjgt7SNG3TxnWaC2mu\nduLGiRXZVnyRnv1DUuI4tiPbulj253VOTvX76qefHvf8pEffu7k7IiIyuwVKHYCIiJSekoGIiCgZ\niIiIkoGIiKBkICIiKBmIiAgwp9QBTIaZaTysiMgkuLuNVl6WyQBA8yPyZ9OmTWzatKnUYYicQ/dm\nfpmNmgcANROJiAhKBiIigpKBAC0tLaUOQWRUujeLx8qx7d3MvBzjFhEpJTMbswNZNQMREVEyEBER\nJQMREUHJQEREKONJZyIyOyWTSTo7O+nrSxAOh2hoaCAYDJY6rLKn0UQiUjZisRjRaCu9vRECgWpS\nqTiVlTFaWpqJRCKlDm/aG280kZKBiJSFZDLJ/fdvJhBYRXV17enyePwEqVQbGzeuVw3hPDS0VETK\nXmdnJ729kbMSAUB1dS29vRE6OztLFNnMoGQgImWhry9BIFA96nOBQDWJxKkiRzSzKBmISFkIh0Ok\nUvFRn0ul4oRC84sc0cyiZCAiZaGhoYHKyhjx+ImzyuPxE1RWxmhoaChRZDODOpBFpGxoNNHUaDSR\niMwY2XkGicQpQqH5mmcwAUoGIiKioaUiIjI+JQMREVEyEBERJQMREUHJQEREUDIQERGUDEREBCUD\nERFByUBERFAyEBERlAxERAQlAxERocDJwMyazOxXZvaymW03s89kymvN7CEz22lmD5pZZNhr7jCz\n3Wa2w8zeXsj4REQkraCrlprZEmCJu28zsyrgOeBm4BNAl7v/rZl9Cah199vNbBXwfWAt0AQ8Alw2\ncolSrVoqIjJxJVu11N2PuPu2zOMeYAfpL/mbge9lTvse8O7M43cBP3T3IXffB+wG1hUyRhERKWKf\ngZldDKwBNgOL3f0opBMGkN2v7gLg4LCXHc6UiYhIARUlGWSaiP4vcFumhjCyjUdtPiIiJTSn0G9g\nZnNIJ4J73P1nmeKjZrbY3Y9m+hU6M+WHgQuHvbwpU3aOTZs2nX7c0tJCS0tLniMXESlv0WiUaDSa\n07kF3/bSzO4Gjrn754aV3QUcd/e7xuhAfhPp5qGHUQeyiEhelGwPZDO7DvgNsJ10U5ADfwZsAX5E\nuhawH/iAu3dnXnMH8IfAIOlmpYdGua6SgYjIBJUsGRSKkoGIyMSVbGipiIiUByUDERFRMhARESUD\nERFByUBERFAyEBERijADWUQkH5LJJJ2dnfT1JQiHQzQ0NBAMBksd1oyheQYiMu3FYjGi0VZ6eyME\nAtWkUnEqK2O0tDQTiUTOfwEBNOlMRMpYMpnk/vs3Ewisorq69nR5PH6CVKqNjRvXq4aQI006E5Gy\n1dnZSW9v5KxEAFBdXUtvb4TOzs4xXikToT6DWUrtr1Iu+voSBALVoz4XCFSTSJwqckQzk5LBLHRu\n+2sXlZV71f4q01I4HCKV6hr1uVQqTihUX+SIZiY1E80yyWSSaLSVQGAVjY3NLFmyjMbGZgKBVUSj\nrSSTyVKHKHKWhoYGKitjxOMnziqPx09QWRmjoaFhjFfKRCgZzDJqf5VyEwwGaWlpJpVqo729lSNH\n9tPe3koq1UZLS7OaN/NEzUSzjNpfpRxFIhE2blxPZ2cnicQpQqF6GhquUCLIIyWDWUbtr1KugsEg\nS5cuLXUYM5aaiWYZtb+KyGg06WwW0mxOkdlJM5DlHNl5Bun21/maZyAyCygZiIiIlqMQEZHxnTcZ\nmNmjuZSJiEj5GnNoqZnNB8LAQjOrBbJViwXABUWITUREimS8eQafBD4LNALPcSYZnAS+UeC4RESk\niM7bgWxmn3b3rxcpnpyoA1lEZOKmPJrIzK4FLmZYTcLd785XgBOlZCAiMnHjJYPzLkdhZvcAlwDb\ngOySlg6ULBmIiEh+5bI20TXAKv0UFxGZuXKZZ9AKLCl0ICIiUjq51AwWAm1mtgXozxa6+7sKFpWI\niBRVLslgU6GDEBGR0tLaRCIis8RURxPFSY8eAqgA5gK97r4gfyGKiEgpnTcZuPvpPRLNzICbgfWF\nDEpERIprUs1EZvaCu19VgHhyfX81E4mITNBUm4neO+wwQHregXZNFxGZQXIZTfQfhz0eAvaRbioS\nESmI7E58fX0JwuGQduIrAo0mEpFpJZc9upUsJmdKC9WZWRPwdeC6TNHjwG3ufiiHN/4OsBE46u5X\nZsruBP4I6Myc9mfu/svMc3cAt5Cugdzm7g+NcV0lA5EZKJlMcv/9mwkEVlFdXUsymaSrq4N9+9pw\nf5UPfegmqqurefzxHeMmCxndVJPBw8C9wD2Zoo8CH3H3t+XwxtcDPcDdI5JB3N2/OuLcKzLvsxZo\nAh4BLhvtW1/JQGRm6ujo4LHHumhsbKanJ8bmzc+wa1ecwcHF9PYeYdmy16iuHmLdunezZMmFp18X\nj58glWpj48b1qiGMY0odyMAid//usON/N7PP5vLG7v6EmS0bLaZRym4GfujuQ8A+M9sNrAOeyeW9\nRKT8jGzuicd7gDCdnYd4+umnOHhwAQsW/C7hcISTJ48TCOxk797fUlNznEWLGk9/8VdX19LeHqGz\ns5OlS5eW9o8qU7kkgy4z+yjwg8zxh4GuKb7vn5jZx4CtwOfdPUZ6K82nh51zGG2vKTJjxWIxHn30\nRTo6kgwOzmHu3CHcO9ixI8GcOSvYuXMePT219PV1ccEFFbj3M39+GLOL6epK0d3dTX19/enrBQLV\nJBIa6DhZuSSDW0j3GXyN9Ezkp4BPTOE9/yfwFXd3M/sr4B+AWyd6kU2bNp1+3NLSQktLyxRCEpFi\nSiaT/OIXm2lrc6CJQKCaoaFudu16iaGhapYvryMUqiSVugizBezbt4+mprnU1NRz8OBxksk59PcP\nnHXNVCpOKFQ/+hvOUtFolGg0mtO5ucxA3g/kbYVSd39t2OG3gZ9nHh8GLhz2XFOmbFTDk4GIlJeO\njg62bu1i4cJ3Eg7XAhCLVREI/A7wCqdO7SGRCJBIhHCfy+BgP0uW1BGJLCQY/A2DgwHmzWs8fb14\n/ASVlTEaGq4o0V80PY38ofzlL395zHPHTAZm9nfAHnf/lxHlnwSWu/vtOcZjDOsjMLMl7n4kc/he\n0vslANwHfN/Mvka6eehSYEuO7yEiZWTfvv2cOFFJKJRgcDBFOLyA48ePMThYSX//QlaurKaq6hi7\ndh2msnIRwWATwaBx6tRJli2bx5w5r5BI1HDkSM9Zo4nUeTx549UMbgS+OEr5t4GXgPMmAzO7F2gB\n6s3sAHAncIOZrQFSpCewfRLA3dvM7EdAGzAI/LGGDInMPLFYjCefbOPw4UUkkwH6+1+jq+t5wuE6\nTpyoIJGAtrb9vPWt1zFv3svs2vUkx44FmTt3DjU1J7nmmnpuuukG+vv7SSROEQrV09BwhRLBFI05\ntNTMWt29eYznXnb31xc0snFoaKlIecrOI+jqivDAAwdYsOAtdHScZHBwLnPmtDMw0MepU+1cf30z\n1dUDrFu3kkOHdvHaa09x/fXNLFt2EUuXLtUX/yRNdmhpwswuc/fdIy52GZDIZ4AiMjt0dnbS2xth\n2bIrWLnyCK+88iInT9YSDi8hHj9FXd0R6uraCQYv5ejROLt3P0NTU4D3ve+DmlBWYOMlg/8KPJAZ\n8fNcpuwa4A4gp3kGIiLD9fUlCASqCQaDrF+/lt7eKD09PYRCAwQCr3HJJUO0tHyEoaF+Dh2Kc/XV\nFaxde41qAkUwZjJw9wfM7N3AnwKfzhS3Au9z9+3FCE5EZpZwOEQqlZ6mVFUVYcOGDaRSr1JZWUVn\n5yDLlzcyNNRPTU0D/f1dLFtWr0RQJOMOLXX3VuDjRYpFRGa4hoYGKiv3Eo+foLq6lvr6eqqr97Fr\n1z7cE3R2LqWzs4tU6nkuuCBGPL4W6NBCdEWQy6QzEZG8CAaDtLQ0E4220t4eAcJ0d+/BvZdly9YR\nCITp64uzf3873d0hoIuBgT3U1Axw881vpq6urtR/woylJaxFpOiyaxLt33+A557r47LLriMej5NI\nJGht3U4wuJJdu16gsbGRSKSRvr4OKiq284UvvE8JYQrGG00UKHYwIiLBYJClS5eyaNEi6utXUFFR\nQX19PaFQkECgkSNHOggErqK6+nXU1S2jqWk9/f1rue++J0kmkySTSTo6Ovjtb/fS0dFBMpks9Z9U\n9nLZ9vJ1wDeBxe7ebGZXAu9y978qeHQiMqMN71AGOHUqQSIxwMBAhIqKEHPnnvmKqqxcQnf3cfbs\n2cMrrxwbtp9BF5WVe7WfwRTlUjP4NunhpIMA7v4S8KFCBiXFoV9XUmrpDuUY8fgJAObPD9Hff4yB\ngTlUVCSoqqoC0vdqLNZBLDbIT3/6a9xX0tjYzJIly2hsbCYQWEU02qp7eApy6UAOu/sWs7OamYYK\nFI8UyblbC+rXlRTfaB3KQ0N7OXWqjzVrfo9gMEgi0UNbWxv9/QlisX6gkqGhnaxePZeqqvS9qv0M\npi6XZHDMzC4hvXw1ZvZ+oKOgUUlBJZNJotFWAoFVNDbWni6Px08QjbZqtygpqkgkwsaN6+ns7CSR\nOMWaNTdw992Pc+zYTkKhBnbubGPevMVcfnkTPT1bCYdbMAvz4outrF9/5l7VfgZTk0sz0aeAfwEu\nN7PDpGcf/5eCRiUFlV0SoLq69qzy6upaenvTv65EiinbobxixXJWrVrFF7/4+7zhDQepqHiGqqoh\nLrooxbx5r7JmzRuAAcLhWhKJCN3dZ+7V9H4G80v3R5S5XPYz2Av8rplVAgF3jxc+LCmk7JIAo9Gv\nK5kO6urq+NjHfo8tW55ly5YhmpoaqKlJr5t54EAbfX1xAoFq+vvT96r2M5i689YMzOw2M1sA9AFf\nM7PnzezthQ9NCiU9gmP0nK5fV1JM4w1iCAaDXHzxMhYtqqG+Pr1SaTAYZM2a5aRS++jqaiMWO0Z7\neyupVJv2M5iinLa9dPd/MrN3APXAx4B7gIcKGpkUzMglAbL060qKKZdBDKPdq1VVVTQ3L6WpaR/X\nXVdLVVWl9jPIg/POQDazl9z9SjP7JyDq7j8xsxfc/arihDhqTJqBPEXnfhDP7Bal0URSaNl9DQKB\nVef8IEml2s4axKB7NX/Gm4GcSzL4LultKJcDq4Eg6aRwdb4DzZWSQX5klwRI7xY1X4uBSdF0dHTw\n2GNdNDaeu39We3srN9xQf9YQUd2r+THZzW2y/hBYA+x19z4zqwc+kc8ApTSyIzhEim2igxh0rxZe\nLqOJUmb2KvA6M1PPoohM2chlKIZLD2KoH/f12ZpCX1+CcDikmkIe5LI20a3AbUATsA1YDzwN3FjY\n0ERkpprKIAbNni+MXPoMtgNrgc3uvsbMLgf+m7u/txgBjhGT+gxEytxkOoYn0vEs55pqn8Epdz9l\nZpjZPHd/xcxW5jlGEZllRi5DEQrVn3eIaHb2/PBlVEBrE+VDLsngkJnVAD8FHjazE8D+woYlIrNB\nrh3D2T6CHTt2cvJkDYsXJ89JGpo9PzW5dCC/J/Nwk5k9BkSAXxY0KhGRjOHNSbFYJa2tHRw7lp6J\nnF3iGnLreJaxjbkchZmtNbN3Di9z91+TXr76DYUOTETk7BV2m3nd695IY+MC+vsjbNv26unlK850\nPDeUOOLyNd7aRHcBbaOUvwz8XWHCkelEm99IqY1cYTcYDLJ6dTPz5x+ko+MQu3dv19pEeTJeM1G1\nu5/TN+Du+81sYQFjkmlAw/dkOhhtclpVVYT169eza9fzXH75SVauvExrE+XBeMmgdpznwvkORKYP\nbX4j08VYk9OCwSCRSIiVK5s0eihPxmsmesTM/tqG7XdpaV8BflX40CSfJtLko81vZLoYuUdylvoI\n8m+8msHngX8F9pjZtkzZamArcGuhA5P8mWiTjza/keli5B7JIyenqYaaP2MmA3fvBT5sZiuA12eK\nX87sfCZlYmBggB//OMrg4EXU19dTU5New2W8Jp+prhsjkk+TmZwmE3fenc7cfa+7/zzzT4mgjMRi\nMe6990G2bq3g0KE6Xnihi82bN9PTExu3yUdVc5lugsEgDQ0NhELz6etL0NnZqdFteZbLDGQpQ9lO\n4IGBS6ivr6OubgkAfX0nePHFVtavXz9mk4+q5jIdDF+ZdGhokLa2oyQSdRrdViBKBjNUthO4vn4x\nBw6cPF0eDtdy7FiE7u7OcZt8VDWXYhhrKerh/VwQZuvW5wiFlrFhw8WnZx1rdFt+5ZQMzCwILB5+\nvrsfKFRQMnXZTuCamhrmz2+nry9OOJzuFA4Eqjl+/CgNDb3jLhWsDUWkkMYa2PDmN1/B44/vOD20\nuaurg3D4SsLhZWzb9iobNqwiGAxqcbo8y2U/g08DdwJHgVSm2IErCxiXTFG2EzgYTK/hsm3bqxw7\nVkUgEKKrq43lywdoaWkZ9xeVNhCRQhlvLst99z3BnDkX0dSULj91Kv3DJhyu5tixKrq7u6mvT9do\nNbotf3KpGdwGrHT30YeXyLQ0cvOQDRtW0d3dzfHjnaxYMcQf/ME7qaioGPP1moEshTTeUtR79oSZ\nP3/gdNn8+WdGtwUCIfr7zzyn0W35c97RRMBBIDaZi5vZd8zsqJm9NKys1sweMrOdZvagmUWGPXeH\nme02sx1m9vbJvKekZTuBU6k22ttbee21Q/T3d9DQcJz3vOct4yaC7HDUo0fDzJtXz6JFTTQ2NhMI\nrCIabdUoDpmy8eayhEL1nDp15rdnTU0DoVCMvr4TpFIJ5s1L37sa3ZZfuex09h1gJfALoD9b7u5f\nPe/Fza4HeoC73f3KTNldQJe7/62ZfQmodffbzWwV8H3Su6o1AY8Al422pZl2Ostdtqkn3Qk8/7xN\nPbFYjJ/85Dc8++wc6utXZX55xVi9upmqqgjt7a3ccEO92mhlSjo6OnjssS4aG5vPee7QoRcZHNzH\nwoVvOT0LvqcnxubNT5NI9HL11VcBfefdFU3ONdWdzg5k/lVk/uXM3Z8ws2Ujim8Gfifz+HtAFLgd\neBfwQ3cfAvaZ2W5gHfDMRN5TzjaRTuCRw1EjkUX09HTT1RXgiSd+w403vkNttJIX4+2BXF3dw5vf\n/GYef7ztrKHNb3xjmNe/fjlz55pGtxVALpvbfBnAzKoyxz1TfM8Gdz+audYRM8vW8S4Anh523uFM\nmRTJ8OGou3cf4fDh1xgYqMKsjpMn9wMPcOmlDYRCF5c6VClz55vLoqHNxZfLaKJm4B6gLnN8DPhP\n7v5ynmJQe880kW3HXbBgAe3tTxMIrCUSWZx5dhU9PQfZu7eV+vqrSxqnzAzn+8LX0ObiyqWZ6FvA\n59z9MQAzawG+DVw7yfc8amaL3f2omS0BsushHAYuHHZeU6ZsVJs2bTr9uKWlhZaWlkmGI1nZ4agn\nT56ksfESjh/vJRbrwGweJ0/+lqVLYcWK1XR1delDKnmhL/zCikajRKPRnM7NpQP5RXdffb6ycV5/\nMfBzd39D5vgu4Li73zVGB/KbSDcPPYw6kItqYGCAe+99kEOHInR1NdDYeAl9fX2cPHmMuXN3ceON\nb+X48Q7WrUuxYsXyUocrIhM01Q7kvWb2l6SbigA+CuS0YJ2Z3Qu0APVmdoD05LW/Af6Pmd0C7Ac+\nAODubWb2I9JbbQ4Cf6xv/OLJzisYGlrI4cN72LPnAO3t+7nwwlrq6gZZvfpaKioqNK5bZIbKpWZQ\nC3wZuD5T9Diwyd1PjP2qwlLNIL+SyST337+ZQGAV1dW1DAwM8OijT9DbW0FVVRc33pieoBaPnyCV\natNaMCJlaryawXmTwXSkZJBfo4357unpYdu2V+noOERz80IikZDGdYuUuUk1E5nZP7r7Z83s54wy\n4sfd35XHGKWERt90vIoNG1axe3cys+l4k4b2icxg4/UZZPsI/r4YgUjpjLfp+IIFc1i58jKN+BCZ\n4cbb9vK5zH9/nS3L9B9c6O4vjfU6KT/jzQZNr/0y9jLXIjIz5NKBHCW9VMQc4DnS8wKedPfPFTy6\nsWNSn0GenbtKaVx9BCIzzJQ6kM3sBXe/ysxuJV0ruNPMXsouPFcKSgaFMdFF7USkvEx1nsEcM1tK\nej7An+c1MplWNBtUZPbKZT+DrwAPAnvc/VkzWwHsLmxYIiJSTJpnICIyS0x2nsEXMxvQfJ3R5xl8\nJo8xiohICY3XZ7Aj89+txQhERERKR81EIiKzxGSbiUZdhiJLy1GIiMwc4zUTaRkKEZFZQs1EIiKz\nxGSbibYzejORAV7KGcgiIpJf4zUTbSxaFCIiUlI5NROZ2WJgbeZwi7t3jnd+oamZSERk4sZrJjrv\nchRm9gFgC/D7pNcnesbM3p/fEEVEpJRyWbX0ReBt2dqAmS0CHnH31UWIb6yYVDMQEZmgKdUMgMCI\nZqGuHF8nIiJlIpclrH9pZg8CP8gcfxB4oHAhiYhIseXagfxe4PrM4ePu/pOCRnX+eNRMJCIyQZPa\n6czMLgUWu/uTI8qvBzrc/bd5jzRHSgYiIhM32T6DfwROjlIeyzwnIiIzxHjJYLG7bx9ZmCm7uGAR\niYhI0Y2XDGrGeS6U70BERKR0xksGW83sj0YWmtmtwHOFC0lERIptvA7kxcBPgAHOfPlfA1QA73H3\nI0WJcPTY1IEsIjJBkxpNNOzFNwDNmcOX3f1XeY5vwpQMREQmbkrJYDpSMhARmbipLkchIiIznJKB\niIgoGYiIiJKBiIigZCAiIigZiIgISgYiIkJum9sUhJntI70CagoYdPd1ZlYL/G9gGbAP+IC7x0oV\no4jIbFHKmkEKaHH3q9x9XabsdtL7K68EfgXcUbLoRERmkVImAxvl/W8Gvpd5/D3g3UWNSERklipl\nMnDgYTN7NrMSKqT3UDgKkFkIr6Fk0YmIzCIl6zMArnP3DjNbBDxkZjtJJ4jhxlyAaNOmTacft7S0\n0NLSUogYRUTKVjQaJRqN5nTutFiozszuBHqAW0n3Ixw1syXAY+5+xSjna6E6EZEJmnYL1ZlZ2Myq\nMo8rgbcD24H7gP+cOe3jwM9KEZ+IyGxTkpqBmS0nvXGOk26q+r67/42Z1QE/Ai4E9pMeWto9yutV\nMxARmSDtZyAiItOvmUhERKYXJQMREVEyEBERJQMREUHJQEREUDIQERGUDEREBCUDERFByUBERFAy\nEBERlAxERAQlAxERQclARERQMhAREZQMREQEJQMREUHJQEREUDIQERGUDEREBCUDERFByUBERFAy\nEBERlAxERAQlAxERQclARERQMhAREZQMREQEJQMREUHJQEREUDIQERGUDEREBCUDERFByUBERFAy\nEBERlAxERAQlAxERQclARESYpsnAzG4ys1fMbJeZfanU8YiIzHTTLhmYWQD4BvAO4PXAh83s8tJG\nNbNFo9FShyAyKt2bxTPtkgGwDtjt7vvdfRD4IXBziWOa0fSBk+lK92bxTMdkcAFwcNjxoUyZiIgU\nyHRMBiIiUmTm7qWO4Sxmth7Y5O43ZY5vB9zd7xp2zvQKWkSkTLi7jVY+HZNBENgJvBXoALYAH3b3\nHSUNTERkBptT6gBGcvekmf0J8BDpZqzvKBGIiBTWtKsZiIhI8akDuUyY2VfN7DPDjn9pZt8advz3\nZvbZSVw3nq8Yh11zmZltz/d1Zfoq5f1pZt81s/dO9NpyNiWD8vEkcC2AmRmwkPSkvKxrgacmcd1C\nVQ1V5Zxdyu3+lBGUDMrHU2Q+bKQ/ZK1A3MwiZlYBXA48b2ZfMLMtZrbNzO7MvtjMPmJmz5jZ82b2\nzcwHlmHPLzSzp8zsnZnjc66T+cXfZmbfMrPWzK+/eZnnrs6c+wLwqYL/35Dpptj35zfMbIeZPQQ0\nDDvvLzPXecnM/jlTtsLMnht2zqXDjyVNyaBMuHsHMGhmTZz5lfUMsAG4BtgO3ABc5u7rgKuAa8zs\n+sxyHh8ErnX3NwIp4CPZa5tZA3A/8Bfu/oCZvW2062ROvxT4urs3AzHgfZnyfwM+5e5XFe7/gkxX\nRb4/35O5zhXAxzmThCB9b77J3a8Ewmb2H9x9L9BtZldmzvkE6ftVhpl2o4lkXE8B15G++f8BaMoc\nx0hX098OvM3MngcMqAQuA1YDVwPPZn5xzQeOZK5ZATxC+ov88UzZWNc5CLzq7tn+gOeAi80sAkTc\n/clM+T3ATfn/82WaK9b9+RbgB5BOQmb2q2ExvNXM/hQIA7Wkayi/AL4DfMLMPk868azN+19f5pQM\nyku2Kt5M+iY/BHye9Iftu0AL8N/d/dvDX5QZqvvv7v7no1xziPSX+k1A9sNmY1xnGdA/rChJ+oOb\nfY3MbsW6P0eVabL8H8Ab3b090wyVvT//H3An8Biw1d1PTOYPnMnUTFRengI2Asc97QRQQ7oq/hTw\nIHCLmVUCmFmjmS0CHgXen3mMmdWa2YWZazpwC3C5mX0xUzbWdWCUL313jwEnzCxbXf/IyHNkVijW\n/fkb4INmFjCzpaSbnyD9xe9Al5lVAe/PBubu/Zn3/ybpxCQjqGZQXrYD9cD/GlEWdvfjwMOZ9ten\nM/1vceCj7r7DzP4CeMjSS4QPkO7kPUh6qQ83sw8DPzOzk+7+z2Z2xcjrkG7LHWt0xy3Av5lZivSE\nQZl9inl/3gi8DBwgM0rJ3WNm9q+Z8uzqBcN9H3g3uj9HpUlnIjIrZPoLFrj7nec9eRZSzUBEZjwz\n+zGwArix1LFMV6oZiIiIOpBFRETJQEREUDIQERGUDEREBCUDERFByUBERID/DwaJbSqvdZkZAAAA\nAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#plt.plot(range(1,32),ordDayCounts)\n", "import random\n", "\n", "ax1 = plt.subplot()\n", "\n", "dayX = [2]*len(ordDayCounts)\n", "wkdX = [1]*len(weekendCounts)\n", "\n", "dayX = [i+(random.random()-0.5)/10 for i in dayX]\n", "wkdX = [i+(random.random()-0.5)/10 for i in wkdX]\n", "\n", "ax1.scatter(dayX+wkdX,ordDayCounts+weekendCounts,s=50,alpha=0.25)\n", "ax1.set_xlim(0.5,2.5)\n", "ax1.set_ylim(0,240)\n", "ax1.set_xticks((1,2))\n", "ax1.set_xticklabels(('Weekend','Weekday'))\n", "ax1.set_ylabel('Collision Case Count')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So there appears to really be minimal difference between weekend & weekday crash counts, at least for this slice of data in January.\n", "\n", "# Modelling collisions as a Poisson process\n", "\n", "A Poisson process is one that occurs randomly over time (or in space). For a process that occurs with expected frequency $\\lambda$ per unit time, the number of events $n$, observed after time $t$ is a random variable Poisson distributed with parameter $\\lambda{t}$.\n", "\n", "The Poisson distribution is defined as follows:\n", "\n", "$\\large P(\\lambda{t},n) = \\frac{\\lambda{t}^{n}e^{-\\lambda{t}}} {n!}$\n", "\n", "If we assume that ${t}$ = 1 day and absorb it into $\\lambda$, then $\\lambda$ becomes the expected number of collisions in one day\n", "\n", "$\\large P(\\lambda,n) = \\frac{\\lambda^{n}e^{-\\lambda}} {n!}$\n", "\n", "We can now determine the joint likelihood of observing the number of collisions each day, under the assumption that the days are iid draws from a poisson distribution. Joint likelihood:\n", "\n", "$\\large Likelihood(\\lambda)={\\displaystyle \\prod_{i \\in Days}^{} P(\\lambda,n_i)}={\\displaystyle \\prod_{i \\in Days}^{} \\frac{\\lambda^{n_i}e^{-\\lambda}}{n_i!}}$\n", "\n", "Since likelihood is only relative, we needn't include the denominator, since the factor of $\\prod_{Days}^{} \\frac{1} {n!}$ will be present at every value of $\\lambda$.\n", "\n", "Numeric stability is better if we calculate this as log-likelihood.\n", "\n", "$RLL = \\sum_{i \\in Days}^{} nlog(\\lambda) - \\lambda$" ] }, { "cell_type": "code", "execution_count": 198, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEPCAYAAAC3NDh4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuUFPWZ//H3A3gBAaMiKBhQEUWNgIarGGnUCDEqXjau\nMbtJvETjalZjzv7UZA2sSVZz/SXxEqOSnMSfhnhJvCSiaHQUhYFBUERBUFEUBVQwKkochuf3x7c6\nNMP0TM/Q1VVd/XmdU8eu6urqpyxmnvnezd0RERFpSaekAxARkfRSkhARkaKUJEREpCglCRERKUpJ\nQkREilKSEBGRomJNEmY21cxWm9nCIu+fYWbPRNsTZnZInPGIiEj7xF2S+C0woZX3XwaOdPehwPeB\nm2KOR0RE2qFLnBd39yfMbEAr79cX7NYD/eKMR0RE2idNbRLnANOTDkJERDaLtSRRKjMbD5wJHJF0\nLCIislniScLMhgA3AhPdfV0r52mSKRGRDnB36+hnK1HdZNG29Rtm/YG7gH9395faupC7Z3abPHly\n4jHo/nR/tXZvtXB/2yrWkoSZ3QbkgN3MbAUwGdgecHe/EbgC2BW43swMaHT3kXHGJCIipYu7d9MZ\nbbz/NeBrccYgIiIdl6beTTUtl8slHUKsdH/VK8v3Btm/v21l5aizqgQz82qJVUQkLcwMT3nDtYiI\nVCklCRERKUpJQkREilKSEBGRopQkRESkKCUJEREpSklCRESKUpIQEZGilCRERKQoJQkRESlKSUJE\nRIpSkhARkaKUJEREpCglCRERKUpJQkREilKSEBGRopQkRESkKCUJEREpSklCRESKUpIQEZGilCRE\nRKQoJQkRESlKSUJERIpSkhARkaKUJEREpCglCRERKSrWJGFmU81stZktbOWcX5rZMjN72syGxRmP\niIi0T9wlid8CE4q9aWafAwa6+yDgPOCGmOMREZF2iDVJuPsTwLpWTpkE/D46dw6ws5n1iTMmEREp\nXdJtEv2A1wr2V0bHREQkBbokHYBI0l58Ea68EgYMgEMOgaFDYf/9wSzpyESSl3SSWAl8smB/r+hY\ni6ZMmfLP17lcjlwuF1dcUiOWLoWjj4azzgr706bBRRfBxRfDpZcmG5tIR9TV1VFXV1e265m7l+1i\nLX6B2d7Afe5+SAvvHQdc4O6fN7PRwM/dfXSR63jcsUptWbIEjjkGvvc9OPPMzcdfew1GjIA77oDP\nfCa5+ETKwcxw9w6Xi2NNEmZ2G5ADdgNWA5OB7QF39xujc64FJgLrgTPdfX6RaylJSNksWRJKEFdd\nBV/+8tbvT58O554LTz0FvXtXPj6Rckl1kignJQkpp2OPheOPh//8z+LnfPvbIUncfz907ly52ETK\naVuTRNK9m0QqrqEBXngBvv711s+78krYsAF+8pPKxCWSRipJSM05+eRQ1XThhW2fu3QpHHEErFgB\nO+4Yf2wi5aaShEg7LFoE9fVw9tmlnb///nDYYXD77fHGJZJWShJSU66+OnRv7dq19M9ccAFcd118\nMYmkmaqbpGa89BKMGgUvvww9e5b+uaYm2G+/UJoYMSK++ETioOomkRL96Edw/vntSxAQejadf75K\nE1KbVJKQmrB+PfTtG6bg2H339n/+7bdDaeLFF6FXr/LHJxIXlSRESjB9Oowe3bEEASExnHwyTJ1a\n3rhE0k5JQmrCXXfBqadu2zUuuACuvx42bSpPTCLVQElCMm/DBnjgATjppG27zvDh0KMHzJ1bnrhE\nqoGShGTejBkwbFh55mA68US4995tv45ItVCSkMwrR1VT3oknwn33ledaItVAvZsk0z7+GPbcExYu\nhH5lWPNw06ZwvdmzYd99t/16InFT7yaRVjzyCAweXJ4EAdCpU5g9VqUJqRVKEpJp5axqylO7hNQS\nVTdJZm3cGAbQzZ0Le+9dvuuuXx+qnFasgE98onzXFYmDqptEinjySejfv7wJAmCnneDII0O3WpGs\nU5KQzHr4YZgwIZ5rq8pJaoWShGTWI4/AUUfFc+3jjw8licbGeK4vkhZKEpJJ778PzzwDhx8ez/X7\n9oWBA2HmzHiuL5IWShKSSU88EdZ+aM/iQu113HHw0EPxXV8kDZQkJJPirGrKy+Wgri7e7xBJmpKE\nZFIlksTo0fDss/DBB/F+j0iSlCQkc9auDYsDjRwZ7/d07Qqf/nToaiuSVUoSkjl1dTB2LGy3Xfzf\npSonyTolCcmcSlQ15SlJSNYpSUjmVDJJqF1Csk5JQjLlzTdh1SoYOrQy35dvl5g1qzLfJ1JpShKS\nKY8+GqqAOneu3HeOG6cqJ8mu2JOEmU00syVmttTMLm3h/Z5mdq+ZPW1mz5rZV+OOSbLrkUdg/PjK\nfqfaJSTLYp0q3Mw6AUuBo4E3gAbgdHdfUnDO5UBPd7/czHoBLwB93H1js2tpqnBp04EHwh/+ENa0\nrpQPPwzrZ69aBd27V+57RUqR9qnCRwLL3P1Vd28EpgGTmp3jQI/odQ/gneYJQqQUa9fCypXwqU9V\n9nu7dYPDDlO7hGRT3EmiH/Bawf7r0bFC1wIHmdkbwDPARTHHJBlVXx/ma+rSpfLfrSonyaoEfpy2\nMgFY4O5HmdlA4CEzG+LuW3UqnDJlyj9f53I5crlcxYKU9Js9G8aMSea7x42DK65I5rtFCtXV1VFX\nxr9Y4m6TGA1McfeJ0f5lgLv7DwvO+Qtwlbs/Ge3/DbjU3ec1u5baJKRVRx8N3/pWmJ210j74APr0\nCVVeO+xQ+e8XKSbtbRINwH5mNsDMtgdOB5qv5/UqcAyAmfUB9gdejjkuyZimJmhoCIPbktC9O+y/\nPyxYkMz3i8Ql1iTh7k3AhcAM4DlgmrsvNrPzzOzc6LTvA4eb2ULgIeD/uPvaOOOS7Fm0KCwEtOuu\nycUwZkyo8hLJklirm8pJ1U3Sml/9KpQkfvOb5GK45Ra47z64/fbkYhBpriLVTVF1Ub5KqKuZ9Wjr\nMyKVlGSjdZ5KEpJFbSYJM/sacCfw6+jQXsDdcQYl0l6zZsW3nnWpBg6EDRvg9deTjUOknEopSVwA\njAXeA3D3ZUDvOIMSaY81a+Dtt8No6ySZhYbz+vpk4xApp1KSxD/c/eP8jpl1IYySFkmF2bPDL+dO\nKZiuUlVOkjWl/Fg9ZmbfBrqa2WeBO4D74g1LpHRpaI/IU5KQrCklSVwGvAU8C5wH3A/8d5xBibTH\n7NnJt0fkjRgBzzwD//hH0pGIlIe6wEpVa2yEXXYJE/vtvHPS0QSHHhq65CY1sE+k0LZ2gW1z7iYz\nW04LbRDuvm9Hv1SkXBYtggED0pMgYHOVk5KEZEEpE/wNL3i9I/AFIMFxrSKbNTSEKp40GTMG/vKX\npKMQKY822yTc/Z2CbaW7/xz4fAViE2nT3LnpSxKjR6vxWrKjlOqmwwp2OxFKFmmYYlyEhgY477yk\no9jSfvvBRx+FdpJ+zVdPEakypfyy/2nB643AK8BpsUQj0g7r18OyZTBkSNKRbMkMRo0KpZyTT046\nGpFt02aScPcKLysvUpoFC+Dgg9O5fsPIkUoSkg1Fk4SZXdLaB939Z+UPR6R0DQ3hl3EajRgBP/lJ\n0lGIbLvWShKa6VVSraEBjj026ShaNmIEzJsHmzalY7oQkY7SYDqpWoMGwd13hyqnNBo4EP76Vxg8\nOOlIpJZVYjDdjsDZwMGEcRIAuPtZHf1SkW21di2sWpXuX8D5dok0xyjSllIKwrcAewATgMcI60m8\nH2dQIm2ZNw8+/Wno3DnpSIrLJwmRalZKktjP3a8A1rv77wgD6UbFG5ZI69I40rq5ESOUJKT6lZIk\nGqP/vmtmnwJ2RosOScLSONK6uUMPheee04ywUt1KSRI3mtkuwBXAvcDzwA9jjUqkDdVQkthpp9C4\nvnBh0pGIdFwpI65/6+5NhPYIzfwqiVu5MkwRvvfeSUfStny7RNoTmkgxpZQklpvZjWZ2tJl1uBuV\nSLnkf+lWw79GtUtItSslSQwGHgYuAF4xs2vN7Ih4wxIpbt48GD687fPSYOTIUDUmUq1KmSr8Q3e/\n3d1PAYYBPQlVTyKJmDeveqpvDj4YVqyAv/896UhEOqakCQPMbJyZXQ88RRhQp1lgJRHum8dIVIMu\nXUIvp6eeSjoSkY5pM0mY2SvAxcBM4BB3P83d74o7MJGWvPIK7Lgj9O2bdCSl06A6qWal9G4a4u7v\nxR6JSAmqqaopb8QIuOOOpKMQ6ZhS2iS2KUGY2UQzW2JmS83s0iLn5MxsgZktMrNHt+X7JNuqqdE6\nTz2cpJrFOomxmXUCriXM+3Qw8EUzG9zsnJ2B64Dj3f1TwBfijEmqWzUmiX33hQ8/DBMSilSbuGe6\nHwksc/dX3b0RmAZManbOGcBd7r4SwN3fjjkmqVKbNoUG4GpptM4zC4lNXWGlGpXSJoGZfZ6tpwq/\nsoSP9gNeK9h/nZA4Cu0PbBdVM3UHfunut5QSl9SWF1+ET3wCdt896UjaL994fcIJSUci0j6lrCdx\nA9ANGA/cDPwLUM4a1i7AYcBRwE7AbDOb7e4vNj9xypQp/3ydy+XI5XJlDEPSrhqrmvJGjIDrr086\nCqkFdXV11NXVle16ba5MZ2YL3X1IwX+7A9Pd/TNtXtxsNDDF3SdG+5cB7u4/LDjnUmBHd/+faP/m\n6Pp3NbuWVqarcZdcAr17w2WXJR1J+735JnzqU/D229UxnYhkx7auTFdKm8RH0X8/NLO+hKnD9yzx\n+g3AfmY2wMy2B04nzCRb6B7gCDPrbGbdCGtVLC7x+lJDqrkkseee0LUrLF+edCQi7VNKm8RfzOwT\nwI+B+YATqp3a5O5NZnYhMIOQkKa6+2IzOy+87Te6+xIzexBYCDQBN7r78x25GcmupiZYsKD6Gq0L\n5bvC7qu5lKWKtFndtMXJZjsQqoYqPhONqptq2/PPw6RJsGxZ0pF03FVXheqmn/406UiklmxrdVPR\nkoSZndLGl/6po18q0l7VXNWUN2IEXFlKn0CRFGmtuinfWa83cDjwSLQ/HpgFKElIxWQhSQwfHqrM\nNm4ME/+JVIOiDdfufqa7nwlsBxzk7qe6+6mE8RLbVSpAEQgD0aq5PQLCGI8994TF6pYhVaSU3k2f\ndPc3C/ZXA/1jikdkK42NYZ3oak8SEKqcNPJaqkkpSeJvZvagmX3VzL4K/JWwUp1IRSxaBAMGQI8e\nSUey7bRSnVSbUmaBvRC4ARgabTe6+zfiDkwkr6Eh/HLNAs0IK9Wm1Oaz2cCmaNPfQVJRDQ3Vt4ZE\nMYceGtokPvooDK4TSbtSVqY7hzBX08mEeZvqzeysuAMTyctSkujaFQYPDr2cRKpBKXM3vQAc7u7v\nRPu7AbPc/YAKxFcYhwbT1aAPP4RevWDdOthhh6SjKY/zz4cDDoCLL046EqkFlZi76R3g/YL996Nj\nIrF7+mk46KDsJAiAUaNgzpykoxApTWsjri+JXr4IzDGzewjzNk0izLMkEru5c7NT1ZQ3apRGXkv1\naK3hOt/h8KVoy7snvnBEttTQAEcfnXQU5XXAAbB2Lbz1VnUuoCS1pWiSyK/vIJKkhga4/PKkoyiv\nTp1C6WjOHDj++KSjEWldKb2bhpvZn81svpktzG+VCE5q27vvhsV6Djww6UjKT+0SUi1KGSdxK/Bf\nwLOEcRIiFTFvXhhX0Llz0pGU36hRcM01SUch0rZSksRb7t58NTmR2GVpfERzo0bBl78MmzaF6ieR\ntColSUyO1p3+G/CP/EGtJyFxa2iA005LOop49O4dZoVdujQMrhNJq1KSxJnAYML04PnqJkfrSUjM\nGhrgxz9OOor45NsllCQkzUpJEiMqPbpa5M03w2jrLK8HnU8SX/lK0pGIFFdKbegsMzso9khECsyZ\nE2Z+tQ5PJpB+6uEk1aCUksRo4GkzW05okzDA3X1IrJFJTauvhzFjko4iXpoRVqpBKUliYuxRiDRT\nXw/f/nbSUcSra9cwBmT+fBg7NuloRFpWyqJDr7r7q8BHhAbr/CYSi40b4amnsrPQUGvGjIHZs5OO\nQqS4UkZcn2hmy4DlwGPAK8D0mOOSGvbss9C/f+gimnVjx8KTTyYdhUhxpTRcf4/QLrHU3fcBjgbq\nY41Katrs2TB6dNJRVMbYsTBrFmipFEmrUpJEY7TgUCcz6+TujwLDY45Lalh9fe0kif79Yfvt4aWX\n2j5XJAmlJIl3zaw78Dhwq5n9Algfb1hSy2qhZ1Ohww9XlZOkVylJYhKh0fqbwAOEtSVOiDMoqV1v\nvw2rV2dz5tdi1C4haVZK76b17t7k7hvd/Xfu/sv8etelMLOJZrbEzJaa2aWtnDfCzBrN7JRSry3Z\nM2dOmNQvizO/FqMkIWnW2vKl79NyV9f8YLqebV3czDoB1xIau98AGszsHndf0sJ5VwMPtiN2yaDZ\ns2urqglg6FBYsQLWrYNddkk6GpEtFS1JuHsPd+/ZwtajlAQRGQksi8ZaNALTCNVXzX0DuBNY0+47\nkEyppUbrvC5dQulJ4yUkjeKeyb4f8FrB/uvRsX8ys77ASe7+K0IpRWpUU1OY+XXUqKQjqTxVOUla\nlTItR9x+DhS2VRRNFFOmTPnn61wuRy6Xiy0oqbznn4c+faBXr6QjqbyxY+Hqq5OOQrKgrq6Ourq6\nsl3PPMZRPGY2Gpji7hOj/csI7Rk/LDjn5fxLoBehe+25zVfDMzOPM1ZJ3k03wcyZ8PvfJx1J5f39\n79CvX2iX2G67pKORLDEz3L3DtTRxVzc1APuZ2QAz2x44Hdjil7+77xtt+xDaJf5Dy6XWplmzaq/R\nOm/nncPaGU8/nXQkIluKNUm4exNwITADeA6Y5u6Lzew8Mzu3pY/EGY+k2+OPw5FHJh1FcjSoTtIo\n1uqmclJ1U7a9/joMGwZr1kCnuMu3KXXLLXDPPXDnnUlHIlmS9uomkZLMnAmf+UztJgiAceNCaUp/\nC0ma1PCPpKRJrVc1QZjsr0cPeO65pCMR2UxJQlJBSSIYPx4efTTpKEQ2U5KQxL39dmiTGDo06UiS\npyQhaaMkIYl74onQ9bVLGoZ2JiyXg8ceg02bko5EJFCSkMSpqmmzfv1gt91g4cKkIxEJlCQkcUoS\nW1KVk6SJkoQk6r33YMmSMAuqBOPHQxmn3hHZJkoSkqhZs2D4cNhhh6QjSY9cLpSumpqSjkRESUIS\npqqmre2xB+y5p+ZxknRQkpBEzZypJNGSXE7tEpIOShKSmPXrw1/LtbYSXSnUeC1poSQhiXnssdAe\n0b170pGkTy4Xxo9s3Jh0JFLrlCQkMTNmwLHHJh1FOu2+O+yzT1jzWyRJShKSmAcfVJJozXHHwfTp\nSUchtU5JQhKxYkWYs+nQQ5OOJL2OOw7uvz/pKKTWKUlIIh56CI45prbXj2jL6NHw6quwcmXSkUgt\n04+oJGLGDJgwIeko0q1LF/jsZ+GBB5KORGqZkoRUXFMTPPxw+AUorVO7hCRNSUIqbv78MKK4X7+k\nI0m/iRNDQm1sTDoSqVVKElJx6vpauj59YNAgePLJpCORWqUkIRWnrq/to15OkiQlCamo996DBQs0\nX1N7fO5zShKSHCUJqahHHgldO7t1SzqS6jFiBKxaFcaWiFSakoRU1N13w6RJSUdRXTp3Dg3Yf/1r\n0pFILVKSkIppbIT77oOTTko6kupz8slw111JRyG1SElCKqauLvTU2WuvpCOpPp/7HMybB2vWJB2J\n1BolCamYP/0JTjkl6SiqU7duoZeTShNSabEnCTObaGZLzGypmV3awvtnmNkz0faEmR0Sd0xSeU1N\noT3i5JOTjqR6nXYa3H570lFIrYk1SZhZJ+BaYAJwMPBFMxvc7LSXgSPdfSjwfeCmOGOSZNTXhzUS\nBg1KOpLqNXFiWMlv1aqkI5FaEndJYiSwzN1fdfdGYBqwRd8Wd693979Hu/WAJmvIIFU1bbsdd4Tj\nj4c770w6EqklcSeJfsBrBfuv03oSOAfQdGYZ464kUS6qcpJKS03DtZmNB84Etmq3kOq2YEHo63+I\nWpu22bHHwqJFWmNCKqdLzNdfCfQv2N8rOrYFMxsC3AhMdPd1xS42ZcqUf77O5XLkcrlyxSkxypci\nzJKOpPrtsAOceGKocrrooqSjkTSqq6ujrq6ubNczdy/bxba6uFln4AXgaOBNYC7wRXdfXHBOf+Bv\nwL+7e9Fl383M44xV4rFpE+y3H/zxj2F6Cdl206fDlVfC7NlJRyLVwMxw9w7/iRZrdZO7NwEXAjOA\n54Bp7r7YzM4zs3Oj064AdgWuN7MFZjY3zpiksurqoHt3GD486Uiy45hjwjxOixYlHYnUglhLEuWk\nkkR1+tKXYORIVY2U2+TJsHYtXHNN0pFI2m1rSUJJQmKzbh3ssw+89BLstlvS0WTLa6/B0KHhvzvt\nlHQ0kmaprm6S2nbbbTBhghJEHD75STjiCJg2LelIJOuUJCQ2U6fC2WcnHUV2ff3rcMMNSUchWack\nIbFYsADeeSc0sko8JkyAt94Ks8OKxEVJQmIxdSqceSZ00r+w2HTuDOeeC7/+ddKRSJap4VrK7qOP\nQp35U0/BgAFJR5Ntq1bBgQfCK6/AzjsnHY2kkRquJXWmToWxY5UgKmGPPcI6E9dfn3QkklUqSUhZ\nffxxGGF9xx0walTS0dSGxYth3LjQ1bhHj6SjkbRRSUJS5dZbYf/9lSAq6cADQweBa69NOhLJIpUk\npGyamsIvrF//GsaPTzqa2rJkCRx5pEoTsrVtLUnEPQus1JA774RevSCXC/8wpbLOOMO59lq4/PKk\nI5EsUUlCysIdhg2D//1f+Pznk46mNqk0IS1Rm4Skwl/+EtaLOO64pCOpXYMHw2c/C7/4RdKRSJao\nJCHbbMMGGDIEfvazsAazJGf58rBux9y5sO++SUcjaaCShCTuxz+Ggw9WgkiDffaBSy8N8zrpbyop\nB5UkZJu89FLo7qrR1emxcWMoTXzrW/Bv/5Z0NJI0rSchiXEPbRC5XPjrVdJj3rxQslu0KPQ4k9ql\n6iZJzJ//HJbR/OY3k45Emhs+HM44Ay65JOlIpNqpJCEdsmpV+EV0661hSghJnw8+gEMPhe98B776\n1aSjkaRoMJ1U3Mcfw6mnhmmqlSDSq3t3uPfe8IwOOADGjEk6IqlGKklIu33967B6Ndx1l9aLqAb3\n3w9f+xrU14cp3KW2qCQhFXXTTfDYYzBnTkgQmn4j/dydiy+Gk06CmTOhW7ekI5JqopKElGzGjNCl\ncubMUH0h1cMdzjoLXn0V7r4bevZMOiKpFPVukoq4446QIP70JyWIamQGN98cpu4YPx7WrEk6IqkW\nShLSpptugosugocegiOOSDoa6ajOneG66+CEE8JzXL486YikGqhNQopqbITvfQ9uuQUefzysOCfV\nzQymTAkD7MaMgWuugS98IemoJM3UJiEtWrQIvvIV2H13+M1voG/fpCOScquvD+Mnhg0Lq9ppZHY2\nqU1CyurDD+EHPwj11uefD9OnK0Fk1ejRsGBBeL5DhoRE8dFHSUclaRN7kjCziWa2xMyWmlmLM/yY\n2S/NbJmZPW1mw+KOSbb27rshOeyzT5isb948OOecUD0h2dW1a5ji/d57Q5vTwIHw05/C++8nHZmk\nRaxJwsw6AdcCE4CDgS+a2eBm53wOGOjug4DzgBvijCmt6urqKv6djY2hW+t554VfDkuXwqOPhh5M\n5Z7RNYn7q6Rqv7/hw+Gee8LAuzlzYK+94LTTwr+FBx+sSzq8WFX7s4tb3CWJkcAyd3/V3RuBacCk\nZudMAn4P4O5zgJ3NrE/McaVOJf6hfvxxKCFcdx186Uuwxx7w3e+GBukFC+B3v4ODDornu7P+g5iV\n+xs2DG6/PUwBf8wx+d5QdRx1FEyeHEobb72VdJTllZVnF5e4ezf1A14r2H+dkDhaO2dldGx1vKFl\nh3uoS37vvVBttGZNmDZj9eoweOrFF2HZMnj55VBiGD06tDlcfbWmaZCW9eoV5uY699wwDfy4cfDE\nE6G327PPhu60Bx4Y/sDo1y9sffvCbrvBrrvCLruEAXtdu2rqlmpXVV1gTzgh6QiKa63jVeF7zc9z\nD9uyZTBrVni9adPmralp89bYuHnbsCFsH30UGpu33z78UPbsCb17Q58+YdtrrzBl9KBB4Qe6e/d4\n7l+yq2vXsG5Ifv1y9/AHyOLFocTxxhvw9NOhqmrt2rCtWxf+aNmwIXy+WzfYYYfN23bbQZcuYevc\nOWydOm3ezDb/F7b8b/NjzV8311a72gsvhHY4aVmsXWDNbDQwxd0nRvuXAe7uPyw45wbgUXf/Y7S/\nBBjn7qubXUv9X0VEOiDNE/w1APuZ2QDgTeB04IvNzrkXuAD4Y5RU3m2eIGDbblJERDom1iTh7k1m\ndiEwg9BIPtXdF5vZeeFtv9Hd7zez48zsRWA9cGacMYmISOmqZsS1iIhUXmr6HZjZVDNbbWYLC479\nyMwWR4Ps7jKzngXvXR4NwFtsZscmE3Xp2nN/ZjbAzD40s/nRdn1ykbetyL1daWbPmNkCM3vAzPYo\neC8Lz67F+6u2Zwct31/Be98ys01mtmvBsap/fgXvbXF/WXl+ZjbZzF4vuI+JBe+17/m5eyo24Ahg\nGLCw4NgxQKfo9dXAVdHrg4AFhOqyvYEXiUpFad3aeX8DCs9L+1bk3roXvP4G8KuMPbti91dVz67Y\n/UXH9wIeAJYDu0bHDszC82vl/jLx/IDJwCUtnNvu55eakoS7PwGsa3bsYXffFO3WEx4qwInANHff\n6O6vAMvYevxFqrTz/gCqpqG+yL19ULC7E5C/z6w8u2L3B1X07KDl+4v8X+C/mh2bRAaeX6Sl+4Ps\nPL+W7qPdzy81SaIEZwH3R6+LDcCrZmcB0wv2946KiY+aWVWu4mBm3zezFcAZwHejw5l5dkXuD7Lx\n7E4EXnP3Z5u9lYnn18r9QQaeX+TCqCr7ZjPbOTrW7udXFUnCzL4DNLr7H5KOJQ4F93dbdOgNoL+7\nHwZ8C7jNzKpuGJy7/7e79wduJVTJZEqR+3uTKn92ZtYV+DahyiJzitxf/q/uTPzsAdcD+7r7MGAV\n8NOOXij81FZoAAAEYUlEQVT1ScLMvgocR/hrLW8lUDihxF7RsarT0v25e6O7r4tezwdeAvZPJMDy\nuA04JXqdmWdX4DbgVAB3/zgDz24gob76GTNbTnhG882sN+FZ9S84txqfX0v395SZ9c7Kz567v+VR\nIwRwE5urlNr985e2JGEU1KNFLfL/BZzo7v8oOO9e4HQz297M9gH2A+ZWNNKOKen+zKyXhRl0MbN9\nCff3coVjba/m91a4jt1JwJLodVaeXfP7Wxwdr8ZnBwX35+6L3H0Pd9/X3fchzLl2qLuvITy/f63m\n59fa/WXh+QEU9iYk/IG2KHrd/p+/pFvmC1rdbyMU9f4BrCAMqlsGvArMj7brC86/nNAyvxg4Nun4\ny3l/BQ91PjAPOC7p+Dtwb3cCzwJPA/cAe2bs2bV4f9X27IrdX7P3Xybq/ZOV51fs/rLy/Agzay+M\n/n3eDfTp6PPTYDoRESkqbdVNIiKSIkoSIiJSlJKEiIgUpSQhIiJFKUmIiEhRShIiIlKUkoRUPTN7\nv0zXmWxml5Rw3m/N7JS2zos7DpFKUJKQLNBgH5GYKElIZpjZTmb2sJnNixYEOjE6PiBaYOW3ZvaC\nmf0/MzvazJ6I9ocXXGaYmc2Kjp9TcO1ro2vMAHoXHL/CzOaY2UIzu6GFmHqa2SsF+93MbIWZdTaz\nc8xsroWFi+4wsx1b+PyjZnZY9Hq3aK4hzKyThUWr5kQzfX6tDP8LRbaiJCFZsgE4yd2HA0ex5cyX\nA4Efu/sBwGDgi+5+BGHurO8UnHcIkAMOB75rZnuY2cnAIHc/EPhK9F7eNe4+yt2HAN3M7POFAbn7\ne8ACMxsXHToeeMDdm4C73H2kux9KmNvq7BLuMV9qOht4191HESZvO9fMBpTweZF2UZKQLDHgKjN7\nBngY6BvNXAqw3N2fj14/B/wtev0sYTWyvHs8zOT6DvAIMAo4EvgDgLu/GR3PO9rM6qOlI8cDB7cQ\n1+3Av0avTwf+GL0eYmaPR589o8hnizkW+LKZLQDmALsCg9rxeZGSdEk6AJEy+hLQizCj56aoaiZf\nhVM4i/Cmgv1NbPlzUNi+YWy54twWzGwH4DrgMHd/w8wmF3xfoXuBH5jZLsCn2ZxkfkuYAXiRmX0F\nGNfCZzey+Y+5wmsb8A13f6hYfCLloJKEZEF+iuSdgTVRghjPliWEUpeknBRNo7wb4Zd2A/A4YXrs\nTma2J6HEAOGXtgPvRAvT/EtLF3T39YQZRX8B3OebZ9XsDqwys+0ICa4lrwD5NpMvFBx/EPgPM+sC\nYGaDosV0RMpKJQnJgvwv3VuB+6LqpnlEazw0O6f56+YWAnXAbsCV7r4K+LOZHUWoploBzAJw97+b\n2c3R8TdpfV7+PxKqnQpLC1dEn1lDqDLq0cLnfgLcHjVM/7Xg+M2EhXPmm5lF1ziple8X6RBNFS4i\nIkWpuklERIpSkhARkaKUJEREpCglCRERKUpJQkREilKSEBGRopQkRESkKCUJEREp6v8DvGwgOvtJ\nLdUAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "MLE = 133.64\n" ] } ], "source": [ "from math import log,exp\n", "import numpy as np\n", "\n", "totDayCounts = ordDayCounts+weekendCounts\n", "lambdavals=np.linspace(120,150,100)\n", "def like(lam,n): return n*log(lam)-lam \n", "\n", "rloglik=[sum([like(lamval,n) for n in totDayCounts]) for lamval in lambdavals]\n", "rloglik = [i-max(rloglik) for i in rloglik]\n", "rlik = [exp(i) for i in rloglik]\n", "rlikeighth = [j for (i,j) in zip(rlik,lambdavals) if i>0.125]\n", "rliksixtnth = [j for (i,j) in zip(rlik,lambdavals) if i>0.0625]\n", "\n", "ax2 = plt.subplot()\n", "ax2.plot(lambdavals,rlik)\n", "ax2.set_ylim(0,1.2)\n", "ax2.set_xlabel('lambda value')\n", "ax2.set_ylabel('relative likelihood')\n", "ax2.plot(rlikeighth,[0.125]*len(rlikeighth), color='k')\n", "ax2.plot(rliksixtnth,[0.0625]*len(rliksixtnth), color='k')\n", "\n", "plt.show()\n", "\n", "print 'MLE = {:03.2f}'.format([j for (i,j) in zip(rlik,lambdavals) if i==1.0][0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Interpreting likelihood\n", "\n", "So the maximum likelihood estimate (MLE) for $\\lambda$ is 133.64. Of course, this was really a roundabout way of calculation if all we really wanted was the point estimate for MLE, since the point estimate is of course nothing but the mean of the observations." ] }, { "cell_type": "code", "execution_count": 201, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "133.74193548387098" ] }, "execution_count": 201, "metadata": {}, "output_type": "execute_result" } ], "source": [ "float(sum(totDayCounts))/len(totDayCounts)\n", "\n", "#The number is slightly off due to the discretization of lambda above" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, the point estimate doesn't provide us with any sense of the uncertainty of our estimate. The full likelihood plot does. The lines for one eigth and one sixteenth likelihoods have been plotted. One eighth is a traditional cutoff on likelihood plots. The interpretation of these parameter values is that within this range there is no parameter value better than eight times better supported by the evidence." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Auto collision lawyer billboards\n", "\n", "Let's say we are an auto collision lawyer looking to place a billboard somewhere in LA county. We'd like to know what street to put our billboard on. Maybe it makes sense to put the billboard where many collisions are occuring. In that case, we'd like to get the street name where the most collisions are occurring.\n", "\n", "There are many ways to get this information about street names. One would be to use a web API to provide information about a location based on lat/lon coordinates. The Google Maps Geocoding API is a very feature-rich example. It is limited to 2500 queries per day. Project OSRM is another no-frills example. \n", "\n", "We can ask the OSRM server which street is closest to the 2nd collision in the database with the following:" ] }, { "cell_type": "code", "execution_count": 112, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "'Gardena Freeway (CA 91)'" ] }, "execution_count": 112, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import urllib2\n", "import ast\n", "\n", "sampX=data.POINT_X[1]\n", "sampY=data.POINT_Y[1]\n", "page = urllib2.urlopen(\"http://router.project-osrm.org/nearest?loc=\"+str(sampY)+','+str(sampX))\n", "stuff = page.next()\n", "stuffdict = ast.literal_eval(stuff)\n", "stuffdict['name']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Don't reinvent the wheel\n", "\n", "Luckily, however, there already exists a column in the database with the streetname, and the nearest cross-street. Let's find the most collision-able streets:" ] }, { "cell_type": "code", "execution_count": 162, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "RT 5 114\n", "RT 10 106\n", "RT 405 105\n", "RT 101 89\n", "RT 110 60\n", "RT 210 54\n", "RT 60 54\n", "RT 605 50\n", "RT 710 36\n", "RT 91 35\n", "CRENSHAW BL 35\n", "WESTERN AV 35\n", "WILSHIRE BL 33\n", "VERMONT AV 31\n", "SANTA MONICA BL 30\n", "OLYMPIC BL 30\n", "FLORENCE AV 30\n", "RT 14 28\n", "SEPULVEDA BL 28\n", "VICTORY BL 28\n", "Name: PRIMARYRD, dtype: int64" ] }, "execution_count": 162, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#let's look at the top 200\n", "rdCounts = data.PRIMARYRD.value_counts()[:20]\n", "rdCounts" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Unsurprisingly, the freeways have the most collisions, followed by Crenshaw, Wilshire & Western." ] } ], "metadata": { "kernelspec": { "display_name": "Python [Root]", "language": "python", "name": "Python [Root]" }, "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 }