{
"cells": [
{
"cell_type": "markdown",
"id": "259c89bc",
"metadata": {},
"source": [
"# Probability Mass Functions\n",
"\n",
"In the previous chapter we represented distributions using a `Hist` object, which contains a set of quantities and their frequencies -- that is, the number of times each one appears.\n",
"In this chapter we'll introduce another way to represent a distribution: a `Pmf` object, which contains a set of quantities and their probabilities.\n",
"\n",
"And we'll use `Pmf` objects to compute the mean and variance of a distribution, and the skewness, which indicates whether it is skewed to the left or right.\n",
"Finally, we will explore how a phenomenon called the \"inspection paradox\" can cause a sample to give a biased view of a distribution."
]
},
{
"cell_type": "markdown",
"id": "11b3a5ef",
"metadata": {
"tags": [
"remove-print"
]
},
"source": [
"[Click here to run this notebook on Colab](https://colab.research.google.com/github/AllenDowney/ThinkStats/blob/v3/nb/chap03.ipynb)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fd972111",
"metadata": {
"tags": [
"remove-print",
"hide-cell"
]
},
"outputs": [],
"source": [
"%load_ext nb_black\n",
"%load_ext autoreload\n",
"%autoreload 2"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "c4578587",
"metadata": {
"tags": [
"remove-print",
"hide-cell"
]
},
"outputs": [],
"source": [
"from os.path import basename, exists\n",
"\n",
"\n",
"def download(url):\n",
" filename = basename(url)\n",
" if not exists(filename):\n",
" from urllib.request import urlretrieve\n",
"\n",
" local, _ = urlretrieve(url, filename)\n",
" print(\"Downloaded \" + local)\n",
"\n",
"\n",
"download(\"https://github.com/AllenDowney/ThinkStats/raw/v3/nb/thinkstats.py\")"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "f66668e0",
"metadata": {
"tags": [
"remove-print",
"hide-cell"
]
},
"outputs": [],
"source": [
"try:\n",
" import empiricaldist\n",
"except ImportError:\n",
" !pip install empiricaldist"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "526c4c90",
"metadata": {},
"outputs": [],
"source": [
"try:\n",
" import statadict\n",
"except ImportError:\n",
" !pip install statadict"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "e47f1afd",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from thinkstats import decorate"
]
},
{
"cell_type": "markdown",
"id": "5f56b6a7",
"metadata": {},
"source": [
"## PMFs\n",
"\n",
"A `Pmf` object is like a `Hist` that contains probabilities instead of frequencies.\n",
"So one way to make a `Pmf` is to start with a `Hist`.\n",
"For example, here's a `Hist` that represents the distribution of values in a short sequence."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "8f418fed",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
freqs
\n",
"
\n",
" \n",
" \n",
"
\n",
"
1
\n",
"
1
\n",
"
\n",
"
\n",
"
2
\n",
"
2
\n",
"
\n",
"
\n",
"
3
\n",
"
1
\n",
"
\n",
"
\n",
"
5
\n",
"
1
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"1 1\n",
"2 2\n",
"3 1\n",
"5 1\n",
"Name: , dtype: int64"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from empiricaldist import Hist\n",
"\n",
"hist = Hist.from_seq([1, 2, 2, 3, 5])\n",
"hist"
]
},
{
"cell_type": "markdown",
"id": "f70e8df4",
"metadata": {},
"source": [
"The sum of the frequencies is the size of the original sequence."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "8e1ee1ca",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = hist.sum()\n",
"n"
]
},
{
"cell_type": "markdown",
"id": "eeed5af8",
"metadata": {},
"source": [
"If we divide the frequencies by `n`, they represent proportions, rather than counts."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "07dab3d2",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
probs
\n",
"
\n",
" \n",
" \n",
"
\n",
"
1
\n",
"
0.2
\n",
"
\n",
"
\n",
"
2
\n",
"
0.4
\n",
"
\n",
"
\n",
"
3
\n",
"
0.2
\n",
"
\n",
"
\n",
"
5
\n",
"
0.2
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"1 0.2\n",
"2 0.4\n",
"3 0.2\n",
"5 0.2\n",
"Name: , dtype: float64"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pmf = hist / n\n",
"pmf"
]
},
{
"cell_type": "markdown",
"id": "1e0a7d3c",
"metadata": {},
"source": [
"This result indicates that 20% of the values in the sequence are 1, 40% are 2, and so on.\n",
"\n",
"We can also think of these proportions as probabilities in the following sense: if we choose a random value from the original sequence, the probability we choose the value 1 is 0.2, the probability we choose the value 2 is 0.4, and so on.\n",
"\n",
"Because we divided through by `n`, the sum of the probabilities is 1, which means that this distribution is **normalized**."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "af296dc9",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.0"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pmf.sum()"
]
},
{
"cell_type": "markdown",
"id": "c81d485a",
"metadata": {},
"source": [
"A normalized `Hist` object represents a **probability mass function** (PMF), so-called because probabilities associated with discrete quantities are also called \"probability masses\".\n",
"\n",
"The `empiricaldist` library provides a `Pmf` object that represents a probability mass function, so instead of creating a `Hist` object and then normalizing it, we can create a `Pmf` object directly."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "3fff30ec",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
probs
\n",
"
\n",
" \n",
" \n",
"
\n",
"
1
\n",
"
0.2
\n",
"
\n",
"
\n",
"
2
\n",
"
0.4
\n",
"
\n",
"
\n",
"
3
\n",
"
0.2
\n",
"
\n",
"
\n",
"
5
\n",
"
0.2
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"1 0.2\n",
"2 0.4\n",
"3 0.2\n",
"5 0.2\n",
"Name: , dtype: float64"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from empiricaldist import Pmf\n",
"\n",
"pmf = Pmf.from_seq([1, 2, 2, 3, 5])\n",
"pmf"
]
},
{
"cell_type": "markdown",
"id": "d8428fe8",
"metadata": {},
"source": [
"The `Pmf` is normalized so the total probability is 1."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "5e01b656",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.0"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pmf.sum()"
]
},
{
"cell_type": "markdown",
"id": "4cc44f90",
"metadata": {},
"source": [
"`Pmf` and `Hist` objects are similar in many ways.\n",
"To look up the probability associated with a value, we can use the bracket operator."
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "41280992",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.4"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pmf[2]"
]
},
{
"cell_type": "markdown",
"id": "ff05014f",
"metadata": {},
"source": [
"Or use parentheses to call the `Pmf` like a function."
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "2ac2919b",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.4"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pmf(2)"
]
},
{
"cell_type": "markdown",
"id": "98af6144",
"metadata": {},
"source": [
"To assign a probability to a quantity, you have to use the bracket operator."
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "b53f5742",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.2"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pmf[2] = 0.2\n",
"pmf(2)"
]
},
{
"cell_type": "markdown",
"id": "69db5cfb",
"metadata": {},
"source": [
"You can modify an existing `Pmf` by incrementing the probability associated with a value:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "25ecc3b0",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.5"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pmf[2] += 0.3\n",
"pmf[2]"
]
},
{
"cell_type": "markdown",
"id": "b4de74dd",
"metadata": {},
"source": [
"Or you can multiply a probability by a factor:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "df663bda",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.25"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pmf[2] *= 0.5\n",
"pmf[2]"
]
},
{
"cell_type": "markdown",
"id": "3abe8d9e",
"metadata": {},
"source": [
"If you modify a `Pmf`, the result may not be normalized -- that is, the probabilities may no longer add up to 1."
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "864a09be",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.8500000000000001"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pmf.sum()"
]
},
{
"cell_type": "markdown",
"id": "bce848f1",
"metadata": {},
"source": [
"The `normalize` method renormalizes the `Pmf` by dividing through by the sum -- and returning the sum. "
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "a06845d8",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.8500000000000001"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pmf.normalize()"
]
},
{
"cell_type": "markdown",
"id": "4ae297e6",
"metadata": {},
"source": [
"`Pmf` objects provide a `copy` method so you can make and modify a copy without affecting the original."
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "7738244e",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
probs
\n",
"
\n",
" \n",
" \n",
"
\n",
"
1
\n",
"
0.235294
\n",
"
\n",
"
\n",
"
2
\n",
"
0.294118
\n",
"
\n",
"
\n",
"
3
\n",
"
0.235294
\n",
"
\n",
"
\n",
"
5
\n",
"
0.235294
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"1 0.235294\n",
"2 0.294118\n",
"3 0.235294\n",
"5 0.235294\n",
"Name: , dtype: float64"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pmf.copy()"
]
},
{
"cell_type": "markdown",
"id": "94f5f531",
"metadata": {},
"source": [
"Like a `Hist` object, a `Pmf` object has a `qs` attribute that accesses the quantities and a `ps` attribute that accesses the probabilities.\n",
"\n",
"It also has a `bar` method that plots the `Pmf` as a bar graph and a `plot` method that plots it as a line graph."
]
},
{
"cell_type": "markdown",
"id": "3e77da45",
"metadata": {},
"source": [
"## Summary Statistics\n",
"\n",
"In Section XXX we computed the mean of a sample by adding up the elements and dividing by `n`.\n",
"Here's a simple example."
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "07579612",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.6"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"seq = [1, 2, 2, 3, 5]\n",
"\n",
"mean = np.sum(seq) / len(seq)\n",
"mean"
]
},
{
"cell_type": "markdown",
"id": "7ca46b48",
"metadata": {},
"source": [
"Now suppose we compute the PMF of the values in the sequence."
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "fa1ca2c9",
"metadata": {},
"outputs": [],
"source": [
"pmf = Pmf.from_seq(seq)"
]
},
{
"cell_type": "markdown",
"id": "ea4a034d",
"metadata": {},
"source": [
"Given the `Pmf`, we can still compute the mean, but the process is different -- we have to multiply the probabilities and quantities and add up the products."
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "f73de148",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.6"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mean = np.sum(pmf.ps * pmf.qs)\n",
"mean"
]
},
{
"cell_type": "markdown",
"id": "b7913b14",
"metadata": {},
"source": [
"The `mean` method does the same thing."
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "b7a352d8",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.6"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pmf.mean()"
]
},
{
"cell_type": "markdown",
"id": "444f74bd",
"metadata": {},
"source": [
"Notice that we *don't* have to divide by `n`, because we already did that when we normalized the `Pmf`.\n",
"\n",
"Give a `Pmf`, we can compute the variance by computing the deviation of each quantity from the mean."
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "06ecd356",
"metadata": {},
"outputs": [],
"source": [
"deviations = pmf.qs - mean"
]
},
{
"cell_type": "markdown",
"id": "b4a01857",
"metadata": {},
"source": [
"Then we multiply the squared deviations by the probabilities and add up the products."
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "bd600815",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.84"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"var = np.sum(pmf.ps * deviations**2)\n",
"var"
]
},
{
"cell_type": "markdown",
"id": "f8a9a9b6",
"metadata": {},
"source": [
"The `var` method does the same thing."
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "a220d0e6",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.84"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pmf.var()"
]
},
{
"cell_type": "markdown",
"id": "9a1c2c57",
"metadata": {},
"source": [
"From the variance, we can compute the standard deviation in the usual way."
]
},
{
"cell_type": "code",
"execution_count": 26,
"id": "5e152cda",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.3564659966250536"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.sqrt(var)"
]
},
{
"cell_type": "markdown",
"id": "f2accd22",
"metadata": {},
"source": [
"Or the `std` method does the same thing."
]
},
{
"cell_type": "code",
"execution_count": 27,
"id": "e7557612",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.3564659966250536"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pmf.std()"
]
},
{
"cell_type": "markdown",
"id": "6645531f",
"metadata": {},
"source": [
"`Pmf` also provides a `mode` method that finds the value with the highest probability."
]
},
{
"cell_type": "code",
"execution_count": 28,
"id": "19405109",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pmf.mode()"
]
},
{
"cell_type": "markdown",
"id": "51fc69ca",
"metadata": {},
"source": [
"## The Class Size Paradox\n",
"\n",
"As an example of what we can do with `Pmf` objects, let's consider a phenomenon I call \"the class size paradox.\"\n",
"\n",
"At many American colleges and universities, the student-to-faculty ratio is about 10:1.\n",
"But students are often surprised to discover that their average class size is bigger than 10. There are two reasons for the discrepancy:\n",
"\n",
"- Students typically take 4 or 5 classes per semester, but professors often teach 1 or 2.\n",
"\n",
"- The number of students who enjoy a small class is small, but the number of students in a large class is large.\n",
"\n",
"The first effect is obvious, at least once it is pointed out; the second is more subtle.\n",
"Let's look at an example.\n",
"Suppose that a college offers 65 classes in a given semester, and we are given the number of classes in each of the following size ranges."
]
},
{
"cell_type": "code",
"execution_count": 29,
"id": "36fa48bc",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
count
\n",
"
\n",
"
\n",
"
class size
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
[5, 10)
\n",
"
8
\n",
"
\n",
"
\n",
"
[10, 15)
\n",
"
8
\n",
"
\n",
"
\n",
"
[15, 20)
\n",
"
14
\n",
"
\n",
"
\n",
"
[20, 25)
\n",
"
4
\n",
"
\n",
"
\n",
"
[25, 30)
\n",
"
6
\n",
"
\n",
"
\n",
"
[30, 35)
\n",
"
12
\n",
"
\n",
"
\n",
"
[35, 40)
\n",
"
8
\n",
"
\n",
"
\n",
"
[40, 45)
\n",
"
3
\n",
"
\n",
"
\n",
"
[45, 50)
\n",
"
2
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" count\n",
"class size \n",
"[5, 10) 8\n",
"[10, 15) 8\n",
"[15, 20) 14\n",
"[20, 25) 4\n",
"[25, 30) 6\n",
"[30, 35) 12\n",
"[35, 40) 8\n",
"[40, 45) 3\n",
"[45, 50) 2"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ranges = pd.interval_range(start=5, end=50, freq=5, closed=\"left\")\n",
"ranges.name = \"class size\"\n",
"\n",
"data = pd.DataFrame(index=ranges)\n",
"data[\"count\"] = [8, 8, 14, 4, 6, 12, 8, 3, 2]\n",
"data"
]
},
{
"cell_type": "markdown",
"id": "2e43baa6",
"metadata": {},
"source": [
"Since we don't know the sizes of the classes in each range, let's assume that all sizes are at the midpoint of the range."
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "36f208f9",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Index([7, 12, 17, 22, 27, 32, 37, 42, 47], dtype='int64')"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sizes = ranges.left + 2\n",
"sizes"
]
},
{
"cell_type": "markdown",
"id": "7551d814",
"metadata": {},
"source": [
"Now let's make a `Pmf` that represents the distribution of class sizes.\n",
"Because we already know the frequency of each size, we don't have to use `from_seq`; instead, we can create a `Pmf` directly, passing as arguments the counts, sizes, and name.\n",
"When we normalize the new `Pmf`, the result is the sum of the counts."
]
},
{
"cell_type": "code",
"execution_count": 31,
"id": "c1db28bc",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"65"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"counts = data[\"count\"]\n",
"actual_pmf = Pmf(counts, sizes, name=\"actual\")\n",
"actual_pmf.normalize()"
]
},
{
"cell_type": "markdown",
"id": "ca39aa3a",
"metadata": {},
"source": [
"If you ask the college for the average class size, they report the mean of this distribution, which is 23.7."
]
},
{
"cell_type": "code",
"execution_count": 32,
"id": "6f8eb2c2",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"mean 23.692307692307693\n"
]
}
],
"source": [
"print(\"mean\", actual_pmf.mean())"
]
},
{
"cell_type": "markdown",
"id": "0d22bd69",
"metadata": {},
"source": [
"But if you survey a group of students, ask them how many students are in their classes, and compute the mean, the average is bigger.\n",
"Let's see how much bigger.\n",
"\n",
"The following function takes the actual `Pmf` of class sizes and makes a new `Pmf` that represents the class sizes as seen by students.\n",
"The quantities in the two distributions are the same, but the probabilities in the distribution are multiplied by the quantities, because in a class with size `x`, there are `x` students who observe that class.\n",
"So the probability of observing a class is proportional to its size."
]
},
{
"cell_type": "code",
"execution_count": 33,
"id": "3cef6126",
"metadata": {},
"outputs": [],
"source": [
"def bias(pmf, name):\n",
" qs = pmf.qs\n",
" ps = pmf.ps * pmf.qs\n",
"\n",
" new_pmf = Pmf(ps, pmf.qs, name=name)\n",
" new_pmf.normalize()\n",
" return new_pmf"
]
},
{
"cell_type": "markdown",
"id": "f1449478",
"metadata": {},
"source": [
"Now we can compute the biased `Pmf` as observed by students."
]
},
{
"cell_type": "code",
"execution_count": 34,
"id": "6df3cccb",
"metadata": {},
"outputs": [],
"source": [
"observed_pmf = bias(actual_pmf, name=\"observed\")"
]
},
{
"cell_type": "markdown",
"id": "e96df61b",
"metadata": {},
"source": [
"Here's what the two distributions look like."
]
},