import pandas as pd
import numpy as np
from scipy import stats

df = pd.read_csv('superstore_sales.csv')
pop_mean = df['Sales'].mean()
print(f'Real population mean Sales: {pop_mean:.2f}')

sample = df.sample(n=100, random_state=21)['Sales']
n = len(sample)
mean = sample.mean()
se = sample.std(ddof=1) / np.sqrt(n)
print(f'Real sample mean: {mean:.2f}')
print(f'Real standard error: {se:.2f}')

t_crit = stats.t.ppf(0.975, df=n - 1)
margin = t_crit * se
ci_low, ci_high = mean - margin, mean + margin
print(f'Real t critical value (95% CI, df={n - 1}): {t_crit:.4f}')
print(f'Real 95% CI: ({ci_low:.2f}, {ci_high:.2f})')
print(f'Real true population mean: {pop_mean:.2f}')
print(f'Does this real interval contain the true mean? {ci_low <= pop_mean <= ci_high}')

ci_scipy = stats.t.interval(0.95, df=n - 1, loc=mean, scale=se)
print(f'scipy shortcut, same real result: ({ci_scipy[0]:.2f}, {ci_scipy[1]:.2f})')

sample_p = df.sample(n=300, random_state=21)
p_hat = (sample_p['Category'] == 'Technology').mean()
n_p = len(sample_p)
se_p = np.sqrt(p_hat * (1 - p_hat) / n_p)
z_crit = stats.norm.ppf(0.975)
ci_p_low, ci_p_high = p_hat - z_crit * se_p, p_hat + z_crit * se_p
true_p = (df['Category'] == 'Technology').mean()
print(f'Real sample proportion Technology: {p_hat:.4f}')
print(f'Real 95% CI: ({ci_p_low:.4f}, {ci_p_high:.4f})')
print(f'Real true population proportion: {true_p:.4f}')

west = df[df['Region'] == 'West'].sample(n=150, random_state=5)['Sales']
central = df[df['Region'] == 'Central'].sample(n=150, random_state=5)['Sales']
diff = west.mean() - central.mean()
se_diff = np.sqrt(west.var(ddof=1) / 150 + central.var(ddof=1) / 150)
ci_diff_low, ci_diff_high = diff - z_crit * se_diff, diff + z_crit * se_diff
print(f'Real West sample mean: {west.mean():.2f}, Central sample mean: {central.mean():.2f}')
print(f'Real observed difference: {diff:.2f}')
print(f'Real 95% CI for the difference: ({ci_diff_low:.2f}, {ci_diff_high:.2f})')

true_west_mean = df[df['Region'] == 'West']['Sales'].mean()
true_central_mean = df[df['Region'] == 'Central']['Sales'].mean()
true_diff = true_west_mean - true_central_mean
print(f'Real true West population mean: {true_west_mean:.2f}')
print(f'Real true Central population mean: {true_central_mean:.2f}')
print(f'Real true difference: {true_diff:.2f}')
print(f'Does this real interval contain the true difference? {ci_diff_low <= true_diff <= ci_diff_high}')

rng = np.random.default_rng(seed=21)
sample_vals = df.sample(n=100, random_state=21)['Sales'].values
sample_median = np.median(sample_vals)
boot_medians = np.array([np.median(rng.choice(sample_vals, size=len(sample_vals), replace=True)) for _ in range(5000)])
boot_ci_low, boot_ci_high = np.percentile(boot_medians, [2.5, 97.5])
print(f'Real sample median: {sample_median:.2f}')
print(f'Real bootstrap 95% CI: ({boot_ci_low:.2f}, {boot_ci_high:.2f})')

true_median = df['Sales'].median()
print(f'Real true population median: {true_median:.2f}')
print(f'Does the real bootstrap interval contain the true median? {boot_ci_low <= true_median <= boot_ci_high}')

rng2 = np.random.default_rng(seed=100)
sales_values = df['Sales'].values
for sample_n in [30, 100, 500, 1000]:
    hits = 0
    trials = 1000
    for _ in range(trials):
        s = rng2.choice(sales_values, size=sample_n, replace=False)
        m = s.mean()
        se_s = s.std(ddof=1) / np.sqrt(sample_n)
        tcrit_s = stats.t.ppf(0.975, df=sample_n - 1)
        lo, hi = m - tcrit_s * se_s, m + tcrit_s * se_s
        if lo <= pop_mean <= hi:
            hits += 1
    print(f'n={sample_n}: real empirical coverage = {hits / trials:.1%} (target: 95%)')