Sampling Demo

#!/usr/bin/env python3
"""Optional repeated-sampling illustration, separate from both course datasets.

Known teaching process: independent Normal(mean=24, SD=4) observations.
Generate 5,000 independent samples of size 20 with a fixed NumPy seed. The
fraction of 95% t intervals covering the KNOWN simulated mean illustrates
long-run coverage; it is not a probability assigned to an observed parameter.
Run: python tools/sampling_demo.py
"""
from __future__ import annotations

import json
from pathlib import Path
import numpy as np
from scipy import stats

ROOT = Path(__file__).resolve().parents[1]
SEED = 20261016


def simulate(seed=SEED, repetitions=5000, sample_size=20) -> dict:
    if repetitions < 2 or sample_size < 2:
        raise ValueError('Need at least two repetitions and two observations per sample.')
    rng = np.random.default_rng(seed)
    mu, sigma = 24.0, 4.0
    samples = rng.normal(mu, sigma, size=(repetitions, sample_size))
    means = samples.mean(axis=1)
    sample_sds = samples.std(axis=1, ddof=1)
    standard_errors = sample_sds / np.sqrt(sample_size)
    critical = float(stats.t.ppf(.975, df=sample_size - 1))
    low = means - critical * standard_errors
    high = means + critical * standard_errors
    covered = int(np.count_nonzero((low <= mu) & (mu <= high)))
    return {'seed': seed, 'repetitions': repetitions, 'sample_size': sample_size,
            'known_process_mean': mu, 'known_process_sd': sigma,
            'average_sample_mean': float(means.mean()),
            'sd_of_sample_means': float(means.std(ddof=1)),
            'theoretical_sd_of_sample_means': sigma / np.sqrt(sample_size),
            'average_sample_sd': float(sample_sds.mean()),
            'nominal_confidence_level': .95, 'intervals_covering_known_mean': covered,
            'empirical_coverage_fraction': covered / repetitions,
            'first_five_intervals': [{'mean': float(means[i]), 'low': float(low[i]), 'high': float(high[i]),
                                       'covers_known_mean': bool(low[i] <= mu <= high[i])} for i in range(min(5, repetitions))],
            'scope': 'Independent normal teaching process, not shipment or real-workforce evidence. Each repetition creates a NEW sample.',
            'interpretation': 'The method has 95% theoretical coverage under this normal independent model. A finite simulation need not cover in exactly 95% of repetitions. The parameter 24 is fixed throughout.'}


def main() -> None:
    result = simulate()
    target = ROOT / 'outputs/sampling_demo.json'
    target.parent.mkdir(exist_ok=True)
    target.write_text(json.dumps(result, indent=2, allow_nan=False) + '\n', encoding='utf-8')
    print(f"{result['intervals_covering_known_mean']}/{result['repetitions']} intervals covered mean 24 "
          f"({result['empirical_coverage_fraction']:.2%}).")
    print(f"SD of sample means: {result['sd_of_sample_means']:.4f}; theoretical: "
          f"{result['theoretical_sd_of_sample_means']:.4f}.")


if __name__ == '__main__':
    main()