📈 Statistics

Probability goes from a known model to data: "if the coin is fair, how many heads?" Statistics goes the other way: "I saw 7 heads in 10 flips, what is the coin?" Two ideas cover most interview questions: maximum likelihood (how to pick an estimate) and bias & variance (how to judge it).

Maximum likelihood estimation

We observe data $x_1,\dots,x_n$, drawn i.i.d. from a distribution with an unknown parameter $\theta$. Maximum likelihood says: pick the $\theta$ under which the data we actually saw would have been most probable.

$$L(\theta)=\prod_{i=1}^n p(x_i;\theta),\qquad \hat\theta_{\text{MLE}}=\arg\max_\theta L(\theta).$$

The log-likelihood trick

Products are awkward to differentiate, and multiplying hundreds of small probabilities underflows to 0 on a computer. Since $\log$ is increasing, maximizing $\log L$ gives the same $\hat\theta$, and it turns the product into a sum:

$$\ell(\theta)=\log L(\theta)=\sum_{i=1}^n\log p(x_i;\theta).$$

Recipe

1. Write down $\ell(\theta)=\sum_i \log p(x_i;\theta)$ and simplify.
2. Differentiate, set $\ell'(\theta)=0$, and solve for $\hat\theta$.
3. Check it really is a maximum: $\ell''(\hat\theta)\lt 0$, and compare with the endpoints / boundaries of the allowed parameter range. Sometimes the maximum is on the boundary, where the derivative isn't zero at all (see the uniform example below).

Example: Bernoulli (coin flips)

$x_i\in\{0,1\}$ i.i.d. $\Bern(p)$, with PMF $P(X=x;p)=p^x(1-p)^{1-x}$. Let $k=\sum_i x_i$ be the number of 1s.

$$\ell(p)=\sum_{i=1}^n\log\big[p^{x_i}(1-p)^{1-x_i}\big]=\sum_{i=1}^n\big[x_i\log p+(1-x_i)\log(1-p)\big]=k\log p+(n-k)\log(1-p).$$ $$\ell'(p)=\frac kp-\frac{n-k}{1-p}=0\quad\Longrightarrow\quad \hat p=\frac kn.$$

The MLE is just the sample frequency, which is what you would have guessed.

Check it's a maximum

$\ell''(p)=-\dfrac{k}{p^2}-\dfrac{n-k}{(1-p)^2}\lt 0$ for all $p\in(0,1)$, so $\ell$ is concave and the critical point is the global max. Edge case: if $k=0$, then $\ell(p)=n\log(1-p)$ is decreasing and there is no interior critical point; the max is at the boundary $\hat p=0$, which still equals $k/n$.

Example: Exponential (waiting times)

$x_i$ i.i.d. $\Exp(\lambda)$ with density $f(x;\lambda)=\lambda e^{-\lambda x}$.

$$\ell(\lambda)=\sum_{i=1}^n\log\big(\lambda e^{-\lambda x_i}\big)=\sum_{i=1}^n\big(\log\lambda-\lambda x_i\big)=n\log\lambda-\lambda\sum_i x_i.$$ $$\ell'(\lambda)=\frac n\lambda-\sum_i x_i=0\quad\Longrightarrow\quad \hat\lambda_{\text{MLE}}=\frac{n}{\sum_i x_i}=\frac1{\bar x}.$$

