Stats & ML Guide

Python Example

# Two-Sample Hypothesis Testing Workflow
import numpy as np
from scipy import stats

# Step 1: Explore data
control = df[df['group'] == 'control']['metric'].values
treatment = df[df['group'] == 'treatment']['metric'].values
print(f"Control: mean={control.mean():.3f}, std={control.std():.3f}, n={len(control)}")
print(f"Treatment: mean={treatment.mean():.3f}, std={treatment.std():.3f}, n={len(treatment)}")

# Step 2: Check normality
_, p_shapiro = stats.shapiro(control[:500])  # Shapiro sensitive at large n
print(f"Shapiro-Wilk p={p_shapiro:.4f}")

# Step 3: Welch's t-test (default — doesn't assume equal variance)
t_stat, p_value = stats.ttest_ind(treatment, control, equal_var=False)
print(f"Welch's t-test: t={t_stat:.3f}, p={p_value:.4f}")

# Step 4: Non-parametric alternative
u_stat, p_mw = stats.mannwhitneyu(treatment, control, alternative='two-sided')
print(f"Mann-Whitney U: U={u_stat:.0f}, p={p_mw:.4f}")

# Step 5: Effect size (Cohen's d)
pooled_std = np.sqrt((control.std()**2 + treatment.std()**2) / 2)
cohens_d = (treatment.mean() - control.mean()) / pooled_std
print(f"Cohen's d: {cohens_d:.3f}")

# Step 6: Z-test for proportions (binary outcomes)
from statsmodels.stats.proportion import proportions_ztest
z_stat, p_prop = proportions_ztest([conv_treat, conv_ctrl], [n_treat, n_ctrl])

Key Assumptions & Pitfalls

