Estimation and Hypothesis Testing

Given i.i.d. data $X_1, \ldots, X_n \sim P_\theta$, a point estimator $\hat{\theta}_n = T(X_1, \ldots, X_n)$ is any function of the data that produces a single guess for $\theta$. Common examples:

• Sample mean $\bar{X}_n$ estimates the population mean $\mu$.
• Sample variance $s^2 = \frac{1}{n-1}\sum(X_i - \bar{X})^2$ estimates $\sigma^2$.
• Empirical proportion $\hat{p} = k/n$ estimates a Bernoulli parameter.

An estimator is itself a random variable — it changes across different datasets. Two questions define its quality:

1. Bias: is its expectation equal to the true $\theta$?
2. Variance: how much does it fluctuate across repeated samples?

Notice $n-1$ in $s^2$ instead of $n$. That's not a typo — it's the correction that makes $s^2$ unbiased. The naive $1/n$ version systematically underestimates $\sigma^2$.

The bias of $\hat{\theta}$ is $\text{Bias}(\hat{\theta}) = \mathbb{E}[\hat{\theta}] - \theta$. An estimator is unbiased if $\text{Bias} = 0$.

The mean squared error decomposes cleanly: $\text{MSE}(\hat{\theta}) = \mathbb{E}[(\hat{\theta} - \theta)^2] = \underbrace{\text{Var}(\hat{\theta})}_{\text{variance}} + \underbrace{(\text{Bias}(\hat{\theta}))^2}_{\text{bias}^2}.$

This is the bias-variance decomposition for estimation. A biased estimator can have lower MSE if it dramatically reduces variance — the logic behind ridge regression, early stopping, and shrinkage estimators like James-Stein.

Example: The sample mean $\bar{X}_n$ is unbiased ($\text{Bias}=0$) with $\text{Var}(\bar{X}_n) = \sigma^2/n$. As $n$ grows, variance shrinks to zero and $\hat{\theta} \to \theta$ in probability — this is consistency.

Principle: choose $\theta$ that makes the observed data most probable.

$\hat{\theta}_{\text{MLE}} = \arg\max_\theta\, \mathcal{L}(\theta) = \arg\max_\theta \prod_{i=1}^n P(X_i \mid \theta).$

