Cover Photo by Nasonov Aleksandr on Unsplash

## Overview

This is a continuation of my progress through Data Science from Scratch by Joel Grus. We'll use a classic coin-flipping example in this post because it is simple to illustrate with both **concept** and **code**. The goal of this post is to connect the dots between several concepts including the Central Limit Theorem, hypothesis testing, p-Values and confidence intervals, using python to build our intuition.

## Central Limit Theorem

Terms like "null" and "alternative" hypothesis are used quite frequently, so let's set some context. The "null" is the **default** position. The "alternative", alt for short, is something we're *comparing to* the default (null).

The classic coin-flipping exercise is to test the *fairness* off a coin. If a coin is fair, it'll land on heads 50% of the time (and tails 50% of the time). Let's translate into hypothesis testing language:

**Null Hypothesis**: Probability of landing on Heads = 0.5.

**Alt Hypothesis**: Probability of landing on Heads != 0.5.

Each coin flip is a **Bernoulli trial**, which is an experiment with two outcomes - outcome 1, "success", (probability *p*) and outcome 0, "fail" (probability *p - 1*). The reason it's a Bernoulli trial is because there are only two outcome with a coin flip (heads or tails). Read more about Bernoulli here.

Here's the code for a single Bernoulli Trial:

```
def bernoulli_trial(p: float) -> int:
"""Returns 1 with probability p and 0 with probability 1-p"""
return 1 if random.random() < p else 0
```

When you **sum the independent Bernoulli trials**, you get a **Binomial(n,p)** random variable, a variable whose *possible* values have a probability distribution. The **central limit theorem** says as **n** or the *number* of independent Bernoulli trials get large, the Binomial distribution approaches a normal distribution.

Here's the code for when you sum all the Bernoulli Trials to get a Binomial random variable:

```
def binomial(n: int, p: float) -> int:
"""Returns the sum of n bernoulli(p) trials"""
return sum(bernoulli_trial(p) for _ in range(n))
```

**Note**: A single 'success' in a Bernoulli trial is 'x'. Summing up all those x's into X, is a Binomial random variable. Success doesn't imply desirability, nor does "failure" imply undesirability. They're just terms to count the cases we're looking for (i.e., number of heads in multiple coin flips to assess a coin's fairness).

Given that our **null** is (p = 0.5) and **alt** is (p != 0.5), we can run some independent bernoulli trials, then sum them up to get a binomial random variable.

Each `bernoulli_trial`

is an experiment with either 0 or 1 as outcomes. The `binomial`

function sums up **n** bernoulli(0.5) trails. We ran both twice and got different results. Each bernoulli experiment can be a success(1) or faill(0); summing up into a binomial random variable means we're taking the probability p(0.5) *that a coin flips head* and we ran the experiment 1,000 times to get a random binomial variable.

The first 1,000 flips we got 510. The second 1,000 flips we got 495. We can repeat this process many times to get a *distribution*. We can plot this distribution to reinforce our understanding. To this we'll use `binomial_histogram`

function. This function picks points from a Binomial(n,p) random variable and plots their histogram.

```
from collections import Counter
import matplotlib.pyplot as plt
def normal_cdf(x: float, mu: float = 0, sigma: float = 1) -> float:
return (1 + math.erf((x - mu) / math.sqrt(2) / sigma)) / 2
def binomial_histogram(p: float, n: int, num_points: int) -> None:
"""Picks points from a Binomial(n, p) and plots their histogram"""
data = [binomial(n, p) for _ in range(num_points)]
# use a bar chart to show the actual binomial samples
histogram = Counter(data)
plt.bar([x - 0.4 for x in histogram.keys()],
[v / num_points for v in histogram.values()],
0.8,
color='0.75')
mu = p * n
sigma = math.sqrt(n * p * (1 - p))
# use a line chart to show the normal approximation
xs = range(min(data), max(data) + 1)
ys = [normal_cdf(i + 0.5, mu, sigma) -
normal_cdf(i - 0.5, mu, sigma) for i in xs]
plt.plot(xs, ys)
plt.title("Binomial Distribution vs. Normal Approximation")
plt.show()
# call function
binomial_histogram(0.5, 1000, 10000)
```

This plot is then rendered:

What we did was sum up independent `bernoulli_trial`

(s) of 1,000 coin flips, where the probability of head is p = 0.5, to create a `binomial`

random variable. We then repeated this a large number of times (N = 10,000), then plotted a histogram of the distribution of all binomial random variables. And because we did it so many times, it approximates the standard normal distribution (smooth bell shape curve).

Just to demonstrate how this works, we can generate several `binomial`

random variables:

If we do this 10,000 times, we'll generate the above histogram. You'll notice that because we are testing whether the coin is fair, the probability of heads (success) *should* be at 0.5 and, from 1,000 coin flips, the **mean**(`mu`

) should be a 500.

We have another function that can help us calculate `normal_approximation_to_binomial`

:

```
import random
from typing import Tuple
import math
def normal_approximation_to_binomial(n: int, p: float) -> Tuple[float, float]:
"""Return mu and sigma corresponding to a Binomial(n, p)"""
mu = p * n
sigma = math.sqrt(p * (1 - p) * n)
return mu, sigma
# call function
# (500.0, 15.811388300841896)
normal_approximation_to_binomial(1000, 0.5)
```

When calling the function with our parameters, we get a mean `mu`

of 500 (from 1,000 coin flips) and a standard deviation `sigma`

of 15.8114. Which means that 68% of the time, the binomial random variable will be 500 +/- 15.8114 and 95% of the time it'll be 500 +/- 31.6228 (see 68-95-99.7 rule)

## Hypothesis Testing

Now that we have seen the results of our "coin fairness" experiment plotted on a binomial distribution (approximately normal), we will be, for the purpose of testing our hypothesis, be interested in the probability of its realized value (binomial random variable) lies **within or outside a particular interval**.

This means we'll be interested in questions like:

- What's the probability that the binomial(n,p) is below a threshold?
- Above a threshold?
- Between an interval?
- Outside an interval?

First, the `normal_cdf`

(normal cummulative distribution function), which we learned in a previous post, *is* the probability of a variable being *below* a certain threshold.

Here, the probability of X (success or heads for a 'fair coin') is at 0.5 (`mu`

= 500, `sigma`

= 15.8113), and we want to find the probability that X falls below 490, which comes out to roughly 26%

```
def normal_cdf(x: float, mu: float = 0, sigma: float = 1) -> float:
return (1 + math.erf((x - mu) / math.sqrt(2) / sigma)) / 2
normal_probability_below = normal_cdf
# probability that binomal random variable, mu = 500, sigma = 15.8113, is below 490
# 0.26354347477247553
normal_probability_below(490, 500, 15.8113)
```

On the other hand, the `normal_probability_above`

, probability that X falls *above* 490 would be

1 - 0.2635 = 0.7365 or roughly 74%.

```
def normal_probability_above(lo: float,
mu: float = 0,
sigma: float = 1) -> float:
"""The probability that an N(mu, sigma) is greater than lo."""
return 1 - normal_cdf(lo, mu, sigma)
# 0.7364565252275245
normal_probability_above(490, 500, 15.8113)
```

To make sense of this we need to recall the binomal distribution, that approximates the normal distribution, but we'll draw a vertical line at 490.

We're asking, given the binomal distribution with `mu`

500 and `sigma`

at 15.8113, what is the probability that a binomal random variable falls below the threshold (left of the line); the answer is approximately 26% and correspondingly falling above the threshold (right of the line), is approximately 74%.

### Between an Interval

We may also wonder what the probability of a binomial random variable **falling between 490 and 520**:

Here is the function to calculate this probability and it comes out to approximately 63%. *note*: Bear in mind the full area under the curve is 1.0 or 100%.

```
def normal_probability_between(lo: float,
hi: float,
mu: float = 0,
sigma: float = 1) -> float:
"""The probability that an N(mu, sigma) is between lo and hi."""
return normal_cdf(hi, mu, sigma) - normal_cdf(lo, mu, sigma)
# 0.6335061861416337
normal_probability_between(490, 520, 500, 15.8113)
```

Finally, the area outside of the interval should be 1 - 0.6335 = 0.3665:

```
def normal_probability_outside(lo: float,
hi: float,
mu: float = 0,
sigma: float = 1) -> float:
"""The probability that an N(mu, sigma) is not between lo and hi."""
return 1 - normal_probability_between(lo, hi, mu, sigma)
# 0.3664938138583663
normal_probability_outside(490, 520, 500, 15.8113)
```

In addition to the above, we may also be interested in finding (symmetric) intervals around the mean that account for a *certain level of likelihood*, for example, 60% probability centered around the mean.

For this operation we would use the `inverse_normal_cdf`

:

```
def inverse_normal_cdf(p: float,
mu: float = 0,
sigma: float = 1,
tolerance: float = 0.00001) -> float:
"""Find approximate inverse using binary search"""
# if not standard, compute standard and rescale
if mu != 0 or sigma != 1:
return mu + sigma * inverse_normal_cdf(p, tolerance=tolerance)
low_z = -10.0 # normal_cdf(-10) is (very close to) 0
hi_z = 10.0 # normal_cdf(10) is (very close to) 1
while hi_z - low_z > tolerance:
mid_z = (low_z + hi_z) / 2 # Consider the midpoint
mid_p = normal_cdf(mid_z) # and the CDF's value there
if mid_p < p:
low_z = mid_z # Midpoint too low, search above it
else:
hi_z = mid_z # Midpoint too high, search below it
return mid_z
```

First we'd have to find the cutoffs where the upper and lower tails each contain 20% of the probability. We calculate `normal_upper_bound`

and `normal_lower_bound`

and use those to calculate the `normal_two_sided_bounds`

.

```
def normal_upper_bound(probability: float,
mu: float = 0,
sigma: float = 1) -> float:
"""Returns the z for which P(Z <= z) = probability"""
return inverse_normal_cdf(probability, mu, sigma)
def normal_lower_bound(probability: float,
mu: float = 0,
sigma: float = 1) -> float:
"""Returns the z for which P(Z >= z) = probability"""
return inverse_normal_cdf(1 - probability, mu, sigma)
def normal_two_sided_bounds(probability: float,
mu: float = 0,
sigma: float = 1) -> Tuple[float, float]:
"""
Returns the symmetric (about the mean) bounds
that contain the specified probability
"""
tail_probability = (1 - probability) / 2
# upper bound should have tail_probability above it
upper_bound = normal_lower_bound(tail_probability, mu, sigma)
# lower bound should have tail_probability below it
lower_bound = normal_upper_bound(tail_probability, mu, sigma)
return lower_bound, upper_bound
```

So if we wanted to know what the cutoff points were for a 60% probability around the mean and standard deviation (`mu`

= 500, `sigma`

= 15.8113), it would be between **486.69 and 513.31**.

Said differently, this means roughly 60% of the time, we can expect the binomial random variable to fall between 486 and 513.

```
# (486.6927811021805, 513.3072188978196)
normal_two_sided_bounds(0.60, 500, 15.8113)
```

### Significance and Power

Now that we have a handle on the binomial normal distribution, thresholds (left and right of the mean), and cut-off points, we want to make a **decision about significance**. Probably the most important part of *statistical significance* is that it is a decision to be made, not a standard that is externally set.

Significance is a decision about how willing we are to make a *type 1* error (false positive), which we explored in a previous post. The convention is to set it to a 5% or 1% willingness to make a type 1 error. Suppose we say 5%.

We would say that out of 1,000 coin flips, 95% of the time, we'd get between 469 and 531 heads on a "fair coin" and 5% of the time, outside of this 469-531 range.

```
# (469.0104394712448, 530.9895605287552)
normal_two_sided_bounds(0.95, 500, 15.8113)
```

If we recall our hypotheses:

**Null Hypothesis**: Probability of landing on Heads = 0.5 (fair coin)

**Alt Hypothesis**: Probability of landing on Heads != 0.5 (biased coin)

Each binomial distribution (test) that consist of 1,000 bernoulli trials, each *test* where the number of heads falls outside the range of 469-531, we'll **reject the null** that the coin is fair. And we'll be wrong (false positive), 5% of the time. It's a false positive when we **incorrectly reject** the null hypothesis, when it's actually true.

We also want to avoid making a type-2 error (false negative), where we **fail to reject** the null hypothesis, when it's actually false.

**Note**: Its important to keep in mind that terms like *significance* and *power* are used to describe **tests**, in our case, the test of whether a coin is fair or not. Each test is the sum of 1,000 independent bernoulli trials.

For a "test" that has a 95% significance, we'll assume that out of a 1,000 coin flips, it'll land on heads between 469-531 times and we'll determine the coin is fair. For the 5% of the time it lands outside of this range, we'll determine the coin to be "unfair", but we'll be wrong because it actually is fair.

To calculate the power of the test, we'll take the assumed `mu`

and `sigma`

with a 95% bounds (based on the assumption that the probability of the coin landing on heads is 0.5 or 50% - a fair coin). We'll determine the lower and upper bounds:

```
lo, hi = normal_two_sided_bounds(0.95, mu_0, sigma_0)
lo # 469.01026640487555
hi # 530.9897335951244
```

And if the coin was *actually biased*, we should reject the null, but we fail to. Let's suppose the actual probability that the coin lands on heads is 55% ( **biased** towards head):

```
mu_1, sigma_1 = normal_approximation_to_binomial(1000, 0.55)
mu_1 # 550.0
sigma_1 # 15.732132722552274
```

Using the same range 469 - 531, where the coin is assumed 'fair' with `mu`

at 500 and `sigma`

at 15.8113:

If the coin, in fact, had a bias towards head (p = 0.55), the distribution would shift right, but if our 95% significance test remains the same, we get:

The probability of making a type-2 error is 11.345%. This is the probability that we're see that the coin's distribution is within the previous interval 469-531, thinking we should accept the null hypothesis (that the coin is fair), but in actuality, failing to see that the distribution has shifted to the coin having a *bias* towards heads.

```
# 0.11345199870463285
type_2_probability = normal_probability_between(lo, hi, mu_1, sigma_1)
```

The other way to arrive at this is to find the probability, under the *new* `mu`

and `sigma`

(new distribution), that X (number of successes) will fall *below* 531.

```
# 0.11357762975476304
normal_probability_below(531, mu_1, sigma_1)
```

So the probability of making a type-2 error or the probability that the *new* distribution falls below 531 is approximately 11.3%.

The **power to detect** a type-2 error is 1.00 minus the probability of a type-2 error (1 - 0.113 = 0.887), or 88.7%.

```
power = 1 - type_2_probability # 0.8865480012953671
```

Finally, we may be interested in **increasing power** to detect a type-2 error. Instead of using a `normal_two_sided_bounds`

function to find the cut-off points (i.e., 469 and 531), we could use a *one-sided test* that rejects the null hypothesis ('fair coin') when X (number of heads on a coin-flip) is much larger than 500.

Here's the code, using `normal_upper_bound`

:

```
# 526.0073585242053
hi = normal_upper_bound(0.95, mu_0, sigma_0)
```

This means shifting the upper bounds from 531 to 526, providing more probability in the upper tail. This means the probability of a type-2 error goes down from 11.3 to 6.3.

```
# previous probability of type-2 error
# 0.11357762975476304
normal_probability_below(531, mu_1, sigma_1)
# new probability of type-2 error
# 0.06356221447122662
normal_probability_below(526, mu_1, sigma_1)
```

And the new (stronger) **power to detect** type-2 error is 1.0 - 0.064 = 0.936 or 93.6% (up from 88.7% above).

## p-Values

p-Values represent *another way* of deciding whether to accept or reject the Null Hypothesis. Instead of choosing bounds, thresholds or cut-off points, we could compute the probability, assuming the Null Hypothesis is true, that we would see a value *as extreme as* the one we just observed.

Here is the code:

```
def two_sided_p_values(x: float, mu: float = 0, sigma: float = 1) -> float:
"""
How likely are we to see a value at least as extreme as x (in either
direction) if our values are from an N(mu, sigma)?
"""
if x >= mu:
# x is greater than the mean, so the tail is everything greater than x
return 2 * normal_probability_above(x, mu, sigma)
else:
# x is less than the mean, so the tail is everything less than x
return 2 * normal_probability_below(x, mu, sigma)
```

If we wanted to compute, assuming we have a "fair coin" (`mu`

= 500, `sigma`

= 15.8113), what is the probability of seeing a value like 530? (**note**: We use 529.5 instead of 530 below due to continuity correction)

Answer: approximately 6.2%

```
# 0.06207721579598835
two_sided_p_values(529.5, mu_0, sigma_0)
```

The p-value, 6.2% is higher than our (hypothetical) 5% significance, so we don't reject the null. On the other hand, if X was slightly more extreme, 532, the probability of seeing that value would be approximately 4.3%, which is less than 5% significance, so we would reject the null.

```
# 0.04298479507085862
two_sided_p_values(532, mu_0, sigma_0)
```

For one-sided tests, we would use the `normal_probability_above`

and `normal_probability_below`

functions created above:

```
upper_p_value = normal_probability_above
lower_p_value = normal_probability_below
```

Under the `two_sided_p_values`

test, the extreme value of 529.5 had a probability of 6.2% of showing up, but not low enough to reject the null hypothesis.

However, with a one-sided test, `upper_p_value`

for the same threshold is now 3.1% and we would reject the null hypothesis.

```
# 0.031038607897994175
upper_p_value(529.5, mu_0, sigma_0)
```

## Confidence Intervals

A *third* approach to deciding whether to accept or reject the null is to use confidence intervals. We'll use the 530 as we did in the p-Values example.

```
p_hat = 530/1000
mu = p_hat
sigma = math.sqrt(p_hat * (1 - p_hat) / 1000) # 0.015782902141241326
# (0.4990660982192851, 0.560933901780715)
normal_two_sided_bounds(0.95, mu, sigma)
```

The confidence interval for a coin flipping heads 530 (out 1,000) times is (0.4991, 0.5609). Since this interval **contains** the p = 0.5 (probability of heads 50% of the time, assuming a fair coin), we do not reject the null.

If the extreme value were *more* extreme at 540, we would arrive at a different conclusion:

```
p_hat = 540/1000
mu = p_hat
sigma = math.sqrt(p_hat * (1 - p_hat) / 1000)
(0.5091095927295919, 0.5708904072704082)
normal_two_sided_bounds(0.95, mu, sigma)
```

Here we would be 95% confident that the mean of this distribution is contained between 0.5091 and 0.5709 and this **does not** contain 0.500 (albiet by a slim margin), so we reject the null hypothesis that this is a fair coin.

**note**: Confidence intervals are about the *interval* not probability p. We interpret the confidence interval as, if you were to repeat the experiment many times, 95% of the time, the "true" parameter, in our example p = 0.5, would lie within the observed confidence interval.

## Connecting the Dots

We used several python functions to build intuition around statistical hypothesis testing. To highlight this "from scratch" aspect of the book here is a diagram tying together the various python function used in this post:

This post is part of an ongoing series where I document my progress through Data Science from Scratch by Joel Grus.

For more content on data science, machine learning, R, Python, SQL and more, find me on Twitter.

## Top comments (1)

Hello Paul Apivat,

Thank you for your article. It is well written, short and easy to understand. I imagine making a 23 part series is not an easy task.

Nevertheless, there are some issues I will comment here, may it help you to improve your article.

import random" for the first code snippetimport math" for the first plotA helpful reminder at the beginning of the article what imports are required for the codes to work would solve this problem.

"

When calling the function with our parameters, we get a mean mu of 500 (from 1,000 coin flips) and a standard deviation sigma of 15.8114. Which means that 68% of the time, the binomial random variable will be 500 +/- 15.8114"This finally helped me understand why the standard deviation is needed. Thanks for that.

"_...and 95% of the time it'll be 500 +/- 31.6228 (see 68-95-99.7 rule) _" That's a bit too complicated for me. The only thing I can see is that it's twice 15.8114. Thanks for the link. It gives me the opportunity to search deeper for an answer.

"

lo, hi = normal_two_sided_bounds(0.95, mu_0, sigma_0)"lo # 469.01026640487555

hi # 530.9897335951244

Here the code produces "

NameError: name 'mu_0' is not defined". Possibly same issue with "sigma_0".At this point I stopped reading the article as it is tiring to learn and fix things at the same time.