
{
"cell_type": "code",
"execution_count": 35,
"id": "1c2b81d9",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAboAAAD/CAYAAACHFRPuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAAuJAAALiQE3ycutAAAfX0lEQVR4nO3df1RUdf4/8OeAQBqEopbiDJEQJsOPUVhbAoXogBWBMU1IelhJRtecXa20vhqe3T5uq51o3SinVQ8UdHBdQdTFtIzdBTQ9J6kkUFYClQNzsB8CMpaGwtzvHx3vNotL4MwFfPt8nDPnzNx5v++85n3CZ+87976vSpIkCURERIJyGe4CiIiIlMSgIyIioTHoiIhIaAw6IiISGoOOiIiExqAjIiKhMeiIiEhoo4a7AGe74447oFarh7sMIiIaQhaLBVar9brvCRd0arUa9fX1w10GERENoeDg4P/5Hg9dEhGR0Bh0REQkNOEOXRIRjTRcUth5VCrVoPsoGnSVlZUwmUzo7u5GXFwctm7dCldXV/n9mpoamEwmXLhwASqVCkuXLsWKFSsAAAUFBVi1ahU0Gg0AIDo6GmazWclyiYic6urVq2htbUV3d/dwlyIMDw8PaDQauLm5DbiPYkFns9lgNBpRVlaG4OBgpKWloaioCIsWLZLbjBkzBu+88w6mTZsGq9WKiIgIzJkzBzqdDgCQmpqKvLw8pUokIlJUa2srvLy84O/vf0MzEbInSRLa29vR2tqKqVOnDrifYr/RVVdXw9fXVz4TJisrC6WlpXZtgoKCMG3aNAA/XhYwffp0tLa2KlUSEdGQkSQJ3d3dGD9+PFxcXKBSqfhw8OHi4oLx48eju7t7UIeDFQs6i8UiH3YEAD8/v35D7PTp0/j0008RHR0tb9u3bx/Cw8ORmJiI6upqpUolIlIMZ3LONaJ+oxtM2l64cAGPP/44cnNz4ePjAwB47LHHMH/+fIwePRpVVVVITU3FqVOn4OnpadfXbDbb/XbX2dnpnC9ARM61b6Vj/ZNznVMH/aw33ngDixcvxh133HHD+4iLi8Mrr7yCmJgYJ1Z2YxQLOo1GYzeDa2lpue6KJZcuXUJSUhKWLFmCJ598Ut4+YcIE+XlsbCzUajUaGhoQERFh199kMsFkMsmv+7tokIhoOK3dXef0fW7Uhzp9n2+88QYMBoNDQTeSKHboMjIyEhaLRV6lJD8/H3q93q7N1atXodfrkZCQIJ9teU1bW5v8/OTJk2hubh7Uj49ERPQfBoMBkZGRCAkJwXPPPQcAOH/+PObPn4+wsDCEhYVhx44dyMnJQVtbG+bOnQudTocrV67A398fFosFwI8/S/n7+wMALl++jISEBERERECr1SInJ2e4vl6/FJvRubq6Ii8vDwaDAd3d3YiNjUVGRgbKyspQVlaGvLw8FBcXo7y8HF999RX27t0LAFi3bh0MBgM2b96Mv//973Bzc4ObmxsKCgowbtw4pcolIhLa1q1bMX78eNhsNuj1ehw8eBDvvfceQkJCsHPnTgBAR0cHfHx8YDabcfDgwZ9dN9jd3R0lJSUYO3Ysrly5gujoaCQnJ+O+++4biq80YIpeRxcfH99n3cmUlBSkpKQAABYuXIiFCxdet++GDRuwYcMGJcsjIrplbNu2DcXFxejt7cU333yDmJgYHDx4EH/5y1/kNtfOkRgoSZKwfv16/POf/4QkSbBYLDhx4sStFXRERDT8qqqqUFpaikOHDsHLywurVq3CDz/8AGBgZzGOGjUKNpsNAOR+ALB9+3acPn0ax44dg4eHB5544gm790cKrnVJRCS4rq4ujB07Fl5eXmhvb5evaZ47dy5yc/9zNmtHRwcAwMvLy+6WN/fccw8+++wzAMCuXbvs9jthwgR4eHjg7NmzKC8vH4qvM2gMOiIiwT388MPw9PTEtGnTYDAYMGfOHABAbm4uamtrERISgvDwcDmonnnmGaSmpsono/zf//0f1q5di4iICHz//ffyfjMyMtDY2AitVovf/OY3iI2NHZbv93NUkmCrjQYHB/N+dEQj0S12HZ0kSTh16hTuu+8+XjTuRP9rXPv7t58zOiIiEhqDjoiIhMagIyIioTHoiIhIaAw6IiISGoOOiIiExqAjIrpFZWZmoqioaLjLuK6XX34Zr7zyilP2xSXAiIiGiqPXEl7PCLy+sLe3F66ursNdhowzOiKiW8CePXsQHh6O0NBQLFiwABcvXgQAfPzxx5g1axbuvfdebNq0CcCP9wlNTU1FWFgYQkJCkJ2dDeDH2/IsW7YMs2bNQmhoKDZv3gwAaG5uxj333IOlS5ciPDwcr776KubOnSt/9ldffYW7774bNpsNHR0dSE9Px6xZs6DT6eQlxSRJwvPPP4+goCDMnj0bjY2NTvvunNEREQnuq6++wvLly/HJJ5/Az88Pv/3tb+XDgidOnMDHH3+MS5cuISIiAgkJCWhsbMSdd96JPXv2AAA6OzsBABs3bsTMmTOxZcsWdHd3Izo6GvHx8RgzZgyam5sxf/58bNu2DTabDVu3bsXXX3+Nu+66Czt37kRaWhpcXFzw7LPPYvHixUhMTMSFCxfwi1/8AvHx8aiqqsJnn32GkydP4vvvv8fMmTMxffp0p3x/zuiIiAT3ySefICYmBn5+fgCArKwsVFRUAACeeuopuLu7Y+zYsUhJScGhQ4cQHh6Of/zjH3j++efxwQcfyHcaP3DgAHJzc6HT6XD//fejo6MDDQ0NAIBJkybhoYceAgC4uLjgySeflO9z99e//hULFiwAAHzwwQd48cUXodPpEBcXhytXruDMmTM4dOgQ5s+fDzc3N4wdOxbz5s1z2vfnjI6ISHD/vdbmT19fbx3OgIAAfP755ygvL0dhYSE2b96M/fv3Q5IkbN++HTqdzq59c3Mzbr/9drttCxYsgMlkQlJSEr777jvMmDEDAGCz2VBZWYmxY8fatd++fbsD37B/nNEREQlu1qxZOHLkCCwWCwDg3XffRXx8PADgb3/7G65evYquri7s27cPc+bMgcVigZubGwwGAzZt2oRjx44BAB555BHk5uait7cXANDY2Gh3O5+fioiIQFdXF/7whz/Is7lr+/jzn/8svz5+/DgkSUJsbCyKi4vlWsrKypz2/Rl0RESCmzRpEjZv3oykpCSEhobim2++kU8wCQ4ORnR0NCIjI7F8+XKEhoairq4OUVFR0Ol0eOSRR/DWW28BANatWwdPT0+Eh4cjJCQES5YswZUrV/7n5y5YsACFhYV2Qffmm2/i9OnTCA0NhVarxdq1ayFJEubNm4eZM2dCq9Xisccewy9/+UunfX/epoeIhgZv00NOwNv0EBER/RcGHRERCY1BR0REQuPlBeR0a3fX3XDfjfpQJ1ZCNPwkSeJvdE50I6eVMOiIiBSgUqng4eGB9vZ2jB8/nmHnBJIkob29HR4eHoMaTwYdEZFCNBoNWltbcf78+eEuRRgeHh7QaDSD6sOgIyJSiJubG6ZOnXpDh9vo+m5kZsygIyJSGA9bDi+edUlEREJj0BERkdAYdEREJDQGHRERCU3RoKusrIRWq0VgYCCMRqN8a4drampqEB0dDa1Wi5CQELz55pvyez09PXj66acRGBgIrVaLw4cPK1kqEREJSrGgs9lsMBqNKCkpQVNTE6xWK4qKiuzajBkzBu+88w5OnjyJo0eP4q233kJNTQ0AoKCgAN3d3WhqasKOHTuQlZXFU3SJiGjQFAu66upq+Pr6Ijg4GMCPt24vLS21axMUFIRp06YBAO644w5Mnz4dra2tAIDS0lIYjUYAQFhYGMaPH4/jx48rVS4REQlKsaCzWCx2V6/7+fnJIXY9p0+fxqefforo6Ogb6k9ERHQ9il0wPpjDjBcuXMDjjz+O3Nxc+Pj4DKq/2WyG2WyWX3d2dg6uUCIiEppiM7pra7xd09LSArVa3afdpUuXkJSUhCVLluDJJ58cdH+TyYT6+nr5MW7cOCd/EyIiupkpFnSRkZGwWCzyrc3z8/Oh1+vt2ly9ehV6vR4JCQlYsWKF3Xt6vR55eXkAgNraWnz77beYOXOmUuUSEZGgFAs6V1dX5OXlwWAwICAgAJ6ensjIyEBZWZl8kklxcTHKy8uxd+9e6HQ66HQ67Nq1CwCQmZkJNzc3BAYGIj09Hfn5+VwvjoiIBk0lCXbOfnBwsDyLpOHBG6/Sde1b6Vj/5Fzn1EFC6u/ffq6MQkREQuNteoiIAMdmnJxtjmic0RERkdAYdEREJDQGHRERCY1BR0REQmPQERGR0HjWJRENiCPXRwLARjcnFUI0SJzRERGR0Bh0REQkNAYdEREJjUFHRERCY9AREZHQGHRERCQ0Bh0REQmNQUdEREJj0BERkdAYdEREJDQGHRERCY1BR0REQmPQERGR0Bh0REQkNAYdEREJjUFHRERCY9AREZHQGHRERCQ0Bh0REQmNQUdEREJj0BERkdAYdEREJDQGHRERCW2UkjuvrKyEyWRCd3c34uLisHXrVri6utq1SUtLQ0VFBby9vdHU1CRvLygowKpVq6DRaAAA0dHRMJvNSpZLI8G+lY71T851Th1EJAzFZnQ2mw1GoxElJSVoamqC1WpFUVFRn3bLli3DwYMHr7uP1NRU1NTUoKamhiFHREQ3RLGgq66uhq+vL4KDgwEAWVlZKC0t7dMuPj4ePj4+SpVBRES3OMWCzmKxyIcdAcDPzw+tra2D2se+ffsQHh6OxMREVFdXO7tEIiK6BSj2G50kSQ71f+yxxzB//nyMHj0aVVVVSE1NxalTp+Dp6WnXzmw22x3W7OzsdOhziYhILIrN6DQajd0MrqWlBWq1esD9J0yYgNGjRwMAYmNjoVar0dDQ0KedyWRCfX29/Bg3bpzjxRMRkTAUC7rIyEhYLBbU19cDAPLz86HX6wfcv62tTX5+8uRJNDc3Y+rUqU6vk4iIxKZY0Lm6uiIvLw8GgwEBAQHw9PRERkYGysrKYDQa5XZJSUmIiopCc3Mz1Go1Nm7cCADYvHkztFotdDodMjMzUVBQwNkaERENmqLX0cXHx8szumtSUlKQkpIiv96/f/91+27YsAEbNmxQsjwiIroFcGUUIiISWr9Bt2PHDvn5+fPnFS+GiIjI2foNupycHPl5YmKi4sUQERE5W79B99Nr4Ry9Lo6IiGg49HsySk9PD86dOwebzSY//2ng+fr6Kl4gERGRI/oNuosXLyI6OloOtwceeEB+T6VS4cyZM8pWR0RE5KB+g665uXmIyiAiIlJGv0HX0tLSb2c/Pz+nFkNERORs/Qadv78//Pz84O7u3udkFJVKhS+//FLR4oiIiBzVb9A999xzKC8vxwMPPICFCxdi9uzZQ1UXERGRU/R7ecGf/vQnfPHFF0hLS0NBQQHCwsKwZs0a/nZHREQ3jZ9d61KlUiE+Ph5RUVHYuXMnVq9eDY1GA5PJNBT1EREROaTfoLPZbDh48CC2b9+Ouro6JCcno6qqClqtdqjqIyIicki/QTd58mRMmTIFCxYswDPPPAOVSoWuri4cPXoUgP11dURERCNRv0E3ffp0qFQq7N+/v8/tdFQqFf71r38pWhwREZGj+g26ysrKISqDiIhIGf0G3dmzZ7Fq1So0NjZCp9Ph9ddfx1133TVUtQ2btbvrHOq/UR/qpEqub6TXRzQcHP67cHNSITTi9Ht5weLFixEUFIScnBx4eXnh2WefHaKyiIiInKPfGd358+fx6quvAgDmzp2LmTNnDklRREREztLvjM7N7T9zeZVKBZVKpXhBREREztTvjK6mpgbu7u7y656eHnndS5VKhStXriheIBERkSN+9oJxIiKim1m/hy6JiIhudgw6IiISGoOOiIiExqAjIiKhMeiIiEhoDDoiIhIag46IiITGoCMiIqH1e8G4oyorK2EymdDd3Y24uDhs3boVrq6udm3S0tJQUVEBb29vNDU1ydt7enqwZMkSHD58GB4eHtiyZQtmz56tZLlEN799K2+8b3Ku8+ogGkEUm9HZbDYYjUaUlJSgqakJVqsVRUVFfdotW7YMBw8e7LO9oKAA3d3daGpqwo4dO5CVlQVJkpQql4iIBKVY0FVXV8PX1xfBwcEAgKysLJSWlvZpFx8fDx8fnz7bS0tLYTQaAQBhYWEYP348jh8/rlS5REQkKMWCzmKxQKPRyK/9/PzQ2to6ZP2JiIgABX+jc/Qw40D7m81mmM1m+XVnZ6dDn0tERGJRbEan0WjsZmAtLS1Qq9VO728ymVBfXy8/xo0b51jhREQkFMWCLjIyEhaLBfX19QCA/Px86PX6AffX6/XIy8sDANTW1uLbb7/lHc6JiGjQFAs6V1dX5OXlwWAwICAgAJ6ensjIyEBZWZl8kgkAJCUlISoqCs3NzVCr1di4cSMAIDMzE25ubggMDER6ejry8/N5h3MiIho0Ra+ji4+Pl2d016SkpCAlJUV+vX///uv2dXNzQ2FhoZLlERHRLYAroxARkdAYdEREJDRFD10S0eCs3V3nUP+Nbk4qhEggnNEREZHQGHRERCQ0Bh0REQmNQUdEREJj0BERkdAYdEREJDQGHRERCY3X0RERjXT7VjrWPznXOXXcpDijIyIioTHoiIhIaAw6IiISGoOOiIiExqAjIiKhMeiIiEhoDDoiIhIag46IiITGoCMiIqEx6IiISGgMOiIiEhqDjoiIhMagIyIioTHoiIhIaLxNz62It/wgolsIg45uKWt31znUf6M+1EmVEAlkhP/PMw9dEhGR0Bh0REQkNAYdEREJjb/REQ3GCP8tgoj6UnRGV1lZCa1Wi8DAQBiNRvT29vZpU1xcjKCgIAQEBCA7O1ve/vLLL2Py5MnQ6XTQ6XT44x//qGSpREQkKMWCzmazwWg0oqSkBE1NTbBarSgqKrJr09XVhdWrV6OqqgoNDQ2oqKhAVVWV/L7JZEJNTQ1qamrsQpCIiGigFAu66upq+Pr6Ijg4GACQlZWF0tJSuzYffvgh4uLiMHnyZIwaNQqLFi3q04aIiMgRigWdxWKBRqORX/v5+aG1tXVQbbZt24awsDCkpqaisbFRqVKJiEhgigWdJEkOtVm2bBlOnz6N2tpapKenIzk5+brtzGYzgoOD5UdnZ+cN10xEROJRLOg0Go3d7KylpQVqtXrAbSZNmgQ3NzcAwPz583HhwgW0t7f3+RyTyYT6+nr5MW7cOCW+DhER3aQUC7rIyEhYLBbU19cDAPLz86HX6+3aPPzww6ioqMC5c+fQ09ODwsJCuU1bW5vcrqKiAm5ubvDx8VGqXCIiEpRi19G5uroiLy8PBoMB3d3diI2NRUZGBsrKylBWVoa8vDx4e3sjJycHc+bMgc1mQ1paGh588EEAwJo1a3D8+HG4urrCy8sLu3fvhkqlUqpcIiLFOLzGqpuTCrlFKXrBeHx8vDyjuyYlJQUpKSny6/T0dKSnp/fp+9577ylZGhER3SK4BBgREQmNQUdEREJj0BERkdAYdEREJDQGHRERCY1BR0REQuP96JTAe5YREY0YnNEREZHQGHRERCQ0Bh0REQmNQUdERELjyShERLc40Red5oyOiIiExqAjIiKhMeiIiEhoDDoiIhIag46IiITGoCMiIqEx6IiISGgMOiIiEhqDjoiIhMagIyIioTHoiIhIaAw6IiISGoOOiIiExqAjIiKhMeiIiEhoDDoiIhIag46IiITGoCMiIqEx6IiISGiKB11lZSW0Wi0CAwNhNBrR29vbp01xcTGCgoIQEBCA7OxsebvVakVycjLuvfdeREZG4t///rfS5RIRkWAUDTqbzQaj0YiSkhI0NTXBarWiqKjIrk1XVxdWr16NqqoqNDQ0oKKiAlVVVQCA1157DaGhoWhsbMT69euxfPlyJcslIiIBKRp01dXV8PX1RXBwMAAgKysLpaWldm0+/PBDxMXFYfLkyRg1ahQWLVoktyktLYXRaAQAPProo/jyyy/R3t6uZMlERCQYRYPOYrFAo9HIr/38/NDa2jrgNv/9nlqthsViUbJkIiISzCgldy5JklPa9MdsNsNsNsuvz5w5I88gh8vfHd3B/+u//s7OTowbN+6Gd690fY7g2DnGofo4djeOY+cYJ4xff5MgRYNOo9HYzeBaWlqgVqv7tKmtrb1uG7VajdbWVkydOhXAj19kypQpdv1NJhNMJpNSX2FECg4ORn19/XCXcVPi2N04jt2N49g5xtHxU/TQZWRkJCwWi1xgfn4+9Hq9XZuHH34YFRUVOHfuHHp6elBYWCi30ev1yMvLAwAcOHAAgYGBmDBhgpIlExGRYBQNOldXV+Tl5cFgMCAgIACenp7IyMhAWVmZfJKJt7c3cnJyMGfOHEybNg2xsbF48MEHAQAvvvgiamtrce+992LdunV4++23lSyXiIgEpOihSwCIj4/vM+VMSUlBSkqK/Do9PR3p6el9+np7e+P9999XusSbzq12qNaZOHY3jmN34zh2jnF0/FSSo2eDEBERjWBcAoyIiITGoBvhVq5cCbVajVGj7I8yr1mzBoGBgQgKCupzET4Bra2teOihhzB9+nRotVqsXbtWfo9jNzCJiYnQ6XQIDQ2FwWCA1WoFwPEbKJPJZPd3y3EbGH9/f2i1Wuh0Ouh0OtTV1QFwcPwkGtEOHz4snTt3TnJ1dZW3lZeXS7Nnz5Z6enoki8UiaTQa6eLFi8NY5cjT1tYmVVdXS5IkSd3d3VJMTIy0d+9ejt0gXLhwQX6+cuVK6fe//z3Hb4AOHTok/epXv5L/bjluA3f33XdLra2tdtscHT/O6Ea4mJgYTJo0yW5baWkpMjMz4erqiilTpiA6OhofffTRMFU4Mk2ePBmRkZEAAHd3d8yYMQMtLS0cu0Hw9vYG8OOatT/88ANUKhXHbwC6u7uxZs0avP766/I2jptjHB0/Bt1NaCBLq9F/dHR0YO/evUhISODYDVJqairuvPNONDQ0YNWqVRy/AVi/fj2ysrIwceJEeRvHbXCSk5Oh0+mQnZ2Nq1evOjx+DLqbkMQTZQfsypUrMBgMWLlyJe677z6O3SDt2bMHbW1tUKvV2LVrF8fvZ9TW1uKTTz7B008/bbed4zZwhw8fxvHjx3HkyBE0NDTg9ddfd3j8GHQ3oYEsrUZAb28vFixYAJ1Oh1WrVgHg2N0Id3d3pKenY8+ePRy/n3HkyBHU19fjnnvugb+/P3p7e+Hv74+JEydy3Abo2szt9ttvh9FoxNGjRx3/786JvyGSgn56MspHH31k98OsWq2WrFbrMFY3Mi1evFjKzMyUbDabvI1jNzBWq1Vqa2uTJEmSent7paVLl0pr167l+A3Stb9bjtvAfPfdd1JXV5ckSZLU09MjGY1G6aWXXnJ4/BRfGYUc8+tf/xr79+9Hb28v1Go15s2bB7PZjPLycgQFBcHFxQWbNm2Cl5fXcJc6ohw5cgTvvPMOQkJCMGPGDADA4sWLsWLFCo7dAFy8eBHz5s1Dd3c3bDYb7r//fqxbtw5jxozh+N2AhIQEjtsAfP3119Dr9bDZbOjt7UVUVBSys7Md/u+OK6MQEZHQ+BsdEREJjUFHRERCY9AREZHQGHRERCQ0Bh0REQmNQUdEREJj0BENkcuXL+PZZ5/F1KlTMWPGDNx///0oKioCAGRmZsrPh8qWLVuQl5c3pJ9JNBx4wTjREFm6dClsNhtOnjyJ0aNH48KFC9i7d++w1bNs2bJh+2yiocQZHdEQOHPmDPbs2YO3334bo0ePBgCMHTsWmZmZfdq+++67mDVrFmbMmIG4uDicPXsWAHDq1ClERUVhxowZ0Gq1ckj+7ne/g1arRVhYGBISEvrs79KlS0hNTUVYWBhCQkKQnZ0NAHj55ZfxyiuvAAAeffRR+UaXEydOxHPPPQcAqKysRExMDCIiIpCYmMgV9+mmxBkd0RA4efIkAgMD5Xu89SclJUVe/X737t146aWXsGPHDrz99ttYsWIFnnrqKUiSBKvVio6ODuzatQsnTpyAi4sLOjs7++zvww8/xJ133ok9e/YAwHXbHDhwAABw9uxZJCYm4plnnkFHRwdeeuklfPDBB/D29kZJSQlWr16NnTt3OjIUREOOQUc0wjQ0NCA7Oxvnz59Hb28vXFx+PPASExOD9evXo7GxEUlJSYiIiEBvby9uv/12ZGZmYu7cuUhOTu6zv/DwcLzwwgt4/vnnkZCQgMTExOt+7nfffQe9Xg+z2YygoCC8//77aGhoQGxsLIAfb8DK9RnpZsRDl0RDQKvVoqmpCVar9WfbLly4EBs2bEBdXR127tyJH374AQCQlpaGAwcOYPLkyVi+fDlycnLg6uqKo0ePYtGiRfjiiy8QHh6Orq4uu/0FBATg888/xwMPPIDCwkKkpKT0+UxJkrBw4UJkZmbKQShJEmbPno2amhrU1NSgtrYWR44cccJoEA0tBh3REJg6dSrmzZsHk8kkB5fVakVhYWGftlarFVOmTAEAbNu2Td7e1NQEjUaDJUuWYOXKlTh27BguXryI9vZ2PPTQQ3j11Vdx2223wWKx2O3PYrHAzc0NBoMBmzZtwrFjx/p8ZnZ2NiZOnIiVK1fK26KionDs2DGcOHECAHD16lXU1dU5PhhEQ4yHLomGyLZt27BmzRpMnz4dnp6ecHd3l28I+1OvvfYaYmNj4ePjg6SkJHl7cXExtm/fDnd3d9x2223YsmULurq68MQTT+Dy5cuw2WyYN28etFqt3f7q6uqwZs0aqFQqSJKEt956q89nbty4EcHBwdDpdAB+nFW+8MIL2LFjB4xGIy5fvoyenh4sX74coaGhzh0YIoXxNj1ERCQ0HrokIiKhMeiIiEhoDDoiIhIag46IiITGoCMiIqEx6IiISGgMOiIiEhqDjoiIhPb/AYJeZas2LRGxAAAAAElFTkSuQmCC",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from thinkstats import two_bar_plots\n",
"\n",
"two_bar_plots(actual_pmf, observed_pmf, width=2)\n",
"decorate(xlabel=\"Class size\", ylabel=\"PMF\")"
]
},
{
"cell_type": "markdown",
"id": "2842ea59",
"metadata": {},
"source": [
"In the observed distribution there are fewer small classes and more large ones.\n",
"And the mean is 29.1, almost 25% higher than the actual mean."
]
},
{
"cell_type": "code",
"execution_count": 36,
"id": "61623eb5",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"29.123376623376622"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"observed_pmf.mean()"
]
},
{
"cell_type": "markdown",
"id": "d6820e00",
"metadata": {},
"source": [
"It is also possible to invert this operation.\n",
"Suppose you want to find the distribution of class sizes at a college, but you can't get reliable data.\n",
"One option is to choose a random sample of students and ask how many students are in their classes.\n",
"\n",
"The result would be biased for the reasons we've just seen, but you can use it to estimate the actual distribution.\n",
"Here's the function that unbiases a `Pmf` by dividing the probabilities by the sizes."
]
},
{
"cell_type": "code",
"execution_count": 37,
"id": "4c52202b",
"metadata": {},
"outputs": [],
"source": [
"def unbias(pmf, name):\n",
" qs = pmf.qs\n",
" ps = pmf.ps / pmf.qs\n",
"\n",
" new_pmf = Pmf(ps, pmf.qs, name=name)\n",
" new_pmf.normalize()\n",
" return new_pmf"
]
},
{
"cell_type": "markdown",
"id": "9728cc41",
"metadata": {},
"source": [
"And here's the result."
]
},
{
"cell_type": "code",
"execution_count": 38,
"id": "e69c6e88",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"23.692307692307693"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"debiased_pmf = unbias(observed_pmf, \"debiased\")\n",
"debiased_pmf.mean()"
]
},
{
"cell_type": "markdown",
"id": "d9b23744",
"metadata": {},
"source": [
"The mean of the debiased `Pmf` is the same as the mean of the actual distribution we started with.\n",
"If you think this example is interesting, you might like Chapter 2 of *Probably Overthinking It*, which includes this and several other examples of what's called the \"inspection paradox\"."
]
},
{
"cell_type": "markdown",
"id": "089ffa11",
"metadata": {},
"source": [
"## NSFG Data\n",
"\n",
"In the previous chapter, we plotting histograms for the distributions of pregnancy lengths for first babies and others.\n",
"But the sizes of the groups are not the same, so we can't compare the histograms directly.\n",
"Because PMFs are normalized so they add up to 1, we can compare them.\n",
"\n",
"TODO: Intro NSFG data\n",
"Instructions for downloading the data are in the notebook for this chapter."
]
},
{
"cell_type": "markdown",
"id": "ead45565",
"metadata": {
"tags": [
"remove-print"
]
},
"source": [
"The following cells download the data files and install `statadict`, which we need to read the data."
]
},
{
"cell_type": "code",
"execution_count": 39,
"id": "31216ce5",
"metadata": {
"tags": [
"remove-print"
]
},
"outputs": [],
"source": [
"try:\n",
" import statadict\n",
"except ImportError:\n",
" !pip install statadict"
]
},
{
"cell_type": "code",
"execution_count": 40,
"id": "7deaf238",
"metadata": {
"tags": [
"remove-print"
]
},
"outputs": [],
"source": [
"download(\"https://github.com/AllenDowney/ThinkStats/raw/v3/nb/nsfg.py\")\n",
"download(\"https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemPreg.dct\")\n",
"download(\"https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemPreg.dat.gz\")"
]
},
{
"cell_type": "markdown",
"id": "550534f3",
"metadata": {},
"source": [
"Again, we can use `read_fem_preg` to read the pregnancy file and use the `query` method to select live births."
]
},
{
"cell_type": "markdown",
"id": "1b94e467",
"metadata": {},
"source": [
"We can use `birthord` to separate first babies and others into two groups."
]
},
{
"cell_type": "code",
"execution_count": 41,
"id": "945aea2d",
"metadata": {},
"outputs": [],
"source": [
"from nsfg import read_fem_preg\n",
"\n",
"\n",
"def get_nsfg_groups():\n",
" \"\"\"Read the NSFG pregnancy file and split into groups.\n",
" \n",
" returns: all live births, first babies, other babies\n",
" \"\"\"\n",
" preg = read_fem_preg()\n",
" live = preg.query(\"outcome == 1\")\n",
" firsts = live.query(\"birthord == 1\")\n",
" others = live.query(\"birthord != 1\")\n",
" return live, firsts, others"
]
},
{
"cell_type": "code",
"execution_count": 42,
"id": "4424061e",
"metadata": {},
"outputs": [],
"source": [
"live, firsts, others = get_nsfg_groups()"
]
},
{
"cell_type": "markdown",
"id": "db6ed39a",
"metadata": {},
"source": [
"And make a `Pmf` for the pregnancy lengths in each group."
]
},
{
"cell_type": "code",
"execution_count": 43,
"id": "900062cf",
"metadata": {},
"outputs": [],
"source": [
"first_pmf = Pmf.from_seq(firsts[\"prglngth\"], name=\"firsts\")\n",
"other_pmf = Pmf.from_seq(others[\"prglngth\"], name=\"others\")"
]
},
{
"cell_type": "markdown",
"id": "26d3185a",
"metadata": {},
"source": [
"Here are the PMFs of pregnancy length for first babies and others, plotted as bar graphs."
]
},