08_distributions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Distributions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "[Click here to run this notebook on Colab](https://colab.research.google.com/github/AllenDowney/ElementsOfDataScience/blob/master/08_distributions.ipynb) or\n",
    "[click here to download it](https://github.com/AllenDowney/ElementsOfDataScience/raw/master/08_distributions.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this chapter we'll see three ways to describe a set of values:\n",
    "\n",
    "* A probability mass function (PMF), which represents a set of values and the number of times each one appears in a dataset.\n",
    "\n",
    "* A cumulative distribution function (CDF), which contains the same information as a PMF in a form that makes it easier to visualize, make comparisons, and perform some computations.\n",
    "\n",
    "* A kernel density estimate (KDE), which is like a smooth, continuous version of a histogram.\n",
    "\n",
    "For examples, we'll use data from the General Social Survey (GSS) to look at distributions of age and income, and to explore the relationship between income and education.\n",
    "\n",
    "But we'll start with one of the most important ideas in statistics, the distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Distributions\n",
    "\n",
    "A distribution is a set of values and their corresponding probabilities.  For example, if you roll a six-sided die, there are six possible outcomes, the numbers `1` through `6`, and they all have the same probability, `1/6`.\n",
    "\n",
    "We can represent this distribution of outcomes with a table, like this:\n",
    "\n",
    "\n",
    "| Value   |      Probability |\n",
    "| -----   | -----    |\n",
    "| 1       |      1/6 |\n",
    "| 2       |      1/6 |\n",
    "| 3       |      1/6 |\n",
    "| 4       |      1/6 |\n",
    "| 5       |      1/6 |\n",
    "| 6       |      1/6 |\n",
    "\n",
    "More generally, a distribution can have any number of values, the values can be any type, and the probabilities do not have to be equal.\n",
    "\n",
    "To represent distributions in Python, we will use a library called `empiricaldist`, for \"empirical distribution\", where \"empirical\" means it is based on data rather than a mathematical formula."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "The following cell checks whether `empiricaldist` is installed and installs it if necessary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    import empiricaldist\n",
    "except ImportError:\n",
    "    !pip install empiricaldist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`empiricaldist` provides an object called `Pmf`, which stands for \"probability mass function\".\n",
    "A `Pmf` object contains a set of possible outcomes and their probabilities.\n",
    "\n",
    "For example, here's a `Pmf` that represents the outcome of rolling a six-sided die:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from empiricaldist import Pmf\n",
    "\n",
    "outcomes = [1,2,3,4,5,6]\n",
    "die = Pmf(1/6, outcomes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first argument is the probability of each outcome; the second argument is the list of outcomes.\n",
    "We can display the result like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "die"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A `Pmf` object is a specialized version of a Pandas `Series`, so it provides all of the attributes and methods of a `Series`, plus some additional methods we'll see soon."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The General Social Survey\n",
    "\n",
    "Now we'll use `Pmf` objects to represent distributions of values from a new dataset, the General Social Survey (GSS).\n",
    "The GSS surveys a representative sample of adult residents of the U.S. and asks questions about demographics, personal history, and beliefs about social and political issues.\n",
    "It is widely used by politicians, policy makers, and researchers.\n",
    "\n",
    "The GSS dataset contains hundreds of columns; using an online tool call [GSS Explorer](https://gssdataexplorer.norc.org/) I've selected just a few and created a subset of the data, called an **extract**.\n",
    "\n",
    "Like the NSFG data we used in the previous chapter, the GSS data is stored in a fixed-width format, described by a Stata data dictionary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "dict_file = 'GSS.dct'\n",
    "data_file = 'GSS.dat.gz'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "The following cells download the data extract."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "download('https://github.com/AllenDowney/' +\n",
    "         'ElementsOfDataScience/raw/master/data/' +\n",
    "         dict_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "download('https://github.com/AllenDowney/' +\n",
    "         'ElementsOfDataScience/raw/master/data/' +\n",
    "         data_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use the `statadict` library to read the data dictionary."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "The following cell checks whether you have `statadict` installed and installs it if necessary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    from statadict import parse_stata_dict\n",
    "except ImportError:\n",
    "    !pip install statadict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statadict import parse_stata_dict\n",
    "\n",
    "stata_dict = parse_stata_dict(dict_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data file is compressed, but we can use the `gzip` library to open it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gzip\n",
    "\n",
    "fp = gzip.open(data_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is an object that behaves like a file, so we can pass it as an argument to `read_fwf`, which reads a fixed-width file, as we saw in the previous chapter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "gss = pd.read_fwf(fp, \n",
    "                  names=stata_dict.names, \n",
    "                  colspecs=stata_dict.colspecs)\n",
    "gss.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is a `DataFrame` with one row for each respondent, and one column for each variable.\n",
    "Here are the first few rows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "gss.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I'll explain these variables as we go along, but if you want more information, you can read the online documentation at <https://gssdataexplorer.norc.org/variables/vfilter>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Distribution of Education\n",
    "\n",
    "To get started with this dataset, let's look at the distribution of `EDUC`, which records the number of years of education for each respondent.  First we'll select a column from the `DataFrame` and use `value_counts` to see what values are in it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "gss['EDUC'].value_counts().sort_index()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result from `value_counts` is a set of possible values and the number of times each one appears, so it is a kind of distribution.\n",
    "\n",
    "The values `98` and `99` are special codes for \"Don't know\" and \"No answer\".\n",
    "We'll use `replace` to replace these codes with `NaN`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "educ = gss['EDUC'].replace([98, 99], np.nan)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've already seen one way to visualize a distribution, a histogram.  Here's the histogram of education level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "educ.hist(grid=False)\n",
    "plt.xlabel('Years of education')\n",
    "plt.ylabel('Number of respondents')\n",
    "plt.title('Histogram of education level');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the histogram, we can see the general shape of the distribution and the central tendency -- it looks like the peak is near 12 years of education.\n",
    "But a histogram is not the best way to visualize this distribution because it obscures some important details.\n",
    "\n",
    "An alternative is to use a `Pmf`.\n",
    "The function `Pmf.from_seq` takes any kind of sequence -- like a list, tuple, or Pandas `Series` -- and computes the distribution of the values in the sequence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmf_educ = Pmf.from_seq(educ, normalize=False)\n",
    "type(pmf_educ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The keyword argument `normalize=False` indicates that we don't want to normalize this PMF.  I'll explain what that means soon.\n",
    "\n",
    "Here's what the first few rows look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmf_educ.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this dataset, there are `165` respondents who report that they have had no formal education, and `47` who have only one year.\n",
    "Here the last few rows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmf_educ.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are `1439` respondents who report that they have 20 or more years of formal education, which probably means they attended college and graduate school.\n",
    "\n",
    "You can use the bracket operator to look up a value in a `Pmf` and get the corresponding count:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmf_educ[20]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Usually when we make a PMF, we want to know the *fraction* of respondents with each value, rather than the counts.  We can do that by setting `normalize=True`; then we get a **normalized** PMF, that is, a PMF where the values in the second column add up to 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmf_educ_norm = Pmf.from_seq(educ, normalize=True)\n",
    "pmf_educ_norm.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now if we use the bracket operator to look up a value, the result is a fraction rather than a count.\n",
    "For example, the fraction of people with 12 years of education is about 30%:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmf_educ_norm[12]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Pmf` provides a `bar` method that plots the values and their probabilities as a bar chart."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmf_educ_norm.bar(label='EDUC')\n",
    "\n",
    "plt.xlabel('Years of education')\n",
    "plt.xticks(range(0, 21, 4))\n",
    "plt.ylabel('PMF')\n",
    "plt.title('Distribution of years of education')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this figure, we can see that the most common value is 12 years, but there are also peaks at 14 and 16, which correspond to two and four years of college.\n",
    "\n",
    "For this data, the PMF is probably a better choice than the histogram.  The PMF shows all unique values, so we can see where the peaks are.\n",
    "Because the histogram puts values into bins, it obscures these details.  With this dataset, and the default number of bins, we couldn't see the peaks at 14 and 16 years.\n",
    "But PMFs have limitations, too, as we'll see.\n",
    "\n",
    "First, here's an exercise where you can practice with PMFs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Let's look at the `YEAR` column in the `DataFrame`, which represents the year each respondent was interviewed.\n",
    "\n",
    "Make an unnormalized `Pmf` for `YEAR` and display the result.\n",
    "How many respondents were interviewed in 2018?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cumulative Distribution Functions\n",
    "\n",
    "Now we'll see another way to represent a distribution, the cumulative distribution function (CDF).\n",
    "`empiricaldist` provides a `Cdf` object that represents a CDF.\n",
    "We can import it like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "from empiricaldist import Cdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example, suppose we have a sequence of five values: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "values = 1, 2, 2, 3, 5  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the `Pmf` of these values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "Pmf.from_seq(values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you draw a random value from `values`, the `Pmf` tells you the chance of getting `x`, for any value of `x`.\n",
    "So the probability of the value `1` is `1/5`; the probability of the value `2` is `2/5`; and the probabilities for `3` and `5` are `1/5` each."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A CDF is similar to a PMF in the sense that it contains values and their probabilities; the difference is that the probabilities in the CDF are the cumulative sum of the probabilities in the PMF.\n",
    "\n",
    "Here's a `Cdf` object for the same five values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "Cdf.from_seq(values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you draw a random value from `values`, `Cdf` tells you the chance of getting a value *less than or equal to* `x`, for any given `x`.\n",
    "\n",
    "So the `Cdf` of `1` is `1/5` because one of the five values in the sequence is less than or equal to 1.\n",
    "\n",
    "The `Cdf` of 2 is `3/5` because three of the five values are less than or equal to 2.\n",
    "\n",
    "And the `Cdf` of 5 is `5/5` because all of the values are less than or equal to 5."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CDF of Age\n",
    "\n",
    "Now let's look at a more substantial `Cdf`, the distribution of ages for respondents in the General Social Survey.\n",
    "\n",
    "The variable we'll use is `'AGE'`.\n",
    "According to the codebook, the range of the values is from `18` to `89`, where `89` means \"89 or older\".\n",
    "The special codes `98` and `99` mean \"Don't know\" and \"Didn't answer\".\n",
    "See <https://gssdataexplorer.norc.org/variables/53/vshow>.\n",
    "\n",
    "We can use `replace` to replace the special codes with `NaN`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "age = gss['AGE'].replace([98, 99], np.nan)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute the `Cdf` of these values like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "cdf_age = Cdf.from_seq(age)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Cdf` provides a method called `plot` that plots the CDF as a line.  Here's what it looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "cdf_age.plot()\n",
    "\n",
    "plt.xlabel('Age (years)')\n",
    "plt.ylabel('CDF')\n",
    "plt.title('Distribution of age');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $x$-axis is the ages, from 18 to 89.  The $y$-axis is the cumulative probabilities, from 0 to 1.\n",
    "\n",
    "`cdf_age` can be used as a function, so if you give it an age, it returns the corresponding probability (in a NumPy array)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = 51\n",
    "p = cdf_age(q)\n",
    "p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`q` stands for \"quantity\", which is what we are looking up.  `p` stands for probability, which is the result.\n",
    "In this example, the quantity is age 51, and the corresponding probability is about `0.63`.\n",
    "That means that about 63% of the respondents are 51 years old or younger.\n",
    "\n",
    "The arrow in the following figure shows how you could read this value from the CDF, at least approximately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def draw_line(p, q, x):\n",
    "    xs = [q, q, x]\n",
    "    ys = [0, p, p]\n",
    "    plt.plot(xs, ys, ':', color='gray')\n",
    "\n",
    "def draw_arrow_left(p, q, x):\n",
    "    dx = 3\n",
    "    dy = 0.025\n",
    "    xs = [x+dx, x, x+dx]\n",
    "    ys = [p-dy, p, p+dy]\n",
    "    plt.plot(xs, ys, ':', color='gray')\n",
    "    \n",
    "def draw_arrow_down(p, q, y):\n",
    "    dx = 1.25\n",
    "    dy = 0.045\n",
    "    xs = [q-dx, q, q+dx]\n",
    "    ys = [y+dy, y, y+dy]\n",
    "    plt.plot(xs, ys, ':', color='gray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cdf_age.plot()\n",
    "\n",
    "x = 17\n",
    "draw_line(p, q, x)\n",
    "draw_arrow_left(p, q, x)\n",
    "\n",
    "plt.xlabel('Age (years)')\n",
    "plt.xlim(x-1, 91)\n",
    "plt.ylabel('CDF')\n",
    "plt.title('Distribution of age');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The CDF is an invertible function, which means that if you have a probability, `p`, you can look up the corresponding quantity, `q`.\n",
    "`Cdf` provides a method called `inverse` that computes the inverse of the cumulative distribution function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1 = 0.25\n",
    "q1 = cdf_age.inverse(p1)\n",
    "q1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we look up the probability `0.25` and the result is `31`.  \n",
    "That means that 25% of the respondents are age 31 or less.  Another way to say the same thing is \"age 31 is the 25th percentile of this distribution\".\n",
    "\n",
    "If we look up probability `0.75`, it returns `59`, so 75% of the respondents are 59 or younger."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "p2 = 0.75\n",
    "q2 = cdf_age.inverse(p2)\n",
    "q2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following figure, the arrows show how you could read these values from the CDF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "cdf_age.plot()\n",
    "\n",
    "x = 17\n",
    "draw_line(p1, q1, x)\n",
    "draw_arrow_down(p1, q1, 0)\n",
    "\n",
    "draw_line(p2, q2, x)\n",
    "draw_arrow_down(p2, q2, 0)\n",
    "\n",
    "plt.xlabel('Age (years)')\n",
    "plt.xlim(x-1, 91)\n",
    "plt.ylabel('CDF')\n",
    "plt.title('Distribution of age');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The distance from the 25th to the 75th percentile is called the **interquartile range**, or IQR.  It measures the spread of the distribution, so it is similar to standard deviation or variance.\n",
    "\n",
    "Because it is based on percentiles, it doesn't get thrown off by extreme values or outliers, the way standard deviation does.  So IQR is more **robust** than variance, which means it works well even if there are errors in the data or extreme values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  Using `cdf_age`, compute the fraction of the respondents in the GSS dataset that are *older* than 65."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** The distribution of income in almost every country is long-tailed, which means there are a small number of people with very high incomes.\n",
    "In the GSS dataset, the column `REALINC` represents total household income, converted to 1986 dollars.  We can get a sense of the shape of this distribution by plotting the CDF.\n",
    "\n",
    "Select `REALINC` from the `gss` dataset, make a `Cdf` called `cdf_income`, and plot it.  Remember to label the axes!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing Distributions\n",
    "\n",
    "So far we've seen two ways to represent distributions, PMFs and CDFs.  Now we'll use PMFs and CDFs to compare distributions, and we'll see the pros and cons of each.\n",
    "\n",
    "One way to compare distributions is to plot multiple PMFs on the same axes.  For example, suppose we want to compare the distribution of age for male and female respondents.\n",
    "\n",
    "First we'll create a Boolean Series that's true for male respondents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "male = (gss['SEX'] == 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And another that's true for female respondents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "female = (gss['SEX'] == 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can select ages for the male and female respondents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "male_age = age[male]\n",
    "female_age = age[female]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And plot a PMF for each."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmf_male_age = Pmf.from_seq(male_age)\n",
    "pmf_male_age.plot(label='Male')\n",
    "\n",
    "pmf_female_age = Pmf.from_seq(female_age)\n",
    "pmf_female_age.plot(label='Female')\n",
    "\n",
    "plt.xlabel('Age (years)') \n",
    "plt.ylabel('PMF')\n",
    "plt.title('Distribution of age by sex')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A plot like this, which is highly variable, is often described as *noisy*.\n",
    "\n",
    "If we ignore the noise, it looks like the PMF is higher for men between ages 40 and 50, \n",
    "and higher for women between ages 70 and 80.\n",
    "But both of those differences might be due to random variation.\n",
    "\n",
    "Now let's do the same thing with CDFs; everything is the same except we replace `Pmf` with `Cdf`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "cdf_male_age = Cdf.from_seq(male_age)\n",
    "cdf_male_age.plot(label='Male')\n",
    "\n",
    "cdf_female_age = Cdf.from_seq(female_age)\n",
    "cdf_female_age.plot(label='Female')\n",
    "\n",
    "plt.xlabel('Age (years)') \n",
    "plt.ylabel('CDF')\n",
    "plt.title('Distribution of age by sex')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, CDFs are smoother than PMFs.  Because they smooth out randomness, we can often get a better view of real differences between distributions.\n",
    "In this case, the lines are close together until age 40; after that, the CDF is higher for men than women.\n",
    "So what does that mean?\n",
    "\n",
    "One way to interpret the difference is that the fraction of men below a given age is generally more than the fraction of women below the same age.\n",
    "For example, about 79% of men are 60 or less, compared to 76% of women."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "cdf_male_age(60), cdf_female_age(60)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Going the other way, we could also compare percentiles.\n",
    "For example, the median age woman is older than the median age man, by about one year."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "cdf_male_age.inverse(0.5), cdf_female_age.inverse(0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** What fraction of men are over 80?  What fraction of women?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing Incomes\n",
    "\n",
    "As another example, let's look at household income and compare the distribution before and after 1995 (I chose 1995 because it's roughly the midpoint of the survey).\n",
    "The variable `REALINC` represents household income in 1986 dollars.\n",
    "\n",
    "We'll make a Boolean `Series` to select respondents interviewed before and after 1995."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "pre95 = (gss['YEAR'] < 1995)\n",
    "post95 = (gss['YEAR'] >= 1995)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can plot the PMFs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "income = gss['REALINC'].replace(0, np.nan)\n",
    "\n",
    "Pmf.from_seq(income[pre95]).plot(label='Before 1995')\n",
    "Pmf.from_seq(income[post95]).plot(label='After 1995')\n",
    "\n",
    "plt.xlabel('Income (1986 USD)')\n",
    "plt.ylabel('PMF')\n",
    "plt.title('Distribution of income')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a lot of unique values in this distribution, and none of them appear very often.  As a result, the PMF is so noisy and we can't really see the shape of the distribution.\n",
    "\n",
    "It's also hard to compare the distributions.\n",
    "It looks like there are more people with high incomes after 1995, but it's hard to tell.  We can get a clearer picture with a CDF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "Cdf.from_seq(income[pre95]).plot(label='Before 1995')\n",
    "Cdf.from_seq(income[post95]).plot(label='After 1995')\n",
    "\n",
    "plt.xlabel('Income (1986 USD)')\n",
    "plt.ylabel('CDF')\n",
    "plt.title('Distribution of income')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below $30,000 the CDFs are almost identical; above that, we can see that the post-1995 distribution is shifted to the right.\n",
    "In other words, the fraction of people with high incomes is about the same, but the income of high earners has increased.\n",
    "\n",
    "In general, I recommend CDFs for exploratory analysis.  They give you a clear view of the distribution, without too much noise, and they are good for comparing distributions, especially if you have more than two."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** In the previous figure, the dollar amounts are big enough that the labels on the `x` axis are crowded.  Improve the figure by expressing income in 1000s of dollars (and update the `x` label accordingly)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Let's compare incomes for different levels of education in the GSS dataset\n",
    "\n",
    "To do that we'll create Boolean `Series` to identify respondents with different levels of education.\n",
    "\n",
    "* In the U.S, 12 years of education usually means the respondent has completed high school (secondary education).  \n",
    "\n",
    "* A respondent with 14 years of education has probably completed an associate degree (two years of college)\n",
    "\n",
    "* Someone with 16 years has probably completed a bachelor's degree (four years of college or university).\n",
    "\n",
    "Define Boolean `Series` named `high`, `assc`, and `bach` that are true for respondents with\n",
    "\n",
    "* 12 or fewer years of education,\n",
    "\n",
    "* 13, 14, or 15 years, and\n",
    "\n",
    "* 16 or more.\n",
    "\n",
    "Compute and plot the distribution of income for each group.  Remember to label the CDFs, display a legend, and label the axes.\n",
    "Write a few sentences that describe and interpret the results."
   ]
  },
  {
   "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": [
    "## Modeling Distributions\n",
    "\n",
    "Some distributions have names.  For example, you might be familiar with the normal distribution, also called the Gaussian distribution or the bell curve.  And you might have heard of others like the exponential distribution, binomial distribution, or maybe Poisson distribution.\n",
    "\n",
    "These \"distributions with names\" are called **analytic** because they are described by analytic mathematical functions, as contrasted with empirical distributions, which are based on data.\n",
    "\n",
    "It turns out that many things we measure in the world have distributions that are well approximated by analytic distributions, so these distributions are sometimes good models for the real world.  \n",
    "In this context, what I mean by a \"model\" is a simplified description of the world that is accurate enough for its intended purpose.\n",
    "\n",
    "In this section, we'll compute the CDF of a normal distribution and compare it to an empirical distribution of data.\n",
    "But before we get to real data, we'll start with fake data.\n",
    "\n",
    "The following statement uses NumPy's `random` library to generate 1000 values from a normal distribution with mean `0` and standard deviation `1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "np.random.seed(17)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample = np.random.normal(size=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the empirical distribution of the sample looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "cdf_sample = Cdf.from_seq(sample)\n",
    "cdf_sample.plot(label='Random sample')\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('CDF')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we did not know that this sample was drawn from a normal distribution, and we wanted to check, we could compare the CDF of the data to the CDF of an ideal normal distribution, which we can use the SciPy library to compute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "xs = np.linspace(-3, 3)\n",
    "ys = norm(0, 1).cdf(xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we import `norm` from `scipy.stats`, which is a collection of functions related to statistics.\n",
    "\n",
    "Then we use `linspace()` to create an array of equally-spaced points from -3 to 3; those are the `x` values where we will evaluate the normal CDF.\n",
    "\n",
    "Next, `norm(0, 1)` creates an object that represents a normal distribution with mean `0` and standard deviation `1`.\n",
    "\n",
    "Finally, `cdf` computes the CDF of the normal distribution, evaluated at each of the `xs`.\n",
    "\n",
    "I'll plot the normal CDF with a gray line and then plot the CDF of the data again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(xs, ys, color='gray', label='Normal CDF')\n",
    "cdf_sample.plot(label='Random sample')\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('CDF')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The CDF of the random sample agrees with the normal model.\n",
    "And that's not surprising because the data were actually sampled from a normal distribution.\n",
    "When we collect data in the real world, we do not expect it to fit a normal distribution as well as this.  In the next exercise, we'll try it and see. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Is the normal distribution a good model for the distribution of ages in the U.S. population?\n",
    "\n",
    "To answer this question:\n",
    "\n",
    "* Compute the mean and standard deviation of ages in the GSS dataset.\n",
    "\n",
    "* Use `linspace` to create an array of equally spaced values between 18 and 89.\n",
    "\n",
    "* Use `norm` to create a normal distribution with the same mean and standard deviation as the data, then use it to compute the normal CDF for each value in the array.\n",
    "\n",
    "* Plot the normal CDF with a gray line.\n",
    "\n",
    "* Plot the CDF of the ages in the GSS.\n",
    "\n",
    "How well do the plotted CDFs agree?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** In many datasets, the distribution of income is approximately **lognormal**, which means that the logarithms of the incomes fit a normal distribution.  Let's see whether that's true for the GSS data. \n",
    "\n",
    "* Extract `REALINC` from `gss` and compute its logarithm using `np.log10()`.  Hint: Replace the value `0` with `NaN` before computing logarithms.\n",
    "\n",
    "* Compute the mean and standard deviation of the log-transformed incomes.\n",
    "\n",
    "* Use `norm` to make a normal distribution with the same mean and standard deviation as the log-transformed incomes.\n",
    "\n",
    "* Plot the CDF of the normal distribution.\n",
    "\n",
    "* Compute and plot the CDF of the log-transformed incomes.\n",
    "\n",
    "How similar are the CDFs of the log-transformed incomes and the normal distribution?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Kernel Density Estimation\n",
    "\n",
    "We have seen two ways to represent distributions, PMFs and CDFs.  Now we'll learn another way: a probability density function, or PDF.\n",
    "The `norm` function, which we used to compute the normal CDF, can also compute the normal PDF:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs = np.linspace(-3, 3)\n",
    "ys = norm(0,1).pdf(xs)\n",
    "plt.plot(xs, ys, color='gray', label='Normal PDF')\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('PDF')\n",
    "plt.title('Normal density function')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The normal PDF is the classic \"bell curve\".  \n",
    "\n",
    "It is tempting to compare the PMF of the data to the PDF of the normal distribution, but that doesn't work.\n",
    "Let's see what happens if we try:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(xs, ys, color='gray', label='Normal PDF')\n",
    "\n",
    "pmf_sample = Pmf.from_seq(sample)\n",
    "pmf_sample.plot(label='Random sample')\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('PDF')\n",
    "plt.title('Normal density function')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The PMF of the sample is a flat line across the bottom.  In the random sample, every value is unique, so they all have the same probability, one in 1000.\n",
    "\n",
    "However, we can use the points in the sample to estimate the PDF of the distribution they came from.\n",
    "This process is called **kernel density estimation**, or KDE.  It's a way of getting from a PMF, a probability mass function, to a PDF, a probability density function.\n",
    "\n",
    "To generate a KDE plot, we'll use the Seaborn library, imported as `sns`.\n",
    "Seaborn provides `kdeplot`, which takes the sample, estimates the PDF, and plots it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "\n",
    "sns.kdeplot(sample, label='Estimated sample PDF')\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('PDF')\n",
    "plt.title('Normal density function')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can compare the KDE plot and the normal PDF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(xs, ys, color='gray', label='Normal PDF')\n",
    "sns.kdeplot(sample, label='Estimated sample PDF')\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('PDF')\n",
    "plt.title('Normal density function')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The KDE plot matches the normal PDF pretty well, although the differences look bigger when we compare PDFs than they did with the CDFs.\n",
    "That means that the PDF is a more sensitive way to look for differences, but often it is too sensitive.  \n",
    "It's hard to tell whether apparent differences mean anything, or if they are just random, as in this case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** In a previous exercise, we asked \"Is the normal distribution a good model for the distribution of ages in the U.S. population?\"  To answer this question, we plotted the CDF of the data and compared it to the CDF of a normal distribution with the same mean and standard deviation.\n",
    "\n",
    "Now we'll compare the estimated density of the data with the normal PDF.\n",
    "\n",
    "* Again, compute the mean and standard deviation of ages in the GSS dataset.\n",
    "\n",
    "* Use `linspace` to create an array of values between 18 and 89.\n",
    "\n",
    "* Use `norm` to create a normal distribution with the same mean and standard deviation as the data, then use it to compute the normal PDF for each value in the array.\n",
    "\n",
    "* Plot the normal PDF with a gray line.\n",
    "\n",
    "* Use `sns.kdeplot` to estimate and plot the density of the ages in the GSS.  \n",
    "\n",
    "Note: Seaborn can't handle NaNs, but you can use `dropna` to drop them before calling `kdeplot`.\n",
    "\n",
    "How well do the PDF and KDE plots agree?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** In a previous exercise, we used CDFs to see if the distribution of income fits a lognormal distribution.  We can make the same comparison using a PDF and KDE.\n",
    "\n",
    "* Again, extract `REALINC` from `gss` and compute its logarithm using `np.log10()`.\n",
    "\n",
    "* Compute the mean and standard deviation of the log-transformed incomes.\n",
    "\n",
    "* Use `norm` to make a normal distribution with the same mean and standard deviation as the log-transformed incomes.\n",
    "\n",
    "* Plot the PDF of the normal distribution.\n",
    "\n",
    "* Use `sns.kdeplot()` to estimate and plot the density of the log-transformed incomes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "In this chapter, we've seen four ways to visualize distributions, PMFs, CDFs, and KDE plots.\n",
    "\n",
    "In general, I use CDFs when I am exploring data.\n",
    "That way, I get the best view of what's going on without getting distracted by noise.\n",
    "\n",
    "Then, if I am presenting results to an audience unfamiliar with CDFs, I might use a PMF if the dataset contains a small number of unique values, or KDE if there are many unique values."
   ]
  },
  {
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}