• Independent samples — each observation in one group only, no cross-contamination between users
• Approximately normal distributions, or n > 30 for CLT to apply
• Welch's t-test handles unequal variance; use it by default over Student's t
• Mann-Whitney tests for stochastic dominance (ranks), not means — use when normality is violated
• A significant p-value does NOT mean the effect is large — always report effect size (Cohen's d)
• Default to two-tailed tests unless you have a strong directional prior stated before seeing data

Python Example

# Confidence Intervals — Means vs Proportions
import numpy as np
from scipy import stats
from statsmodels.stats.proportion import proportion_confint

# CI for continuous data (t-based)
mean_diff = treatment.mean() - control.mean()
se_diff = np.sqrt(treatment.var()/len(treatment) + control.var()/len(control))
ci_t = (mean_diff - 1.96 * se_diff, mean_diff + 1.96 * se_diff)

# CI for proportions (normal approximation)
p1, n1 = 0.732, 5000
se1 = np.sqrt(p1 * (1 - p1) / n1)
ci_normal = (p1 - 1.96 * se1, p1 + 1.96 * se1)

# Wilson score interval (better for small n or extreme p)
ci_wilson = proportion_confint(int(p1 * n1), n1, alpha=0.05, method='wilson')

# Bootstrap CI (works for any statistic)
boot_means = [np.random.choice(data, size=len(data), replace=True).mean()
              for _ in range(10_000)]
ci_boot = np.percentile(boot_means, [2.5, 97.5])

Key Assumptions & Pitfalls

• WRONG: "95% probability the true value is in this range" — treats CI as random, parameter as fixed (backwards)
• RIGHT: "If repeated many times, 95% of intervals would contain the true parameter"
• Overlapping CIs do NOT mean no significant difference — always test the difference directly
• Normal approximation breaks for p < 0.01 or p > 0.99, or n < 30 — use Wilson or bootstrap
• Bootstrap CIs are the universal fallback — no distributional assumptions needed
• Wider CI = more uncertainty. Narrow CI = precise estimate. Width ∝ 1/√n

Python Example

# Power Analysis & Sample Size Calculation
from scipy import stats
import numpy as np

def required_sample_size(alpha, power, effect_size):
    """Calculate required n per group for a two-sample t-test."""
    z_alpha = stats.norm.ppf(1 - alpha / 2)
    z_beta = stats.norm.ppf(power)
    n = 2 * ((z_alpha + z_beta) / effect_size) ** 2
    return int(np.ceil(n))

# Example: detect 2pp lift in filing completion
baseline_rate = 0.45
mde = 0.02
sigma = np.sqrt(baseline_rate * (1 - baseline_rate))
cohens_d = mde / sigma
n_required = required_sample_size(0.05, 0.80, cohens_d)
print(f"Cohen's d: {cohens_d:.3f}")
print(f"Required n per group: {n_required:,}")

# Using statsmodels (more precise)
from statsmodels.stats.power import TTestIndPower
analysis = TTestIndPower()
n_precise = analysis.solve_power(effect_size=cohens_d, alpha=0.05, power=0.80)
print(f"Statsmodels n per group: {int(np.ceil(n_precise)):,}")

Key Assumptions & Pitfalls

• Power analysis must be done BEFORE the experiment — doing it after is circular reasoning
• MDE comes from the business, not statistics: "What's the smallest change worth caring about?"
• Underpowered tests are wasteful — you run the experiment but can't detect the effect
• Halving the MDE quadruples the required sample size (n ∝ 1/d²)
• Variance reduction (CUPED) effectively increases power without collecting more data
• Account for clustering (design effect) in marketplace experiments — ignoring it leads to underpowered tests

Python Example

# Multiple Testing Correction Workflow
from scipy import stats
from statsmodels.stats.multitest import multipletests
import numpy as np

# 4 variants × 3 metrics = 12 pairwise tests
p_values = [0.003, 0.012, 0.032, 0.048, 0.051, 0.089,
            0.12, 0.23, 0.34, 0.41, 0.67, 0.91]

# Bonferroni correction (FWER control)
reject_bonf, pvals_bonf, _, _ = multipletests(p_values, method='bonferroni')
print("Bonferroni rejections:", sum(reject_bonf))

# Benjamini-Hochberg (FDR control)
reject_bh, pvals_bh, _, _ = multipletests(p_values, method='fdr_bh')
print("BH rejections:", sum(reject_bh))

# Chi-squared omnibus → pairwise z-tests
chi2, p_omnibus, dof, expected = stats.chi2_contingency(observed_table)
if p_omnibus < 0.05:
    # Pairwise z-tests with BH correction
    from statsmodels.stats.proportion import proportions_ztest
    pairwise_p = []
    for variant in variants:
        z, p = proportions_ztest([variant.conv, ctrl.conv], [variant.n, ctrl.n])
        pairwise_p.append(p)
    reject, adj_p, _, _ = multipletests(pairwise_p, method='fdr_bh')

Key Assumptions & Pitfalls

• Always correct across ALL tests you ran, not just the ones with small p-values (this IS p-hacking)
• Chi-squared is the omnibus test for proportions (3+ groups); ANOVA for continuous metrics
• After omnibus rejection, do pairwise tests with correction — never skip the omnibus step
• Adjust CIs to match: if Bonferroni at α/3 = 0.017, report 98.3% CIs for consistency
• BH controls FDR under independence or positive dependence — safe for most A/B scenarios
• Bonferroni is over-conservative with many tests — you'll miss real effects

Python Example

# Sequential Testing: Why Peeking Inflates FPR
import numpy as np
from scipy import stats

def simulate_peeking_fpr(n_peeks, n_sims=10000, alpha=0.05):
    """Simulate false positive rate when peeking at results."""
    false_pos = 0
    n_per_peek = 100
    z_crit = stats.norm.ppf(1 - alpha / 2)
    for _ in range(n_sims):
        a = np.random.normal(0, 1, n_per_peek * n_peeks)
        b = np.random.normal(0, 1, n_per_peek * n_peeks)
        for peek in range(1, n_peeks + 1):
            n = peek * n_per_peek
            z = abs((a[:n].mean() - b[:n].mean()) / np.sqrt(2/n))
            if z > z_crit:
                false_pos += 1
                break
    return false_pos / n_sims

for peeks in [1, 2, 5, 10, 20, 50]:
    fpr = simulate_peeking_fpr(peeks)
    print(f"Peeks: {peeks:2d} → Actual FPR: {fpr:.1%}")
# Peeks:  1 → 5.0%   Peeks:  5 → 14.2%
# Peeks: 10 → 19.3%  Peeks: 50 → 32.1%

Key Assumptions & Pitfalls

• Traditional A/B tests with fixed-horizon analysis are invalid if you peek and stop early
• mSPRT / always-valid p-values allow continuous monitoring without FPR inflation
• O'Brien-Fleming is the most popular group sequential design — conservative early, lenient late
• Sequential tests trade power for flexibility — typically need 20-30% more total samples
• Never use a significant interim result as final if you didn't pre-specify interim analyses
• Uber's experimentation platform uses sequential testing by default for marketplace experiments

Python Example

# CUPED & Delta Method
import numpy as np

# --- CUPED ---
def cuped_adjust(y_post, x_pre):
    """Adjust outcomes using pre-experiment covariate."""
    theta = np.cov(y_post, x_pre)[0,1] / np.var(x_pre)
    return y_post - theta * (x_pre - np.mean(x_pre))

np.random.seed(42)
n = 10000
x_pre = np.random.normal(65, 20, n)
treatment = np.random.binomial(1, 0.5, n)
y_post = 0.7 * x_pre + treatment * 2.0 + np.random.normal(0, 15, n)

y_adj = cuped_adjust(y_post, x_pre)
print(f"Raw SE:   {y_post.std()/np.sqrt(n/2):.3f}")
print(f"CUPED SE: {y_adj.std()/np.sqrt(n/2):.3f}")
print(f"Variance reduction: {1 - y_adj.var()/y_post.var():.1%}")

# --- Delta Method for Ratio Metrics ---
def delta_method_ratio_se(y, x):
    """SE of mean(y)/mean(x) via delta method."""
    n = len(y)
    mu_y, mu_x = y.mean(), x.mean()
    var_y, var_x = y.var(), x.var()
    cov_yx = np.cov(y, x)[0,1]
    return np.sqrt((var_y - 2*(mu_y/mu_x)*cov_yx + (mu_y/mu_x)**2*var_x) / (mu_x**2*n))

Key Assumptions & Pitfalls

• CUPED requires pre-experiment covariate X measured BEFORE randomization
• Higher ρ(X, Y) = more variance reduction — stale pre-period data may have low ρ
• Delta method assumes approximate normality of the estimator (CLT for large n)
• For heavily skewed ratio metrics, block bootstrap may be more reliable than delta method
• CUPED + delta method can be combined for ratio metrics with pre-period covariates
• CUPED is equivalent to regression adjustment — the treatment effect estimate is identical

Python Example

# Bayesian A/B Testing & Thompson Sampling
import numpy as np

# --- Bayesian A/B Test ---
alpha_A, beta_A = 1 + 450, 1 + 5000 - 450   # Control: 9.0%
alpha_B, beta_B = 1 + 520, 1 + 5000 - 520   # Treatment: 10.4%

samples_A = np.random.beta(alpha_A, beta_A, 100000)
samples_B = np.random.beta(alpha_B, beta_B, 100000)
p_b_better = (samples_B > samples_A).mean()
loss_B = np.maximum(samples_A - samples_B, 0).mean()
print(f"P(B > A) = {p_b_better:.3f}")
print(f"Expected loss of choosing B: {loss_B:.5f}")

# --- Thompson Sampling ---
class ThompsonBandit:
    def __init__(self, n_arms):
        self.alpha = np.ones(n_arms)
        self.beta = np.ones(n_arms)
    def select_arm(self):
        samples = [np.random.beta(a, b) for a, b in zip(self.alpha, self.beta)]
        return np.argmax(samples)
    def update(self, arm, reward):
        if reward: self.alpha[arm] += 1
        else: self.beta[arm] += 1

true_probs = [0.05, 0.07, 0.10]
bandit = ThompsonBandit(3)
for _ in range(10000):
    arm = bandit.select_arm()
    reward = np.random.binomial(1, true_probs[arm])
    bandit.update(arm, reward)

Key Assumptions & Pitfalls

• Bayesian results depend on the prior — with enough data, prior washes out. Beta(1,1) is standard.
• Thompson sampling converges to the optimal arm but initial exploration means some "wasted" traffic
• Bandits sacrifice clean causal inference — can't compute unbiased ATEs from bandit data
• Use A/B tests for causal inference; use bandits for optimizing cumulative reward
• Contextual bandits require well-chosen features — bad features = worse than non-contextual
• Expected loss is often more useful than P(best) for business decisions

Python Example

# Distribution Fundamentals
import numpy as np
from scipy import stats

# Normal distribution
x = np.linspace(-4, 4, 200)
pdf = stats.norm.pdf(x, 0, 1)  # 68-95-99.7 rule

# Binomial: n trials, p success probability
pmf = stats.binom.pmf(np.arange(21), n=20, p=0.3)
print(f"Binomial(20,0.3): mean={20*0.3}, var={20*0.3*0.7:.1f}")

# Poisson: λ events per unit time
pmf = stats.poisson.pmf(np.arange(20), mu=5)
print(f"Poisson(5): mean=5, var=5")  # mean = variance!

# MLE for normal
data = np.random.normal(5, 2, 1000)
mu_mle, sigma_mle = stats.norm.fit(data)

# CLT demonstration
means = [np.random.exponential(2, 50).mean() for _ in range(10000)]
# means are approximately normal even though exponential is skewed!

# Correlation and R²
r = np.corrcoef(x, y)[0,1]
print(f"Correlation: {r:.3f}, R²: {r**2:.3f}")

Key Assumptions & Pitfalls

• Normal: 68-95-99.7 rule. CLT makes it the default for sample means.
• Binomial requires: fixed n, independent trials, constant p. Approximate with normal when np ≥ 5 and n(1-p) ≥ 5.
• Poisson: events independent, constant rate. Var = Mean — check this in your data!
• Exponential: memoryless — P(T > t+s | T > t) = P(T > s). Only distribution with this property.
• Correlation measures LINEAR relationships only — zero correlation ≠ independence
• MLE is consistent and efficient, but can overfit with small samples — consider regularization

Python Example

# Regression Diagnostics Pipeline
from sklearn.linear_model import LinearRegression, LogisticRegression, Ridge, Lasso
from sklearn.model_selection import cross_val_score
from sklearn.metrics import roc_auc_score, classification_report
from scipy import stats
import numpy as np

# Linear Regression with diagnostics
model = LinearRegression().fit(X_train, y_train)
residuals = y_test - model.predict(X_test)
_, p_normal = stats.shapiro(residuals[:500])  # Normality check

# Multicollinearity (VIF > 10 = problem)
from statsmodels.stats.outliers_influence import variance_inflation_factor
vif = [variance_inflation_factor(X_train, i) for i in range(X_train.shape[1])]

# Logistic Regression with L1
log_model = LogisticRegression(penalty='l1', C=1.0, solver='saga')
log_model.fit(X_train, y_train)
auc = roc_auc_score(y_test, log_model.predict_proba(X_test)[:,1])

# Regularization comparison
for name, m in [("Ridge", Ridge(1.0)), ("Lasso", Lasso(0.1))]:
    scores = cross_val_score(m, X_train, y_train, cv=5, scoring='r2')
    m.fit(X_train, y_train)
    n_zero = (np.abs(m.coef_) < 1e-6).sum()
    print(f"{name}: R²={scores.mean():.3f}, zero coefs={n_zero}")

Key Assumptions & Pitfalls

• Linear regression: linearity, independence, homoscedasticity, normal residuals — check all four
• Logistic regression: no multicollinearity, linear relationship between features and log-odds
• Ridge when all features matter; Lasso for feature selection; Elastic Net when unsure
• Regularization λ must be tuned via cross-validation, not guessed
• AUC is the primary metric for imbalanced classification; F1 when balancing precision/recall
• R² can be misleading — always check residual plots for patterns

Model Comparison

Python Example

# Tree-Based Models Comparison
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier
from sklearn.tree import DecisionTreeClassifier
from sklearn.model_selection import cross_val_score
from sklearn.metrics import roc_auc_score

models = {
    "Decision Tree": DecisionTreeClassifier(max_depth=5),
    "Random Forest": RandomForestClassifier(n_estimators=200, max_depth=10),
    "Gradient Boosting": GradientBoostingClassifier(n_estimators=200, lr=0.1),
}
for name, model in models.items():
    scores = cross_val_score(model, X_train, y_train, cv=5, scoring='roc_auc')
    print(f"{name}: AUC={scores.mean():.3f} ± {scores.std():.3f}")

# Feature importance
rf = RandomForestClassifier(n_estimators=200).fit(X_train, y_train)
for name, imp in sorted(zip(features, rf.feature_importances_), key=lambda x: -x[1])[:5]:
    print(f"  {name}: {imp:.3f}")

Key Assumptions & Pitfalls

• Decision trees are greedy (locally optimal) — may miss globally optimal partitions
• Random Forest reduces variance via averaging; does NOT reduce bias of individual trees
• Gradient Boosting reduces both bias and variance but overfits without regularization
• XGBoost's scale_pos_weight handles class imbalance — set to (neg count / pos count)
• Learning rate and n_estimators are inversely related — lower rate + more trees = better generalization
• Feature importance from trees measures predictive contribution, NOT causal effect

Python Example

# Clustering, PCA & Special Topics
from sklearn.cluster import KMeans, DBSCAN
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_score

# PCA
X_scaled = StandardScaler().fit_transform(X)
pca = PCA(n_components=0.95)  # Keep 95% of variance
X_pca = pca.fit_transform(X_scaled)
print(f"Reduced {X.shape[1]} → {X_pca.shape[1]} components")

# K-Means with elbow
for k in range(2, 11):
    km = KMeans(n_clusters=k, n_init=10).fit(X_pca)
    print(f"k={k}: silhouette={silhouette_score(X_pca, km.labels_):.3f}")

# DBSCAN
db = DBSCAN(eps=0.5, min_samples=5).fit_predict(X_pca)
print(f"Clusters: {len(set(db))-1}, Outliers: {(db==-1).sum()}")

# Imbalanced dataset handling
from imblearn.over_sampling import SMOTE
X_res, y_res = SMOTE().fit_resample(X_train, y_train)
# Never evaluate on resampled data!

Key Assumptions & Pitfalls

• K-means assumes spherical clusters of similar size — use DBSCAN for arbitrary shapes
• Always scale features before PCA or clustering (StandardScaler)
• DBSCAN eps: use k-distance plot (sort k-NN distances, look for elbow)
• Simpson's paradox: ALWAYS segment by key confounders before drawing aggregate conclusions
• DAGs: conditioning on a collider CREATES bias; conditioning on a confounder REMOVES bias
• Never evaluate models on resampled (SMOTE) data — only resample the training set