Flagship demo — What a 95% confidence interval really means¶

import sys, pathlib
for _p in [pathlib.Path.cwd(), pathlib.Path.cwd().parent]:
    if (_p / 'lib' / 'statslab.py').exists():
        sys.path.insert(0, str(_p / 'lib')); break
import statslab as sl

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
sl.use_course_style()

AMES_URL = 'https://www.openintro.org/data/csv/ames.csv'
ames = sl.load_csv('ames.csv', url=AMES_URL)

population = ames['area'].dropna().to_numpy()    # treat the WHOLE file as the population
MU = population.mean()                            # the TRUE mean (we get to know it!)
SIGMA = population.std(ddof=0)                    # the TRUE population standard deviation
N_POP = len(population)
print(f'Population: {N_POP:,} homes')
print(f'TRUE population mean  MU    = {MU:.1f} sq ft   <- the answer we will try to catch')
print(f'TRUE population stdev SIGMA = {SIGMA:.1f} sq ft')

fig, ax = plt.subplots()
ax.hist(population, bins=40, color='#9aa7c7', edgecolor='white')
ax.axvline(MU, color='black', lw=2, ls='--', label=f'true mean = {MU:.0f} sq ft')
ax.set_xlabel('living area (sq ft)'); ax.set_ylabel('number of homes')
ax.set_title('The population: all 2,930 homes (note the right skew)')
ax.legend(); plt.show()

rng = np.random.default_rng(7)   # change this seed to draw a different sample
n = 50

def confidence_interval(sample, level=0.95):
    n = len(sample)
    xbar = sample.mean()
    s = sample.std(ddof=1)
    se = s / np.sqrt(n)
    tcrit = stats.t.ppf(1 - (1 - level) / 2, df=n - 1)
    return xbar - tcrit * se, xbar + tcrit * se, xbar

sample = rng.choice(population, size=n, replace=False)
lo, hi, xbar = confidence_interval(sample)
print(f'Sample mean  = {xbar:.1f} sq ft')
print(f'95% CI       = ({lo:.1f}, {hi:.1f})')
print(f'True mean MU = {MU:.1f} sq ft')
print('Did THIS interval catch the true mean? ->', 'YES' if lo <= MU <= hi else 'NO')

n_intervals = 100
results = []
for _ in range(n_intervals):
    s = rng.choice(population, size=n, replace=False)
    lo, hi, xbar = confidence_interval(s)
    results.append({'lo': lo, 'hi': hi, 'xbar': xbar, 'hit': lo <= MU <= hi})

hits = sum(r['hit'] for r in results)
print(f'{hits} out of {n_intervals} intervals contained the true mean ({MU:.1f} sq ft).')
print(f'That is {hits}% coverage - the theory predicts about 95%.')

fig, ax = plt.subplots(figsize=(7, 10))
for i, res in enumerate(results):
    hit = res['hit']
    style = sl.hit_miss_style(hit)          # solid circle = a hit, dashed cross = a miss
    color = sl.CAPTURE if hit else sl.MISS  # colourblind-safe teal vs. vermillion
    ax.plot([res['lo'], res['hi']], [i, i], color=color, lw=1.4, ls=style['linestyle'])
    ax.plot(res['xbar'], i, color=color, ms=3.5, marker=style['marker'], ls='none')
ax.axvline(MU, color='black', ls='--', lw=2, label=f'true mean = {MU:.0f} sq ft')
ax.set_xlabel('living area (sq ft)'); ax.set_ylabel('interval number (1-100)')
ax.set_title(f'100 confidence intervals: {hits} caught the true mean, {n_intervals - hits} missed')
ax.legend(loc='lower right'); plt.show()

def coverage(level, n=50, trials=10000):
    hits = 0
    for _ in range(trials):
        s = rng.choice(population, size=n, replace=False)
        lo, hi, _ = confidence_interval(s, level=level)
        hits += (lo <= MU <= hi)
    return hits / trials

print(f'Coverage of the 95% method over 10,000 samples: {coverage(0.95) * 100:.1f}%')

levels = [0.80, 0.90, 0.95, 0.99]
rows = []
for lev in levels:
    widths, hits, trials = [], 0, 5000
    for _ in range(trials):
        s = rng.choice(population, size=50, replace=False)
        lo, hi, _ = confidence_interval(s, level=lev)
        widths.append(hi - lo); hits += (lo <= MU <= hi)
    rows.append({'level': f'{int(lev*100)}%', 'coverage_%': hits / trials * 100, 'avg_width_sqft': np.mean(widths)})
summary = pd.DataFrame(rows)
print(summary.to_string(index=False))

fig, (a1, a2) = plt.subplots(1, 2, figsize=(11, 4))
a1.bar(summary['level'], summary['coverage_%'], color='#4c72b0'); a1.axhline(95, color='grey', ls=':')
a1.set_title('Actual coverage'); a1.set_ylabel('% of intervals catching the mean'); a1.set_ylim(70, 100)
a2.bar(summary['level'], summary['avg_width_sqft'], color='#E69F00')
a2.set_title('Average interval width'); a2.set_ylabel('sq ft (smaller = more precise)')
plt.tight_layout(); plt.show()

many_means = np.array([rng.choice(population, size=50, replace=False).mean() for _ in range(5000)])
fig, ax = plt.subplots()
ax.hist(many_means, bins=40, density=True, color='#9aa7c7', edgecolor='white', label='sample means (n=50)')
xs = np.linspace(many_means.min(), many_means.max(), 200)
ax.plot(xs, stats.norm.pdf(xs, MU, SIGMA / np.sqrt(50)), 'r-', lw=2, label='CLT normal curve')
ax.axvline(MU, color='black', ls='--', label=f'true mean = {MU:.0f}')
ax.set_xlabel('sample mean of living area (sq ft)'); ax.set_ylabel('density')
ax.set_title('Skewed population, but the sample means are normal (CLT)')
ax.legend(); plt.show()
print(f'Std of the sample means   = {many_means.std():.1f} sq ft')
print(f'CLT prediction sigma/sqrt(n) = {SIGMA / np.sqrt(50):.1f} sq ft   (they should match)')