Statistical Inference and Hypothesis Testing Questions

For independent, identically distributed random variables with finite variance, the Central Limit Theorem says the (properly scaled) sample mean converges in distribution to a normal distribution as the sample size grows, regardless of the shape of the original distribution. In practice this means confidence intervals and z/t-based tests built on a sample mean are approximately valid once n is "large enough," even when the underlying data itself is skewed or non-normal.

Structured elaboration

Precise statement

Let X1, X2, ..., Xn be IID with mean mu and finite variance sigma^2 less than infinity. Let X_bar_n be the sample mean. Then as n approaches infinity:

$\sqrt{n}\,(\bar X_n - \mu) \;\xrightarrow{d}\; N(0, \sigma^2)$

Equivalently, for large n:

$\bar X_n \approx N\!\left(\mu, \frac{\sigma^2}{n}\right)$

so the standard error of the mean is sigma divided by the square root of n, shrinking as n grows regardless of the shape of the original distribution X.

Implications for practice

• The sampling distribution of a mean becomes approximately bell-shaped even when individual observations are heavily skewed (latency, revenue, task time).
• This is what licenses z-intervals, t-tests, and z-tests on averages and proportions, even for non-normal raw data.
• Convergence is not instantaneous or uniform: the more skewed or heavy-tailed the underlying distribution, the larger n needs to be before the normal approximation is trustworthy.

Rule of thumb

• Roughly symmetric, light-tailed data: n of 30 or more is often enough.
• Noticeably skewed data: n of 50 to 100.
• Heavy-tailed or extremely skewed data (or when tail behavior itself matters): treat the CLT-based interval with caution regardless of n, and prefer bootstrap or robust methods.

Worked example

Simulate the sampling distribution of the mean of an Exponential(1) distribution (population skewness = 2, a strongly right-skewed case) for growing n, and track how quickly its skewness and standard error approach the CLT prediction:

import random, math

random.seed(11)

def sample_exp(rate=1.0):
    u = random.random()
    return -math.log(1 - u) / rate

def skewness(xs):
    n = len(xs)
    m = sum(xs) / n
    m2 = sum((x - m) ** 2 for x in xs) / n
    m3 = sum((x - m) ** 3 for x in xs) / n
    return m3 / (m2 ** 1.5)

R = 5000
for n in [5, 10, 30, 100, 500]:
    means = [sum(sample_exp() for _ in range(n)) / n for _ in range(R)]
    mean_of_means = sum(means) / R
    sd_of_means = math.sqrt(sum((x - mean_of_means) ** 2 for x in means) / R)
    print(n, sd_of_means, skewness(means))

Running this (seed 11, R = 5,000 replicates per n) gives:

Both the standard error and the skewness of the sample-mean distribution shrink toward the CLT prediction as n grows, even starting from a distribution with skewness 2. By n = 500 the skewness of the sampling distribution of the mean (about 0.11) is far closer to a normal distribution's skewness of 0 than the original Exponential's skewness of 2.

Trade-offs & pitfalls

• The CLT describes the sampling distribution of the mean, not the distribution of individual observations; individual data points from a skewed distribution stay skewed no matter how large n gets.
• "n of 30 or more" is a rough convention, not a law; a heavily skewed or heavy-tailed distribution (or ratio metrics, or rare-event counts) can need a much larger n before the normal approximation is trustworthy.
• Sums also converge under the CLT, not just means, but variance-based reasoning must use the right scaling (variance of a sum grows with n; variance of a mean shrinks with n).
• CLT-based intervals are approximate. When in doubt, especially for skewed metrics with a business-critical decision riding on the interval, a bootstrap confidence interval is a useful cross-check.

Direct answer

At the scale of thousands of experiments and hundreds of metrics each, the design has to separate metrics into tiers by decision importance (a handful of primary, ship/no-ship metrics per experiment vs. many secondary and exploratory ones), apply a stricter false discovery control to the important tier and a looser, FDR-based control to the exploratory tier, and back all of it with an immutable experiment registry so results are reproducible and auditable. Interpretability comes from tiering and clear reporting conventions, not from picking one universal statistical procedure for every metric.

