Central Limit Theorem

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

np.random.seed(42)

# True distribution: Exponential (heavily skewed, not Normal)
lam = 2.0
true_mean = 1 / lam  # 0.5
true_std = 1 / lam   # 0.5

# Simulate CLT: draw samples of increasing size
sample_sizes = [1, 5, 15, 50, 100]
n_experiments = 10000

fig, axes = plt.subplots(1, 5, figsize=(20, 4))

for ax, n in zip(axes, sample_sizes):
    # Draw n_experiments samples, each of size n
    samples = np.random.exponential(1/lam, (n_experiments, n))
    means = samples.mean(axis=1)

    # CLT prediction: N(true_mean, true_std^2 / n)
    theoretical_mean = true_mean
    theoretical_std = true_std / np.sqrt(n)

    # Histogram of sample means
    ax.hist(means, bins=50, density=True, alpha=0.7,
            edgecolor='black', linewidth=0.5)

    # Overlay theoretical Normal
    x = np.linspace(means.min(), means.max(), 200)
    ax.plot(x, stats.norm.pdf(x, theoretical_mean, theoretical_std),
            'r-', linewidth=2, label='CLT prediction')

    # Shapiro-Wilk test for normality
    stat, p_val = stats.shapiro(means[:1000])
    ax.set_title(f'n = {n}\nShapiro p = {p_val:.4f}')
    ax.legend(fontsize=8)

plt.suptitle('CLT: Exponential -> Normal via Averaging', fontsize=14)
plt.tight_layout()
plt.savefig('clt_exponential.png', dpi=150)
plt.show()

import numpy as np
from scipy import stats

np.random.seed(42)

# True conversion rates
p_control = 0.10   # 10% conversion
p_treatment = 0.12  # 12% conversion

n_per_group = 5000
n_simulations = 10000

# Simulate A/B test many times
diffs = []
for _ in range(n_simulations):
    control = np.random.binomial(1, p_control, n_per_group)
    treatment = np.random.binomial(1, p_treatment, n_per_group)
    diffs.append(treatment.mean() - control.mean())

diffs = np.array(diffs)

# CLT prediction for the difference
se_theory = np.sqrt(p_control*(1-p_control)/n_per_group +
                     p_treatment*(1-p_treatment)/n_per_group)
print(f"Observed SE:       {diffs.std():.6f}")
print(f"Theoretical SE:    {se_theory:.6f}")
print(f"True difference:   {p_treatment - p_control:.4f}")
print(f"Mean of diffs:     {diffs.mean():.6f}")

# Normality test
stat, p_val = stats.shapiro(diffs[:1000])
print(f"Shapiro-Wilk p:    {p_val:.4f}")

# 95% CI should contain 0 about 95% of the time under null (p_c = p_t)
# Here, true diff = 0.02, so CI should mostly NOT contain 0
from scipy.stats import norm
z_crit = norm.ppf(0.975)
contains_zero = np.mean(np.abs(diffs) < z_crit * se_theory)
print(f"Fraction within CLT 95% CI of 0: {1 - contains_zero:.1%}")

import numpy as np
from scipy import stats

np.random.seed(42)

def berry_esseen_error(samples, mu, sigma):
    """Empirical max |F_n(x) - Phi(x)|"""
    standardized = (samples - mu) / sigma
    n = len(standardized)
    sorted_data = np.sort(standardized)
    empirical_cdf = np.arange(1, n+1) / n
    normal_cdf = stats.norm.cdf(sorted_data)
    return np.max(np.abs(empirical_cdf - normal_cdf))

distributions = {
    'Uniform(0,1)': {'rvs': np.random.uniform, 'args': (0, 1),
                     'mu': 0.5, 'sigma': 1/np.sqrt(12)},
    'Exponential(1)': {'rvs': np.random.exponential, 'args': (1,),
                       'mu': 1.0, 'sigma': 1.0},
    'Chi-Square(3)': {'rvs': np.random.chisquare, 'args': (3,),
                      'mu': 3.0, 'sigma': np.sqrt(6)},
}

n_values = [10, 30, 50, 100, 500]
n_sims = 5000

for name, dist in distributions.items():
    print(f"\n{name}:")
    for n in n_values:
        errors = []
        for _ in range(n_sims):
            samples = dist['rvs'](*dist['args'], n)
            means = samples.mean(axis=0) if samples.ndim > 1 else samples
            err = berry_esseen_error(
                means, dist['mu'], dist['sigma'] / np.sqrt(n))
            errors.append(err)
        avg_err = np.mean(errors)
        bound = 0.4748 * stats.moment(samples, 3)**0.33 / (dist['sigma']**3 * np.sqrt(n))
        print(f"  n={n:>4d}: empirical={avg_err:.4f}")