
{
"cell_type": "code",
"execution_count": 44,
"id": "62c91cfa",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAboAAAD/CAYAAACHFRPuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAAuJAAALiQE3ycutAAAeJ0lEQVR4nO3de3BU5f3H8c9mDUgJWGRQCdlwSeSyIRIggBhKMIBFA2uJkSqIoFkpQyp0rO1wExVBar1McYhVGnTU0FpDTBpbFRhLoMAME0UECSiRQrJSsUAgKBCS3fP7w7I/Yy4usCebnH2/ZnYme/bZk+/DM/Dhec7NZhiGIQAALCoi1AUAAGAmgg4AYGkEHQDA0gg6AIClEXQAAEsj6AAAlkbQAQAs7YpQF3ApOnfurJiYmFCXAQBoQR6PR9XV1Rf9vTYZdDExMSorKwt1GQCAFuR0Oi/peyxdAgAsjaADAFham1y6BAAr4ZbD9dlstqDuj6ADgBCpra1VZWWlampqQl1Kq9K+fXs5HA5FRkYGZX8EHQCESGVlpTp16qRevXoFfRbTVhmGoePHj6uyslJ9+vQJyj4JOgAIAcMwVFNTo169eikigtMlLrDZbOratauOHTsmwzCC8h8A/nQBIISYyTXEMToA1vH2vMa3T1rZsnVAq1atUk5OjqKjo3X99dfrxRdfDPi7u3btUkVFhVwul4kVXjqCDgBakQVv7Qnq/lZkJAbUbtWqVfr73/+u+Pj4Rj83DEOGYTS6zLpr1y5t3bq11QYdS5cAEObuu+8+HTx4ULfffrtSUlI0btw4SVJJSYlGjhypzMxMDRw4UJWVlZo1a5YSEhJ0ww03aObMmTp9+rSWLFmiwsJCJSUl6fnnn9f+/fs1cuRIDR48WAkJCSoqKgpp/5jRAUCYe+WVV7Rp0yatX79e5eXlWrZsmf+znTt36uWXX9aAAQO0c+dOHT58WHv37pUkVVVVqVOnTlq6dKm2bt2q3NxcSdLcuXM1d+5c3X333TIM45LuTxlMzOgAAE0aMmSIBgwYIEnq06ePDh8+rNmzZ+utt95Shw4dGv3OqFGjtHz5ci1dulQ7d+7UVVdd1ZIlN0DQAQCa1LFjR//PP/7xj/XRRx9p4sSJ2rhxo0aMGCGfz9fgO1OmTNE777yj7t27a86cOXr66adbsuQGWLoEAATkv//9ryIjIzVx4kSNHTtW0dHROnPmjDp16lRvebK8vFxxcXF64IEH1LFjRxUWFoawaoIOABCgyspKud1ueb1e+Xw+PfbYY4qKilJaWpqeeuopJSUl6f7779fXX3+ttWvXql27drryyisv6lIFM9iMNng3UafTyfPoACsI4+voDMPQ/v371b9/fy4a/56m/mwu9d9+jtEBACyNoAMAWBpBBwCwNIIOAGBppgZdSUmJEhISFB8f7z9T5/tsNpuSkpL8r+PHj5tZEgAgzJgWdD6fT263W/n5+SovL1d1dbXy8vIatLPb7dq1a5f/1bVrV7NKAgCEIdOCrrS0VNHR0XI6nZKkrKwsFRQUmPXrAABBVlRUpN27d/vfP/bYY/Xug9lWmHbBuMfjkcPh8L+PjY1VZWVlg3Y+n0/Dhg2Tz+fTtGnT9NBDD5lVEgC0fk1dW3ipLuOaxKKiIo0bN0433HBDUEqpq6vTFVe0/H1KTJvRBXod+uHDh1VaWqoNGzaosLBQf/nLXxq0ycnJkdPp9L+qqqqCXS4AhLXCwkINGjRIiYmJmjp1qt5++20VFxdr0aJFSkpK0vbt2yVJBw4c0NixYxUfH6+FCxf6v797926lpaVp6NChGjVqlPbs+fa5eo899pimTZum0aNHa/z48Tp69KjS0tKUlJSkhIQEvfDCC6b3zbRodTgc9WZwFRUViomJabSdJHXt2lXTpk3T9u3bdffdd9drk52drezsbP/7C8uhAIDL9+WXX2rOnDnasWOHYmNj9eCDD/ofpDpu3Djdc889kqQNGzZoz5492rZtm3w+n+Lj4zVnzhxde+21mjVrlgoKCtSjRw+VlpbK7XZrx44dkr59MOuOHTsUFRWl5557TuPHj9eCBQskqUUmLqYFXXJysjwej8rKyuR0OrVmzRplZGTUa1NVVaUOHTroyiuv1Llz51RcXKzJkyebVRIAoBE7duzQqFGjFBsbK+nbcypmzZrV6KQiPT3d/3ie/v3769ChQzp58qT27t2r9PR0f7sTJ074f3a5XIqKipIkjRgxQvfdd59Onz6tCRMmaPTo0WZ2TZKJS5d2u125ubnKzMxUXFycoqKiNH36dBUXF8vtdkuS9u/fr2HDhmnQoEEaOnSoEhMTlZWVZVZJAIBGfP9em83de7N9+/b+n+12u+rq6mQYhuLi4uqdQV9RUeFv991H/aSkpGjbtm0aMGCAli9fXm+1ziymHhVMS0trcANOl8sll8slSRo5cqR/HRcAEBrDhw/X7Nmz5fF4FBMTo1deeUVpaWn65ptvAno6eP/+/XX69Gm9//77Gjt2rAzD0K5duzR48OAGbQ8dOqQePXpo+vTpiouL07x5QT75phE8pgcAwtx1112nVatWKT09XT6fT4mJiXrppZf0ySefKCsrS6tXr272pJHIyEgVFRVp7ty5euihh1RbW6uMjIxGg27Tpk169tlnFRkZqYiIiBZ5KCuP6QEQOjymh8f0NCLYj+lhRgfAdAveavwQxYrIFi4EYYmbOgMALI2gAwBYGkEHACHUBk+TMF2w/0w4RgcAIWCz2dS+fXsdP35cXbt25YSU/zEMQ8ePH1f79u2D9mdC0AFAiFy4VeKxY8dCXUqr0r59+3oPBbhcBB0AhEhkZKT69OnD8uX3BHt2S9ABQIixbGkuTkYBAFgaQQcAsDSCDgBgaQQdAMDSCDoAgKURdAAASyPoAACWRtABACyNoAMAWBpBBwCwNIIOAGBpBB0AwNIIOgCApRF0AABLI+gAAJZG0AEALI2gAwBYmqlBV1JSooSEBMXHx8vtdsvr9TbZNj09XfHx8WaWAwAIQ6YFnc/nk9vtVn5+vsrLy1VdXa28vLxG265du1ZXX321WaUAAMKYaUFXWlqq6OhoOZ1OSVJWVpYKCgoatDt27JhycnK0aNEis0oBAIQx04LO4/HI4XD438fGxqqysrJBu1/96ldatmyZrrzySrNKAQCEMdOCzjCMH2zz7rvvym63Ky0trdl2OTk5cjqd/ldVVVWwygQAWNwVZu3Y4XDUm8FVVFQoJiamXpstW7bo/fffV69evVRXV6ejR4/qhhtu0O7du+u1y87OVnZ2tv/9heVQAAB+iGkzuuTkZHk8HpWVlUmS1qxZo4yMjHptVqxYIY/Ho0OHDmnr1q3q2bNng5ADAOBymBZ0drtdubm5yszMVFxcnKKiojR9+nQVFxfL7Xab9WsBAKjHtKVLSUpLS/PP6C5wuVxyuVwN2vbq1Uvl5eVmlgMACEPcGQUAYGkEHQDA0gg6AIClEXQAAEsj6AAAlkbQAQAsjaADAFgaQQcAsDSCDgBgaQQdAMDSCDoAgKUFFHSPP/64PB6P2bUAABB0AQVdu3btlJaWpttuu02FhYXyer1m1wUAQFAEFHQLFizQZ599pocfflj5+fmKj4/X/PnzdfDgQbPrAwDgslzUMbqYmBj16NFDNptNX375pSZOnKhHHnnErNoAALhsAT2P7vXXX1dubq7OnDmjX/ziF9qzZ486duyo2tpa9e3bV0888YTZdQIAcEkCCrotW7bo2WefVXJycr3tkZGRev31100pDACAYAho6bJbt24NQm7hwoWSpFGjRgW/KgAAgiSgoHvvvfcabHvnnXeCXgwAAMHW7NLl66+/rtdee02ff/65brnlFv/26upqXXPNNaYXBwDA5Wo26EaPHi2Hw6GHHnpIixYt8m/v1KmTBg0aZHpxAABcrmaDrmfPnurZs6d27tzZUvUAABBUzQbdnDlz9MILL2j8+PGy2WwNPt+wYYNphQEAEAzNBl1WVpYkafHixS1SDAAAwdZs0A0dOlSSlJqa2iLFAAAQbM0GXVNLlhewdAkAaO2aDTqWLAEAbV2zQceSJQCgrTP1rMuSkhJlZ2erpqZGY8aM0UsvvSS73e7//JtvvtHNN9+s8+fPq7a2VikpKXrhhRd0xRUB3YITAIAfZNpZlz6fT263W8XFxXI6nZoyZYry8vI0Y8YMf5sOHTron//8p6KiouTz+ZSZmam8vDzNnDnzon8fAACNCfisy7q6Oh04cECSdP311//grKu0tFTR0dFyOp2Svg3NnJycekEXERGhqKgoSVJdXZ1qamqaPfkFAICLFdBNnUtKStSnTx/dc889mjZtmuLj47V58+Zmv+PxeORwOPzvY2NjVVlZ2WjbESNGqFu3burcubOmTZt2EeUDANC8gILul7/8pYqLi/Xhhx9q586d+tvf/qbs7Oxmv2MYRsBF7NixQ1988YVOnDihkpKSBp/n5OTI6XT6X1VVVQHvGwAQ3gIKusjISCUlJfnfDxo0SJGRkc1+x+Fw1JvBVVRUKCYmpsn2UVFRcrlcevvttxt8lp2drbKyMv+rS5cugZQNAEDzQXfkyBEdOXJEP/3pTzV//nzt27dP+/fv18KFC3Xrrbc2u+Pk5GR5PB6VlZVJktasWaOMjIx6bb766iudPHlSklRTU6N33nlHCQkJl9EdAADqa/aMkpSUFNlsNv8y5F//+lf/ZzabTU8++WST37Xb7crNzVVmZqZqamqUmpqq6dOnq7i4WMXFxcrNzdWRI0d03333yev1yuv16rbbbpPb7Q5S1wAAkGzGxRxMayWcTqd/pgig9Vvw1p5Gt6+IzG38C5NWmlgN2qpL/bc/4Cuzq6urVV5ernPnzvm33XTTTRf9CwEAaEkBBd3atWv1+OOP68svv1Tfvn318ccfa/jw4dq2bZvZ9QEAcFkCOuvy97//vT788EP16dNHH3zwgbZu3arevXubXRsAAJct4MsLOnXqJJ/PJ+nbC7w/+eQTUwsDACAYAlq6jIqK0rlz5zRs2DA9+OCDuu666+rdnBkAgNYqoBldXl6ebDabVq5cqW7duunEiRNat26d2bUBAHDZAprRxcTEqK6uTgcPHtSdd94Z0E2dAQBoDQJKq5KSEt17773q1q2bDMPQiRMn9Oqrr/JgVgBAqxdQ0F24qfOF+11+/PHHmjZtGiekAABaPdNu6gwAQGvQ7IzuyJEjkuS/qfOMGTNks9n02muv/eBNnQEAaA1Mu6kzAACtQbNB9+9//7ul6gAAwBQBXyOwfv16bdq0SZI0duxYjR8/3rSiAAAIloBORnniiSe0ZMkSde/eXdHR0VqyZImWL19udm0AAFy2gGZ0+fn5Ki0tVfv27SVJs2bN0vDhw7Vo0SJTiwMA4HIFNKMzDEMREf/fNCIiQm3wea0AgDAU0Ixu8uTJuvnmmzVt2jRJ0p///GfdcccdphYGAEAwBBR0S5cuVXJyskpKSmSz2fTb3/5WkyZNMrs2AAAu2w8GndfrVWJiosrKyuRyuVqiJgAAguYHj9HZ7XbFxsbq2LFjLVEPAABBFdDSZWRkpBISEjRu3Dh17NjRv3316tWmFQYAQDAEFHSZmZnKzMw0uxYAAIIuoKCbMWOGJOns2bOSpA4dOphXEQAAQRTQdXR79+7VsGHD5HA45HA4NHz4cO3du9fs2gAAuGwBBd2MGTP0yCOP6NixYzp27JiWLFnin+UBANCaBRR058+fr3dpwcSJE1VbW2taUQAABEtAQTdp0iS99NJLOnv2rM6ePavVq1cHdE1dSUmJEhISFB8fL7fbLa/XW+/zXbt2KSUlRQkJCRo4cKCef/75S+sFAABNsBkB3LQyMjKyQUjZ7fZvd2Cz6fz58w2+4/P51LdvXxUXF8vpdGrKlClKT0+vt+T52WefyTAM9evXT9XV1Ro6dKjy8/OVlJTUbD1Op1NlZWWB9A9AK7DgrT2Nbl8Rmdv4FyatNLEatFWX+m9/QDO62tpa+Xy+eq/a2lrV1tY2GnKSVFpaqujoaDmdTklSVlaWCgoK6rXp27ev+vXrJ0nq3LmzBgwYoMrKyovuBAAATQko6C6Fx+ORw+Hwv4+NjW02xD7//HN98MEHSklJMaskAEAYCvgJ4xfrYh7jc/LkSf3sZz/TypUrdfXVVzf4PCcnRzk5Of73VVVVQakRAGB9ps3oHA5HvRlcRUWFYmJiGrQ7c+aM0tPT9cADD+jOO+9sdF/Z2dkqKyvzv7p06WJW2QAAizEt6JKTk+XxePwHDtesWaOMjIx6bWpra5WRkaHx48dr7ty5ZpUCAAhjpgWd3W5Xbm6uMjMzFRcXp6ioKE2fPl3FxcVyu92SpDfffFMbN25UUVGRkpKSlJSUpHXr1plVEgAgDAV0eUFrw+UFQNvC5QUIBlMvLwAAoK0i6AAAlkbQAQAsjaADAFgaQQcAsDSCDgBgaQQdAMDSCDoAgKURdAAASyPoAACWRtABACyNoAMAWBpBBwCwNIIOAGBpBB0AwNIIOgCApRF0AABLI+gAAJZG0AEALI2gAwBYGkEHALA0gg4AYGkEHQDA0q4IdQEAcFHentf49kkrW7YOtBkEHYBWacFbexrdviKyhQtBm8fSJQDA0gg6AIClEXQAAEszNehKSkqUkJCg+Ph4ud1ueb3eBm2mTJmibt26KT4+3sxSAABhyrSg8/l8crvdys/PV3l5uaqrq5WXl9eg3ezZs7V+/XqzygAAhDnTgq60tFTR0dFyOp2SpKysLBUUFDRol5aWpquvvtqsMgAAYc60oPN4PHI4HP73sbGxqqysNOvXAQDQKNOuozMMI2j7ysnJUU5Ojv99VVVV0PYNALA202Z0Doej3gyuoqJCMTExl7Sv7OxslZWV+V9dunQJVpkAAIszLeiSk5Pl8XhUVlYmSVqzZo0yMjLM+nUAADTKtKCz2+3Kzc1VZmam4uLiFBUVpenTp6u4uFhut9vfLj09XSNHjtShQ4cUExOjFStWmFUSACAMmXqvy7S0NP+M7gKXyyWXy+V//49//MPMEgAAYY47owAALI2nFwCwjCafeJCR2MKVoDVhRgcAsDSCDgBgaQQdAMDSOEYHIDjentf49kkrW7YO4HuY0QEALI2gAwBYGkEHALA0jtEBuGiNXa+2IjIEhQABYEYHALA0gg4AYGksXQJoiEsFYCHM6AAAlsaMDghzYXFiSWMz1P/NThvtPzeBthRmdAAASyPoAACWxtIlYCE8jw1oiKADWilCCwgOgg5oazj1H7goBB0QDpo56xCwOk5GAQBYGjM6APg+locthRkdAMDSmNEBwEXgbNi2h6BDq8U/KACCgaBD2xPk4ydNBmpkbtB+B8IAx/VaLVOP0ZWUlCghIUHx8fFyu93yer0N2rz55pvq27ev4uLitGjRIjPLAQCEIdOCzufzye12Kz8/X+Xl5aqurlZeXl69NqdOndLDDz+szZs369NPP9WmTZu0efNms0oCAIQh05YuS0tLFR0dLafTKUnKyspSTk6OZsyY4W/z3nvvacyYMerevbskacaMGSooKFBqaqpZZSFIwv34WeOPtmlkqVOSJq28uOXR/30H1sH4h5ZpQefxeORwOPzvY2NjVVlZ+YNt3nvvPbNKCrmWCIeQ/4W62Od+BbmukPcfQKtjMwzDMGPH69atU2FhodauXStJ2rdvn6ZOnaqPPvrI3+aZZ55RVVWVli9fLkl69913tXr1ahUWFtbbV05OjnJycvzvDxw4oOuvv96MstuEqqoqdenSJdRlhAz9p//0Pzz7f/DgQZ07d+6iv2fajM7hcNSbwVVUVCgmJqZBm927dzfbRpKys7OVnZ3tf+90OlVWVmZC1W0D/af/9J/+h6MLh8IulmknoyQnJ8vj8fgHZM2aNcrIyKjXZsKECdq0aZP+85//qK6uTq+++mqDNgAAXA7Tgs5utys3N1eZmZmKi4tTVFSUpk+fruLiYrndbknSVVddpaefflqjR49Wv379lJqaqptvvtmskgAAYcjUC8bT0tIaTLFdLpdcLpf//V133aW77rrrovb73WXMcET/6X84o//h2/9L7btpJ6MAANAa8PQCAICltfqgq6ys1NixYzVgwAAlJCRowYIF/s/mz5+v+Ph49e3bVwUFBSGs0hxN9b2kpESdOnVSUlKSkpKSNHny5BBXap5bbrlFSUlJSkxMVGZmpqqrqyVZf+wvaKz/4TT+F2RnZ+uKK/7/SEu4jP8F3+1/OI1/r169lJCQ4O/rnj3fXid70eNvtHJHjhwxSktLDcMwjJqaGmPUqFFGUVGRsXHjRuMnP/mJUVdXZ3g8HsPhcBinT58OcbXB1VTfN23aZIwdOzbE1bWMkydP+n+eN2+e8eijj4bF2F/QWP/DafwNwzC2bNli3HvvvYbdbjcMwwir8TeMhv0Pp/Hv2bOnUVlZWW/bpYx/q5/Rde/eXcnJyZKkdu3aafDgwaqoqFBBQYFmzpwpu92uHj16KCUlRRs2bAhxtcHVVN/DyVVXXSXp23unnjt3TjabLSzG/oLG+h9OampqNH/+fD3zzDP+beE0/o31P9xdyvi3+qD7rhMnTqioqEjjx48P6BZjVvLdvkvShx9+qMGDB2v06NFav359iKsz1+TJk3XNNdfo008/1a9//euwG/vv918Kn/FfunSpsrKy1K1bN/+2cBr/xvovhc/4S9KkSZOUlJSkRYsWqba29pLGv808j+78+fPKzMzUvHnz1L9/fxlhdLLo9/seHR2tw4cPq3Pnztq7d68mTJigLVu2qHfv3qEu1RSFhYU6f/68srKytG7durAae6lh/zMyMsJi/Hfv3q0dO3Zo2bJl9baHy/g31f8hQ4aExfhL0r/+9S85HA598803mjFjhp555plLGv82MaPzer2aOnWqkpKS/P+jDeQWY1bQWN87d+6szp07S5ISEhKUkpKinTt3hrJM07Vr10533XWXCgsLw2bsv+u7/Q+X8d+2bZvKysrUu3dv9erVS16vV7169VK3bt3CYvyb6r+ksBh/Sf6ZW8eOHeV2u7V9+/ZL+/tv1kHEYLr//vuNmTNnGj6fz79tw4YN9Q5IxsTEGNXV1SGs0hyN9f3IkSP+9x6Px4iNjTX27dsXqhJNU11dbRw5csQwDMPwer3GrFmzjAULFoTN2DfV/3AZ/++7cDJGuIz/913of7iM/9dff22cOnXKMAzDqKurM9xut7Fw4cJLGv9Wv3S5bds2vfzyyxo4cKAGDx4sSbr//vs1d+5cbdy4UX379lVERISee+45derUKcTVBldTfY+IiNAf//hHRUZGSpKWLVum/v37h7JUU5w+fVq33367ampq5PP5NGLECC1evFg/+tGPLD/2UtP9f/nll8Ni/Jsyfvz4sBj/phQUFITF+B89elQZGRny+Xzyer0aOXKkFi1adEl//7kzCgDA0trEMToAAC4VQQcAsDSCDgBgaQQdAMDSCDoAgKURdAAASyPoAAsYM2aMtm7dGuoygFaJoAMAWBpBB7SwP/zhD3r88cf9P19zzTWSpHPnzik2NlZnz57V7NmzNXz4cCUmJmrVqlX+7+bn52vEiBEaPHiw7rjjDp06darB/hcvXqypU6eqrq5OL774ogYOHKhBgwZpyJAhOnfuXMt0EmhFWv0twACrSU1N9d+ge/Pmzerdu7fKysr01Vdfafjw4VqxYoWGDBmiF198UTU1NUpJSVFaWprsdrv+9Kc/acuWLWrfvr2efvppPfnkk3rqqackffvMujlz5sgwDOXl5SkiIkLLly/XZ599pg4dOujUqVNq165dKLsOhAS3AANamM/nU0xMjA4ePKgbb7xRc+bMkdfr1dGjR9W1a1e9+uqrOnv2rP9ehtXV1Xr22Wf1xRdfaPny5br22mslSbW1tUpMTNQbb7yhMWPG6NSpUxo9erRWrlzp/10ul0t2u13p6elyuVz+2SMQTli6BFpYRESEkpOT/TfsTk1N1ebNm1VSUqLU1FQZhqG1a9dq165d2rVrlw4ePKjJkyfLMAz9/Oc/92/fu3ev3njjDf9+b7rpJm3fvl1VVVX+bUVFRfrNb36jL774QsnJyTp48GAougyEFEEHhEBqaqp+97vfacyYMerXr5/279+v8vJyJSYm6tZbb9XKlSvl9XolSQcOHFB1dbXGjRunwsJCeTweSdKZM2e0f/9+/z7vvvtuPfjgg5owYYJOnTqluro6HTp0SDfddJMeffRRDRgwQPv27QtJf4FQIuiAEEhNTVVlZaXGjBkjSerfv7+GDRsmm82mxYsXKyoqSoMGDdLAgQP1wAMP6Pz58xowYICee+45uVwuDRo0SDfeeKP27t1bb7/33nuvZs2apVtvvdX/VObExEQlJiaqR48euuWWW0LQWyC0OEYHALA0ZnQAAEsj6AAAlkbQAQAsjaADAFgaQQcAsDSCDgBgaQQdAMDSCDoAgKX9H6gsH5KA/fKaAAAAAElFTkSuQmCC",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"two_bar_plots(first_pmf, other_pmf)\n",
"decorate(xlabel=\"weeks\", ylabel=\"probability\", xlim=[20, 50])"
]
},