Structured elaboration

Pipeline flow

flowchart LR
  A[Experiment registry: pre-registered plan, metrics, hypotheses] --> B[Metric pipeline: batch-compute per experiment]
  B --> C{Metric tier}
  C -->|Primary, ~1-3 per experiment| D[Per-experiment alpha, no correction needed within this small set]
  C -->|Secondary/exploratory, hundreds| E[Benjamini-Hochberg within the metric family]
  D --> F[Results store: adjusted + unadjusted, effect size, confidence interval]
  E --> F
  F --> G[Dashboard: verdict, q-value, effect size by segment]
  F --> H[Continuous monitoring: realized FDR, power, registry-plan mismatches]

Statistical procedures

• Primary metrics (the one to three metrics that actually decide ship/no-ship for a given experiment) are tested at the experiment's own pre-registered $\alpha$, since there are only a few of them per experiment and they were chosen before the data was seen; multiplicity correction here is unnecessary if the primary set stays small and is genuinely fixed in advance.
• Secondary and exploratory metrics (the hundreds of additional metrics reported per experiment) are controlled with Benjamini-Hochberg (BH) within that family, since the goal there is to surface promising leads for follow-up, not to make a single high-stakes decision, and FDR control keeps a bounded, known fraction of those leads spurious.
• Cross-experiment aggregation (e.g. "how many experiments this quarter showed a significant lift on metric X") is its own separate family and needs its own correction; don't let a metric's significance in one experiment's local family imply anything about a cross-experiment claim without re-testing at that level.

Data infrastructure

• An experiment registry that stores the hypothesis, the declared primary and secondary metrics, the randomization scheme, and the analysis plan before the experiment starts, and treats that record as immutable once registered.
• Raw event logging (not just pre-aggregated metrics), so any later question about what actually happened can be answered by recomputation rather than trusting a snapshot.
• Reproducible, versioned pipelines (a workflow orchestrator with pinned code, data snapshots, and seeds) so a reported result can be regenerated exactly.

Reporting conventions

• One decision-relevant verdict per primary metric per experiment (ship / hold / inconclusive), always paired with its effect size and confidence interval, not the p-value alone.
• Secondary and exploratory metrics reported with their BH-adjusted q-value alongside the raw p-value, explicitly labeled as hypothesis-generating rather than confirmatory.
• A short plain-language note translating each significant primary result into business terms, so the interpretability requirement is met at the reporting layer, not just the statistical one.

Worked example

A single experiment reports 12 exploratory secondary metrics with unadjusted p-values: $0.0009, 0.006, 0.041, 0.09, 0.14, 0.19, 0.27, 0.33, 0.41, 0.58, 0.71, 0.85$. Comparing Bonferroni (family $\alpha=0.05$) against BH (target FDR $q=0.10$):

(both computed directly from the sorted p-values and the two threshold rules above)

At platform scale, this difference compounds: across thousands of experiments each reporting hundreds of exploratory metrics, a flat Bonferroni-style correction across everything would leave the exploratory tier almost powerless, while a per-experiment BH pass on just that tier keeps a known, bounded false-discovery share without needing an ever-shrinking per-test threshold as the platform grows.

Computational scalability