The MLE is the reciprocal of the sample mean. That makes sense because $\Exp(\lambda)$ has mean $1/\lambda$: if buses arrive every 10 minutes on average, the rate is 0.1 per minute. ($\ell''=-n/\lambda^2\lt 0$, so it's a max.)

Example: Uniform(0, θ), where calculus fails

$x_i$ i.i.d. $\Unif(0,\theta)$ with density $1/\theta$ on $[0,\theta]$. If $\theta$ is smaller than any data point, that point would have probability 0, so

$$L(\theta)=\begin{cases}\dfrac{1}{\theta^n} & \text{if } \theta\ge\max_i x_i\\[4pt] 0&\text{otherwise.}\end{cases}$$

$1/\theta^n$ is decreasing, so we want $\theta$ as small as allowed: $\hat\theta=\max_i x_i$. Setting the derivative to zero would get you nowhere ($-n\theta^{-n-1}\ne0$); the maximum is at the boundary of the feasible region. This is exactly why the recipe says to check boundaries.

But the MLE here is biased: the max of the sample is always below $\theta$, so it underestimates on average.

$$\E[\max_i x_i]=\frac{n}{n+1}\,\theta \qquad\Longrightarrow\qquad \frac{n+1}{n}\max_i x_i \text{ is unbiased.}$$

Proof of $\E[\max]=\frac{n}{n+1}\theta$

$P(\max\le x)=P(\text{all }x_i\le x)=(x/\theta)^n$, so the density is $n x^{n-1}/\theta^n$ on $[0,\theta]$ and

$$\E[\max]=\int_0^\theta x\cdot\frac{n x^{n-1}}{\theta^n}\,dx=\frac{n}{\theta^n}\cdot\frac{\theta^{n+1}}{n+1}=\frac{n}{n+1}\theta.$$

Shortcut: $n$ uniform points cut $[0,\theta]$ into $n+1$ gaps with equal expected length $\theta/(n+1)$, and the max sits one gap below $\theta$ (see uniform).

MLE = minimizing cross-entropy loss

Neural nets and logistic regression are trained to minimize the cross-entropy (negative log-likelihood) loss:

$$\text{CE}(\theta)=-\frac1N\sum_{i=1}^N\log p_\theta(x_i)\qquad\Longleftrightarrow\qquad \text{maximize }\ \ell(\theta)=\sum_{i=1}^N\log p_\theta(x_i).$$

Minimizing one is the same as maximizing the other (the $-\frac1N$ is just a negative constant). In words: models are trained to make the observed data as probable as possible. For classification the "data" are labels, $p_\theta(y_i\mid x_i)$. Similarly, minimizing squared error is MLE under Gaussian noise, since $-\log\N(y;\hat y,\sigma^2)=\frac{(y-\hat y)^2}{2\sigma^2}+\text{const}$.

MAP: MLE plus a prior

MLE uses only the data. If you have a prior belief $p(\theta)$, Bayes' rule gives the posterior $p(\theta\mid x)\propto p(x\mid\theta)\,p(\theta)$, and the maximum a posteriori (MAP) estimate maximizes it:

$$\hat\theta_{\text{MAP}}=\arg\max_\theta\Big[\underbrace{\ell(\theta)}_{\text{data fit}}+\underbrace{\log p(\theta)}_{\text{prior = regularizer}}\Big].$$

• Coin with a $\text{Beta}(\alpha,\beta)$ prior: $\hat p_{\text{MAP}}=\dfrac{k+\alpha-1}{n+\alpha+\beta-2}$. A $\text{Beta}(2,2)$ prior acts like one extra head and one extra tail ("add-one smoothing").
• Gaussian prior on weights → $\log p(w)=-\frac{\|w\|^2}{2\tau^2}+\text{const}$ → L2 regularization (weight decay). A Laplace prior gives L1.
• A flat prior gives MAP = MLE. As $n\to\infty$, the data term dominates and MAP → MLE.

Bias & variance

An estimator $\hat\theta$ is a function of the data. The data are random, so $\hat\theta$ is a random variable: rerun the experiment and you get a different estimate. We judge an estimator by the distribution of its estimates around the true $\theta$.

Proof of the decomposition

Add and subtract $\E[\hat\theta]$ inside the square:

$$\begin{aligned} \E[(\hat\theta-\theta)^2]&=\E\Big[\big((\hat\theta-\E\hat\theta)+(\E\hat\theta-\theta)\big)^2\Big]\\ &=\underbrace{\E[(\hat\theta-\E\hat\theta)^2]}_{\Var(\hat\theta)}+2(\E\hat\theta-\theta)\underbrace{\E[\hat\theta-\E\hat\theta]}_{=0}+\underbrace{(\E\hat\theta-\theta)^2}_{\operatorname{Bias}^2}. \end{aligned}$$

The cross term vanishes because $\E\hat\theta-\theta$ is a constant and $\hat\theta$ minus its own mean has mean zero. It's the same identity as $\E[X^2]=\Var(X)+\E[X]^2$, applied to $X=\hat\theta-\theta$.

Worked example: three ways to estimate θ in Uniform(0, θ)

Estimator	Bias	Variance	MSE
$\max_i x_i$ (MLE)	$-\dfrac{\theta}{n+1}$	$\dfrac{n\theta^2}{(n+1)^2(n+2)}$	$\dfrac{2\theta^2}{(n+1)(n+2)}$
$\frac{n+1}{n}\max_i x_i$	$0$	$\dfrac{\theta^2}{n(n+2)}$	$\dfrac{\theta^2}{n(n+2)}$
$2\bar x$ (method of moments)	$0$	$\dfrac{\theta^2}{3n}$	$\dfrac{\theta^2}{3n}$

Both corrected estimators are unbiased, but the max-based one has variance shrinking like $1/n^2$, versus $1/n$ for $2\bar x$. The max uses the most informative statistic. (Fun fact: the MSE-optimal multiple of the max is $\frac{n+2}{n+1}\max$, which is slightly biased.)

Where the variances come from

$\E[\max^2]=\int_0^\theta x^2\,\frac{nx^{n-1}}{\theta^n}dx=\frac{n}{n+2}\theta^2$, so $\Var(\max)=\frac{n}{n+2}\theta^2-\frac{n^2}{(n+1)^2}\theta^2=\frac{n\theta^2}{(n+1)^2(n+2)}$. Scaling by $\frac{n+1}{n}$ multiplies the variance by $\frac{(n+1)^2}{n^2}$.

For the mean: $\Var(x_i)=\theta^2/12$, so $\Var(2\bar x)=4\cdot\frac{\theta^2}{12n}=\frac{\theta^2}{3n}$.