In [1]:

```
# Set up packages for lecture. Don't worry about understanding this code,
# but make sure to run it if you're following along.
import numpy as np
import babypandas as bpd
import pandas as pd
from matplotlib_inline.backend_inline import set_matplotlib_formats
import matplotlib.pyplot as plt
set_matplotlib_formats("png")
plt.style.use('ggplot')
np.set_printoptions(threshold=20, precision=2, suppress=True)
pd.set_option("display.max_rows", 7)
pd.set_option("display.max_columns", 8)
pd.set_option("display.precision", 2)
def draw_cutoff(results):
ax = bpd.DataFrame().assign(results=results).plot(kind='hist', bins=np.arange(160, 240, 4),
density=True, figsize=(10, 5),
title='Empirical Distribution of the Number of Heads in 400 Flips of a Fair Coin')
for bar in ax.containers[0]:
x = bar.get_x() + 0.5 * bar.get_width()
if x < 184:
bar.set_color('#796fb3')
plt.annotate('likely biased\ntowards tails', (160, 0.007), size=16, color='#796fb3', weight='bold')
plt.annotate('likely fair', (225, 0.008), size=16, color='#e24a33', weight='bold');
```

- Homework 5 is due
**tomorrow at 11:59PM**. - Lab 6 is due
**Monday 11/27 at 11:59PM**. - Happy Thanksgiving! 🦃
- There will be lecture and discussion section on Wednesday as usual, but there is no quiz.
- There will be no office hours on Thursday, Friday, or Saturday.

- In the project, you'll explore data from NASA about meteorites. You'll look at where meteorites landed and where they were seen falling, and explore whether there is any relationship between these locations and the locations where people tend to live. For example, here are the places where meteorites have been seen falling.
- If you want to work with a partner, follow these guidelines and get started soon!

- Example: Is our coin fair?
- p-values and conventions for consistency.

- Example: Midterm scores.
- Does this sample look like it was drawn from this population?

- Example: Jury selection in Alameda County.
- Total variation distance.

Last time we looked at an example where we found a coin on the ground, flipped it 400 times, and used the results to determine whether the coin was likely fair.

In [2]:

```
# The results of our 400 flips.
flips_400 = bpd.read_csv('data/flips-400.csv')
flips_400.groupby('outcome').count()
```

Out[2]:

flip | |
---|---|

outcome | |

Heads | 188 |

Tails | 212 |

**Question**: Does our coin look like a fair coin, or not?

How "weird" is it to flip a fair coin 400 times and see only 188 heads?

- The
**hypotheses**describe two views of how our data was generated.- Remember, the null hypothesis needs to be a well-defined probability model about how the data was generated, so that we can use it for simulation.

- One pair of hypotheses is:
**Null Hypothesis:**The coin is fair.**Alternative Hypothesis:**The coin is not fair.

- A different pair is:
**Null Hypothesis:**The coin is fair.**Alternative Hypothesis:**The coin is biased towards tails.

- We plan to simulate 400 flips of a fair coin, many times, and keep track of a single number, called a
**test statistic**, each time.

- Our test statistic should be chosen such that
**high observed values lean towards one hypothesis and low observed values lean towards the other**.

The purpose of a test statistic is to help you distinguish between the two hypotheses.

The test statistic you choose depends heavily on the pair of hypotheses you are testing. That said, these general guidelines may help.

- Once we have stated our hypotheses and chosen an appropriate test statistic, we repeatedly generate samples under the assumption that the null hypothesis is true and record the value of the test statistic for each generated sample.

Here, we'll repeatedly flip a fair coin 400 times.

**Suppose we're choosing between "the coin is fair" and "the coin is biased towards tails".**

Here, our test statistic is the number of heads in 400 flips.

In [3]:

```
# Computes a single simulated test statistic.
np.random.multinomial(400, [0.5, 0.5])[0]
```

Out[3]:

194

In [4]:

```
# Computes 10,000 simulated test statistics.
results = np.array([])
for i in np.arange(10000):
result = np.random.multinomial(400, [0.5, 0.5])[0]
results = np.append(results, result)
results
```

Out[4]:

array([206., 205., 191., ..., 201., 211., 200.])

Let's visualize the empirical distribution of the test statistic $\text{number of heads}$.

In [5]:

```
bpd.DataFrame().assign(results=results).plot(kind='hist', bins=np.arange(160, 240, 4),
density=True, ec='w', figsize=(10, 5),
title='Empirical Distribution of the Number of Heads in 400 Flips of a Fair Coin');
plt.legend();
```

