test_scenario_sim.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bayesian interpretation of medical tests\n",
    "-----------------------------------------\n",
    "\n",
    "This notebooks explores several problems related to interpreting the results of medical tests.\n",
    "\n",
    "Copyright 2016 Allen Downey\n",
    "\n",
    "MIT License: http://opensource.org/licenses/MIT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "\n",
    "from thinkbayes2 import Pmf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from random import random\n",
    "\n",
    "def flip(p):\n",
    "    return random() < p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def run_single_simulation(func, iters=1000000):\n",
    "    pmf_t = Pmf([0.2, 0.4])\n",
    "    p = 0.1\n",
    "    s = 0.9\n",
    "\n",
    "    outcomes = Pmf()\n",
    "    post_t = Pmf()\n",
    "    for i in range(iters):\n",
    "        test, sick, t = func(p, s, pmf_t)\n",
    "        if test:\n",
    "            outcomes[sick] += 1\n",
    "            post_t[t] += 1\n",
    "\n",
    "    outcomes.Normalize()\n",
    "    post_t.Normalize()\n",
    "    return outcomes, post_t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scenario A:  Choose `t` for each patient, yield all patients regardless of test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False 0.75061282845\n",
      "True 0.24938717155\n",
      "0.2 0.375486862013\n",
      "0.4 0.624513137987\n"
     ]
    }
   ],
   "source": [
    "def generate_patient_A(p, s, pmf_t):\n",
    "    while True:\n",
    "        t = pmf_t.Random()\n",
    "        sick = flip(p)\n",
    "        test = flip(s) if sick else flip(t)\n",
    "        return test, sick, t\n",
    "                \n",
    "outcomes, post_t = run_single_simulation(generate_patient_A)\n",
    "outcomes.Print()\n",
    "post_t.Print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scenario B:  Choose `t` before generating patients, yield all patients regardless of test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False 0.751086814612\n",
      "True 0.248913185388\n",
      "0.2 0.375297571924\n",
      "0.4 0.624702428076\n"
     ]
    }
   ],
   "source": [
    "def generate_patient_B(p, s, pmf_t):\n",
    "    t = pmf_t.Random()\n",
    "    while True:\n",
    "        sick = flip(p)\n",
    "        test = flip(s) if sick else flip(t)\n",
    "        return test, sick, t\n",
    "                \n",
    "outcomes, post_t = run_single_simulation(generate_patient_B)\n",
    "outcomes.Print()\n",
    "post_t.Print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scenario C:  Choose `t` for each patient, only yield patients who test positive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False 0.750118\n",
      "True 0.249882\n",
      "0.2 0.374743\n",
      "0.4 0.625257\n"
     ]
    }
   ],
   "source": [
    "def generate_patient_C(p, s, pmf_t):\n",
    "    while True:\n",
    "        t = pmf_t.Random()\n",
    "        sick = flip(p)\n",
    "        test = flip(s) if sick else flip(t)\n",
    "        if test:\n",
    "            return test, sick, t\n",
    "                \n",
    "outcomes, post_t = run_single_simulation(generate_patient_C)\n",
    "outcomes.Print()\n",
    "post_t.Print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scenario D:  Choose `t` before generating patients, only yield patients who test positive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False 0.73333\n",
      "True 0.26667\n",
      "0.2 0.498557\n",
      "0.4 0.501443\n"
     ]
    }
   ],
   "source": [
    "def generate_patient_D(p, s, pmf_t):\n",
    "    t = pmf_t.Random()\n",
    "    while True:\n",
    "        sick = flip(p)\n",
    "        test = flip(s) if sick else flip(t)\n",
    "        if test:\n",
    "            return test, sick, t\n",
    "                \n",
    "outcomes, post_t = run_single_simulation(generate_patient_D)\n",
    "outcomes.Print()\n",
    "post_t.Print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a variation of the Scenario D where we only consider cases where patient 1 is the first to test positive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dist of t\n",
      "0.2 0.502046233824\n",
      "0.4 0.497953766176\n",
      "Dist of status\n",
      "False 0.733476787564\n",
      "True 0.266523212436\n"
     ]
    }
   ],
   "source": [
    "from random import choice\n",
    "import numpy as np \n",
    "N = 100\n",
    "patients = range(N)\n",
    "\n",
    "p = 0.1\n",
    "s = 0.9\n",
    "num_sick = 0\n",
    "\n",
    "pmf_t = Pmf()\n",
    "pmf_sick = Pmf()\n",
    "\n",
    "for i in range(10000000):\n",
    "    # decide what the value of t is\n",
    "    t = choice([0.2, 0.4])\n",
    "    np.random.shuffle(patients)\n",
    "    \n",
    "    # generate patients until we get a positive test\n",
    "    for patient in patients:\n",
    "        sick = flip(p)\n",
    "        test = flip(s) if sick else flip(t)\n",
    "        if test:\n",
    "            if patient==1:\n",
    "                #print(patient, sick, t)\n",
    "                pmf_t[t] += 1\n",
    "                pmf_sick[sick] += 1\n",
    "            break\n",
    "            \n",
    "pmf_t.Normalize()\n",
    "pmf_sick.Normalize()\n",
    "\n",
    "print('Dist of t')\n",
    "pmf_t.Print()\n",
    "print('Dist of status')\n",
    "pmf_sick.Print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num_sick"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}