# --- Setup: load our helpers and the data --------------------------------
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        # the one-sample t-test we'll check against
sl.use_course_style()

rng = np.random.default_rng(0)   # ONE seeded generator drives every random draw below

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

# Treat the WHOLE file as the population, exactly like the core lesson.
population = ames['area'].dropna().to_numpy()    # above-grade living area, sq ft
MU = population.mean()                            # the TRUE mean — we get to peek at it
N_POP = len(population)
print(f'Loaded a population of {N_POP:,} Ames homes')
print(f'TRUE mean MU = {MU:.1f} sq ft   <- normally hidden; here we can look')

Loaded a population of 2,930 Ames homes
TRUE mean MU = 1499.7 sq ft   <- normally hidden; here we can look

mu0 = 1400.0          # H0: the true mean is 1,400 sq ft (the realtor's claim)
n = 40                # our sample size

sample = rng.choice(population, size=n, replace=True)   # one sample, drawn with replacement
xbar = sample.mean()
s = sample.std(ddof=1)                                  # ddof=1 -> divide by n-1 (a SAMPLE sd)
se = s / np.sqrt(n)                                      # standard error of the mean
t_obs = (xbar - mu0) / se                                # the OBSERVED test statistic

print(f'sample mean  xbar = {xbar:.1f} sq ft')
print(f'sample stdev s    = {s:.1f} sq ft   (ddof=1)')
print(f'standard error    = {se:.1f} sq ft')
print(f'observed t        = {t_obs:.3f}   (our sample sits {t_obs:.2f} standard errors from the claim)')

sample mean  xbar = 1468.9 sq ft
sample stdev s    = 599.6 sq ft   (ddof=1)
standard error    = 94.8 sq ft
observed t        = 0.727   (our sample sits 0.73 standard errors from the claim)

# A population whose true mean IS the null value (same shape, just shifted).
null_population = population - MU + mu0
print(f'null world mean = {null_population.mean():.1f} sq ft  (exactly mu0, by construction)')

# Draw many samples FROM THE NULL WORLD and record the t statistic each time.
reps = 20_000
null_ts = np.empty(reps)
for i in range(reps):
    s_null = rng.choice(null_population, size=n, replace=True)
    null_ts[i] = (s_null.mean() - mu0) / (s_null.std(ddof=1) / np.sqrt(n))

print(f'Built a null distribution of {reps:,} t statistics.')
print(f'  mean   ~ {null_ts.mean():+.3f}  (centered on 0, as the null demands)')
print(f'  stdev  ~ {null_ts.std(ddof=1):.3f}')

null world mean = 1400.0 sq ft  (exactly mu0, by construction)

Built a null distribution of 20,000 t statistics.
  mean   ~ -0.090  (centered on 0, as the null demands)
  stdev  ~ 1.036

extreme = np.abs(null_ts) >= np.abs(t_obs)   # a True/False mask, one entry per simulated sample
p_sim = extreme.mean()                        # mean of a boolean array = the fraction that are True

print(f'observed |t|                = {abs(t_obs):.3f}')
print(f'simulated two-sided p-value = {p_sim:.4f}')
print(f'(the share of {reps:,} make-believe samples with |t| >= {abs(t_obs):.2f})')
print(f'\nPlain English: a sample this far from the claim turns up by pure chance about')
print(f'{p_sim*100:.0f} times in 100, even when the claim is exactly right.')

observed |t|                = 0.727
simulated two-sided p-value = 0.4704
(the share of 20,000 make-believe samples with |t| >= 0.73)

Plain English: a sample this far from the claim turns up by pure chance about
47 times in 100, even when the claim is exactly right.

fig, ax = plt.subplots()
ax.hist(null_ts, bins=60, color='#9aa7c7', edgecolor='white', density=True,
        label='null distribution of t (simulated)')
for sign in (-1, 1):
    ax.axvline(sign * abs(t_obs), color='#D55E00', lw=2,
               label='observed |t|' if sign == 1 else None)
xs = np.linspace(null_ts.min(), null_ts.max(), 300)
tail = xs[np.abs(xs) >= abs(t_obs)]
ax.fill_between(tail, stats.t.pdf(tail, df=n - 1), color='#D55E00', alpha=0.30)
ax.set_xlabel('t statistic in the null world'); ax.set_ylabel('density')
ax.set_title(f'Null distribution of t — shaded tails (the p-value) = {p_sim:.3f}')
ax.legend(); plt.show()

result = stats.ttest_1samp(sample, popmean=mu0, alternative='two-sided')
print(f'simulated p-value (null world) = {p_sim:.4f}')
print(f'scipy ttest_1samp p-value      = {result.pvalue:.4f}')
print(f'scipy t statistic              = {result.statistic:.3f}   (matches ours: {t_obs:.3f})')
print('\nmatch?', np.isclose(p_sim, result.pvalue, atol=0.01))
print('The formula and the simulation agree — two views of the same idea.')

simulated p-value (null world) = 0.4704
scipy ttest_1samp p-value      = 0.4719
scipy t statistic              = 0.727   (matches ours: 0.727)

match? True
The formula and the simulation agree — two views of the same idea.

# The loop way (how you might first think of it) — already done above as null_ts.