- Our hypothesis test boils down to checking
**whether our observed statistic is a "typical value" in the distribution of our test statistic.**- If we see a small number of heads, we'll
**reject the null hypothesis**in favor of the alternative hypothesis, "the coin is biased towards tails." - Otherwise we'll
**fail to reject the null hypothesis**because it's plausible that "the coin is fair."

- If we see a small number of heads, we'll

**Question:**Where do we draw the line between rejecting the null and failing to reject the null? How many heads, exactly, is considered a small number of heads?

In [6]:

```
bpd.DataFrame().assign(results=results).plot(kind='hist', bins=np.arange(160, 240, 4),
density=True, ec='w', figsize=(10, 5),
title='Empirical Distribution of the Number of Heads in 400 Flips of a Fair Coin');
plt.legend();
```

- If we were to observe, say, 170 heads, even though such a result
*could*happen with a fair coin, seeing 170 heads from a fair coin is sufficiently rare that we'd think the coin was biased towards tails.

- Let's say an outcome is
**sufficiently rare**if it falls in the lowest**five percent**of outcomes in our simulation under the assumptions of the null hypothesis.- Five percent could be something else - five is just a nice round number.

In [7]:

```
np.percentile(results, 5)
```

Out[7]:

184.0

In [8]:

```
draw_cutoff(results)
```

- Now we have a clear cutoff that tells us which hypothesis to side with.
- If we find a coin, flip it 400 times, and observe less than 184 heads, we'll side with "the coin is biased towards tails."
- If we find a coin, flip it 400 times, and observe 184 heads or more, we'll side with "the coin is fair."

- To quantify how rare our observed statistic is, under the null hypothesis, we can compute what’s called a
**p-value**.

- The p-value is defined as the probability, under the null hypothesis, that the test statistic
**is equal**to the value that was observed in the data**or is even further in the direction of the alternative**.

Its formal name is the observed significance level.

- p-values correspond to the "tail areas" of a histogram, starting at the observed statistic.

In [9]:

```
bpd.DataFrame().assign(results=results).plot(kind='hist', bins=np.arange(160, 240, 4),
density=True, ec='w', figsize=(10, 5),
title='Empirical Distribution of the Number of Heads in 400 Flips of a Fair Coin');
plt.legend();
```

- For example, it's extremely rare to flip a fair coin 400 times and see 170 heads or fewer.

In [10]:

```
np.count_nonzero(results <= 170) / len(results)
```

Out[10]:

0.0014

- But it's much less rare to flip a fair coin 400 times and see 195 heads or fewer.

In [11]:

```
np.count_nonzero(results <= 195) / len(results)
```

Out[11]:

0.3275

- The larger the p-value, the more likely our observation is, if the null hypothesis is true. Therefore,
- Larger p-values mean we should fail to reject the null.
- Smaller p-values mean we should reject the null.

- If the p-value is sufficiently large, we say the data is
**consistent**with the null hypothesis and so we "**fail to reject the null hypothesis**".- We never say that we "accept" the null hypothesis! The null hypothesis may be plausible, but there are many other possible explanations for our data.

- If the p-value is below some cutoff, we say the data is
**inconsistent**with the null hypothesis, and we**"reject the null hypothesis"**.- If a p-value is less than 0.05, the result is said to be "statistically significant".
- If a p-value is less than 0.01, the result is said to be "highly statistically significant".
- These conventions are historical and completely arbitrary! (And controversial.)

- Recall that we found a coin, flipped it 400 times, and saw 188 heads. This is the
**observed statistic**.

In [12]:

```
bpd.DataFrame().assign(results=results).plot(kind='hist', bins=np.arange(160, 240, 4),
density=True, ec='w', figsize=(10, 5),
title='Empirical Distribution of the Number of Heads in 400 Flips of a Fair Coin');
plt.axvline(188, color='black', linewidth=4, label='observed statistic (188)')
plt.legend();
```

In [13]:

```
np.count_nonzero(results <= 188) / len(results)
```

Out[13]:

0.1258

- It happens about 12% of the time. It's not so rare. Since the p-value is at least 0.05, we
**fail to reject**the null hypothesis at the standard 0.05**significance level**and conclude that it's plausible that our coin is fair.

- Unfortunately, we'll never know. 🤷♂️

- In this way, we can interpret our p-value cutoff as a probability of error.

- In this example, we calculated the p-value as the area in the
**left**tail of the histogram.

- But other times, when large values of the observed statistic indicate the alternative hypothesis, the p-value corresponds to the area in the
**right**tail.

- We always calculate a p-value with
`>=`

or`<=`

because it includes the probability that the test statistic exactly equals the observed statistic.

- In Fall 2022, there were four sections of DSC 10 – A, B, C, and D.

- One of the four sections – the one taught by Suraj – had a much lower average than the other three sections.
- All midterms were graded by the same people, with the same rubrics.

In [14]:

```
# Midterm scores from DSC 10, Fall 2022, slightly perturbed for anonymity.
scores = bpd.read_csv('data/fa22-midterm-scores.csv')
scores
```

Out[14]:

Section | Score | |
---|---|---|

0 | A | 54.5 |

1 | D | 62.0 |

2 | B | 23.5 |

... | ... | ... |

390 | C | 62.5 |

391 | B | 47.5 |

392 | D | 72.5 |

393 rows × 2 columns

In [15]:

```
scores.plot(kind='hist', density=True, figsize=(10, 5), ec='w', bins=np.arange(10, 85, 5), title='Distribution of Midterm Exam Scores in DSC 10');
```

In [16]:

```
# Total number of students who took the exam.
scores.shape[0]
```

Out[16]:

393

In [17]:

```
# Calculate the number of students in each section.
scores.groupby('Section').count()
```

Out[17]:

Score | |
---|---|

Section | |

A | 117 |

B | 115 |

C | 108 |

D | 53 |

In [18]:

```
# Calculate the average midterm score in each section.
scores.groupby('Section').mean()
```

Out[18]:

Score | |
---|---|

Section | |

A | 51.54 |

B | 49.48 |

C | 46.17 |

D | 50.88 |

- Suppose we
*randomly*place all 393 students into one of four sections and compute the average midterm score within each section.

- One of the sections would
*have to*have a lower average score than the others (unless multiple are tied for the lowest).- In any set of four numbers, one of them has to be the minimum!

- But is Section C's average
*lower*than we'd expect due to chance? Let's perform a hypothesis test!

**Null Hypothesis:**Section C's scores are drawn randomly from the distribution of scores in the course overall. The observed difference between the average score in Section C and the average score in the course overall is due to random chance.

**Alternative Hypothesis:**Section C's average score is too low to be explained by chance alone.

- Even the fair coin example can be posed in this way.
- In this case, the population distribution is the uniform distribution (50% heads, 50% tails), and when we flip a coin, we are sampling from this distribution with replacement.
- We wanted to know if the sample of coin flips we observed looked like it came from this population, which represents a fair coin.

In [19]:

```
section_size = scores.groupby('Section').count().get('Score').loc['C']
observed_avg = scores.groupby('Section').mean().get('Score').loc['C']
print(f'Section C had {section_size} students and an average midterm score of {observed_avg}.')
```

Section C had 108 students and an average midterm score of 46.1712962962963.

- Model: There is no significant difference between the exam scores in different sections. Section C had a lower average purely due to chance.
- To simulate: sample 108 students uniformly at random without replacement from the class.

- Test statistic: the average midterm score of a section.
- The observed statistic is the average midterm score of Section C (about 46.17).

In [20]:

```
# Sample 108 students from the class, independent of section,
# and compute the average score.
scores.sample(int(section_size), replace=False).get('Score').mean()
```

Out[20]:

49.03703703703704

In [21]:

```
averages = np.array([])
repetitions = 1000
for i in np.arange(repetitions):
random_sample = scores.sample(int(section_size), replace=False)
new_average = random_sample.get('Score').mean()
averages = np.append(averages, new_average)
averages
```

Out[21]:

array([47.92, 49.72, 48.28, ..., 49.95, 51.04, 48.64])

In [22]:

```
observed_avg
```

Out[22]:

46.1712962962963

In [23]:

```
bpd.DataFrame().assign(RandomSampleAverage=averages).plot(kind='hist', bins=20,
density=True, ec='w', figsize=(10, 5),
title='Empirical Distribution of Midterm Averages for 108 Randomly Selected Students')
plt.axvline(observed_avg, color='black', linewidth=4, label='observed statistic')
plt.legend();
```

In [24]:

```
# p-value
np.count_nonzero(averages <= observed_avg) / repetitions
```

Out[24]:

0.003

**reject the null hypothesis.** It's not looking good for you, Suraj!

Recall from Lecture 19:

Section 197 of California's Code of Civil Procedure says,

"All persons selected for jury service shall be selected at random, from a source or sources inclusive of a representative cross section of the population of the area served by the court."

- 1453 people reported for jury duty in total (we will call them "panelists").

In [25]:

