{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"tags": [
"remove-cell"
]
},
"outputs": [],
"source": [
"# Reference: https://jupyterbook.org/interactive/hiding.html\n",
"# Use {hide, remove}-{input, output, cell} tags to hiding content\n",
"\n",
"import numpy as np\n",
"import scipy as sp\n",
"import matplotlib.pyplot as plt\n",
"import pandas as pd\n",
"import seaborn as sns\n",
"%matplotlib inline\n",
"import ipywidgets as widgets\n",
"from ipywidgets import interact, interactive, fixed, interact_manual\n",
"from IPython.display import display\n",
"\n",
"sns.set()\n",
"sns.set_context('talk')\n",
"np.set_printoptions(threshold=20, precision=2, suppress=True)\n",
"pd.set_option('display.max_rows', 7)\n",
"pd.set_option('display.max_columns', 8)\n",
"pd.set_option('precision', 2)\n",
"# This option stops scientific notation for pandas\n",
"# pd.set_option('display.float_format', '{:.2f}'.format)\n",
"\n",
"def display_df(df, rows=pd.options.display.max_rows,\n",
" cols=pd.options.display.max_columns):\n",
" with pd.option_context('display.max_rows', rows,\n",
" 'display.max_columns', cols):\n",
" display(df)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(sec:theory_urn)=\n",
"# The Urn Model\n",
"\n",
"The urn model is a simple abstraction of the chance mechanism for drawing indistinguishable marbles from a container, an urn. The randomness in the selection process of drawing marbles from an urn can be extended to many chance processes in real-life examples, and we can simulate this random behavior and use our findings to better understand the accuracy of our data. To explain the urn model, we use a small example with seven marbles. The urn is small enough that we can list all possible outcomes that might result from drawing marbles from the urn.\n",
"\n",
"Let's use the SpaceX Starship prototypes to make up a small example. The protoypyes are called $SN1$, $SN2$, ..., where $SN$ stands for \"serial number\", and in the first half of 2020, seven of these prototypes were built. Before deploying them a few were pressure tested. Suppose we want to select three of the seven Starship prototypes for pressure testing. (While this example is made up, the context is based on the actual SpaceX program; pressure tests were made on the Starship prototypes, $SN1$, $SN2$, ..., $SN7$.) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We set up the urn model as follows: write a unique label on each of seven marbles, place all the marbles in the urn, mix them well, and draw three without looking and without replacing marbles between draws. The urn is small enough that we can list all possible samples of three marbles that can be drawn: \n",
"\n",
"$$ABC ~~ ABD ~~ ABE ~~ ABF ~~ ABG ~~ ACD ~~ ACE \\\\ ACF ~~ ACG ~~ ADE ~~ ADF ~~ ADG ~~ AEF ~~ AEG \\\\ AFG ~~ BCD ~~ BCE ~~ BCF ~~ BCG ~~ BDE ~~ BDF \\\\ BDG ~~ BEF ~~ BEG ~~BFG ~~CDE ~~ CDF ~~ CDG \\\\ CEF ~~ CEG ~~ CFG ~~ DEF ~~ DEG ~~ DFG ~~ EFG $$\n",
"\n",
"We use the labels $A$, $B$, etc. rather than $SN1$, $SN2$, etc. because they are shorter and easier to distinguish. Our list shows that we could wind up with any one of the $35$ unique sets of three from the seven marbles.\n",
"\n",
"We draw an analogy to data scope from {numref}'Chapter %s ': a set of marbles drawn from the urn is a *sample*, and the collection of all marbles placed in the urn is the *population*. This particular urn model prescribes a particular selection method, called the *Simple Random Sample (SRS)*. We describe the SRS and other sampling techniques based on the SRS in the next section. \n",
"\n",
"\n",
"## Sampling Designs\n",
"\n",
"The urn model for the SpaceX prototypes reduced to a few basics. We specified: the number of marbles in the urn; what is written on each marble; the number of marbles drawn from the urn; whether or not they were replaced between draws. This process is equivalent to a Simple Random Sample. Our example is a SRS of three draws from a population of seven, and by design, each of the $35$ samples is equally likely to be chosen because the marbles are indistinguishable and well mixed. This means the chance of any one particular sample must be $1/35$,\n",
"\n",
"$${\\mathbb{P}}(ABC) = {\\mathbb{P}}(\\textrm{ABD}) = \\cdots = {\\mathbb{P}}(\\textrm{EFG}) = \\frac{1}{35}.$$\n",
"\n",
"We use the special symbol ${\\mathbb{P}}$ to stand for \"probability\" or \"chance\", and we read the statement ${\\mathbb{P}}(ABC)$ as the \"chance the sample contains the marbles labeled A, B, and C.\"\n",
"\n",
"We now have a more formal definition of \"representative data\" that is very useful:\n",
"__In a *Simple Random Sample*, every sample has the same chance of being selected.__\n",
"\n",
":::{note}\n",
"\n",
"Many people mistakingly think that the defining property of a SRS is that every unit has an equal chance of being in the sample. However, this is not the case. A SRS of $n$ units from a population of $N$ means that every possible subset of $n$ units has the same chance of being selected.\n",
"\n",
":::\n",
"\n",
"We can use the enumeration of all of the possible samples from the urn to answer additional questions about this chance process. For example, to find the chance that marble $A$ is in the sample, we can add up the chance of all samples that contain $A$. There are 15 of them so the chance is:\n",
"\n",
"$${\\mathbb{P}}(\\textrm{A is in the sample}) = \\frac{15}{35} = \\frac{3}{7}.$$\n",
"\n",
"When it's too difficult to list and count all of the possible samples, we can use simulation to help understand this chance process. \n",
"\n",
"The SRS (and its corresponding urn) is the main building block for more complex survey designs. We briefly describe two of the more widely used designs.\n",
"\n",
"+ *Stratified Sampling* Divide the population into non-overlapping groups, called *strata* (one group is called a *stratum* and more than one are strata), and then take a simple random sample from each. This is like having a separate urn for each stratum and drawing marbles from each urn, independently. The strata do not have to be the same size, and we need not take the same number of units from each.\n",
"\n",
"+ *Cluster Sampling* Divide the population into non-overlapping subgroups (these tend to be smaller than strata), take a simple random sample of the clusters, and include all of the units in a cluster in the sample. We can think on this as a SRS from one urn that contains large marbles that are themselves containers of small marbles. When opened, the sample of large marbles turns into the sample of small marbles.\n",
"\n",
"Often, we are interested in a summary of the sample; that is, some *statistic*. For any sample, we can calculate the statistic, and the urn model helps us find the distribution of possible values that statistic may take on. In the next section, we examine the distribution of a statistic for our example."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Sampling Distribution of a Statistic\n",
"\n",
"Suppose we are interested in whether or not the prototypes can pass a pressure test. It's expensive to carry out the pressure test, so we first test only a sample of prototypes. We can use the urn model to choose the protoypes to be pressure tested, and then, we can summarize our test results by, say, the proportion of prototypes that fail the test. The urn model provides us the knowledge that each of the 35 possible samples has the same chance of being selected and so the pressure test results are representative of the population. \n",
"\n",
"For concreteness, suppose prototypes $A, B, D, and F$ would fail the pressure test, if chosen. For each sample of three marbles, we can find the proportion of failures according to how many of these four defective prototypes are in the sample. Below are a few examples of this calculation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"| | | | | | \n",
"| :--- | :---- | :--- | :--- | :--- | \n",
"| Sample | ABC | BCE | BDF | CEG | \n",
"| Proportion | 2/3 | 1/3 | 1 | 0 | "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since we are drawing three marbles from the urn, the only possible samples proportions are $0$, $1/3$, $2/3$ and $1$, and for each triple, we can calculate its corresponding proprotion. For example, there are 4 samples that give us all failed tests (a sample proportion of 1). These are: $ABD$ , $ABF$, $ADF$, $BDF$, so the chance of observing a sample proportion of $1$ is $4/35$. Below we have summarized the distribution of values for the sample proportion into a table."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"| Proportion of Fails | No. of Samples | Fraction of Samples |\n",
"| :---: | :---: | :---: |\n",
"| 1 | 4 | 4/35 |\n",
"| 2/3 | 18 | 18/35 |\n",
"|1/3 | 12 | 12/35 |\n",
"| 0 | 1 | 1/35 |\n",
"|Total | 35 | 1 |"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"While these calculations are relatively straight forward, we can approximate them through a simulation study. To do this, we take samples of three from our population over and over, say 10,000 times. For each sample, we calculate the proportion of failures. That gives us 10,000 simulated sample proportions. \n",
"The table of the simulated proportions should match the distribution table above. We confirm this with a simulation study. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simulating the Sampling Distribution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our original urn had seven marbles marked $A$ through $G$. However, since we care only whether the prototype fails or passes the test, we can re-label each marble as 'fail' or 'pass'. We create this revised urn as an array. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"urn = ['fail', 'fail', 'fail', 'fail', 'pass', 'pass', 'pass']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We simulate the draw of three marbles from our urn without replacement between draws using numpy's 'random.choice' method as follows. "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array(['fail', 'pass', 'fail'], dtype='"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fracs = counts_els/10000\n",
"\n",
"fig = plt.figure()\n",
"ax = fig.add_subplot(111)\n",
"ax.plot(unique_els, fracs, 'bo')\n",
"ax.vlines(unique_els, 0, fracs, lw=2)\n",
"ax.set_xlabel('Sample Proportion')\n",
"ax.set_ylabel('Fraction of 10,000 simulations')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The simulation results closely match the table. \n",
"\n",
"This simulation study does not *prove*, say, that we expect 18/35 samples to have two fails, but it does give us excellent approximations to our earlier calculations, which is reassuring. More importantly, when we have a more complex setting where it might be difficult to list all possibilities, a simulation study can offer valuable insights. \n",
"\n",
":::{note}\n",
"\n",
"A simulation study repeats a random process many many times. A summary of the patterns that result from the simulation can approximate the theoretical properties of the chance process. This summary is not the same as proving the theoretical properties, but often the guidance we get from the simulation is adequate for our purposes.\n",
"\n",
":::\n",
"\n",
"Drawing marbles from an urn with 0s and 1s is such a popular framework for understanding randomness that this chance process has been given a formal name, the hypergeometric. And, most software provide the functionality to rapidly carry out simulations of the hypergeomteric. We redo our simulation using the hypergeometric to complete this section. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The Hypergeometric\n",
"\n",
"The version of the urn model where we count the number of marbles of a certain type (in our case 'fail' marbles) is so common that there is a random chance process named for it: the hypergeometric. \n",
"Instead of using `random.choice`, we can use numpy's `random.hypergeometric` to simulate drawing marbles from the urn and counting the number of `fails`. The `random.hypergeometric` method is optimzed for the 0-1 urn and allows us to ask for 10,000 simulations in one call. For completeness, we repeat our simulation study and calculate the empirical proportions. "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"simulations_fast = np.random.hypergeometric(ngood=4, nbad=3, nsample=3, size=10000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note: we don't think that a pass is \"bad\"; it's just a naming convention to call the type you want to count \"good\" and the other \"bad\". "
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"unique_els, counts_els = np.unique(np.array( simulations_fast ), return_counts=True)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[0. , 1. , 2. , 3. ],\n",
" [0.03, 0.34, 0.52, 0.11]])"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.array((unique_els, counts_els/10000))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You might have asked yourself already - since the hypergeometric is so popular, why not provide the exact distribution of the possible values. In fact, these are available, and we show how to calculate them below. "
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.03, 0.34, 0.51, 0.11])"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from scipy.stats import hypergeom\n",
"\n",
"x = np.arange(0, 4)\n",
"hypergeom.pmf(x, 7, 4, 3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Perhaps the two most common chance processes are those that arise from counting the number of 1s drawn from a 0-1 urn: drawing without replacement is the *hypergeometric* and drawing with replacement is the *binomial*.\n",
"\n",
":::{note}\n",
"\n",
"Whenever possible, it's a good idea to use the functionality provided in a third party package for simulating from a named distribution, rather than writing your own function, such as the random number generators offered in numpy. It's best to take advanatge of efficient and accurate code that others have devloped.\n",
"\n",
":::\n",
"\n",
"While this simulation was simple, so simple that we could have used `hypergeom.pmf` to complete our distribution, we wanted to demonstrate the intuition that a simulation study can reveal. The approach we take in this book is to develop understanding about chance processes based on simulation studies. However, we do formalize the notion of a probability distribution of a statistics (like the proportion of fails in a sample) in\n",
"{numref}`Section %s `. \n",
"\n",
"Now that we have simulation as a tool for understanding accuracy, we can revisit the election example from {numref}`Chapter %s ` and carry out a post-election study of what might have gone wrong with the voter polls. This simulation study imitates drawing more than a thousand marbles (voters who participate in the poll) from an urn of six million. We can examine potential sources of bias and the variation in the polling results, and carry out a what-if analysis, where we examine how the predictions might have gone if even a larger number of draws from the urn were taken. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"celltoolbar": "Tags",
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.4"
}
},
"nbformat": 4,
"nbformat_minor": 4
}