11_inference.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "*Elements of Data Science*\n",
    "\n",
    "Copyright 2021 [Allen B. Downey](https://allendowney.com)\n",
    "\n",
    "License: [Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-nc-sa/4.0/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Click here to run this notebook on Colab](https://colab.research.google.com/github/AllenDowney/ElementsOfDataScience/blob/master/11_inference.ipynb) or\n",
    "[click here to download it](https://github.com/AllenDowney/ElementsOfDataScience/raw/master/11_inference.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This chapter introduces **statistical inference**, which is the process of using a sample to make inferences about a population.  \n",
    "\n",
    "The **population** is the group we are interested in.  Sometimes it is a group of people, but in general it can be any kind of group.\n",
    "\n",
    "If we can observe the entire group, we might not need statistical inference.  If not, sometimes we can observe a **sample** (or subset) of the population and use those observations to make claims about the population.\n",
    "\n",
    "If you have studied statistics before, you might have encountered some of these ideas before: hypothesis testing, p-values, estimation, standard error, and confidence intervals.\n",
    "\n",
    "In this chapter we'll approach these topics using computation and simulation, as opposed to mathematical analysis.  I hope this approach makes the ideas clearer.\n",
    "\n",
    "If you have not seen these ideas before, don't worry.  That might even be better.\n",
    "\n",
    "We'll look at three examples:\n",
    "\n",
    "* Testing whether a coin is \"fair\".\n",
    "\n",
    "* Testing whether first babies are more like to be born early (or late).\n",
    "\n",
    "* Estimating the average height of men in the U.S., and quantifying the precision of the estimate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Euro problem\n",
    "\n",
    "In *Information Theory, Inference, and Learning Algorithms*, David MacKay writes, \"A statistical statement appeared in *The Guardian* on Friday January 4, 2002:\n",
    "\n",
    "> When spun on edge 250 times, a Belgian one-euro coin came\n",
    "up heads 140 times and tails 110. ‘It looks very suspicious\n",
    "to me’, said Barry Blight, a statistics lecturer at the London\n",
    "School of Economics. ‘If the coin were unbiased the chance of\n",
    "getting a result as extreme as that would be less than 7%’.*\n",
    "\n",
    "But do these data give evidence that the coin is biased rather than fair?\"\n",
    "\n",
    "Before we answer MacKay's question, let's unpack what Dr. Blight said:\n",
    "\n",
    "\"If the coin were unbiased the chance of getting a result as extreme as that would be less than 7%\".\n",
    "\n",
    "To see where that comes from, I'll simulate the result of spinning an \"unbiased\" coin, meaning that the probability of heads is 50%.\n",
    "\n",
    "Here's an example with 10 spins:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "spins = np.random.random(10) < 0.5\n",
    "spins"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`np.random.random` returns numbers between 0 and 1, uniformly distributed.  So the probability of being less than 0.5 is 50%.  \n",
    "\n",
    "The sum of the array is the number of `True` elements, that is, the number of heads:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(spins)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can wrap that in a function that simulates `n` spins with probability `p`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def spin(n, p):\n",
    "    return np.sum(np.random.random(n) < p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's an example with the actual sample size (250) and hypothetical probability (50%)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "heads, tails = 140, 110\n",
    "sample_size = heads + tails"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "hypo_prob = 0.5\n",
    "spin(sample_size, hypo_prob)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we are generating random numbers, we expect to see different values if we run the experiment more than once.\n",
    "\n",
    "Here's a loop that runs `spin` 10 times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 250\n",
    "p = 0.5\n",
    "\n",
    "for i in range(10):\n",
    "    print(spin(n, p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, the results vary from one run to the next.\n",
    "\n",
    "Now let's run the simulated experiment 10000 times and store the results in a NumPy array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "outcomes = np.empty(1000)\n",
    "\n",
    "for i in range(len(outcomes)):\n",
    "    outcomes[i] = spin(n, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`np.empty` creates an empty array with the given length.  Each time through the loop, we run `spin` and assign the result to an element of `outcomes`.\n",
    "\n",
    "The result is an array of 10000 integers, each representing the number of heads in a simulated experiment.  \n",
    "\n",
    "The mean of `outcomes` is about 125:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(outcomes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which makes sense.  On average, the expected number of heads is the product of the hypothetical probability and the sample size:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "expected = hypo_prob * sample_size\n",
    "expected"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's see how much the values in `outcomes` differ from the expected value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "diffs = outcomes - expected"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`diffs` is an array that contains the deviation of each experiment from the expected value, 125.\n",
    "\n",
    "Here's the mean of the absolute deviations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(abs(diffs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So a typical experiment deviates from the mean by about 6.\n",
    "\n",
    "To see the whole distribution of deviations, we can plot a histogram.\n",
    "The following function uses Matplotlib to plot a histogram and adjust some of the settings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def plot_hist(values):\n",
    "    xs, ys, patches = plt.hist(values,\n",
    "                               density=True,\n",
    "                               histtype='step',\n",
    "                               linewidth=2,\n",
    "                               alpha=0.5)\n",
    "    \n",
    "    \n",
    "    plt.ylabel('Density')\n",
    "    plt.tight_layout()\n",
    "    return patches[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the distribution of deviations looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_hist(diffs)\n",
    "\n",
    "plt.title('Sampling distribution (n=250)')\n",
    "plt.xlabel('Deviation from expected number of heads');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the \"sampling distribution\" of deviations.  It shows how much variation we should expect between experiments with this sample size (n = 250)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## P-values\n",
    "\n",
    "Getting get back to the Euro example, Dr. Bright reported:\n",
    "\n",
    "\"If the coin were unbiased the chance of getting a result as extreme as that would be less than 7%\".\n",
    "\n",
    "The article doesn't say so explicitly, but this is a \"p-value\".  To understand what that means, let's count how many times, in 10000 attempts, the outcome is \"as extreme as\" the observed outcome, 140 heads.\n",
    "\n",
    "The observed deviation is the difference between the observed and expected number of heads:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "observed_diff = heads - expected\n",
    "observed_diff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how many times the simulated `diffs` exceed the observed deviation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(diffs >= observed_diff)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's around 3%.  But Dr. Blight said 7%.  Where did that come from?\n",
    "\n",
    "So far, we only counted the cases where the outcome is *more* heads than expected.  We might also want to count the cases where the outcome is *fewer* than expected.\n",
    "\n",
    "Here's the probability of falling below the expected number by 15 or more."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(diffs <= -observed_diff)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the total probability of a result \"as extreme as that\", we can use the absolute value of the simulated differences:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(abs(diffs) >= observed_diff)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So that's consistent with what Dr. Blight reported.\n",
    "\n",
    "To show what that looks like graphically, I'll use the following function, which fills in the histogram between `low` and `high`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fill_hist(low, high, patch):\n",
    "    fill = plt.axvspan(low, high, \n",
    "                       clip_path=patch,\n",
    "                       alpha=0.5, \n",
    "                       color='C0')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following plot shows the sampling distribution of `diffs` with two regions shaded.  These regions represent the probability that an unbiased coin yields a deviation from the expected as extreme as 15."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "patch = plot_hist(diffs)\n",
    "\n",
    "# fill the right tail of the hist\n",
    "low = observed_diff\n",
    "high = diffs.max()\n",
    "fill_hist(low, high, patch)\n",
    "\n",
    "# fill the left tail of the hist\n",
    "low = diffs.min()\n",
    "high = -observed_diff\n",
    "fill_hist(low, high, patch)\n",
    "\n",
    "plt.title('Sampling distribution (n=250)')\n",
    "plt.xlabel('Deviation from expected number of heads')\n",
    "plt.ylabel('Density');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These results show that there is a non-negligible chance of getting a result as extreme as 140 heads, even if the coin is actually fair.\n",
    "\n",
    "So even if the results are \"suspicious\" they don't provide compelling evidence that the coin is biased."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** There are a few ways to make \"crooked\" dice, that is, dice that are more likely to land on one side than the others.  Suppose you run a casino and you suspect that a patron is using a die that comes up 3 more often than it should.  \n",
    "\n",
    "You confiscate the die, roll it 300 times, and 63 times it comes up 3.  Does that support your suspicions?\n",
    "\n",
    "- To answer this question, use `spin` to simulate the experiment, assuming that the die is fair.\n",
    "\n",
    "- Use a for loop to run `spin` 1000 times and store the results in a NumPy array.\n",
    "\n",
    "- What is the mean of the results from the simulated experiments?\n",
    "\n",
    "- What is the expected number of 3s if the die is fair?\n",
    "\n",
    "- Use `plot_hist` to plot the results.  The histogram you plot approximates the sampling distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Continuing the previous exercise, compute the probability of seeing a deviation from the expected value that is \"as extreme\" as the observed difference.\n",
    "\n",
    "For this context, what do you think is the best definition of \"as extreme\"?\n",
    "\n",
    "Plot the histogram of the random deviations again, and use `fill_hist` to fill the region of the histogram that corresponds to the p-value you computed. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Are first babies more likely to be late?\n",
    "\n",
    "The examples so far have been based on coins and dice, which are relatively simple.  In this section we'll look at an example that's based on real-world data.\n",
    "\n",
    "Here's the motivation for it: When my wife and I were expecting our first baby, we heard that first babies are more likely to be late.  We also hear that first babies are more likely to be early.  Neither claim was supported by evidence.\n",
    "\n",
    "Fortunately, I am a data scientist!  Also fortunately, the CDC runs the National Survey of Family Growth (NSFG), which \"gathers information on family life, marriage and divorce, pregnancy, infertility, use of contraception, and men’s and women’s health.\"\n",
    "\n",
    "I got the data from their web page, https://www.cdc.gov/nchs/nsfg/index.htm, and wrote some code to get it into a Pandas Dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from os.path import basename, exists\n",
    "\n",
    "def download(url):\n",
    "    filename = basename(url)\n",
    "    if not exists(filename):\n",
    "        from urllib.request import urlretrieve\n",
    "        local, _ = urlretrieve(url, filename)\n",
    "        print('Downloaded ' + local)\n",
    "    \n",
    "download('https://github.com/AllenDowney/' +\n",
    "         'ElementsOfDataScience/raw/master/nsfg.hdf5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "nsfg = pd.read_hdf('nsfg.hdf5')\n",
    "nsfg.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `nsfg` DataFrame contains 9358 rows, one for each recorded pregnancy, and 11 columns, one of each of the variables I selected.\n",
    "\n",
    "Here are the first few lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsfg.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variables we need are `birthord`, which indicates birth order, and `prglength`, which is pregnancy length in weeks.\n",
    "\n",
    "I'll make two boolean Series, one for first babies and one for others."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "firsts = (nsfg.birthord == 1)\n",
    "others = (nsfg.birthord > 1)\n",
    "\n",
    "np.sum(firsts), np.sum(others)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the boolean Series to select pregnancy lengths for the two groups and compute their means:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "prglngth = nsfg['prglngth']\n",
    "\n",
    "mean_first = prglngth[firsts].mean() \n",
    "mean_other = prglngth[others].mean()\n",
    "\n",
    "mean_first, mean_other"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the difference in means, in weeks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "diff = mean_first - mean_other\n",
    "diff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here it is converted to days:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "diff * 7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like first babies are born 1.4 days later than other babies, on average. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hypothesis testing\n",
    "\n",
    "The apparent difference between these groups is based on a random sample that is much smaller than the actual population.  So we can't be sure that the difference we see in the sample reflects a real difference in the population.  There are two other possibilities we should keep in mind:\n",
    "\n",
    "*  Systematic errors: The sample might be more likely to include some pregancies, and less likely to include others, in a way that causes an apparent difference in the sample, even if there is no such difference in the population.\n",
    "\n",
    "*  Sampling errors: Even if every pregnancy is equally likely to appear in the sample, it is still possible to see a difference in the sample that is not in the population, just because of random variation.\n",
    "\n",
    "We can never rule out the possibility of systematic errors, but we *can* test whether an apparent effect could be explained by random sampling.\n",
    "\n",
    "Here's how:\n",
    "\n",
    "1. First we choose a \"test statistic\" that measures the size of the effect; the test statistic in this example is the difference in mean pregnancy length.\n",
    "\n",
    "2. Next we define a model of the population under the assumption that there is actually no difference between the groups.  This assumption is called the \"null hypothesis\".\n",
    "\n",
    "3. Then we use the model to compute the distribution of the test statistic under the null hypothesis.\n",
    "\n",
    "We have already done step 1, but to make it easier to repeat, I'll wrap it in a function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_stat(group1, group2):\n",
    "    \"\"\"Difference in means.\n",
    "    \n",
    "    group1: sequence of values\n",
    "    group2: sequence of values\n",
    "    \n",
    "    returns: float difference in means\n",
    "    \"\"\"\n",
    "    diff = np.mean(group1) - np.mean(group2)\n",
    "    return diff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`test_stat` takes two sequences and computes the difference in their means.\n",
    "\n",
    "Here's how we use it to compute the actual difference in the sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "group1 = prglngth[firsts]\n",
    "group2 = prglngth[others]\n",
    "\n",
    "actual_diff = test_stat(group1, group2)\n",
    "actual_diff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have to define the null hypothesis, which is a model of the world where there is no difference in pregnancy length between first babies and others.\n",
    "\n",
    "One way to do that is to put the two groups together and then divide them up again at random.  That way the distribution of pregnancy lengths is the same for both groups.\n",
    "\n",
    "I'll use `concatenate` to pool the groups. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(group1), len(group2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "pool = np.concatenate([group1, group2])\n",
    "pool.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I'll use `shuffle` to reorder them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.shuffle(pool)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then I'll use `split` to make two simulated groups, the same size as the originals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = len(group1)\n",
    "sim_group1, sim_group2 = np.split(pool, [n])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(sim_group1), len(sim_group2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can compute the test statistic for the simulated data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_stat(sim_group1, sim_group2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the simulated data, the distribution of pregnancy lengths is the same for both groups, so the difference is usually close to 0.\n",
    "\n",
    "But because it is based on a random shuffle of the groups, we get a different value each time we run it.\n",
    "\n",
    "To see what the whole distribution looks like, we can run the simulation many times and store the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "diffs = np.empty(1000)\n",
    "\n",
    "for i in range(len(diffs)):\n",
    "    np.random.shuffle(pool)\n",
    "    sim_group1, sim_group2 = np.split(pool, [n])\n",
    "    diffs[i] = test_stat(sim_group1, sim_group2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is the \"sampling distribution of the test statistic under the null hypothesis\".\n",
    "\n",
    "The mean of this distribution should close to zero, because it is based on the assumption that there is actually no difference between the groups."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(diffs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's what the whole distribution looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_hist(diffs)\n",
    "\n",
    "plt.xlabel('Difference in mean (weeks)')\n",
    "plt.title('Distribution of test statistic under null hypothesis');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If there were actually no difference between the groups, we would expect to see a difference as big as 0.15 weeks by chance, at least occasionally.  But a difference as big as 0.2 would be rare.\n",
    "\n",
    "To quantify that surprise, we can estimate the probability that the test statistic, under the null hypothesis, exceeds the observed differences in the means.\n",
    "\n",
    "The result is a \"p-value\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_value = np.mean(diffs >= actual_diff)\n",
    "p_value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example the result is 0, which is to say that in 1000 simulations of the null hypothesis, we never saw a difference as big as 0.2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpreting p-values\n",
    "\n",
    "To interpret this result, remember that we started with three possible explanations for the observed difference between the groups:\n",
    "\n",
    "1. The observed difference might be \"real\"; that is, there might be an actual difference in pregnancy length between first babies and others.\n",
    "\n",
    "2. There might be no real difference between the groups, and the observed difference might be because of a systematic error in the sampling process or the data collection process.  For example, maybe reported pregnancy lengths are less accurate for first time mothers.\n",
    "\n",
    "3. There might be no real difference between the groups, and the observed difference might be due to random variation in the sampling process.\n",
    "\n",
    "By computing a p-value, we have established that it would be rare to see a difference as big as 0.2 due to sampling alone.  So we can conclude that the third explanation is unlikely.\n",
    "\n",
    "That makes it more likely that the difference is real, but we still can't rule out the second possibility."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** The test statistic we chose is the difference in means between the two groups.\n",
    "\n",
    "But suppose we would like to know whether first babies are more unpredictable than other babies.  In that case the test statistic we choose might be the standard deviation of pregnancy length, which is one way to quantify unpredictability.\n",
    "\n",
    "As an exercise:\n",
    "\n",
    "1. Write a version of `test_stat` that computes the difference in standard deviation between the groups.\n",
    "\n",
    "2. Write a loop that estimates the distribution of this test statistic under the null hypothesis.\n",
    "\n",
    "3. Compute a p-value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimation\n",
    "\n",
    "Suppose we want to estimate the average height of men in the U.S.\n",
    "\n",
    "We can use data from the [BRFSS](https://www.cdc.gov/brfss/index.html):\n",
    "\n",
    "\"The Behavioral Risk Factor Surveillance System (BRFSS) is the nation's premier system of health-related telephone surveys that collect state data about U.S. residents regarding their health-related risk behaviors, chronic health conditions, and use of preventive services.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "from os.path import basename, exists\n",
    "\n",
    "def download(url):\n",
    "    filename = basename(url)\n",
    "    if not exists(filename):\n",
    "        from urllib.request import urlretrieve\n",
    "        local, _ = urlretrieve(url, filename)\n",
    "        print('Downloaded ' + local)\n",
    "    \n",
    "download('https://github.com/AllenDowney/' +\n",
    "         'ElementsOfDataScience/raw/master/brfss.hdf5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "brfss = pd.read_hdf('brfss.hdf5', 'brfss')\n",
    "brfss.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use `SEX` to select male respondents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "male = (brfss.SEX == 1)\n",
    "np.mean(male)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we select height data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "heights = brfss['HTM4']\n",
    "data = heights[male]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use `isnan` to check for NaN values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(np.isnan(data)) * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "About 4% of the values are missing.\n",
    "\n",
    "Here are the mean and standard deviation, ignoring missing data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Mean male height in cm =', np.nanmean(data))\n",
    "print('Std male height in cm =', np.nanstd(data))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantifying precision\n",
    "\n",
    "At this point we have an estimate of the average adult male height.  We'd like to know how accurate this estimate is, and how precise.  In the context of estimation, these words have a [technical distinction](https://en.wikipedia.org/wiki/Accuracy_and_precision):\n",
    "\n",
    ">Given a set of data points from repeated measurements of the same quantity, the set can be said to be precise if the values are close to each other, while the set can be said to be accurate if their average is close to the true value of the quantity being measured.\n",
    "\n",
    "Usually accuracy is what we really care about, but it's hard to measure accuracy unless you know the true value.  And if you know the true value, you don't have to estimate it.\n",
    "\n",
    "Quantifying precision is not as useful, but it is much easier.  Here's one way to do it:\n",
    "\n",
    "1.  Use the data you have to make a model of the population.\n",
    "\n",
    "2.  Use the model to simulate the data collection process.\n",
    "\n",
    "3.  Use the simulated data to compute an estimate.\n",
    "\n",
    "By repeating these steps, we can quantify the variability of the estimate due to random sampling."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To model the population, I'll use **resampling**; that is, I will treat the observed measurements as if they were taken from the entire population, and I will draw random samples from them.\n",
    "\n",
    "We can use `np.random.choice` to resample the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "size = len(data)\n",
    "sim_data = np.random.choice(data, size, replace=True)\n",
    "sim_data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With `replace=True`, we sample with replacement, which means that some measurements might be chosen more than once, and some might not be chosen at all.\n",
    "\n",
    "(If we sample *without* replacement, the resampled data is always identical to the original, so that's no good.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use `nanmean` to compute the mean of the simulated data, ignoring missing values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.nanmean(sim_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we repeat this process 1000 times, we can see how much the results vary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "outcomes = np.empty(1000)\n",
    "size = len(data)\n",
    "\n",
    "for i in range(len(outcomes)):\n",
    "    sim_data = np.random.choice(data, size, replace=True)\n",
    "    outcomes[i] = np.nanmean(sim_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is the \"sampling distribution\", which shows how much the results of the experiment would vary if we ran it many times.  Here's what it looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_hist(outcomes)\n",
    "plt.title('Sampling distribution of the mean')\n",
    "plt.xlabel('Mean adult male height, U.S.');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The width of this distribution shows how much the results vary from one experiment to the next.\n",
    "\n",
    "We can quantify this variability by computing the standard deviation of the sampling distribution, which is called \"standard error\".\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "std_err = np.std(outcomes)\n",
    "std_err"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also summarize the sampling distribution with a \"confidence interval\", which is a range that contains a specified fraction, like 90%, of the values in `sampling_dist_mean`.\n",
    "\n",
    "The central 90% confidence interval is between the 5th and 95th percentiles of the sampling distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "ci_90 = np.percentile(outcomes, [5, 95])\n",
    "ci_90"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function plots a histogram and shades the 90% confidence interval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_sampling_dist(outcomes):\n",
    "    \"\"\"Plot sampling distribution.\n",
    "    \n",
    "    outcomes: sequence of values\n",
    "    \"\"\"\n",
    "    patch = plot_hist(outcomes)\n",
    "    low, high = np.percentile(outcomes, [5, 95])\n",
    "    fill_hist(low, high, patch)\n",
    "    print('Mean = ', np.mean(outcomes))\n",
    "    print('Std error = ', np.std(outcomes))\n",
    "    print('90% CI = ', (low, high))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what it looks like for the sampling distribution of mean adult height:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_sampling_dist(outcomes)\n",
    "plt.xlabel('Mean adult male height, U.S. (%)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For an experiment like this, we can compute the standard error analytically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "size = len(data)\n",
    "analytic_std_err = np.std(data) / np.sqrt(size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is close to what we observed computationally."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "analytic_std_err, std_err"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This result indicates that our estimate of the mean is *precise*; that is, if we ran this experiment many times, the results would fall in a narrow range.\n",
    "\n",
    "But this range reflects only variability due to random sampling.  If there are systematic errors in the sampling process, or in the measurement process, the result would not be *accurate*.\n",
    "\n",
    "Computing a standard error or confidence interval can be useful, but it only quantifies variability due to random sampling, not other sources of error."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** One nice thing about using resampling is that it is easy to compute the sampling distribution for other statistics.\n",
    "\n",
    "For example, suppose we want to estimate the coefficient of variation (standard deviation as a fraction of the mean) for adult male height.  Here's how we can compute it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv = np.nanstd(data) / np.nanmean(data)\n",
    "cv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the standard deviation is about 5% of the mean. \n",
    "\n",
    "Write a loop that uses resampling to estimate the sampling distribution of `cv`; store the results in an array named `outcomes`.\n",
    "\n",
    "Then use `plot_sampling_dist` to plot the sampling distribution of the coefficient of variation.\n",
    "\n",
    "What is the standard error of the estimated coefficient of variation?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "This chapter presents computational methods for computing p-values, standard errors, and confidence intervals.  The two processes are similar, but they answer different questions.\n",
    "\n",
    "The following diagram outlines the hypothesis testing process:\n",
    "\n",
    "<img src=\"https://github.com/AllenDowney/ElementsOfDataScience/raw/master/figs/hypothesis_testing.png\" width=\"400\">\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, the key steps are\n",
    "\n",
    "1. Choose a test statistic that quantifies the observed effect.\n",
    "\n",
    "2. Define a model of the null hypothesis and use it to generate simulated data.\n",
    "\n",
    "3. Compute the distribution of the test statistic under the null hypothesis.\n",
    "\n",
    "4. Compute a p-value, which is probability, under the null hypothesis, of seeing an effect as extreme as what you saw."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following figure shows the similar process for computing standard errors and confidence intervals.\n",
    "\n",
    "<img src=\"https://github.com/AllenDowney/ElementsOfDataScience/raw/master/figs/resampling.png\" width=\"400\">\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The essential steps are:\n",
    "\n",
    "1. Choose a sample statistic that quantifies the thing you want to estimate.\n",
    "\n",
    "2. Use the data to make a model of the population, assuming that the estimate is accurate.\n",
    "\n",
    "3. Use the model to simulate the sampling process and generate simulated data.\n",
    "\n",
    "4. Compute the sampling distribution of the estimate.\n",
    "\n",
    "5. Use the sampling distribution to compute the standard error, confidence interval, or both.\n",
    "\n",
    "Finally, remember that both processes only account for variability due to random sampling.  They don't tell us anything about systematic errors in the sampling process, measurement error, or other sources of error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}