```
jury = bpd.DataFrame().assign(
Ethnicity=['Asian', 'Black', 'Latino', 'White', 'Other'],
Eligible=[0.15, 0.18, 0.12, 0.54, 0.01],
Panels=[0.26, 0.08, 0.08, 0.54, 0.04]
)
jury
```

Out[25]:

Ethnicity | Eligible | Panels | |
---|---|---|---|

0 | Asian | 0.15 | 0.26 |

1 | Black | 0.18 | 0.08 |

2 | Latino | 0.12 | 0.08 |

3 | White | 0.54 | 0.54 |

4 | Other | 0.01 | 0.04 |

What do you notice? 👀

**Null Hypothesis:**Panelists were selected at random from the eligible population.

**Alternative Hypothesis:**Panelists were*not*selected at random from the eligible population.

- Observation: 1453 panelists and the distribution of their ethnicities.

- Test statistic: ???
- How do we deal with multiple categories?

In [26]:

```
jury.plot(kind='barh', x='Ethnicity', figsize=(10, 5));
```

- Panelists are categorized into one of 5 ethnicities. In other words, ethnicity is a
**categorical**variable.

- To see whether the the distribution of ethnicities for the panelists is similar to that of the eligible population, we have to measure the distance between two categorical distributions.
- We've done this for distributions with just two categories – heads and tails, for instance – but not when there are more than two categories.

- Let's start by considering the difference in proportions for each category.

In [27]:

```
with_diffs = jury.assign(Difference=(jury.get('Panels') - jury.get('Eligible')))
with_diffs
```

Out[27]:

Ethnicity | Eligible | Panels | Difference | |
---|---|---|---|---|

0 | Asian | 0.15 | 0.26 | 0.11 |

1 | Black | 0.18 | 0.08 | -0.10 |

2 | Latino | 0.12 | 0.08 | -0.04 |

3 | White | 0.54 | 0.54 | 0.00 |

4 | Other | 0.01 | 0.04 | 0.03 |

- Note that if we sum these differences, the result is 0 (you'll see the proof in DSC 40A).
- To avoid cancellation of positive and negative differences, we can take the absolute value of these differences.

In [28]:

```
with_abs_diffs = with_diffs.assign(AbsoluteDifference=np.abs(with_diffs.get('Difference')))
with_abs_diffs
```

Out[28]:

Ethnicity | Eligible | Panels | Difference | AbsoluteDifference | |
---|---|---|---|---|---|

0 | Asian | 0.15 | 0.26 | 0.11 | 0.11 |

1 | Black | 0.18 | 0.08 | -0.10 | 0.10 |

2 | Latino | 0.12 | 0.08 | -0.04 | 0.04 |

3 | White | 0.54 | 0.54 | 0.00 | 0.00 |

4 | Other | 0.01 | 0.04 | 0.03 | 0.03 |

The **Total Variation Distance (TVD)** of two categorical distributions is **the sum of the absolute differences of their proportions, all divided by 2**.

- We divide by 2 so that, for example, the distribution [0.51, 0.49] is 0.01 away from [0.50, 0.50].

- It would also be valid not to divide by 2. We just wouldn't call that statistic TVD anymore.

In [29]:

```
with_abs_diffs
```

Out[29]:

Ethnicity | Eligible | Panels | Difference | AbsoluteDifference | |
---|---|---|---|---|---|

0 | Asian | 0.15 | 0.26 | 0.11 | 0.11 |

1 | Black | 0.18 | 0.08 | -0.10 | 0.10 |

2 | Latino | 0.12 | 0.08 | -0.04 | 0.04 |

3 | White | 0.54 | 0.54 | 0.00 | 0.00 |

4 | Other | 0.01 | 0.04 | 0.03 | 0.03 |

In [30]:

```
with_abs_diffs.get('AbsoluteDifference').sum() / 2
```

Out[30]:

0.14

In [31]:

```
def total_variation_distance(dist1, dist2):
'''Computes the TVD between two categorical distributions,
assuming the categories appear in the same order.'''
return np.abs((dist1 - dist2)).sum() / 2
jury.plot(kind='barh', x='Ethnicity', figsize=(10, 5))
plt.annotate('If you add up the total amount by which the blue bars\n are longer than the red bars, you get 0.14.', (0.08, 3.9), bbox=dict(boxstyle="larrow,pad=0.3", fc="#e5e5e5", ec="black", lw=2));
plt.annotate('If you add up the total amount by which the red bars\n are longer than the blue bars, you also get 0.14!', (0.23, 0.9), bbox=dict(boxstyle="larrow,pad=0.3", fc="#e5e5e5", ec="black", lw=2));
```