• Metric computation is embarrassingly parallel per experiment and per metric; pre-aggregate with a batch or streaming pipeline, and only run the multiplicity correction step (cheap, since it's just a sort and threshold comparison) once all p-values for a family are in.
• The correction step itself is $O(m\log m)$ per family (dominated by sorting the p-values), trivial compared to the cost of computing the underlying statistics; scalability bottlenecks are almost always in the metric computation and data pipeline, not the multiplicity correction.
• Reserve heavier machinery (empirical-Bayes shrinkage, hierarchical pooling across similar experiments) for cases where the added interpretability or power is worth the added modeling and validation cost; default to plain BH within each experiment's exploratory family as the baseline that always works.

Monitoring strategies

• Backtest the pipeline on historical experiments with known outcomes (including synthetic A/A tests, where the two arms are identical) to confirm the realized false-positive rate matches what the correction claims.
• Continuously track the realized discovery rate over time; a sudden jump is a signal of either a real shift in how many effects are real, or a broken pipeline inflating false positives.
• Alert on registry-plan mismatches (a primary metric changed after the experiment started, or a metric reported that wasn't pre-registered), since that's the most common way multiplicity control gets silently defeated.

Trade-offs & pitfalls

• Applying a single strict correction (Bonferroni across everything, primary and secondary combined) is the safest failure mode but starves the exploratory tier of power at scale; applying no correction anywhere invites exactly the selective-reporting problem multiplicity control exists to prevent. The tiered approach is a deliberate middle ground, not a shortcut.
• Letting a metric drift from "exploratory" to "primary" after seeing a promising result defeats the whole design; the tier has to be fixed by the pre-registration, not re-assigned post hoc.
• BH's guarantee weakens under strong negative dependence between metrics; correlated metrics (which are common when many metrics are derived from the same few underlying events) may need a dependence-aware variant or an empirical-Bayes model rather than plain BH.
• None of this replaces basic experiment hygiene: a broken randomization or a biased sample invalidates every downstream statistical control, no matter how carefully the multiplicity correction is designed.

Direct answer

A Binomial$(n, p)$ distribution approximates a Poisson$(\lambda)$ distribution when $n$ is large, $p$ is small, and their product $\lambda = np$ stays moderate and roughly fixed as $n$ grows. Practically: $n$ in the hundreds or more, $p$ under about 0.05, with $np$ somewhere in the single digits to low tens works well; this is exactly the regime of "rare events counted over many independent opportunities," like defects per kilometer of cable or errors across many server requests.

Structured elaboration

The mathematical rule

Intuition for why: with $p$ small, the chance of two or more "successes" landing on the same trial-count window becomes negligible, so the combinatorial factor $\binom{n}{k}$ and the $(1-p)^{n-k}$ term simplify. In particular, $(1-p)^n \approx e^{-np} = e^{-\lambda}$ for small $p$, which is the Poisson kernel's leading term.

Practical thresholds

The mean check ($np$ moderate) is the most reliable single diagnostic, since $n$ and $p$ can each individually look "extreme" while their product is not in the rare-event regime.

Worked example

Rare server errors: 10,000 requests in an hour, each independently erroring with probability $p = 0.0003$.

Compare the probability of zero errors under each model:

import math
n, p = 10000, 0.0003
lam = n*p                      # 3.0 (up to floating-point noise)

p0_binom = (1-p)**n            # 0.04976466474747146
p0_pois = math.exp(-lam)       # 0.049787068367863965
abs(p0_binom - p0_pois)        # 2.24e-05

Both models give $\approx 0.0498$ for the probability of a clean hour, agreeing to four significant figures despite one requiring the exact binomial formula and the other only needing $\lambda = 3$. For an applied model of rare server errors or defect counts, this means you can drop the exact binomial machinery (which needs both $n$ and $p$ tracked separately) and work directly with the single-parameter Poisson$(\lambda)$, which is simpler to fit, simpler to reason about (mean = variance = $\lambda$), and simpler to compose across time windows (Poisson counts over disjoint intervals just add their $\lambda$'s).

Trade-offs & pitfalls

• The approximation breaks down when $p$ isn't actually small, even if $n$ is large - e.g. a coin-flip-like process with $p=0.4$ over 1000 trials is nowhere near Poisson-shaped regardless of how large $n$ is; check $p$, not just $n$.
• It also breaks down at small $n$ even with small $p$, if that makes $\lambda = np$ too small to be a meaningful count (e.g., $n=20, p=0.01 \to \lambda=0.2$) - the approximation quality depends on the full regime, not any single threshold in isolation.
• Real event counts are often overdispersed relative to Poisson. Poisson forces mean = variance; if observed variance is substantially larger than the mean (common when the underlying rate itself varies, e.g., error rates spike during incidents rather than staying constant), a negative binomial model fits better and a plain Poisson will understate uncertainty.
• Always validate numerically for decisions that matter, the way the worked example does above, rather than trusting the heuristic thresholds blindly - the thresholds are guidelines, not guarantees, and the actual approximation error is cheap to check directly.

Quick answer

Each user receiving 10 independent emails, each clicked with probability $p = 0.3$, is a Binomial$(n=10, p=0.3)$ model. The probability of at least 3 clicks is easiest to get via the complement of 0, 1, or 2 clicks, which comes out to about 61.7%. The expected number of clicks is $np = 3$, with variance $np(1-p) = 2.1$.

Setup

$X \sim \text{Binomial}(n=10,\ p=0.3)$, with PMF:

(a) $P(X \ge 3)$

Computing the complement is less arithmetic than summing 8 terms directly:

from math import comb

n, p = 10, 0.3
p0 = comb(n, 0) * p**0 * (1 - p)**10
p1 = comb(n, 1) * p**1 * (1 - p)**9
p2 = comb(n, 2) * p**2 * (1 - p)**8
p_ge3 = 1 - (p0 + p1 + p2)
print(f"P(X=0)={p0:.6f}, P(X=1)={p1:.6f}, P(X=2)={p2:.6f}")
print(f"P(X>=3)={p_ge3:.4f}")
# P(X=0)=0.028248, P(X=1)=0.121061, P(X=2)=0.233474
# P(X>=3)=0.6172

$P(X \ge 3) \approx 0.6172$: roughly 61.7% of users click at least 3 of the 10 emails.

(b) Expected value and variance

For a Binomial random variable:

Standard deviation $= \sqrt{2.1} \approx 1.449$.

Interpretation

An average user clicks 3 of the 10 emails, with a standard deviation of about 1.45 clicks, so a typical user's click count lands somewhere between about 1.5 and 4.5. The fact that over 60% of users clear the "at least 3" bar despite the average also being 3 reflects the right-skew-free, moderately spread-out shape of this particular Binomial (n=10, p=0.3 isn't close enough to either 0 or 1 to be heavily skewed).

Trade-offs and pitfalls

• Independence across the 10 emails is the load-bearing assumption. If a user who clicks one email is more likely to click the next (engagement momentum) or less likely (fatigue), the true distribution is overdispersed or underdispersed relative to Binomial, and the variance formula above will be wrong.
• $p=0.3$ is itself an estimate, not a known constant. In practice $p$ is estimated from historical click data, and the Binomial calculation above doesn't account for uncertainty in that estimate; a more careful treatment would put a Beta prior on $p$ and integrate over it (a Beta-Binomial), widening the resulting probability estimates.

Direct answer

When daily user metrics are autocorrelated, a test that assumes independent observations understates the true variance of the sample mean and inflates Type I error, because consecutive days carry redundant, not fresh, information. The fix is to either restore independence by aggregating to a unit that genuinely is independent (e.g. one summary per user, not per user-day), or to keep the daily structure but use a variance estimator built to handle dependence (block bootstrap, Newey-West/HAC standard errors, or an explicit time-series model).

Structured elaboration

Why naive variance is wrong. For i.i.d. observations, $\mathrm{Var}(\bar x)=\sigma^2/n$. Autocorrelated observations carry positive covariance terms that the i.i.d. formula ignores, so the true variance of the mean is larger, and a naive t-test's denominator is too small - it rejects the null more often than its nominal $\alpha$ would suggest.

Four approaches, with their cost to effective sample size and power:

In plain terms, Newey-West's "kernel-weighted sum of autocovariances" means: compute how correlated the series is with a copy of itself shifted by 1 day, by 2 days, and so on (each of these is an autocovariance), then add a weighted sum of them into the SE formula, with the weights (the kernel) shrinking smoothly toward zero at longer lags rather than cutting off abruptly at a fixed lag.

Effective sample size. For a stationary AR(1) process (one whose average level and variability do not drift over time, so a single autocorrelation value describes the whole series) with lag-1 autocorrelation $\rho$, the effective sample size behind the mean of $n$ observations is approximately:

which is far below $n$ whenever $\rho$ is sizable - the whole reason naive variance formulas fail here.

Worked example

Simulating a 30-day AR(1) metric with $\phi=0.6$ day-to-day autocorrelation, comparing the naive i.i.d. SE formula to the true (empirical) SE of the 30-day mean:

import numpy as np

def simulate_ar1(rng, n, phi, sigma_innov):
    x = np.empty(n)
    x[0] = rng.normal(0, sigma_innov/np.sqrt(1-phi**2))
    for t in range(1, n):
        x[t] = phi*x[t-1] + rng.normal(0, sigma_innov)
    return x

n_days, phi, sigma_innov = 30, 0.6, 1.0
ess_theory = n_days * (1-phi)/(1+phi)

R, master = 20000, np.random.default_rng(31)
means = np.array([simulate_ar1(master, n_days, phi, sigma_innov).mean() for _ in range(R)])

marginal_var = sigma_innov**2 / (1-phi**2)
naive_se = np.sqrt(marginal_var / n_days)
empirical_se = means.std(ddof=1)
understatement_factor = naive_se / empirical_se
implied_ess = n_days * understatement_factor**2

print(f"AR(1) phi={phi}, n={n_days} days")
print(f"theoretical effective sample size (large-n formula)   = {ess_theory:.2f}   (vs naive n={n_days})")
print()
print(f"naive iid SE of the {n_days}-day mean                        = {naive_se:.4f}")
print(f"empirical SE of the {n_days}-day mean (Monte Carlo, R={R}) = {empirical_se:.4f}")
print(f"understatement factor (naive/empirical)                 = {understatement_factor:.3f}")
print(f"implied ESS from empirical SE: n*(naive_se/empirical_se)^2 = {implied_ess:.2f}")

Output:

AR(1) phi=0.6, n=30 days
theoretical effective sample size (large-n formula)   = 7.50   (vs naive n=30)

naive iid SE of the 30-day mean                        = 0.2282
empirical SE of the 30-day mean (Monte Carlo, R=20000) = 0.4437
understatement factor (naive/empirical)                 = 0.514
implied ESS from empirical SE: n*(naive_se/empirical_se)^2 = 7.94

The naive formula understates the true SE by about a factor of two, matching the theoretical prediction that only about 7-8 of the 30 "daily" observations are behaving like independent information (7.50 theoretical vs 7.94 empirically implied - close agreement, with the small gap explained by the theoretical formula being a large-$n$ asymptotic approximation applied at a modest $n=30$). A naive t-test on this data would use a denominator roughly twice too small, directly inflating the false-positive rate.

Trade-offs & pitfalls

Per-user aggregation is the safest default when there are enough independent users, but it silently assumes the treatment effect is constant over the observation window - if the effect ramps up or decays over time, collapsing to one number per user discards exactly the information needed to see that. Block bootstrap and Newey-West both require choosing a lag/block length; too short under-corrects (leaves residual autocorrelation unaccounted for), too long over-corrects (wastes information and can bias small-sample SEs upward more than necessary) - there's no free parameter-free choice, and it should be reported alongside a sensitivity check across a few lag lengths. A full time-series model is the most information-efficient option when it fits, but a misspecified ARIMA order can understate uncertainty just as badly as ignoring autocorrelation altogether, so out-of-sample residual diagnostics (checking the residuals are no longer autocorrelated) are not optional.