{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "%matplotlib inline\n", "%config InlineBackend.figure_format = 'svg'\n", "\n", "import numpy as np\n", "import pandas as pd\n", "from scipy.stats import beta, gamma, binom\n", "\n", "import matplotlib.pyplot as plt\n", "plt.style.use('ggplot')\n", "plt.rcParams['figure.figsize'] = 9, 6" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Simulation exercise\n", "## Simulate a queuing system\n", "\n", "\n", "A clinic has three doctors. Patients come into the clinic at random, starting at 9 a.m., according to a Poisson process with time parameter 10 minutes: that is, the time after opening at which the first patient appears follows an exponential distribution with expectation 10 minutes and then, after each patient arrives, the waiting time until the next patient is independently exponentially distributed, also with expectation 10 minutes. When a patient arrives, he or she waits until a doctor is available. The amount of time spent by each doctor with each patient is a random variable, uniformly distributed between 5 and 25 minutes. The office stops admitting new patients at 4 p.m. and closes when the last patient is through with the doctor.\n", "\n", "1. Simulate this process once. How many patients came to the office? How many had to wait for a doctor? What was the average wait time of those who did have to wait? What was the longest wait? How many minutes past 4 p.m. did the office close?\n", "\n", "2. Simulate the process 2000 times and estimate the quartiles for each of the sum- maries in part (a). Summarize your results in a 5 × 3 table with the median in the middle column.\n", "\n", "3. Use your 2000 simulations to estimate the probability, on a particular day, that\n", " 1. not a single patient will have to wait to see a doctor;\n", " 2. the office will close on time;\n", " 3. the office will close more than 30 minutes late." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "n_minutes = 420\n", "\n", "class Patient:\n", " \n", " def __init__(self):\n", " self.wait_time = 0\n", " self.waited = False\n", " \n", " def __str__(self):\n", " return ''.format(self.wait_time)\n", " \n", " def __repr__(self):\n", " return self.__str__()\n", " \n", "class Doctor:\n", " \n", " def __init__(self):\n", " self.time_left = None\n", " self.is_free = True\n", " \n", " def take_appt(self):\n", " self.time_left = np.ceil(np.random.uniform(5,25)).astype(int)\n", " self.is_free = False\n", " \n", " def __str__(self):\n", " if self.is_free: \n", " return ''\n", " return ''.format(self.time_left)\n", " \n", " def __repr__(self):\n", " return self.__str__()\n", " \n", " def __nonzero__(self):\n", " return self.is_free\n", " \n", " def __bool__(self):\n", " return self.is_free\n", " \n", " \n", "def get_arrival_times(rate=10):\n", " times = []\n", " t = int(np.random.exponential(rate))\n", " while t < n_minutes:\n", " times.append(t)\n", " t = t + int(np.random.exponential(rate))\n", " return frozenset(times)\n", " \n", " \n", "def run_simulation(time_param=10, verbose=False):\n", " \n", " def get_open_doc():\n", " for doc in doctors:\n", " if doc.is_free:\n", " return doc\n", " \n", " arrival_times = get_arrival_times(time_param) \n", " doctors = [Doctor() for _ in range(3)]\n", " patients = []\n", " wait_times = []\n", " \n", " stats = {\n", " 'n_patients': 0,\n", " 'n_waiters': 0\n", " }\n", " \n", " # run the simulation\n", " t = 0 \n", " while t < n_minutes or not all(doctors): # in case any doctors are working\n", " if verbose:\n", " print(t)\n", " \n", " # add a patient upon arrival\n", " if t in arrival_times and t < n_minutes:\n", " if verbose:\n", " print('patient arrived!')\n", " stats['n_patients'] += 1\n", " patients.append(Patient())\n", " \n", " # decrement counters, free up doctors whom are done\n", " for doc in doctors:\n", " if not doc.is_free:\n", " doc.time_left -= 1\n", " if not doc.is_free and doc.time_left <= 0:\n", " doc.is_free = True\n", " doc.time_left = None\n", " \n", " # send patient to doctor if one is open\n", " while any(doctors) and patients:\n", " p = patients.pop(0)\n", " wait_times.append(p.wait_time)\n", " if p.waited:\n", " stats['n_waiters'] += 1\n", " \n", " d = get_open_doc()\n", " d.take_appt()\n", " \n", " # increment wait times\n", " for patient in patients:\n", " patient.wait_time += 1\n", " if not patient.waited:\n", " patient.waited = True\n", " \n", " if verbose:\n", " print('patients:', patients)\n", " print('doctors: ', doctors)\n", " print()\n", "\n", " t += 1\n", " \n", " stats['avg_wait'] = np.mean([x for x in wait_times if x != 0]) if any(wait_times) else 0\n", " stats['longest_wait'] = max(wait_times)\n", " stats['mins_past_close'] = t - 420\n", " stats['wait_times'] = wait_times\n", " \n", " return stats" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0\n", "patients: []\n", "doctors: [, , ]\n", "\n", "1\n", "patients: []\n", "doctors: [, , ]\n", "\n", "2\n", "patients: []\n", "doctors: [, , ]\n", "\n", "3\n", "patients: []\n", "doctors: [, , ]\n", "\n", "4\n", "patients: []\n", "doctors: [, , ]\n", "\n", "5\n", "patients: []\n", "doctors: [, , ]\n", "\n", "6\n", "patients: []\n", "doctors: [, , ]\n", "\n", "7\n", "patients: []\n", "doctors: [, , ]\n", "\n", "8\n", "patients: []\n", "doctors: [, , ]\n", "\n", "9\n", "patients: []\n", "doctors: [, , ]\n", "\n", "10\n", "patients: []\n", "doctors: [, , ]\n", "\n", "11\n", "patients: []\n", "doctors: [, , ]\n", "\n", "12\n", "patients: []\n", "doctors: [, , ]\n", "\n", "13\n", "patients: []\n", "doctors: [, , ]\n", "\n", "14\n", "patients: []\n", "doctors: [, , ]\n", "\n", "15\n", "patients: []\n", "doctors: [, , ]\n", "\n", "16\n", "patients: []\n", "doctors: [, , ]\n", "\n", "17\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "18\n", "patients: []\n", "doctors: [, , ]\n", "\n", "19\n", "patients: []\n", "doctors: [, , ]\n", "\n", "20\n", "patients: []\n", "doctors: [, , ]\n", "\n", "21\n", "patients: []\n", "doctors: [, , ]\n", "\n", "22\n", "patients: []\n", "doctors: [, , ]\n", "\n", "23\n", "patients: []\n", "doctors: [, , ]\n", "\n", "24\n", "patients: []\n", "doctors: [, , ]\n", "\n", "25\n", "patients: []\n", "doctors: [, , ]\n", "\n", "26\n", "patients: []\n", "doctors: [, , ]\n", "\n", "27\n", "patients: []\n", "doctors: [, , ]\n", "\n", "28\n", "patients: []\n", "doctors: [, , ]\n", "\n", "29\n", "patients: []\n", "doctors: [, , ]\n", "\n", "30\n", "patients: []\n", "doctors: [, , ]\n", "\n", "31\n", "patients: []\n", "doctors: [, , ]\n", "\n", "32\n", "patients: []\n", "doctors: [, , ]\n", "\n", "33\n", "patients: []\n", "doctors: [, , ]\n", "\n", "34\n", "patients: []\n", "doctors: [, , ]\n", "\n", "35\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "36\n", "patients: []\n", "doctors: [, , ]\n", "\n", "37\n", "patients: []\n", "doctors: [, , ]\n", "\n", "38\n", "patients: []\n", "doctors: [, , ]\n", "\n", "39\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "40\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "41\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "42\n", "patients: []\n", "doctors: [, , ]\n", "\n", "43\n", "patients: []\n", "doctors: [, , ]\n", "\n", "44\n", "patients: []\n", "doctors: [, , ]\n", "\n", "45\n", "patients: []\n", "doctors: [, , ]\n", "\n", "46\n", "patients: []\n", "doctors: [, , ]\n", "\n", "47\n", "patients: []\n", "doctors: [, , ]\n", "\n", "48\n", "patients: []\n", "doctors: [, , ]\n", "\n", "49\n", "patients: []\n", "doctors: [, , ]\n", "\n", "50\n", "patients: []\n", "doctors: [, , ]\n", "\n", "51\n", "patients: []\n", "doctors: [, , ]\n", "\n", "52\n", "patients: []\n", "doctors: [, , ]\n", "\n", "53\n", "patients: []\n", "doctors: [, , ]\n", "\n", "54\n", "patients: []\n", "doctors: [, , ]\n", "\n", "55\n", "patients: []\n", "doctors: [, , ]\n", "\n", "56\n", "patients: []\n", "doctors: [, , ]\n", "\n", "57\n", "patients: []\n", "doctors: [, , ]\n", "\n", "58\n", "patients: []\n", "doctors: [, , ]\n", "\n", "59\n", "patients: []\n", "doctors: [, , ]\n", "\n", "60\n", "patients: []\n", "doctors: [, , ]\n", "\n", "61\n", "patients: []\n", "doctors: [, , ]\n", "\n", "62\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "63\n", "patients: []\n", "doctors: [, , ]\n", "\n", "64\n", "patients: []\n", "doctors: [, , ]\n", "\n", "65\n", "patients: []\n", "doctors: [, , ]\n", "\n", "66\n", "patients: []\n", "doctors: [, , ]\n", "\n", "67\n", "patients: []\n", "doctors: [, , ]\n", "\n", "68\n", "patients: []\n", "doctors: [, , ]\n", "\n", "69\n", "patients: []\n", "doctors: [, , ]\n", "\n", "70\n", "patients: []\n", "doctors: [, , ]\n", "\n", "71\n", "patients: []\n", "doctors: [, , ]\n", "\n", "72\n", "patients: []\n", "doctors: [, , ]\n", "\n", "73\n", "patients: []\n", "doctors: [, , ]\n", "\n", "74\n", "patients: []\n", "doctors: [, , ]\n", "\n", "75\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "76\n", "patients: []\n", "doctors: [, , ]\n", "\n", "77\n", "patients: []\n", "doctors: [, , ]\n", "\n", "78\n", "patients: []\n", "doctors: [, , ]\n", "\n", "79\n", "patients: []\n", "doctors: [, , ]\n", "\n", "80\n", "patients: []\n", "doctors: [, , ]\n", "\n", "81\n", "patients: []\n", "doctors: [, , ]\n", "\n", "82\n", "patients: []\n", "doctors: [, , ]\n", "\n", "83\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "84\n", "patients: []\n", "doctors: [, , ]\n", "\n", "85\n", "patients: []\n", "doctors: [, , ]\n", "\n", "86\n", "patients: []\n", "doctors: [, , ]\n", "\n", "87\n", "patients: []\n", "doctors: [, , ]\n", "\n", "88\n", "patients: []\n", "doctors: [, , ]\n", "\n", "89\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "90\n", "patients: []\n", "doctors: [, , ]\n", "\n", "91\n", "patients: []\n", "doctors: [, , ]\n", "\n", "92\n", "patients: []\n", "doctors: [, , ]\n", "\n", "93\n", "patients: []\n", "doctors: [, , ]\n", "\n", "94\n", "patients: []\n", "doctors: [, , ]\n", "\n", "95\n", "patients: []\n", "doctors: [, , ]\n", "\n", "96\n", "patients: []\n", "doctors: [, , ]\n", "\n", "97\n", "patients: []\n", "doctors: [, , ]\n", "\n", "98\n", "patients: []\n", "doctors: [, , ]\n", "\n", "99\n", "patients: []\n", "doctors: [, , ]\n", "\n", "100\n", "patients: []\n", "doctors: [, , ]\n", "\n", "101\n", "patients: []\n", "doctors: [, , ]\n", "\n", "102\n", "patients: []\n", "doctors: [, , ]\n", "\n", "103\n", "patients: []\n", "doctors: [, , ]\n", "\n", "104\n", "patients: []\n", "doctors: [, , ]\n", "\n", "105\n", "patients: []\n", "doctors: [, , ]\n", "\n", "106\n", "patients: []\n", "doctors: [, , ]\n", "\n", "107\n", "patients: []\n", "doctors: [, , ]\n", "\n", "108\n", "patients: []\n", "doctors: [, , ]\n", "\n", "109\n", "patients: []\n", "doctors: [, , ]\n", "\n", "110\n", "patients: []\n", "doctors: [, , ]\n", "\n", "111\n", "patients: []\n", "doctors: [, , ]\n", "\n", "112\n", "patients: []\n", "doctors: [, , ]\n", "\n", "113\n", "patients: []\n", "doctors: [, , ]\n", "\n", "114\n", "patients: []\n", "doctors: [, , ]\n", "\n", "115\n", "patients: []\n", "doctors: [, , ]\n", "\n", "116\n", "patients: []\n", "doctors: [, , ]\n", "\n", "117\n", "patients: []\n", "doctors: [, , ]\n", "\n", "118\n", "patients: []\n", "doctors: [, , ]\n", "\n", "119\n", "patients: []\n", "doctors: [, , ]\n", "\n", "120\n", "patients: []\n", "doctors: [, , ]\n", "\n", "121\n", "patients: []\n", "doctors: [, , ]\n", "\n", "122\n", "patients: []\n", "doctors: [, , ]\n", "\n", "123\n", "patients: []\n", "doctors: [, , ]\n", "\n", "124\n", "patients: []\n", "doctors: [, , ]\n", "\n", "125\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "126\n", "patients: []\n", "doctors: [, , ]\n", "\n", "127\n", "patients: []\n", "doctors: [, , ]\n", "\n", "128\n", "patients: []\n", "doctors: [, , ]\n", "\n", "129\n", "patients: []\n", "doctors: [, , ]\n", "\n", "130\n", "patients: []\n", "doctors: [, , ]\n", "\n", "131\n", "patients: []\n", "doctors: [, , ]\n", "\n", "132\n", "patients: []\n", "doctors: [, , ]\n", "\n", "133\n", "patients: []\n", "doctors: [, , ]\n", "\n", "134\n", "patients: []\n", "doctors: [, , ]\n", "\n", "135\n", "patients: []\n", "doctors: [, , ]\n", "\n", "136\n", "patients: []\n", "doctors: [, , ]\n", "\n", "137\n", "patients: []\n", "doctors: [, , ]\n", "\n", "138\n", "patients: []\n", "doctors: [, , ]\n", "\n", "139\n", "patients: []\n", "doctors: [, , ]\n", "\n", "140\n", "patients: []\n", "doctors: [, , ]\n", "\n", "141\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "142\n", "patients: []\n", "doctors: [, , ]\n", "\n", "143\n", "patients: []\n", "doctors: [, , ]\n", "\n", "144\n", "patients: []\n", "doctors: [, , ]\n", "\n", "145\n", "patients: []\n", "doctors: [, , ]\n", "\n", "146\n", "patients: []\n", "doctors: [, , ]\n", "\n", "147\n", "patients: []\n", "doctors: [, , ]\n", "\n", "148\n", "patients: []\n", "doctors: [, , ]\n", "\n", "149\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "150\n", "patients: []\n", "doctors: [, , ]\n", "\n", "151\n", "patients: []\n", "doctors: [, , ]\n", "\n", "152\n", "patients: []\n", "doctors: [, , ]\n", "\n", "153\n", "patients: []\n", "doctors: [, , ]\n", "\n", "154\n", "patients: []\n", "doctors: [, , ]\n", "\n", "155\n", "patients: []\n", "doctors: [, , ]\n", "\n", "156\n", "patients: []\n", "doctors: [, , ]\n", "\n", "157\n", "patients: []\n", "doctors: [, , ]\n", "\n", "158\n", "patients: []\n", "doctors: [, , ]\n", "\n", "159\n", "patients: []\n", "doctors: [, , ]\n", "\n", "160\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "161\n", "patients: []\n", "doctors: [, , ]\n", "\n", "162\n", "patients: []\n", "doctors: [, , ]\n", "\n", "163\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "164\n", "patients: []\n", "doctors: [, , ]\n", "\n", "165\n", "patients: []\n", "doctors: [, , ]\n", "\n", "166\n", "patients: []\n", "doctors: [, , ]\n", "\n", "167\n", "patients: []\n", "doctors: [, , ]\n", "\n", "168\n", "patients: []\n", "doctors: [, , ]\n", "\n", "169\n", "patients: []\n", "doctors: [, , ]\n", "\n", "170\n", "patients: []\n", "doctors: [, , ]\n", "\n", "171\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "172\n", "patients: []\n", "doctors: [, , ]\n", "\n", "173\n", "patients: []\n", "doctors: [, , ]\n", "\n", "174\n", "patients: []\n", "doctors: [, , ]\n", "\n", "175\n", "patients: []\n", "doctors: [, , ]\n", "\n", "176\n", "patients: []\n", "doctors: [, , ]\n", "\n", "177\n", "patients: []\n", "doctors: [, , ]\n", "\n", "178\n", "patients: []\n", "doctors: [, , ]\n", "\n", "179\n", "patients: []\n", "doctors: [, , ]\n", "\n", "180\n", "patients: []\n", "doctors: [, , ]\n", "\n", "181\n", "patients: []\n", "doctors: [, , ]\n", "\n", "182\n", "patients: []\n", "doctors: [, , ]\n", "\n", "183\n", "patients: []\n", "doctors: [, , ]\n", "\n", "184\n", "patients: []\n", "doctors: [, , ]\n", "\n", "185\n", "patients: []\n", "doctors: [, , ]\n", "\n", "186\n", "patients: []\n", "doctors: [, , ]\n", "\n", "187\n", "patients: []\n", "doctors: [, , ]\n", "\n", "188\n", "patients: []\n", "doctors: [, , ]\n", "\n", "189\n", "patients: []\n", "doctors: [, , ]\n", "\n", "190\n", "patients: []\n", "doctors: [, , ]\n", "\n", "191\n", "patients: []\n", "doctors: [, , ]\n", "\n", "192\n", "patients: []\n", "doctors: [, , ]\n", "\n", "193\n", "patients: []\n", "doctors: [, , ]\n", "\n", "194\n", "patients: []\n", "doctors: [, , ]\n", "\n", "195\n", "patients: []\n", "doctors: [, , ]\n", "\n", "196\n", "patients: []\n", "doctors: [, , ]\n", "\n", "197\n", "patients: []\n", "doctors: [, , ]\n", "\n", "198\n", "patients: []\n", "doctors: [, , ]\n", "\n", "199\n", "patients: []\n", "doctors: [, , ]\n", "\n", "200\n", "patients: []\n", "doctors: [, , ]\n", "\n", "201\n", "patients: []\n", "doctors: [, , ]\n", "\n", "202\n", "patients: []\n", "doctors: [, , ]\n", "\n", "203\n", "patients: []\n", "doctors: [, , ]\n", "\n", "204\n", "patients: []\n", "doctors: [, , ]\n", "\n", "205\n", "patients: []\n", "doctors: [, , ]\n", "\n", "206\n", "patients: []\n", "doctors: [, , ]\n", "\n", "207\n", "patients: []\n", "doctors: [, , ]\n", "\n", "208\n", "patients: []\n", "doctors: [, , ]\n", "\n", "209\n", "patients: []\n", "doctors: [, , ]\n", "\n", "210\n", "patients: []\n", "doctors: [, , ]\n", "\n", "211\n", "patients: []\n", "doctors: [, , ]\n", "\n", "212\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "213\n", "patients: []\n", "doctors: [, , ]\n", "\n", "214\n", "patients: []\n", "doctors: [, , ]\n", "\n", "215\n", "patients: []\n", "doctors: [, , ]\n", "\n", "216\n", "patients: []\n", "doctors: [, , ]\n", "\n", "217\n", "patients: []\n", "doctors: [, , ]\n", "\n", "218\n", "patients: []\n", "doctors: [, , ]\n", "\n", "219\n", "patients: []\n", "doctors: [, , ]\n", "\n", "220\n", "patients: []\n", "doctors: [, , ]\n", "\n", "221\n", "patients: []\n", "doctors: [, , ]\n", "\n", "222\n", "patients: []\n", "doctors: [, , ]\n", "\n", "223\n", "patients: []\n", "doctors: [, , ]\n", "\n", "224\n", "patients: []\n", "doctors: [, , ]\n", "\n", "225\n", "patients: []\n", "doctors: [, , ]\n", "\n", "226\n", "patients: []\n", "doctors: [, , ]\n", "\n", "227\n", "patients: []\n", "doctors: [, , ]\n", "\n", "228\n", "patients: []\n", "doctors: [, , ]\n", "\n", "229\n", "patients: []\n", "doctors: [, , ]\n", "\n", "230\n", "patients: []\n", "doctors: [, , ]\n", "\n", "231\n", "patients: []\n", "doctors: [, , ]\n", "\n", "232\n", "patients: []\n", "doctors: [, , ]\n", "\n", "233\n", "patients: []\n", "doctors: [, , ]\n", "\n", "234\n", "patients: []\n", "doctors: [, , ]\n", "\n", "235\n", "patients: []\n", "doctors: [, , ]\n", "\n", "236\n", "patients: []\n", "doctors: [, , ]\n", "\n", "237\n", "patients: []\n", "doctors: [, , ]\n", "\n", "238\n", "patients: []\n", "doctors: [, , ]\n", "\n", "239\n", "patients: []\n", "doctors: [, , ]\n", "\n", "240\n", "patients: []\n", "doctors: [, , ]\n", "\n", "241\n", "patients: []\n", "doctors: [, , ]\n", "\n", "242\n", "patients: []\n", "doctors: [, , ]\n", "\n", "243\n", "patients: []\n", "doctors: [, , ]\n", "\n", "244\n", "patients: []\n", "doctors: [, , ]\n", "\n", "245\n", "patients: []\n", "doctors: [, , ]\n", "\n", "246\n", "patients: []\n", "doctors: [, , ]\n", "\n", "247\n", "patients: []\n", "doctors: [, , ]\n", "\n", "248\n", "patients: []\n", "doctors: [, , ]\n", "\n", "249\n", "patients: []\n", "doctors: [, , ]\n", "\n", "250\n", "patients: []\n", "doctors: [, , ]\n", "\n", "251\n", "patients: []\n", "doctors: [, , ]\n", "\n", "252\n", "patients: []\n", "doctors: [, , ]\n", "\n", "253\n", "patients: []\n", "doctors: [, , ]\n", "\n", "254\n", "patients: []\n", "doctors: [, , ]\n", "\n", "255\n", "patients: []\n", "doctors: [, , ]\n", "\n", "256\n", "patients: []\n", "doctors: [, , ]\n", "\n", "257\n", "patients: []\n", "doctors: [, , ]\n", "\n", "258\n", "patients: []\n", "doctors: [, , ]\n", "\n", "259\n", "patients: []\n", "doctors: [, , ]\n", "\n", "260\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "261\n", "patients: []\n", "doctors: [, , ]\n", "\n", "262\n", "patients: []\n", "doctors: [, , ]\n", "\n", "263\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "264\n", "patients: []\n", "doctors: [, , ]\n", "\n", "265\n", "patients: []\n", "doctors: [, , ]\n", "\n", "266\n", "patients: []\n", "doctors: [, , ]\n", "\n", "267\n", "patients: []\n", "doctors: [, , ]\n", "\n", "268\n", "patients: []\n", "doctors: [, , ]\n", "\n", "269\n", "patients: []\n", "doctors: [, , ]\n", "\n", "270\n", "patients: []\n", "doctors: [, , ]\n", "\n", "271\n", "patients: []\n", "doctors: [, , ]\n", "\n", "272\n", "patients: []\n", "doctors: [, , ]\n", "\n", "273\n", "patients: []\n", "doctors: [, , ]\n", "\n", "274\n", "patients: []\n", "doctors: [, , ]\n", "\n", "275\n", "patients: []\n", "doctors: [, , ]\n", "\n", "276\n", "patients: []\n", "doctors: [, , ]\n", "\n", "277\n", "patients: []\n", "doctors: [, , ]\n", "\n", "278\n", "patients: []\n", "doctors: [, , ]\n", "\n", "279\n", "patients: []\n", "doctors: [, , ]\n", "\n", "280\n", "patients: []\n", "doctors: [, , ]\n", "\n", "281\n", "patients: []\n", "doctors: [, , ]\n", "\n", "282\n", "patients: []\n", "doctors: [, , ]\n", "\n", "283\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "284\n", "patients: []\n", "doctors: [, , ]\n", "\n", "285\n", "patients: []\n", "doctors: [, , ]\n", "\n", "286\n", "patients: []\n", "doctors: [, , ]\n", "\n", "287\n", "patients: []\n", "doctors: [, , ]\n", "\n", "288\n", "patients: []\n", "doctors: [, , ]\n", "\n", "289\n", "patients: []\n", "doctors: [, , ]\n", "\n", "290\n", "patients: []\n", "doctors: [, , ]\n", "\n", "291\n", "patients: []\n", "doctors: [, , ]\n", "\n", "292\n", "patients: []\n", "doctors: [, , ]\n", "\n", "293\n", "patients: []\n", "doctors: [, , ]\n", "\n", "294\n", "patients: []\n", "doctors: [, , ]\n", "\n", "295\n", "patients: []\n", "doctors: [, , ]\n", "\n", "296\n", "patients: []\n", "doctors: [, , ]\n", "\n", "297\n", "patients: []\n", "doctors: [, , ]\n", "\n", "298\n", "patients: []\n", "doctors: [, , ]\n", "\n", "299\n", "patients: []\n", "doctors: [, , ]\n", "\n", "300\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "301\n", "patients: []\n", "doctors: [, , ]\n", "\n", "302\n", "patients: []\n", "doctors: [, , ]\n", "\n", "303\n", "patients: []\n", "doctors: [, , ]\n", "\n", "304\n", "patients: []\n", "doctors: [, , ]\n", "\n", "305\n", "patients: []\n", "doctors: [, , ]\n", "\n", "306\n", "patients: []\n", "doctors: [, , ]\n", "\n", "307\n", "patients: []\n", "doctors: [, , ]\n", "\n", "308\n", "patients: []\n", "doctors: [, , ]\n", "\n", "309\n", "patients: []\n", "doctors: [, , ]\n", "\n", "310\n", "patients: []\n", "doctors: [, , ]\n", "\n", "311\n", "patients: []\n", "doctors: [, , ]\n", "\n", "312\n", "patients: []\n", "doctors: [, , ]\n", "\n", "313\n", "patients: []\n", "doctors: [, , ]\n", "\n", "314\n", "patients: []\n", "doctors: [, , ]\n", "\n", "315\n", "patients: []\n", "doctors: [, , ]\n", "\n", "316\n", "patients: []\n", "doctors: [, , ]\n", "\n", "317\n", "patients: []\n", "doctors: [, , ]\n", "\n", "318\n", "patients: []\n", "doctors: [, , ]\n", "\n", "319\n", "patients: []\n", "doctors: [, , ]\n", "\n", "320\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "321\n", "patients: []\n", "doctors: [, , ]\n", "\n", "322\n", "patients: []\n", "doctors: [, , ]\n", "\n", "323\n", "patients: []\n", "doctors: [, , ]\n", "\n", "324\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "325\n", "patients: []\n", "doctors: [, , ]\n", "\n", "326\n", "patients: []\n", "doctors: [, , ]\n", "\n", "327\n", "patients: []\n", "doctors: [, , ]\n", "\n", "328\n", "patients: []\n", "doctors: [, , ]\n", "\n", "329\n", "patients: []\n", "doctors: [, , ]\n", "\n", "330\n", "patients: []\n", "doctors: [, , ]\n", "\n", "331\n", "patients: []\n", "doctors: [, , ]\n", "\n", "332\n", "patients: []\n", "doctors: [, , ]\n", "\n", "333\n", "patients: []\n", "doctors: [, , ]\n", "\n", "334\n", "patients: []\n", "doctors: [, , ]\n", "\n", "335\n", "patients: []\n", "doctors: [, , ]\n", "\n", "336\n", "patients: []\n", "doctors: [, , ]\n", "\n", "337\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "338\n", "patients: []\n", "doctors: [, , ]\n", "\n", "339\n", "patients: []\n", "doctors: [, , ]\n", "\n", "340\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "341\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "342\n", "patients: []\n", "doctors: [, , ]\n", "\n", "343\n", "patients: []\n", "doctors: [, , ]\n", "\n", "344\n", "patients: []\n", "doctors: [, , ]\n", "\n", "345\n", "patients: []\n", "doctors: [, , ]\n", "\n", "346\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "347\n", "patients: []\n", "doctors: [, , ]\n", "\n", "348\n", "patients: []\n", "doctors: [, , ]\n", "\n", "349\n", "patients: []\n", "doctors: [, , ]\n", "\n", "350\n", "patients: []\n", "doctors: [, , ]\n", "\n", "351\n", "patients: []\n", "doctors: [, , ]\n", "\n", "352\n", "patients: []\n", "doctors: [, , ]\n", "\n", "353\n", "patients: []\n", "doctors: [, , ]\n", "\n", "354\n", "patients: []\n", "doctors: [, , ]\n", "\n", "355\n", "patients: []\n", "doctors: [, , ]\n", "\n", "356\n", "patients: []\n", "doctors: [, , ]\n", "\n", "357\n", "patients: []\n", "doctors: [, , ]\n", "\n", "358\n", "patients: []\n", "doctors: [, , ]\n", "\n", "359\n", "patients: []\n", "doctors: [, , ]\n", "\n", "360\n", "patients: []\n", "doctors: [, , ]\n", "\n", "361\n", "patients: []\n", "doctors: [, , ]\n", "\n", "362\n", "patients: []\n", "doctors: [, , ]\n", "\n", "363\n", "patients: []\n", "doctors: [, , ]\n", "\n", "364\n", "patients: []\n", "doctors: [, , ]\n", "\n", "365\n", "patients: []\n", "doctors: [, , ]\n", "\n", "366\n", "patients: []\n", "doctors: [, , ]\n", "\n", "367\n", "patients: []\n", "doctors: [, , ]\n", "\n", "368\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "369\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "370\n", "patients: []\n", "doctors: [, , ]\n", "\n", "371\n", "patients: []\n", "doctors: [, , ]\n", "\n", "372\n", "patients: []\n", "doctors: [, , ]\n", "\n", "373\n", "patients: []\n", "doctors: [, , ]\n", "\n", "374\n", "patients: []\n", "doctors: [, , ]\n", "\n", "375\n", "patients: []\n", "doctors: [, , ]\n", "\n", "376\n", "patients: []\n", "doctors: [, , ]\n", "\n", "377\n", "patients: []\n", "doctors: [, , ]\n", "\n", "378\n", "patients: []\n", "doctors: [, , ]\n", "\n", "379\n", "patients: []\n", "doctors: [, , ]\n", "\n", "380\n", "patients: []\n", "doctors: [, , ]\n", "\n", "381\n", "patients: []\n", "doctors: [, , ]\n", "\n", "382\n", "patients: []\n", "doctors: [, , ]\n", "\n", "383\n", "patients: []\n", "doctors: [, , ]\n", "\n", "384\n", "patients: []\n", "doctors: [, , ]\n", "\n", "385\n", "patients: []\n", "doctors: [, , ]\n", "\n", "386\n", "patients: []\n", "doctors: [, , ]\n", "\n", "387\n", "patients: []\n", "doctors: [, , ]\n", "\n", "388\n", "patients: []\n", "doctors: [, , ]\n", "\n", "389\n", "patients: []\n", "doctors: [, , ]\n", "\n", "390\n", "patients: []\n", "doctors: [, , ]\n", "\n", "391\n", "patients: []\n", "doctors: [, , ]\n", "\n", "392\n", "patients: []\n", "doctors: [, , ]\n", "\n", "393\n", "patients: []\n", "doctors: [, , ]\n", "\n", "394\n", "patients: []\n", "doctors: [, , ]\n", "\n", "395\n", "patients: []\n", "doctors: [, , ]\n", "\n", "396\n", "patients: []\n", "doctors: [, , ]\n", "\n", "397\n", "patients: []\n", "doctors: [, , ]\n", "\n", "398\n", "patients: []\n", "doctors: [, , ]\n", "\n", "399\n", "patients: []\n", "doctors: [, , ]\n", "\n", "400\n", "patients: []\n", "doctors: [, , ]\n", "\n", "401\n", "patients: []\n", "doctors: [, , ]\n", "\n", "402\n", "patients: []\n", "doctors: [, , ]\n", "\n", "403\n", "patients: []\n", "doctors: [, , ]\n", "\n", "404\n", "patients: []\n", "doctors: [, , ]\n", "\n", "405\n", "patients: []\n", "doctors: [, , ]\n", "\n", "406\n", "patients: []\n", "doctors: [, , ]\n", "\n", "407\n", "patients: []\n", "doctors: [, , ]\n", "\n", "408\n", "patients: []\n", "doctors: [, , ]\n", "\n", "409\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "410\n", "patients: []\n", "doctors: [, , ]\n", "\n", "411\n", "patients: []\n", "doctors: [, , ]\n", "\n", "412\n", "patients: []\n", "doctors: [, , ]\n", "\n", "413\n", "patients: []\n", "doctors: [, , ]\n", "\n", "414\n", "patients: []\n", "doctors: [, , ]\n", "\n", "415\n", "patients: []\n", "doctors: [, , ]\n", "\n", "416\n", "patients: []\n", "doctors: [, , ]\n", "\n", "417\n", "patients: []\n", "doctors: [, , ]\n", "\n", "418\n", "patients: []\n", "doctors: [, , ]\n", "\n", "419\n", "patient arrived!\n", "patients: []\n", "doctors: [, , ]\n", "\n", "420\n", "patients: []\n", "doctors: [, , ]\n", "\n", "421\n", "patients: []\n", "doctors: [, , ]\n", "\n", "422\n", "patients: []\n", "doctors: [, , ]\n", "\n", "423\n", "patients: []\n", "doctors: [, , ]\n", "\n", "424\n", "patients: []\n", "doctors: [, , ]\n", "\n", "425\n", "patients: []\n", "doctors: [, , ]\n", "\n", "426\n", "patients: []\n", "doctors: [, , ]\n", "\n", "427\n", "patients: []\n", "doctors: [, , ]\n", "\n", "428\n", "patients: []\n", "doctors: [, , ]\n", "\n", "429\n", "patients: []\n", "doctors: [, , ]\n", "\n", "430\n", "patients: []\n", "doctors: [, , ]\n", "\n", "431\n", "patients: []\n", "doctors: [, , ]\n", "\n", "432\n", "patients: []\n", "doctors: [, , ]\n", "\n", "433\n", "patients: []\n", "doctors: [, , ]\n", "\n", "434\n", "patients: []\n", "doctors: [, , ]\n", "\n", "435\n", "patients: []\n", "doctors: [, , ]\n", "\n", "436\n", "patients: []\n", "doctors: [, , ]\n", "\n", "437\n", "patients: []\n", "doctors: [, , ]\n", "\n", "438\n", "patients: []\n", "doctors: [, , ]\n", "\n", "439\n", "patients: []\n", "doctors: [, , ]\n", "\n", "440\n", "patients: []\n", "doctors: [, , ]\n", "\n", "441\n", "patients: []\n", "doctors: [, , ]\n", "\n", "442\n", "patients: []\n", "doctors: [, , ]\n", "\n", "443\n", "patients: []\n", "doctors: [, , ]\n", "\n", "444\n", "patients: []\n", "doctors: [, , ]\n", "\n" ] }, { "data": { "text/plain": [ "{'avg_wait': 5.666666666666667,\n", " 'longest_wait': 15,\n", " 'mins_past_close': 25,\n", " 'n_patients': 30,\n", " 'n_waiters': 3,\n", " 'wait_times': [0,\n", " 0,\n", " 0,\n", " 0,\n", " 1,\n", " 0,\n", " 0,\n", " 0,\n", " 0,\n", " 0,\n", " 0,\n", " 0,\n", " 0,\n", " 0,\n", " 0,\n", " 0,\n", " 0,\n", " 0,\n", " 0,\n", " 0,\n", " 0,\n", " 0,\n", " 0,\n", " 0,\n", " 1,\n", " 15,\n", " 0,\n", " 0,\n", " 0,\n", " 0]}" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "run_simulation(verbose=True)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "data = [run_simulation() for _ in range(2000)]\n", "\n", "patient_counts = [x['n_patients'] for x in data]\n", "avg_waits = [x.get('avg_wait', 0) for x in data]\n", "longest_waits = [x['longest_wait'] for x in data]\n", "waiter_counts = [x['n_waiters'] for x in data]\n", "overtimes = [x['mins_past_close'] for x in data if x['mins_past_close'] > 0]" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 36. 40. 44.]\n", "[ 3.5 5. 6.5]\n", "[ 6. 9. 13.]\n", "[ 3. 6. 10.]\n", "[ 7. 12. 17.25]\n" ] } ], "source": [ "def get_percentiles(data):\n", " return np.percentile(data, [25, 50, 75])\n", "\n", "for stat in [patient_counts, avg_waits, longest_waits,\n", " waiter_counts, overtimes]:\n", " print(get_percentiles(stat))" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/svg+xml": [ "\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", " \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", " \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", " \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", " \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", " \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", " \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", " \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", " \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", " \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", " \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", " \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", " \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", " \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", " \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", " \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" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.rcParams['figure.figsize'] = 7, 9\n", "plt.subplot(3, 2, 1)\n", "plt.hist(patient_counts, color='black')\n", "plt.xlabel('number of patients')\n", "plt.title('Patient Counts')\n", "\n", "plt.subplot(3, 2, 2)\n", "plt.hist(avg_waits, color='black')\n", "plt.xlabel('minutes')\n", "plt.title('Average Waits')\n", "\n", "plt.subplot(3, 2, 3)\n", "plt.hist(longest_waits, color='black')\n", "plt.xlabel('minutes')\n", "plt.title('Longest Waits')\n", "\n", "plt.subplot(3, 2, 4)\n", "plt.hist(waiter_counts, color='black')\n", "plt.xlabel('number of waiters')\n", "plt.title('Number of people who waited')\n", "\n", "plt.subplot(3, 2, 5)\n", "plt.hist(overtimes, color='black')\n", "plt.xlabel('minutes')\n", "plt.title('Number of minutes Overtime')\n", "\n", "plt.tight_layout();\n", "# plt.savefig('stats_2.png', format='png', dpi=400)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Posterior inference and predictive distributions.\n", "1. Suppose you have a Beta(3, 3) prior distribution on the probability θ that a coin will yield ‘heads’ when spun in a specified manner. The coin is independently spun ten times, and ‘heads’ appear fewer than three times. You are not told how may heads were seen, only that the number is less than 3. Calculate your exact posterior density (up to a proportionality constant) for θ and sketch it.\n", "


\n", "2. Consider two coins, C1 and C2, with the following characteristics: Pr(heads|C1) = 0.4 and Pr(heads|C2) = 0.6. Choose one of the coins at random and imagine spinning it repeatedly.\n", " 1. Given that the first two spins from the chosen coin are tails, what is the expected number of additional spins until a heads shows up?\n", " 2. How about if the first three spins are all tails?" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": true }, "outputs": [], "source": [ "x = np.linspace(0, 1, 100)\n", "prior = beta.pdf(x, 3, 3)\n", "likelihood = binom.cdf(3, 10, x)\n", "posterior = likelihood * prior" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/svg+xml": [ "\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", " \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", " \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" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "lw = 4\n", "color = 'black'\n", "plt.rcParams['figure.figsize'] = 8, 6\n", "plt.style.use('ggplot')\n", "\n", "plt.subplot(2, 1, 1)\n", "plt.plot(x, prior, linewidth=lw, color=color)\n", "plt.title('Prior')\n", "plt.xlabel(r'$\\theta$', labelpad=.01)\n", "plt.ylabel(r'$p(\\theta)$')\n", "plt.yticks([])\n", "plt.grid(b=False)\n", "\n", "plt.subplot(2, 1, 2)\n", "plt.plot(x, posterior, linewidth=lw, color=color)\n", "plt.title('Posterior')\n", "plt.xlabel(r'$\\theta$')\n", "plt.ylabel(r'$p(\\theta|y)$')\n", "plt.yticks([])\n", "plt.grid(b=False)\n", "\n", "plt.tight_layout();\n", "plt.savefig('coin_flip_densities.png', format='png', dpi=400)" ] }, { "cell_type": "code", "execution_count": 14, "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", "
accidentsexposures
year
1976243864
1977254300
1978315026
1979315482
1980225815
1981216034
1982265876
1983206224
1984167434
1985227106
\n", "
" ], "text/plain": [ " accidents exposures\n", "year \n", "1976 24 3864\n", "1977 25 4300\n", "1978 31 5026\n", "1979 31 5482\n", "1980 22 5815\n", "1981 21 6034\n", "1982 26 5876\n", "1983 20 6224\n", "1984 16 7434\n", "1985 22 7106" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df = pd.read_csv('data/accident_data.csv', index_col='year')\n", "df" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/svg+xml": [ "\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", " \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", " \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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# check if the accident rate is increasing or decreasing\n", "rates = df['accidents'] / df['exposures']\n", "x = df.index.tolist()\n", "\n", "# plot it\n", "plt.scatter(x, rates, s=70, color='black')\n", "plt.xlim([1975, 1986])\n", "plt.ylim([0, 0.01])\n", "plt.xlabel('Year')\n", "plt.ylabel('Accident Rate')\n", "plt.title('Accident Rates by Year')\n", "plt.show();" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# Multilevel Modeling\n", "\n", "The following table gives the number of fatal accidents on scheduled airline flights, and exposure (in 100s of millions of passenger miles), per year, over a ten-year period.\n", "\n", "- Assume that the numbers of fatal accidents in each year are independent with a Poisson(θ) distribution. Assign a $Gamma(\\alpha,\\beta)$ prior for $\\theta$, selecting reasonable values for $\\alpha$ and $\\beta$.\n", " Using simulation, compute a 95% predictive interval for the number of fatal acci- dents in 1986.\n", " \n", "- For t = 1,...,10, let yt denote the number of fatal accidents in year 1975 + t. Assume the yt are independent $Poisson(xt\\theta)$, where xt denotes the exposure (in 100s of millions of passenger miles) in year 1975 + t. Assign a $Gamma(\\alpha,\\beta)$ prior for $\\theta$; you will probably want different hyperparameter values than in part (1).\n", " Give a 95% predictive interval for the number of fatal accidents in 1986 under the assumption that 7.6 × 1011 passenger miles are flown that year.\n", "\n", "\n", "$\\theta$ ~ $Gamma(\\alpha, \\beta)$\n", "\n", "$y_{i}$ ~ $Poisson(\\theta)$" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [], "source": [ "class BayesianPoissonModel:\n", " '''Spits out data simulated values from a bayesian model'''\n", " \n", " def __init__(self, initial_alpha, initial_beta, n_exposures=1):\n", " self.alpha = initial_alpha\n", " self.beta = initial_beta\n", " self.n_exposures = n_exposures\n", " \n", " \n", " def update(self, data):\n", " '''\n", " updates p(theta) to be p(theta|data)\n", " '''\n", " for y, n in data:\n", " self.alpha += y\n", " self.beta += n\n", " \n", " def simulate(self, n):\n", " theta = np.random.gamma(self.alpha, self.beta**-1, n)\n", " return np.random.poisson(theta * self.n_exposures)\n", " \n", " def __repr__(self):\n", " return ''.format(self.alpha, self.beta)\n", " \n", " def __str__(self):\n", " return self.__repr__()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll disregard number of exposures per year and only model the number of accidents\n", "\n", "$\\beta$ will be 1, as in a period of 1 year" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "95% predictive interval: [14, 34]\n" ] }, { "data": { "image/svg+xml": [ "\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", " \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", " \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", " \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", " \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" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "accidents_data = np.c_[df['accidents'].values, np.ones(df.shape[0])]\n", "\n", "model = BayesianPoissonModel(df['accidents'].mean(), 1)\n", "model.update(accidents_data)\n", "simulated = model.simulate(n=100000)\n", "\n", "print(model)\n", "print('95% predictive interval: [{}, {}]'.format(*np.percentile(simulated, (2.5, 97.5)).astype(int)))\n", "\n", "# plot it\n", "plt.rcParams['figure.figsize'] = 7, 4\n", "plt.hist(simulated, bins=20, color='black')\n", "plt.xlabel('Number of Accidents')\n", "plt.ylabel('Count')\n", "plt.title('Simulated Number of Fatal Accidents')\n", "plt.show();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Incorporate number of exposures into the models to \"weight\" each accident count appropriately" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "95% predictive interval: [21, 44]\n" ] }, { "data": { "image/svg+xml": [ "\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", " \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", " \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", " \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" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "n_exposures = 7600\n", "# use the means as a prior\n", "a = df['accidents'].mean()\n", "b = df['exposures'].mean() \n", "\n", "model_full = BayesianPoissonModel(a, b, n_exposures=n_exposures)\n", "model_full.update(df[['accidents', 'exposures']].values)\n", "simulated = model_full.simulate(n=100000)\n", "\n", "print(model_full)\n", "print('95% predictive interval: [{}, {}]'.format(*np.percentile(simulated, (2.5, 97.5)).astype(int)))\n", "\n", "# plot it\n", "plt.rcParams['figure.figsize'] = 7, 4\n", "plt.hist(simulated, bins=20, color='black')\n", "plt.xlabel('Number of Accidents')\n", "plt.ylabel('count')\n", "plt.title('Simulated Number of Fatal Accidents')\n", "plt.show;" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# Simulation\n", "\n", "An experiment was performed to estimate the effect of beta blockers on mortality of cardiac patients. A group of patients were randomly assigned to treatment and control groups: out of 282 patients receiving the control, 37 died; out of 278 receiving the treatment, 25 died. Assume that the outcomes are independent and binomially distributed, with probabilities of death $p_0$ and $p_1$ under the control and treatment, respectively.\n", "\n", "Set up a noninformative prior distribution on ($p_0$, $p_1$) and obtain posterior simulations. Summarize the posterior distribution for the odds ratio,\n", "\n", "\n", "
$\\rho$ = $\\cfrac{\\frac{p_1}{1 - p_1}}{\\frac{p_0}{1 - p_0}}$
" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "95% predictive interval: [0.403, 1.113]\n" ] }, { "data": { "image/svg+xml": [ "\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", " \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", " \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", " \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" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# draw 1000 p_0's and p_1's from their respective posteriors\n", "n_samples = 1000\n", "p_0 = np.random.beta(37 + 1, 282 - 37 + 1, n_samples) # control\n", "p_1 = np.random.beta(25 + 1, 278 - 25 + 1, n_samples) # treatment\n", "\n", "# compute each simulated odds ratio\n", "odds_samples = (p_1 / (1 - p_1)) / (p_0 / (1 - p_0))\n", "\n", "lower, upper = np.percentile(odds_samples, [2.5, 97.5])\n", "print('95% predictive interval: [{:.3f}, {:.3f}]'.format(lower, upper))\n", "\n", "# plot it\n", "plt.hist(odds_samples, bins=15, color='black')\n", "plt.xlabel('Odds')\n", "plt.ylabel('frequency')\n", "plt.title('Simulated odds ratios of $p_1$ and $p_0$', y=1.03)\n", "# plt.savefig('p_7_plot.png', format='png', dpi=400)\n", "plt.show()" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python [conda env:py3k]", "language": "python", "name": "conda-env-py3k-py" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.5.2" } }, "nbformat": 4, "nbformat_minor": 1 }