
{
"cell_type": "markdown",
"id": "cc505018",
"metadata": {},
"source": [
"By plotting the PMF instead of the histogram, we can compare the two distributions without being misled by the difference in sample size.\n",
"Based on this figure, first babies seem to be less likely than others to arrive on time (week 39) and more likely to be late (weeks 41 and 42)."
]
},
{
"cell_type": "markdown",
"id": "8f2e5b2f",
"metadata": {},
"source": [
"## Other Visualizations\n",
"\n",
"Histograms and PMFs are useful while you are exploring data and trying to identify patterns and relationships.\n",
"Once you have an idea what is going on, a good next step is to design a visualization that makes the patterns you have identified as clear as possible.\n",
"\n",
"In the NSFG data, the biggest differences in the distributions are near the mode.\n",
"So it makes sense to zoom in on that part of the graph, and select data from weeks 35 to 46.\n",
"\n",
"When we call a `Pmf` object like a function, we can look up a sequence of quantities and get a sequence of probabilities."
]
},
{
"cell_type": "code",
"execution_count": 45,
"id": "9f134fc5",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.03602991, 0.03897575, 0.04713347, 0.06163608, 0.4790392 ,\n",
" 0.12145932, 0.08157716, 0.04645366, 0.01971448, 0.00521187,\n",
" 0.00135962])"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"weeks = range(35, 46)\n",
"first_pmf(weeks)"
]
},
{
"cell_type": "code",
"execution_count": 46,
"id": "e185339f",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.03210137, 0.03146779, 0.05216473, 0.07074974, 0.54466737,\n",
" 0.12249208, 0.04794087, 0.02597677, 0.01288279, 0.00485744,\n",
" 0.00084477])"
]
},
"execution_count": 46,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"other_pmf(weeks)"
]
},
{
"cell_type": "markdown",
"id": "332a7946",
"metadata": {},
"source": [
"So we can compute the differences in the probabilities like this."
]
},
{
"cell_type": "code",
"execution_count": 47,
"id": "36d0ba81",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 0.00392854, 0.00750796, -0.00503126, -0.00911366, -0.06562817,\n",
" -0.00103276, 0.03363629, 0.02047689, 0.00683169, 0.00035443,\n",
" 0.00051485])"
]
},
"execution_count": 47,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"diffs = first_pmf(weeks) - other_pmf(weeks)\n",
"diffs"
]
},
{
"cell_type": "markdown",
"id": "cdad3427",
"metadata": {},
"source": [
"Here's what they look like, multiplied by 100 to express the differences in percentage points."
]
},