Taking logs (monotone transform, doesn't change the argmax): $\hat{\theta}_{\text{MLE}} = \arg\max_\theta \sum_{i=1}^n \log P(X_i \mid \theta).$

Log-likelihood is almost always easier to maximize. Key properties of MLE:

• Consistent: $\hat{\theta}_{\text{MLE}} \to \theta^\star$ as $n \to \infty$ (under regularity).
• Asymptotically normal: $\sqrt{n}(\hat{\theta}_{\text{MLE}} - \theta^\star) \xrightarrow{d} \mathcal{N}(0, I(\theta^\star)^{-1})$ where $I(\theta)$ is the Fisher information.
• Asymptotically efficient: achieves the Cramér-Rao lower bound.
• Invariant: if $\hat{\theta}$ is MLE of $\theta$, then $g(\hat{\theta})$ is MLE of $g(\theta)$.

Bernoulli$(p)$: observe $k$ successes in $n$ trials. Log-likelihood: $\ell(p) = k \log p + (n-k) \log(1-p).$ Set $d\ell/dp = k/p - (n-k)/(1-p) = 0$. Solve: $\hat{p}_{\text{MLE}} = k/n$. The sample proportion — obvious in hindsight, derived rigorously here.

Gaussian $\mathcal{N}(\mu, \sigma^2)$: log-likelihood is $\ell(\mu, \sigma^2) = -\frac{n}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^n (X_i - \mu)^2.$ Maximizing over $\mu$: $\hat{\mu} = \bar{X}_n$. Maximizing over $\sigma^2$: $\hat{\sigma}^2 = \frac{1}{n}\sum_i(X_i - \bar{X})^2$.

Note: $\hat{\sigma}^2_{\text{MLE}}$ uses $1/n$, making it biased (underestimates by factor $(n-1)/n$). The unbiased fix is $s^2 = \frac{n}{n-1}\hat{\sigma}^2$.

import numpy as np
np.random.seed(0)
data = np.random.normal(loc=5.0, scale=2.0, size=100)
print(f"MLE mu    = {data.mean():.3f}")
print(f"MLE sigma2 (biased) = {data.var(ddof=0):.3f}")
print(f"Unbiased s2         = {data.var(ddof=1):.3f}")

A $(1-\alpha)$ confidence interval $[L, U]$ satisfies $P(L \le \theta \le U) = 1 - \alpha$ — but $\theta$ is fixed, so this is a statement about the procedure: if you repeated the experiment many times, $(1-\alpha)$ fraction of the intervals you construct would contain $\theta$.

Common mistake: "There is a 95% probability that $\theta$ is in this specific interval." Wrong. Once you have the data and compute $[L, U]$, $\theta$ is either in it or it isn't — no probability involved. The confidence is in the method, not the result.

For a Gaussian mean with known $\sigma$, the 95% CI is: $\bar{X}_n \pm 1.96 \cdot \frac{\sigma}{\sqrt{n}}.$ With unknown $\sigma$ (the real case), replace $1.96$ with $t_{0.975, n-1}$ and $\sigma$ with $s$:

from scipy import stats
import numpy as np
data = np.random.normal(5, 2, 30)
ci = stats.t.interval(0.95, df=len(data)-1,
                      loc=data.mean(), scale=stats.sem(data))
print(f"95% CI: ({ci[0]:.2f}, {ci[1]:.2f})")

A hypothesis test formalizes the question "is this effect real or noise?"

1. Null hypothesis $H_0$: the default (no effect, no difference).
2. Alternative $H_1$: what you're trying to detect.
3. Test statistic $T$: a function of the data designed to be large when $H_1$ is true.
4. Rejection region: reject $H_0$ if $T > c$ (or $|T| > c$ for two-sided tests).

Two error types:

• $\alpha$ = significance level = tolerated false-positive rate (usually 0.05).
• Power $= 1 - \beta$ = probability of detecting a real effect.
• Reducing $\alpha$ increases $\beta$ — you can't minimize both without more data.

The p-value is: $p = P(T \ge t_{\text{obs}} \mid H_0 \text{ true})$ the probability of seeing data at least as extreme as what you observed, assuming the null is true.

You reject $H_0$ if $p < \alpha$.

What p-values are NOT:

• $p$ is not the probability that $H_0$ is true.
• $p$ is not the probability that the result was due to chance.
• $p < 0.05$ does not imply the effect is large or practically important.
• $p > 0.05$ does not mean $H_0$ is true — it means the data is consistent with $H_0$.

A p-value of 0.03 with $n = 1{,}000{,}000$ may correspond to a completely negligible effect size. Always pair hypothesis tests with effect sizes (Cohen's $d$, relative risk, etc.) and confidence intervals.

from scipy import stats
group_a = [0.52, 0.61, 0.58, 0.55, 0.60]
group_b = [0.48, 0.51, 0.49, 0.53, 0.47]
t_stat, p_val = stats.ttest_ind(group_a, group_b)
print(f"t={t_stat:.3f}, p={p_val:.3f}")

One-sample z-test (known $\sigma$): test $H_0: \mu = \mu_0$ using $Z = \frac{\bar{X}_n - \mu_0}{\sigma / \sqrt{n}} \sim \mathcal{N}(0,1) \text{ under } H_0.$

One-sample t-test (unknown $\sigma$): replace $\sigma$ with $s$: $T = \frac{\bar{X}_n - \mu_0}{s / \sqrt{n}} \sim t_{n-1} \text{ under } H_0.$ The $t$-distribution has heavier tails than the standard normal, reflecting uncertainty about $\sigma$. With $n > 30$ the difference is negligible in practice.

Two-sample t-test: compares means of two independent groups. Under $H_0: \mu_1 = \mu_2$, $T = \frac{\bar{X}_1 - \bar{X}_2}{\sqrt{s_1^2/n_1 + s_2^2/n_2}} \sim t_{\nu} \text{ (Welch approximation for df }\nu\text{)}.$

In ML: use two-sample $t$-tests to compare two models' accuracy across $k$-fold cross-validation runs. Paired t-test (one fold = one pair) is preferred when variance across folds is large.

Run 20 independent tests each at $\alpha = 0.05$. Under the null, the probability of at least one false positive is $1 - 0.95^{20} \approx 64\%$. This is the multiple testing problem — a major source of irreproducibility in science and A/B testing.

Corrections:

• Bonferroni: use threshold $\alpha/m$ for $m$ tests. Conservative but simple.
• Benjamini-Hochberg (BH): controls the false discovery rate (FDR) = expected fraction of rejected nulls that are false positives. Less conservative than Bonferroni.

In ML specifically: hyperparameter sweeps, feature selection, neural architecture search, and leaderboard evaluation are all multiple-comparison problems. Always report whether corrections were applied.

from statsmodels.stats.multitest import multipletests
import numpy as np
np.random.seed(1)
pvals = np.random.uniform(size=50)  # all null
reject, pvals_adj, _, _ = multipletests(pvals, alpha=0.05, method='fdr_bh')
print(f"Rejected: {reject.sum()} out of 50")
# Should be ~0-3 (FDR controls the fraction, not the count)