# The vectorized way: one (reps, n) array, then reduce each ROW to a t statistic.
samples = rng.choice(null_population, size=(reps, n), replace=True)   # reps rows, n columns
means = samples.mean(axis=1)                                          # one mean per row
sds   = samples.std(axis=1, ddof=1)                                  # one sd  per row (ddof=1!)
null_ts_vec = (means - mu0) / (sds / np.sqrt(n))                      # reps t's, no loop

p_sim_vec = np.mean(np.abs(null_ts_vec) >= np.abs(t_obs))
print(f'loop p-value       : {p_sim:.4f}')
print(f'vectorized p-value : {p_sim_vec:.4f}   <- same idea, one block, no Python loop')

import time
t0 = time.perf_counter()
for i in range(reps):
    s_null = rng.choice(null_population, size=n, replace=True)
    _ = (s_null.mean() - mu0) / (s_null.std(ddof=1) / np.sqrt(n))
t_loop = time.perf_counter() - t0

t0 = time.perf_counter()
samp = rng.choice(null_population, size=(reps, n), replace=True)
_ = (samp.mean(axis=1) - mu0) / (samp.std(axis=1, ddof=1) / np.sqrt(n))
t_vec = time.perf_counter() - t0

print(f'\n{reps:,} simulated samples:')
print(f'  Python loop : {t_loop*1000:7.1f} ms')
print(f'  vectorized  : {t_vec*1000:7.1f} ms')
print(f'  vectorized is roughly {t_loop/max(t_vec,1e-9):.0f}x faster — and it is a few short lines.')

loop p-value       : 0.4704
vectorized p-value : 0.4713   <- same idea, one block, no Python loop

20,000 simulated samples:
  Python loop :   409.2 ms
  vectorized  :     7.1 ms
  vectorized is roughly 57x faster — and it is a few short lines.

def sim_p_value(observed_stat, simulate_null_stat, reps=20_000):
    '''Two-sided p-value by simulation.

    observed_stat      : the statistic computed from the real sample.
    simulate_null_stat : a function taking no args that returns ONE statistic
                         drawn from the null world.
    Returns the fraction of null-world statistics with |stat| >= |observed_stat|.
    '''
    null_stats = np.array([simulate_null_stat() for _ in range(reps)])
    return np.mean(np.abs(null_stats) >= np.abs(observed_stat))


def one_t_under_null():
    '''Draw one sample from the null world and return its t statistic.'''
    s = rng.choice(null_population, size=n, replace=True)
    return (s.mean() - mu0) / (s.std(ddof=1) / np.sqrt(n))


p_refactored = sim_p_value(t_obs, one_t_under_null, reps=20_000)
print(f'sim_p_value(...) p-value  = {p_refactored:.4f}')
print(f'scipy ttest_1samp p-value = {result.pvalue:.4f}')
print('match?', np.isclose(p_refactored, result.pvalue, atol=0.01))

sim_p_value(...) p-value  = 0.4695
scipy ttest_1samp p-value = 0.4719
match? True

def reject_rate(test_population, h0, n=50, alpha=0.05, reps=4000, seed=1):
    '''Fraction of one-sample t-tests of H0: mu = h0 that REJECT, over many fresh samples.

    We sample WITH replacement, treating the file as an effectively infinite super-population,
    so n is not capped at 2,930 and there is no finite-population correction to distort the SE.
    (Sampling WITHOUT replacement here would shrink the SE and reverse the big-n demo in Sec. 6.)
    '''
    r = np.random.default_rng(seed)
    s = r.choice(test_population, size=(reps, n), replace=True)   # reps tests at once
    t = (s.mean(axis=1) - h0) / (s.std(axis=1, ddof=1) / np.sqrt(n))
    p = 2 * stats.t.sf(np.abs(t), df=n - 1)                       # two-sided p for every row
    return np.mean(p <= alpha)

# Sanity check: testing the TRUE mean should reject ~alpha of the time (see Section 6).
print('reject_rate at the true mean:', round(reject_rate(population, h0=MU, n=50), 3))

reject_rate at the true mean: 0.057

for e in [0, 25, 50, 75, 100, 150, 200, 250]:
    print(f'effect {e:3d} sq ft -> power {reject_rate(population, h0=MU - e, n=50, seed=2):.3f}')

Tool	What it is	What we use it for
NumPy	fast math on whole arrays at once	drawing samples, forming the statistic, the boolean mask
`np.random.default_rng`	a seeded random-number generator	reproducible "chance" — same seed, same numbers every run
scipy.stats	textbook statistical tests	the `ttest_1samp` we check our simulation against
matplotlib	plotting	the null-distribution histogram and the power curve

The code behind it: build a test by simulating the null¶

1) The goal¶

2) The data in code: a claim, a sample, a statistic¶

3) Build it step by step: a p-value by simulating the null¶

The p-value: the mean of an "at least as extreme" mask¶

See it, then check it against `scipy`¶

4) The idiom — vectorized null simulation vs. a loop¶

5) Refactor into a reusable function¶

6) Now you code¶

What you learned¶

The code behind it: build a test by simulating the null¶

1) The goal¶

2) The data in code: a claim, a sample, a statistic¶

3) Build it step by step: a p-value by simulating the null¶

The p-value: the mean of an "at least as extreme" mask¶

See it, then check it against scipy¶

4) The idiom — vectorized null simulation vs. a loop¶

5) Refactor into a reusable function¶

6) Now you code¶

What you learned¶

See it, then check it against `scipy`¶