
{
"cell_type": "code",
"execution_count": 48,
"id": "3f40b276",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAboAAAD/CAYAAACHFRPuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAAuJAAALiQE3ycutAAAeYklEQVR4nO3de3RNZ+LG8eckkVCCSDVFUIuxSFySkaajkYprgrbiOi1FEJ0GnXbGVJVRQzMzpqoWneigLh2qWpRxqakYVFxjNFERdakSSaR6CSGTRC7n94fl/Chi1zk7Jzm+n7Wy2r1zst8ne2kf775arFarVQAAuCg3ZwcAAMBMFB0AwKVRdAAAl0bRAQBcGkUHAHBpFB0AwKVRdAAAl+bh7AA/V+3ateXv7+/sGACASiQzM1N5eXm3/V6VKzp/f3+lp6c7OwYAoBIJCAi44/c4dAkAcGkUHQDApVF0AACXRtEBAFwaRQcAcGkUHQDApVF0AACXVuXuowNgjkcmba6wsc7M7FNhYwHM6AAALo2iAwC4NIoOAODSKDoAgEuj6AAALo2iAwC4NIoOAODSKDoAgEuj6AAALo2iAwC4NIoOAODSKDoAgEuj6AAALo2iAwC4NIoOAODSDL2P7urVqzp06JCys7NVo0YNtWnTRk2aNDE7GwAAdiu36E6fPq34+Hht27ZNrVu3Vv369VVYWKgTJ07Iw8NDL774omJiYmSxWCoqLwAAP0u5Rfe73/1OL7/8shYvXnxLmX377bdaunSpFi9erNjYWFNDAgBwr8otun/96193/J6fn58mTZrk8EAAADiSoYtRNmzYoMuXL0uSZs+erV//+tc6cuSIqcEAAHAEQ0U3ZcoUeXt764svvtDq1asVHR2tF154wexsAADYzVDReXp6SpI2b96sMWPG6Nlnn1VBQYHdg587d07dunVT69atFRgYqNdee83ubQIAcCNDRefl5aWPPvpIH3/8sSIjIyVdu+XAXh4eHvrb3/6mY8eOKSUlRbt37y73vCAAAD+XoaJbunSp9u7dq9dff13+/v76+uuvNXToULsHb9CggUJCQiRdmzUGBwcrIyPD7u0CAHCdoRvGExMTNXfuXNty8+bN5e3t7dAgP/74o9avX6+tW7c6dLsAqo5HJm2usLHOzOxTYWPBuQzN6JYsWWJo3b26evWqBg4cqJdeekmtWrW66XsJCQkKCAiwfeXm5jpsXACA6yt3Rrdx40Zt2LBB586d0/PPP29bn5eXp2rVqjkkQGlpqYYMGaKgoCBNmDDhlu+PGzdO48aNsy0HBAQ4ZFwAwP2h3KJr0qSJOnXqpK1btyosLMy23tvbW127dnVIgOeff17e3t6aPXu2Q7YHAMCNyi269u3bq3379howYIBq1arl8MH37NmjJUuWqE2bNgoODpYkjRo1Sr/97W8dPhYA4P5k6GKUnJwczZo1S998841KSkps67dv327X4GFhYbJarXZtAwCA8hgquoEDB2rYsGEaOnSo3N3dzc4EAIDDGCo6i8Vy2wtFAACo7AzdXtC/f3+tXr1apaWlZucBAMChDM3oZsyYodLSUrm5ucnNzU1Wq1UWi8UhjwEDAMBMhoquuLjY7BwAAJjC0KFLAACqqnKLLjQ0VJJUrVo1eXp6qlq1arav66/uAQCgMiv30GVycrIkDl0CAKouQ+foJOnYsWP6/PPPJUkRERG3PHwZAIDKyNA5usWLF+vJJ5/U8ePHdeLECT311FNaunSp2dkAALCboRnd3Llz9d///lc+Pj6SpKlTp6pz584aOXKkqeEAALCXoRmd1Wq1lZwk+fj48IxKAECVYGhGFx4ersGDB9tmcO+//77Cw8NNDQYAgCMYPnQ5f/58LVq0SBaLRREREXrhhRfMzgYAgN0MFV21atX00ksv2d4yXqNGDVNDAQDgKIbO0aWlpenRRx9V48aN1bhxY4WGhiotLc3sbAAA2M1Q0cXExGjq1Kn6/vvv9f333+v1119XTEyMydEAALCfoaK7evWqnn76advyk08+ydNSAABVgqGie+qpp7RgwQIVFBSooKBACxcuvKn4AACorAxdjPLmm2+qtLRUcXFxtnXu7u6aOXMm76WzwyOTNlfIOGdm9qmQcQCgMuJ9dABwA/4C6np4Hx0AwKVRdAAAl0bRAQBcmuGiu3jxor744gszswAA4HCGiu6TTz5RaGioBgwYIElKTU1V7969TQ0GAIAjGLrqMj4+XgcPHlRERIQkKSgoSBkZGWbmMh1XVgHA/cHQjM7T01N16tS5+QfdOL0HAKj8DM3omjVrpm3btslisSgvL09vv/222rdvb3Y2AADsZmhaNn/+fC1fvlyZmZlq3ry5Tp8+rXnz5jkkwM6dOxUYGKgWLVooNjZWpaWlDtkuAACSwaLz8fHR+++/rwsXLui7777TP//5T/n4+Ng9eFlZmWJjY7V69WqdOnVKeXl5WrFihd3bBQDgOkOHLv/yl7/csq5OnToKDQ3Vo48+es+DHzx4UA0bNlRAQIAkafTo0UpISNCIESPueZsAANzIUNEdPXpUBw4c0NNPPy2LxaINGzYoKChIS5Ys0aBBgzRp0qR7GjwzM1ONGze2LTdp0kTnzp27p20BABzHla5Mt1itVuvdPvTEE09o06ZNql27tiQpLy9PTz31lLZs2aKQkBClp6ff0+Br1qzRunXr9MEHH0iSjh07piFDhiglJcX2mYSEBCUkJNiWc3Nzdf78+XsaD7eqDH+YyVB5MgBVVUBAwB27yNA5ugsXLqhWrVq25Zo1a+rbb7/VAw88IC8vr3sO1rhx45tmcBkZGfL397/pM+PGjVN6errtyxHnBgEA9w9Dhy579eql7t27a/DgwZKuzcR69eql/Px8+fn53fPgISEhyszMVHp6ugICArR48WL179//nrcHVGXMtABzGCq6OXPmaN26dUpKSpIkjR8/XtHR0ZKkf//73/c8uLu7u9577z0NHDhQRUVF6ty5s4YNG3bP2wMA4KcMFZ0k9evXT/369XN4gK5du97zOT4AAO7GUNGdPXtWr732mtLS0lRYWGhbf+LECdOCAQDgCIYuRhk5cqQGDRqk0tJSffLJJ4qMjNSQIUPMzgYAgN0MFd2lS5fUr18/ubu7q02bNnrnnXe0eXPFXAoNAIA9DB269PC49rEHH3xQycnJatiwoX744QdTgwEA4AiGii4uLk65ubn685//rGeeeUaXL1/Wm2++aXY2AADsZqjounbtKh8fH3Xs2FGnT5+WpCr/4lVcw71bAFydoXN01++Zu9s6AAAqm3JndJmZmTp79qzy8/O1b98+XX8s5qVLl/S///2vQgICAGCPcovuP//5j5YtW6bs7GxNnjzZtt7b21szZ840PRwAAPYqt+hGjBihESNGaOXKldw3BwCokgxdjDJkyBDt27dP33zzjUpKSmzrhw8fblowAAAcwVDRxcTE6PDhwwoODpa7u7skyWKxUHQAgErPUNHt27dP6enptpIDAKCqMFR0LVu21OXLl1W3bl2T4+B+xL18AMxkqOjc3d3Vtm1bde3a9aY3ii9cuNC0YAAAOIKhojPrXXQAAJjNUNGNGDFCkvTdd9+pfv36pgYCAMCRDD0CLCkpSc2bN1dISIgkKSUlRbGxsaYGAwDAEQwV3YQJE7Rz507Vq1dPkhQcHKz9+/ebGgwAAEcwVHSS1Lhx45uWr7+jDgCAysxQ0fn5+enYsWOyWCySpKVLl6pJkyamBgMAwBEMTcsSEhI0atQopaeny9fXV4GBgVqxYoXZ2QAAsJuhomvSpIm2bdum/Px8lZWVydvb2+xcAAA4hKFDl3PmzFFubq5q1qwpb29v/fjjj5o7d67Z2QAAsJuhonv//ffl4+NjW65Xr56WLVtmViYAABzGUNGVlZXdtGy1WlVcXGxKIAAAHMlQ0YWFhek3v/mNjh49qqNHjyouLk6dOnUyOxsAAHYzfI6ufv36iomJUUxMjOrXr685c+aYnQ0AALvd9arL0tJSjRkzRsuXL1d8fHxFZAIAwGHuOqNzd3fX119/fct5OnslJiYqJCREbdq0Ufv27fXRRx85dPsAAEgG76Nr166dIiMjNXjwYNWsWdO2fsiQIfc88IMPPqj169fL399f2dnZCg4OVteuXXk7AgDAoQwVXVFRkfz9/bV3717bOovFYlfRBQcH2/69YcOG8vPzU05ODkUHAHAoQ0W3dOlSU0Ps2bNH+fn5at26tanjAADuP4aKLisrSxMmTFB2drZ27dqltLQ0JSUlKS4urtyfi4yMVFZW1i3ro6OjbRe2ZGRkaNiwYVq+fPlt34iQkJCghIQE23Jubq6RyAAASJIsVqvVercPRUVFKTY2VvHx8UpNTVVJSYmCg4N15MgRuwa/cOGCIiIiFB8fr/79+xv6mYCAAKWnp9s1LvBTj0zaXCHjnJnZp0LGAe435XWDofvofvjhBw0cOND2mh4PDw+730eXl5enqKgovfLKK4ZLDgCAn8tQ0VWvXl2XL1+2Fd3Ro0dVvXp1uwZ+5513dOzYMc2dO1dBQUEKCgq66WIXAAAcwdC0LD4+XlFRUcrIyFD//v2VnJysDz/80K6Bp0yZoilTpti1DQAA7sZQ0XXu3FmbNm3S3r17ZbVatWjRIvn6+pqdDQAAuxk6dClJ2dnZysjI0Llz5/Ttt9+amQkAAIcx/FDnvn376tSpUzp58qT69u3Li1cBAFWCoUOX//jHP3To0CHVqVNHkjRt2jSFhobqpZdeMjUcAAD2MjSjq1u3rmrVqmVbrlmzpurWrWtWJgAAHMbQjO6Xv/ylunTpomeffVYWi0WrVq1Shw4dtHLlSkn2PdwZAAAzGSq6wsJCNW/eXMnJyZKkZs2aqaCgQImJiXY/3BkAADNVioc6AwBglnLP0a1bt67cH87JydGBAwccGggAAEcqd0aXkpKi6dOnKzIyUiEhIXrooYdUWFiokydPatu2bbp06ZLmzZtXUVkBAPjZyi26GTNm6JVXXtGaNWu0ZcsWZWdnq0aNGgoMDNTkyZMVGhpaUTkBALgndz1H5+3trZEjR2rkyJEVkQcAAIcy/AgwAACqIooOAODSKDoAgEszVHSlpaVasGCBXn31VUnSmTNntGfPHlODAQDgCIaKLi4uTmlpadq4caOka8++5IHOAICqwNCTUQ4ePKiUlBQFBwdLulZ0xcXFpgYDAMARDM3oPD09JUkWi0WSdOnSJfMSAQDgQIaKbtiwYRo+fLhyc3M1b948RUREaOzYsWZnAwDAboYOXY4fP147duyQn5+fMjIy9NZbb6lbt25mZwMAwG6Giu7ixYsKDw9Xly5dJEnFxcW6ePEiL18FAFR6hg5dRkZGqqioyLZcVFSkXr16mRYKAABHMVR0RUVFqlmzpm25Vq1aKigoMC0UAACOYqjovLy8dPLkSdvyiRMnVK1aNdNCAQDgKIbO0f31r39Vly5d1LFjR1mtVh04cEDLly83OxsAAHYzVHRdu3bV4cOHtX//flmtVi1YsEC+vr5mZwMAwG6Gik669jSU4OBglZSUKD8/X/n5+WrSpImZ2QAAsJuhops/f76mTJmievXqyc3t2mk9i8WiEydOmBoOAAB7GSq62bNnKz09XQ0aNDA7DwAADmXoqsuGDRuaVnLFxcUKCgpS9+7dTdk+AOD+ZmhGFxQUpKefflr9+vWTl5eXbf2QIUPsDvDmm2+qbdu2On/+vN3bAgDgpwzN6K5cuSJfX1/t2rVLiYmJSkxM1LZt2+we/Pjx49q1a5dGjRpl97YAALgdQzO6pUuXOnxgq9WqsWPHKiEhQTk5OQ7fPgAAksGiy8rK0oQJE5Sdna1du3YpLS1NSUlJiouLK/fnIiMjlZWVdcv66OhoNWrUSGFhYWrVqlW5RZeQkKCEhATbcm5urpHIAABIkixWq9V6tw9FRUUpNjZW8fHxSk1NVUlJiYKDg3XkyJF7Hnjo0KFKSkqSm5ubCgsLdenSJXXp0kWffvppuT8XEBCg9PT0ex4XuJ1HJm2ukHHOzOxTIeMA95vyusHQOboffvhBAwcOtL1h3MPDQx4ehu81v60PPvhAGRkZOnPmjFatWqWwsLC7lhwAAD+XoaKrXr26Ll++bCu6o0ePqnr16qYGAwDAEQxNy+Lj4xUVFaWMjAz1799fycnJ+vDDDx0WIiIiQhEREQ7bHgAA19216MrKyuTl5aVNmzZp7969slqtWrRoEQ91BgBUCXctOjc3N8XFxSklJUV9+nAiHQBQtRg6RxcaGqqdO3eaHAUAAMczdI4uMTFRixYtUoMGDVSzZk1ZrVbeXgAAqBIMFd2OHTvMzgEAgCkMHbps2rSpCgoKlJqaqqZNm8rX11fe3t5mZwMAwG6Gim7+/PkaNmyYJkyYIEnKycnRwIEDTQ0GAIAjGCq6BQsWaM+ePbZZXIsWLfTdd9+ZGgwAAEcwVHQ1atSQp6enbbmsrMy0QAAAOJKhogsODtZ7772nkpISffnll4qNjVWXLl3MzgYAgN0MFd2cOXN07tw5Va9eXaNHj1ajRo00a9Yss7MBAGC3cotu6tSpkqTPPvtM06dP18GDB3Xw4EG98cYb8vLyqpCAAADYo9yi27BhgyRp+vTpFRIGAABHK/eG8Xbt2qlx48b6/vvv1bJlS9t6nowCAKgqyi265cuXKycnR926deOlqACAKqncouvVq5e2bNmivn37qmnTphWVCQAAhyn3HN3Zs2d1+vRpbdy4UefPn1d2dvZNXwAAVHblzujGjh2rPn366PTp0woLC5PVarV9z2Kx6PTp06YHBADAHuXO6MaPH69jx45p6NChOn36tL755hvbFyUHAKgKDN0wvmTJErNzAABginKLLjQ0VJJUrVo1eXp63vJPAAAqu3LP0SUnJ0uSiouLKyQMAACOZugN40ePHlV6erosFosCAwPVunVrs3MBAOAQ5RbdlStX1L9/f+3bt08tWrSQ1WrV119/rU6dOmnt2rV64IEHKionAAD3pNxzdNOnT1e9evWUlZWllJQUpaamKisrS/Xq1eP5lwCAKqHcovv000+1aNEi1a5d27audu3aevfdd7Vp0ybTwwEAYK9yi85iscjb2/uW9bVr15bFYjEtFAAAjlJu0ZV3CwG3FwAAqoJyL0ZJTU29baFZrVaVlZWZFgoAAEcpd0ZXVlamq1ev3vJVXFys0tJSuwc/ceKEunbtqoCAALVq1UqHDh2ye5sAANzI0H10ZigrK9PgwYOVkJCgsLAwFRYWqqioyFlxAAAuymlFl5iYqGbNmiksLEySVL16dVWvXt1ZcQAALsrQQ53N8NVXX6lGjRp68sknFRwcrBdffFGFhYXOigMAcFGmzugiIyOVlZV1y/ro6Gj5+Pho+/btOnTokB5++GHFxsZq1qxZmjp16k2fTUhIUEJCgm05NzfXzMgAABdjatF99tlnd/zexx9/rLCwMDVq1EiSNGjQIC1cuPCWz40bN07jxo2zLQcEBDg+KADAZTnt0GVUVJS++uorXbp0SdK1c3aBgYHOigMAcFFOuxildu3amjFjhjp16iRJCgwM1KJFi5wVBwDgopxWdJI0YMAADRgwwJkRAAAuzmmHLgEAqAgUHQDApVF0AACXRtEBAFwaRQcAcGkUHQDApTn19gKgsjgzs4+zIwAwCTM6AIBLo+gAAC6NogMAuDSKDgDg0ig6AIBLo+gAAC6NogMAuDSKDgDg0ixWq9Xq7BA/R+3ateXv7+/sGE6Tm5srHx8fZ8dwOvbDNewH9sF19/t+yMzMVF5e3m2/V+WK7n4XEBCg9PR0Z8dwOvbDNewH9sF17Ic749AlAMClUXQAAJdG0VUx48aNc3aESoH9cA37gX1wHfvhzjhHBwBwaczoAAAujaKrxHr27KmgoCC1bdtWAwcOtF06e+DAAYWGhiowMFCBgYHKzs52clJz3W4/FBUVafjw4Wrbtq1at26t8ePHq6yszNlRK8S4cePk4fH/r5KcNGmSWrRooZYtW2rt2rVOTFZxbtwHH3zwgdq3b6927dopJCRE27dvd3K6ivPTPwuSdPHiRTVq1EixsbFOSlX58OLVSmz16tWqU6eOJOnll1/W22+/rQkTJmj48OFav369WrdurUuXLsnT09PJSc11u/3QqFEjlZSU6MiRIyotLVXnzp21ZcsW9enj2i9QTUpK0pUrV2zL27Zt0969e3X8+HHl5OSoY8eOioyMVK1atZyY0lw/3QfNmjXT9u3b5evrqyNHjqh79+46f/683Nxc++/xP90P102cOFFdu3Z1QqLKy7X/JFRx1//nXlZWpsLCQlksFq1cuVK9e/dW69atbZ+pUaOGM2Oa7nb7QZLy8/NVUlKioqIiFRUV6eGHH3ZmTNMVFRVp0qRJeuutt2zr1q5dq5iYGLm7u6tRo0YKCwvT1q1bnZjSXLfbB48//rh8fX0lSW3atFFRUZHy8/OdFbFC3G4/SNKOHTtUVFSkbt26OSlZ5UTRVXL9+vXTQw89pOPHj2vChAn66quvVFRUpO7duys4OFjTpk3T/XA90U/3w/Dhw+Xp6Sk/Pz/5+fmpZ8+e6tChg7NjmmrGjBkaPXq06tevb1uXmZmpxo0b25abNGmic+fOOSNehbjdPrjRqlWr1LZtW3l7e1dwsop1u/1QUFCg11577ZbyA0VX6a1bt07Z2dny9/fXmjVrVFJSoqSkJK1Zs0b79u3T/v37tWLFCmfHNN1P90NycrJKSkqUnZ2tc+fOafPmzdq7d6+zY5rmyy+/1IEDBzRy5Mib1t8Pf8m57k774LqUlBRNmjRJS5YsqeBkFetO++FPf/qTnn/++Tv+JeB+xjm6KsDT01PPPPOMFi5cqLCwMPXs2VN169aVJPXt21eHDh3SsGHDnBuyAty4H/z9/dW7d295eXnJy8tLPXv21P79+/X44487O6Yp9uzZo/T0dDVr1kySVFpaqkceeURdunS5aQaXkZGh0NBQZ8U01Z32wZdffqmcnBwNGDBAq1at0i9+8QsnJzXXnfaDj4+PPvroI82YMUNXrlxRUVGR3NzctHDhQicndj7uo6ukLl++rCtXrqhBgwYqKytTXFycfH19NWrUKD333HP6/PPPVa1aNQ0aNEhRUVEaM2aMsyOb4k77oW7duvriiy+0cuVKFRcXKyIiQtOmTVNUVJSzI1cIDw8PlZSUKDExUW+88YZ27NihnJwc/epXv1J6errLH7qT/n8fZGZmqkuXLvr73/+uyMhIZ8eqcNf3w42WLVum3bt367333nNSqsqFGV0ldfnyZfXt21dFRUUqKyvTY489pj/+8Y964IEHNGzYMAUHB8vNzU1PPPGERo0a5ey4prnTfrBarRo9erTatGkjq9Wq6Ojo+6bkbtSjRw8lJiaqZcuWcnNz09tvv31flNyNpk+frgsXLujVV1/Vq6++Kklav369HnnkEecGQ6XBjA4A4NK4GAUA4NIoOgCAS6PoAAAujaIDALg0ig4A4NIoOgCAS6PoACeZOHGiZsyYYVv+/e9/r3bt2tmWd+3adU/P74yJibkvHgsHGEXRAU4SHh6u3bt325b37dunatWq6eLFi5Kk3bt3Kzw83EnpANdB0QFOEhYWpgMHDqi0tFT5+fkqKytTjx49bA+nvl50U6ZMUWhoqNq1a6fJkyfbfn7nzp3q1KmTOnTooJ49e972rQULFixQjx49dOXKFa1bt07t2rVTUFCQ2rVrp7Nnz1bY7wo4E0UHOEm9evXUuHFjHT58WPv379djjz2msLAw7d69W1arVfv27dOpU6ckScnJyUpNTVVaWpq2bNmiH3/8UZMnT9bmzZt16NAhjRkzRn/4wx9u2v7MmTP16aefatOmTapVq5amTZumzz77TKmpqTpw4ID8/Pyc8WsDFY5nXQJO1KlTJyUlJSk3N1dhYWF6/PHH9dZbbyktLU0PPfSQDh48qMOHD2vz5s2Srr1s9uTJkyotLdXx48fVuXNnSddeSnvjMy5nz56tZs2aae3atfLwuPafeUREhJ577jlFR0erb9++atKkScX/woATUHSAE4WHh2v9+vW6ePGiYmNj5evrq7y8PG3fvl3h4eHKzc3VrFmzFB0dfdPPbdy40fazt9OhQwft2bNHWVlZatq0qSRp3rx5SklJUWJiojp37qwVK1YoLCzM5N8QcD4OXQJOdH1Gd/78efn7+0u6VlLvvPOOwsPD1atXL7377rsqLCyUJGVnZysnJ0cdO3ZUcnKy0tLSJEnFxcU6cuSIbbsRERGaN2+eevfubTt3d+LECQUHB2vixInq0aOHUlNTK/aXBZyEGR3gRE2bNpWnp6fat29vWxcWFqbFixcrPDxczZo1U2ZmpkJCQmSxWFSzZk0tW7ZMrVq10ocffqjY2FgVFBSopKREY8eOVdu2bW3b6dGjh2bPnq1evXpp69atmjhxok6dOiUPDw81bdpUzz33nDN+ZaDC8ZoeAIBL49AlAMClUXQAAJdG0QEAXBpFBwBwaRQdAMClUXQAAJdG0QEAXBpFBwBwaf8HuEx2BucqE0YAAAAASUVORK5CYII=",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.bar(weeks, diffs * 100)\n",
"decorate(xlabel=\"Weeks\", ylabel=\"Difference (percentage points)\")"
]
},