lectures/Lecture_2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Introduction to Probability & Statistics: I\n",
    "\n",
    "*S. R. Taylor (2021)*\n",
    "\n",
    "Material in this lecture and notebook is based upon the Basic Stats portion of G. Richards' \"Astrostatistics\" class at Drexel University (PHYS 440/540, https://github.com/gtrichards/PHYS_440_540), the Introduction to Probability & Statistics portion of A. Connolly's & Ž. Ivezić's \"Astrostatistics & Machine Learning\" class at the University of Washington (ASTR 598, https://github.com/dirac-institute/uw-astr598-w18), and J. Bovy's mini-course on \"Statistics & Inference in Astrophysics\" at the University of Toronto (http://astro.utoronto.ca/~bovy/teaching.html). \n",
    "\n",
    "##### Reading:\n",
    "\n",
    "- [Textbook](http://press.princeton.edu/titles/10159.html) Chapter 3. \n",
    "- [David Hogg: \"Data analysis recipes: Probability calculus for inference\"](https://arxiv.org/abs/1205.4446)\n",
    "\n",
    "***Exercises required for class participation are in <font color='red'>red</font>.***\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "* [Preliminaries and notation](#one)\n",
    "* [Probability](#two)\n",
    "* [Bayes' Theorem](#three)\n",
    "* [Bayes' Theorem Example: Legos](#four)\n",
    "* [Bayes' Theorem Example: Monty Hall Problem (or \"Deal Or No Deal\")](#five)\n",
    "* [Bayes' Theorem Example: Contingency Table](#six)\n",
    "* [Transformations of random variables](#seven)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Preliminaries and notation <a class=\"anchor\" id=\"one\"></a>\n",
    "\n",
    "* **\"Astrostatistics\"** = extracting knowledge from astronomical data.\n",
    "* **\"Knowledge\"** = summary (physical or phenomenological) of data behavior.\n",
    "* **\"Data\"** = result of measurements.\n",
    "\n",
    "In the textbook, $x$ is a scalar quantity that is measured $N$ times to form a dataset.\n",
    "\n",
    "* $x_i$ is a single measurement with $i=1,...,N$.\n",
    "* $\\{x_i\\}$ refers to the set of all N measurements comprising the dataset. \n",
    "\n",
    "Our data can be real numbers, discrete labels (strings or numbers), or even \"missing values\" (we sometimes pad our datasets with NaNs in this case). \n",
    "\n",
    "**Goal of data mining & statistical inference:**\n",
    "> We are generally trying to *estimate* $h(x)$, the *true* generating distribution from which $\\{x_i\\}$ are drawn. \n",
    "\n",
    "* $h(x)$ is the **probability density (distribution) function** or the **\"pdf\"** and $h(x)dx$ is the propobability of a value lying between $x$ and $x+dx$. This distribution can have several levels-- the population distribution of events (e.g. source redshifts), and a measurement uncertainity distribution that blurs our measured data away from true values.\n",
    "\n",
    "* The \"left to right\" integral of $h(x)$ is the **cumulative distribution function** (**\"cdf\"**), $H(x) = \\int_{-\\infty}^x h(x')dx'$. The inverse function of the cdf is the **quantile function**, e.g. what $x$ value has 90% of the distribution below it?\n",
    "\n",
    "* While $h(x)$ is the \"true\" pdf (or **population** pdf).  What we *measure* from the data is the **empirical** pdf, which is denoted $f(x)$.  So, $f(x)$ is a *model* of $h(x)$.  In principle, with infinite data $f(x) \\rightarrow h(x)$, but in reality the blurring effect of measurement errors keep this from being strictly true. Likewise, the empirical cdf is denoted $F(x)$.\n",
    "\n",
    "* If we are attempting to guess a physical *model* for $h(x)$, then the process is ***parametric***.  With a model solution we can generate new data that should mimic what we measure.  If we are not attempting to guess a model, then the process is ***non-parametric***, i.e. we are just trying to describe the data behavior in a compact practical way."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>**Execute this code cell**</font> (don't worry about warnings or errors)\n",
    "\n",
    "This cell will start with a generating distribution $h(x)$, draw a number of random samples as data $\\{x_i\\}$, and then fit these data with a parametric and non-parametric model $f(x)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# based on Ivezic v2, Figure 6.8; edited by G. T. Richards and S. R. Taylor\n",
    "\n",
    "# Author: Jake VanderPlas\n",
    "# License: BSD\n",
    "#   The figure produced by this code is published in the textbook\n",
    "#   \"Statistics, Data Mining, and Machine Learning in Astronomy\" (2013)\n",
    "#   For more information, see http://astroML.github.com\n",
    "#   To report a bug or issue, use the following forum:\n",
    "#    https://groups.google.com/forum/#!forum/astroml-general\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy import stats\n",
    "\n",
    "from sklearn.neighbors import KernelDensity\n",
    "\n",
    "from astropy.visualization import hist\n",
    "from sklearn.mixture import GaussianMixture\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# This function adjusts matplotlib settings for a uniform feel in the textbook.\n",
    "# Note that with usetex=True, fonts are rendered with LaTeX.  This may\n",
    "# result in an error if LaTeX is not installed on your system.  In that case,\n",
    "# you can set usetex to False.\n",
    "if \"setup_text_plots\" not in globals():\n",
    "    from astroML.plotting import setup_text_plots\n",
    "setup_text_plots(fontsize=14, usetex=True)\n",
    "%config InlineBackend.figure_format='retina' # very useful command for high-res images\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Generate our data: a mix of several Cauchy distributions\n",
    "random_state = np.random.RandomState(seed=0)\n",
    "N = 2000 # number of data points\n",
    "mu_gamma_f = [(5, 1.0, 0.1),\n",
    "              (7, 0.5, 0.5),\n",
    "              (9, 0.1, 0.1),\n",
    "              (12, 0.5, 0.2),\n",
    "              (14, 1.0, 0.1)]\n",
    "hx = lambda x: sum([f * stats.cauchy(mu, gamma).pdf(x)\n",
    "                          for (mu, gamma, f) in mu_gamma_f])\n",
    "x = np.concatenate([stats.cauchy(mu, gamma).rvs(int(f * N), random_state=random_state)\n",
    "                    for (mu, gamma, f) in mu_gamma_f])\n",
    "random_state.shuffle(x)\n",
    "x = x[x > -10]\n",
    "x = x[x < 30]\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# plot the results\n",
    "fig,ax = plt.subplots(figsize=(10, 10))\n",
    "\n",
    "xgrid = np.linspace(-10, 30, 1000)\n",
    "\n",
    "# Compute density with KDE\n",
    "kde = KernelDensity(bandwidth=0.1, kernel='gaussian')\n",
    "kde.fit(x[:, None])\n",
    "dens_kde = np.exp(kde.score_samples(xgrid[:, None]))\n",
    "\n",
    "# Compute density via Gaussian Mixtures using a pre-defined number of clusters (13)\n",
    "gmm = GaussianMixture(n_components=13).fit(x.reshape(-1, 1))\n",
    "logprob = gmm.score_samples(xgrid.reshape(-1, 1))\n",
    "fx = lambda j : np.exp(gmm.score_samples(j.reshape(-1, 1)))\n",
    "\n",
    "# plot the results\n",
    "ax.plot(xgrid, hx(xgrid), ':', color='black', zorder=3,\n",
    "            label=\"$h(x)$, Generating Distribution\")\n",
    "ax.plot(xgrid, fx(np.array(xgrid)), '-', color='gray',\n",
    "            label=\"$f(x)$, parametric (13 Gaussians)\")\n",
    "ax.plot(xgrid, dens_kde, '-', color='black', zorder=3,\n",
    "            label=\"$f(x)$, non-parametric (KDE)\")\n",
    "\n",
    "# label the plot\n",
    "ax.text(0.02, 0.95, \"%i points\" % N, ha='left', va='top',\n",
    "            transform=ax.transAxes)\n",
    "ax.set_ylabel('$p(x)$',fontsize=14)\n",
    "ax.legend(loc='upper right')\n",
    "\n",
    "ax.set_xlabel('$x$',fontsize=14)\n",
    "ax.set_xlim(0, 20)\n",
    "ax.set_ylim(-0.01, 0.4001)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>**Execute this code cell**</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(hx(7.132))             # h(x), the true distribution\n",
    "\n",
    "print(fx(np.array([7.132]))) # f(x) for a parametric model\n",
    "print(np.exp(kde.score_samples(np.atleast_2d(7.132)))) # f(x) for non-parametric model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A note on uncertainties and errors\n",
    "\n",
    "* Technically, errors are systematic biases that we can not mitigate through collecting lots and lots of data. \n",
    "* Statistical uncertainties are the result of random measurement uncertainty. \n",
    "* But \"error\" will be used for both, and denoted as either statistical errors (error bars) or systematic errors (biases).\n",
    "\n",
    "\n",
    "* Statistical error distributions (error bars) that vary from data point to data point are called **heteroscedastic errors**. If they are the same for all points then they are **homoscedastic errors**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##  We could summarize the goal of this course as an attempt to \n",
    "\n",
    "1) estimate $f(x)$ from some real (possibly multi-dimensional) data set, \n",
    "\n",
    "2) find a way to describe $f(x)$ and its uncertainty, \n",
    "\n",
    "3) compare it to models of $h(x)$, and then \n",
    "\n",
    "4) use the knowledge that we have gained to interpret new measurements."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Probability <a class=\"anchor\" id=\"two\"></a>\n",
    "\n",
    "> $p(A)$ = the probability of $A$ (or the probability density at $A$), \n",
    "\n",
    "e.g. the probability that an observed object is a galaxy. This does not mean that the object is in some sort of Schrodinger's cat quantum uncertainity...*the probability reflects our current state of knowledge of the object, and our belief that it is a galaxy*. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Kolmogorov Axioms (knowing that they're called this is good for trivia)\n",
    "\n",
    "1. $p(A)\\geq0 \\quad\\forall\\, A$\n",
    "2. $p(\\Omega) = 1$, where $\\Omega$ is the set of all outcomes, i.e. the sum/integral of all possible outcomes is unity\n",
    "3. $p(\\cup_{i=1}^\\infty A_i) = \\sum_{i=1}^\\infty p(A_i)$ if all events are independent\n",
    "\n",
    "$A \\cup B$ is the *union* of sets $A$ and $B$. **Read as A OR B.**\n",
    "\n",
    "$A \\cap B$ is the *intersection* of sets $A$ and $B$. **Read as A AND B.** Different notations $p(A \\cap B) = p(AB) = p(A,B) = p(A\\,\\mathrm{and}\\,B)$. We will use the comma notation throughout. \n",
    "\n",
    "If we have two events, $A$ and $B$, the possible combinations are illustrated by the following figure:\n",
    "![Figure 3.1](http://www.astroml.org/_images/fig_prob_sum_1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "* The probability that ***either*** $A$ ***or*** $B$ will happen (which could include both) is the *union*, given by\n",
    "$$p(A \\cup B) = p(A) + p(B) - p(A \\cap B)$$\n",
    "The figure makes it clear why the last term is necessary.  Since $A$ and $B$ overlap, we are double-counting the region where *both* $A$ and $B$ happen, so we have to subtract this out.  \n",
    "\n",
    "\n",
    "* The probability that ***both*** $A$ ***and*** $B$ will happen, $p(A \\cap B)$, is \n",
    "$$p(A \\cap B) = p(A|B)p(B) = p(B|A)p(A)$$\n",
    "where p(A|B) is the probability of A *given that* B is true and is called the **conditional probability**.  So the $|$ is short for \"given\".\n",
    "\n",
    "\n",
    "* The **law of total probability** says that (for independent $B_i$)\n",
    "$$p(A) = \\sum_ip(A|B_i)p(B_i)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "It is important to realize that the following is *always* true:\n",
    "\n",
    "$$p(A,B) = p(A|B)p(B) = p(B|A)p(A)$$\n",
    "\n",
    "However, if $A$ and $B$ are independent, then $p(A|B)=p(A)$ and $p(B|A)=p(B)$ and\n",
    "\n",
    "$$p(A,B) = p(A)p(B)$$\n",
    "\n",
    "**EXAMPLE** (classic marbles in bag scenario)\n",
    "\n",
    "If you have a bag with 5 marbles (3 yellow and 2 blue) and you want to know the probability of picking 2 yellow marbles in a row, that would be\n",
    "\n",
    "$$p(Y_1,Y_2) = p(Y_1)p(Y_2|Y_1).$$\n",
    "\n",
    "$p(Y_1) = \\frac{3}{5}$ since you have an equally likely chance of drawing any of the 5 marbles.\n",
    "\n",
    "If you did not put the first marble back in the back after drawing it (sampling *without* \"replacement\"), then the probability\n",
    "\n",
    "$p(Y_2|Y_1) = \\frac{2}{4}$, so that\n",
    "\n",
    "$$p(Y_1,Y_2) = \\frac{3}{5}\\frac{2}{4} = \\frac{3}{10}.$$\n",
    "\n",
    "But if you put the first marble back, then\n",
    "\n",
    "$p(Y_2|Y_1) = \\frac{3}{5} = p(Y_2)$, so that \n",
    "\n",
    "$$p(Y_1,Y_2) = \\frac{3}{5}\\frac{3}{5} = \\frac{9}{25}.$$\n",
    "\n",
    "In the first case $A$ and $B$ (or rather $Y_1$ and $Y_2$) are *not* independent, whereas in the second case they are."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Here is a more complicated example from \n",
    "[Jo Bovy's class at UToronto](http://astro.utoronto.ca/%7Ebovy/teaching.html)\n",
    "![Bovy_L1-StatMiniCourse_page21](figures/JoBovy_L1-StatsMiniCourse.png)\n",
    "\n",
    "As illustrated, \n",
    "\n",
    "$$p(A \\,{\\rm or}\\, B|C) = p(A|C) + p(B|C) - p(A \\, {\\rm and}\\, B|C)$$ \n",
    "\n",
    "Need more help with this?  Try watching some Khan Academy videos and working through the exercises:\n",
    "* [https://www.khanacademy.org/math/probability/probability-geometry](https://www.khanacademy.org/math/probability/probability-geometry)\n",
    "* [https://www.khanacademy.org/math/precalculus/prob-comb](https://www.khanacademy.org/math/precalculus/prob-comb)\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bayes' Theorem <a class=\"anchor\" id=\"three\"></a>\n",
    "\n",
    "In the following figure, we have a 2-D distribution in $x-y$ parameter space.  Here $x$ and $y$ are ***not*** independent as, once you pick a $y$, your values of $x$ are constrained.\n",
    "\n",
    "![http://www.astroml.org/_images/fig_conditional_probability_1.png](http://www.astroml.org/_images/fig_conditional_probability_1.png)\n",
    "\n",
    "We have that \n",
    "$$p(x,y) = p(x|y)p(y) = p(y|x)p(x)$$\n",
    "\n",
    "We can define the ***marginal probability*** as\n",
    "$$p(x) = \\int p(x,y)dy,$$\n",
    "\n",
    "where marginal means projecting on to one axis (integrating over the unwanted variable). The **marginal** distributions are shown on the left and bottom sides of the left panel.  As the equation above says, this is just the integral along the $x$ direction for a given $y$ (left side panel) or the integral along the $y$ direction for a given $x$ (bottom panel).  \n",
    "\n",
    "The three panels on the right show the ***conditional probability*** (of $x$) for three $y$ values: $$p(x|y=y_0)$$  These are just normalized \"slices\" through the 2-D distribution.\n",
    "\n",
    "The marginal probability of $x$ can be re-written as\n",
    "\n",
    "$$p(x) = \\int p(x|y)p(y) dy$$\n",
    "\n",
    "But since $p(x|y)p(y) = p(y|x)p(x)$, we can write\n",
    "\n",
    "> $$p(y|x) = \\frac{p(x|y)p(y)}{p(x)} = \\frac{p(x|y)p(y)}{\\int p(x|y)p(y) dy}$$\n",
    "\n",
    "which in words says that\n",
    "\n",
    "> the (conditional) probability of $y$ given $x$ is just the (conditional) probability of $x$ given $y$ times the (marginal) probability of $y$ divided by the (marginal) probability of $x$, where the latter is just the integral of the numerator.\n",
    "\n",
    "This is **Bayes' Theorem**, which itself is not at all controversial, though its application can be as we'll discuss later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bayes' Theorem Example: Lego's <a class=\"anchor\" id=\"four\"></a>\n",
    "\n",
    "An example with Lego's (it's awesome):\n",
    "[https://www.countbayesie.com/blog/2015/2/18/bayes-theorem-with-lego](https://www.countbayesie.com/blog/2015/2/18/bayes-theorem-with-lego)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bayes' Theorem Example: Monty Hall Problem (or \"Deal Or No Deal\") <a class=\"anchor\" id=\"five\"></a>\n",
    "\n",
    "You are playing a game show and are shown 2 doors.  One has a car behind it, the other a goat.  What are your chances of picking the door with the car?\n",
    "\n",
    "<font color='red'>ANSWER =</font> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "OK, now there are 3 doors: one with a car, two with goats.  The game show host asks you to pick a door, but not to open it yet.  Then the host opens one of the other two doors (that you did not pick), making sure to select one with a goat.  The host offers you the opportunity to switch doors.  Do you?\n",
    " \n",
    "![https://upload.wikimedia.org/wikipedia/commons/thumb/3/3f/Monty_open_door.svg/180px-Monty_open_door.svg.png](https://upload.wikimedia.org/wikipedia/commons/thumb/3/3f/Monty_open_door.svg/180px-Monty_open_door.svg.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Now you are back at the 2 door situation.  But what can you make of your prior information? Let's break this down."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "* Probability of car behind **Door 1** $= 1/3$\n",
    "* Probability of car behind **Doors 2 or 3** $= 2/3$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "But if the host opens **Door 3** to show a goat, you now have prior information to update your probability of the car being behind **Door 2**. The probability of the car behind **Door 1** is still $1/3$, but the total probability must sum to unity, and you only have **Door 2** left. Hence,\n",
    "\n",
    "* With prior knowledge of **Door 3** containing a goat, the probability of car behind **Door 2** $=2/3$.\n",
    "\n",
    "\n",
    "So, ***switching doubles your chances***. You should always switch!! But only because you had prior information.  If someone walked in after the \"bad\" door was opened, then their probability of winning is the expected $1/2$.\n",
    "\n",
    "Try it:\n",
    "https://betterexplained.com/articles/understanding-the-monty-hall-problem/\n",
    "\n",
    "<font color='red'>Add a screen shot showing your results after playing the Monty Hall game above 10 times, either with \"pick and hold\" or \"pick and switch\". You can insert an image in Jupyter by clicking \"Edit\", then \"Insert Image\" down at the bottom. Run the cell to get the image to appear.</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "This actually becomes easier to understand as $N$ increases. For $N$ choices, revealing $N-2$ \"answers\" doesn't change the probability of your original door choice. It is still $\\frac{1}{N}$.  But it *does* change the probability of your knowledge of the *other final* remaining choice by $N-1$ and it is $\\frac{N-1}{N}$. Therefore, by switching, you increase your chance of winning by a factor of $(N-1)$. Shocking, but true. \n",
    "\n",
    "* In the 3-door example, switching doubles your chance of winning (from $1/3$ to $2/3$). \n",
    "* In a 100-door example, switching increases your chance of winning by a factor of $99$.\n",
    "\n",
    "This is an example of the use of *conditional* probability, where we have $p(A|B) \\ne p(A)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bayes' Theorem Example: Contingency Table <a class=\"anchor\" id=\"six\"></a>\n",
    "\n",
    "We can also use Bayes' rule to learn something about false positives and false negatives. **Note below that I am going to give a COVID-19 example below, so if this is a sensitive topic for you then please let me know.**\n",
    "\n",
    "Let's say that we have a test for a disease.  The test can be positive ($T=1$) or negative ($T=0$) and one can either have the disease ($D=1$) or not ($D=0$).  So, there are 4 possible combinations:\n",
    "$$T=0; D=0 \\;\\;\\;  {\\rm true \\; negative}$$\n",
    "$$T=0; D=1 \\;\\;\\; {\\rm false \\; negative}$$\n",
    "$$T=1; D=0 \\;\\;\\; {\\rm false \\; positive}$$\n",
    "$$T=1; D=1 \\;\\;\\; {\\rm true \\; positive}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "All else being equal, you have a 50% chance of being misdiagnosed.  Not good!  But the probability of disease and the accuracy of the test presumably are not random.\n",
    "\n",
    "If the rates of false positive and false negative are:\n",
    "$$p(T=1|D=0) = \\epsilon_{\\rm FP}$$\n",
    "$$p(T=0|D=1) = \\epsilon_{\\rm FN}$$\n",
    "\n",
    "then the true positive and true negative rates are just:\n",
    "$$p(T=0| D=0) = 1-\\epsilon_{\\rm FP}$$\n",
    "$$p(T=1| D=1) = 1-\\epsilon_{\\rm FN}$$\n",
    "\n",
    "Let's assume that $\\epsilon_{\\rm FP}=0.02$ and $\\epsilon_{\\rm FN}=0.001$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In graphical form this $p(T|D)$ matrix is:\n",
    "![http://www.astroml.org/_images/fig_contingency_table_1.png](http://www.astroml.org/_images/fig_contingency_table_1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "If we have a **prior** regarding how likely the disease is, we can take this into account.\n",
    "\n",
    "$$p(D=1)=\\epsilon_D$$\n",
    "\n",
    "and then $p(D=0)=1-\\epsilon_D$. Say, $\\epsilon_D$ = 0.01. \n",
    "\n",
    "Now assume that a person tested positive. What is the probability that this person has the disease? Is it 98% \n",
    "because $\\epsilon_{\\rm FP}=0.02$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "We can't just read $p(D=1|T=1)$ off the table because the table entry is the conditional probability of the *test* given the *data*, $p(T=1|D=1)$. What we want is the conditional probability of the *data* given the *test*, that is, $p(D=1|T=1)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Bayes' Theorem then can be used to help us determine how likely it is that you have the disease if you tested positive:\n",
    "\n",
    "$$p(D=1|T=1) = \\frac{p(T=1|D=1)p(D=1)}{p(T=1)},$$\n",
    "\n",
    "where $$p(T=1) = p(T=1|D=0)p(D=0) + p(T=1|D=1)p(D=1).$$\n",
    "\n",
    "So\n",
    "$$p(D=1|T=1) = \\frac{(1 - \\epsilon_{FN})\\epsilon_D}{\\epsilon_{FP}(1-\\epsilon_D) + (1-\\epsilon_{FN})\\epsilon_D} \\approx \\frac{\\epsilon_D}{\\epsilon_{FP}+\\epsilon_D}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "where in the final approximation we assume that all $\\epsilon$ values are small. That means that to get a reliable diagnosis, we need $\\epsilon_{FP}$ to be quite small.  (Because you *want* the probability to be close to unity if you test positive, otherwise it is a *false* positive).\n",
    "\n",
    "In our example, we have a disease rate of 1% ($\\epsilon_D = 0.01$) and a false positive rate of 2% ($\\epsilon_{\\rm FP}=0.02$).  \n",
    "\n",
    "So we have\n",
    "$$p(D=1|T=1) = \\frac{0.01}{0.02+0.01} = 0.333$$\n",
    "\n",
    "Then in a sample of, e.g.,  1000 people, 10 people will *actually* have the disease $(1000*0.01)$, but another 20 $(1000*0.02)$ will test positive! \n",
    "\n",
    "Therefore, in that sample of 30 people who tested positive, only 1/3 has the disease\n",
    "(not 98%!). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### COVID-19 example\n",
    "\n",
    "Current estimates of the rate of **false positive** swab tests for COVID-19 are between $0.8\\%$ and $4.0\\%$ (for the UK; https://www.thelancet.com/journals/lanres/article/PIIS2213-2600(20)30453-7/fulltext).\n",
    "\n",
    "Current estimates of the rate of **false negative** swab tests for COVID-19 are $\\sim20\\%$ (https://www.health.harvard.edu/blog/which-test-is-best-for-covid-19-2020081020734#:~:text=The%20reported%20rate%20of%20false,over%20just%20a%20few%20months).\n",
    "\n",
    "Current daily average infection count per 100,000 in Davidson County is 78 (https://www.nytimes.com/interactive/2020/us/tennessee-coronavirus-cases.html). Let's assume that translates to a current **disease probability** of $78/100000 = 0.078\\%$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>Plug in some of these numbers to estimate $p(D=1|T=1)$ for COVID-19 in Davidson County.</font>\n",
    "\n",
    "ANSWER here-->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transformations of random variables <a class=\"anchor\" id=\"seven\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $x$ is a random variable then $f(x)$ is also a random variable for any function $f$.\n",
    "\n",
    "To transform probability distributions when taking functions of random variables, we can simply use conservation of dimensionless probability, i.e. \n",
    "\n",
    "$$\\mathrm{Prob}(x, x+dx) = \\mathrm{Prob}(y, y+dy)$$\n",
    "\n",
    "$$p(x)dx = p(y)dy.$$ \n",
    "\n",
    "where $y = f(x)$.\n",
    "\n",
    "Thus, $$p(y) = \\left|\\frac{dx}{dy}\\right| p(x)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**EXAMPLE**\n",
    "\n",
    "* Let $x$ be a random variable drawn from a uniform distribution between $0$ and $1$. So $p(x) = 1/(1-0) = 1$.  \n",
    "* Let's transform to $y = e^x$.\n",
    "* So $p(y) = \\left|dy/dx\\right|^{-1}p(x) = 1/y$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![https://www.astroml.org/_images/fig_transform_distribution_1.png](https://www.astroml.org/_images/fig_transform_distribution_1.png)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:astr8070] *",
   "language": "python",
   "name": "conda-env-astr8070-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.8.5"
  },
  "livereveal": {
   "scroll": true,
   "start_slideshow_at": "selected",
   "theme": "sky"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}