
{
"cell_type": "markdown",
"id": "11615bb9",
"metadata": {},
"source": [
"This figure makes the pattern clearer: first babies are less likely to be born in week 39, and somewhat more likely to be born in weeks 41 and 42.\n",
"\n",
"For now we can't be sure this effect is real -- it might be due to random variation.\n",
"We'll revisit that question in Chapter xxx."
]
},
{
"cell_type": "markdown",
"id": "93d2684e",
"metadata": {},
"source": [
"## Glossary\n",
"\n",
"- **Probability mass function (PMF)**: a representation of a distribution as a function that maps from values to probabilities.\n",
"\n",
"- **probability**: A frequency expressed as a fraction of the sample size.\n",
"\n",
"- **normalization**: The process of dividing a frequency by a sample size to get a probability.\n",
"\n",
"- **index**: In a Pandas `DataFrame`, the index is a special column that contains the row labels."
]
},
{
"cell_type": "markdown",
"id": "e818cde9",
"metadata": {
"collapsed": true
},
"source": [
"## Exercises\n",
"\n",
"For the exercises in this chapter, we'll use the NSFG respondent file, which contains one row for each respondent."
]
},
{
"cell_type": "code",
"execution_count": 49,
"id": "461a9ed7",
"metadata": {
"tags": [
"remove-print"
]
},
"outputs": [],
"source": [
"download(\"https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemResp.dct\")\n",
"download(\"https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemResp.dat.gz\")"
]
},
{
"cell_type": "markdown",
"id": "7c7ce2b7",
"metadata": {
"tags": [
"remove-print"
]
},
"source": [
"The codebook for this dataset is at ."
]
},
{
"cell_type": "markdown",
"id": "1b40d801",
"metadata": {},
"source": [
"The `nsfg.py` module provides a function that reads the respondent file and returns a `DataFrame`."
]
},
{
"cell_type": "code",
"execution_count": 50,
"id": "861e38b1",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(7643, 3092)"
]
},
"execution_count": 50,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from nsfg import read_fem_resp\n",
"\n",
"resp = read_fem_resp()\n",
"resp.shape"
]
},
{
"cell_type": "markdown",
"id": "5495601a",
"metadata": {},
"source": [
"This `DataFrame` contains 7643 rows and 3092 columns."
]
},
{
"cell_type": "markdown",
"id": "8e509710",
"metadata": {},
"source": [
"### Exercise\n",
"\n",
"Select the column `numbabes`, which records the \"number of babies born alive\" to the respondent.\n",
"Make a `Hist` object and display the frequencies of the values in this column.\n",
"Check that they are consistent with the frequencies in the code book.\n",
"Are there any special values that should be replaced with `nan`?\n",
"\n",
"Then make a `Pmf` object and plot it as a bar graph. Is the distribution symmetric, skewed to the left, or skewed to the right?"
]
},
{
"cell_type": "code",
"execution_count": 51,
"id": "e4933509",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"pmf.bar()\n",
"decorate(xlabel=\"Number of babies\", ylabel=\"PMF\")"
]
},
{
"cell_type": "code",
"execution_count": 54,
"id": "64937b8c",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"id": "071ae7e7",
"metadata": {},
"source": [
"### Exercise\n",
"\n",
"In the same way that the mean identifies a central point in a distribution, and variance quantifies its spread, there is another statistic, called **skewness**, that is intended to indicate whether a distribution is skewed to the left or right.\n",
"\n",
"Given a sample, we can compute the skewness by computing the sum of the cubed deviations and dividing by the standard deviation cubed.\n",
"For example, here's how we compute the skewness of `numbabes`.\n"
]
},
{
"cell_type": "code",
"execution_count": 55,
"id": "f67d55f5",
"metadata": {},
"outputs": [],
"source": [
"numbabes = resp[\"numbabes\"].replace(97, np.nan)"
]
},
{
"cell_type": "code",
"execution_count": 56,
"id": "7e1b1e04",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.7018914266755958"
]
},
"execution_count": 56,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"deviations = numbabes - numbabes.mean()\n",
"skewness = np.mean(deviations**3) / numbabes.std(ddof=0) ** 3\n",
"skewness"
]
},
{
"cell_type": "markdown",
"id": "79f68b36",
"metadata": {},
"source": [
"In general, a positive value indicates that a distribution is skewed to the right, and a negative value indicates that it is skewed to the left.\n",
"\n",
"If you are given a `Pmf`, rather than a sequence of values, you can compute skewness like this:\n",
"\n",
"1. Compute the deviation of each quantity in the `Pmf` from the mean.\n",
"\n",
"2. Cube the deviations, multiply by the probabilities in the `Pmf`, and add up the products.\n",
"\n",
"3. Divide the sum by the standard deviation cubed.\n",
"\n",
"Write a function called `pmf_skewness` that takes a `Pmf` object and returns its skewness."
]
},
{
"cell_type": "code",
"execution_count": 57,
"id": "95e9d102",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"id": "bbd5e64f",
"metadata": {},
"source": [
"Use your function and the `Pmf` of `numbabes` to compute skewness, and confirm you get the same result."
]
},
{
"cell_type": "code",
"execution_count": 58,
"id": "4e442b4b",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.7018914266755958"
]
},
"execution_count": 58,
"metadata": {},
"output_type": "execute_result"
}
],
"source": []
},
{
"cell_type": "markdown",
"id": "81e26051",
"metadata": {},
"source": [
"**Exercise:** Something like the class size paradox appears if you survey children and ask how many children are in their family.\n",
"Families with many children are more likely to appear in your sample, and families with no children have no chance to be in the sample.\n",
"\n",
"From `resp`, select `numkdhh`, which records the number of children under 18 in each respondent's household.\n",
"Make a `Pmf` of the values in this column.\n",
"\n",
"Use the `bias` function to compute the distribution we would see if we surveyed the children and asked them how many children under 18 (including themselves) are in their household.\n",
"\n",
"Plot the actual and biased distributions, and compute their means."
]
},
{
"cell_type": "code",
"execution_count": 59,
"id": "8a263d8e",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"