Chapter 4: Probability & Statistics

Before you can assign a probability, you must say what could happen. The sample space $\Omega$ in the formula is the set of all outcomes of an experiment, and an event is any subset of it, such as "the die shows an even number", which is $\{2, 4, 6\}$. Events combine with set operations. The complement $A^c$ is everything outside $A$, the union $A \cup B$ means A or B, and the intersection $A \cap B$ means both. $P(\Omega) = 1$ says some outcome always happens. The usual failure is a badly chosen $\Omega$, such as counting "two heads, one of each, two tails" as three equally likely cases.

Why only three rules? Kolmogorov showed in 1933 that non-negativity, normalization, and additivity over disjoint events are enough to derive everything else. In the formula, $A$ and $B$ are events, $\Omega$ is the sample space, and $\emptyset$ is the empty event. Additivity gives the complement rule $P(A^c) = 1 - P(A)$ and $P(\emptyset) = 0$. Splitting $A \cup B$ into disjoint pieces gives inclusion-exclusion, $P(A \cup B) = P(A) + P(B) - P(A \cap B)$. A softmax layer satisfies the axioms by construction, because exponentials are positive and the outputs are normalized. Independent sigmoid scores do not, since their class scores need not sum to 1.

Learning that $B$ happened shrinks the world to the outcomes inside $B$. The formula renormalizes: $P(A \mid B)$ is the share of $B$'s probability that also lies in $A$, and it is defined only when $P(B) > 0$. Rearranging gives the product rule $P(A \cap B) = P(A \mid B) P(B)$, which the rest of the chapter uses constantly. Every classifier estimates a conditional probability such as $P(\text{spam} \mid \text{words})$, and a language model estimates $P(\text{next token} \mid \text{context})$. The classic failure is swapping $P(A \mid B)$ for $P(B \mid A)$, the prosecutor's fallacy from the spam example in the introduction.

Independence means one event carries no information about another. In the formula, $A \perp B$ holds exactly when the joint probability factors as $P(A) P(B)$, which is equivalent to $P(A \mid B) = P(A)$. Conditional independence applies the same test inside a context $C$: $P(A \cap B \mid C) = P(A \mid C) P(B \mid C)$. Neither kind implies the other. Two symptoms can be dependent overall yet independent given the disease that causes both. Naive Bayes assumes words are conditionally independent given the class, which cuts its parameter count from exponential to linear in the vocabulary. A frequent error is confusing independent with disjoint: disjoint events with positive probability are always dependent.

Some probabilities are hard to compute directly but easy once you know which case you are in. If $A_1, \ldots, A_n$ partition the sample space, meaning they are disjoint and together cover it, the formula writes $P(B)$ as a weighted average of the case probabilities $P(B \mid A_i)$, with weights $P(A_i)$. It follows from the product rule of conditional probability plus additivity. The same move marginalizes a latent variable, $P(x) = \sum_z P(x \mid z) P(z)$, which is how mixture models score data, and it produces the denominator of Bayes' theorem. The failure mode is a set of cases that overlap or leave gaps, so the weights no longer sum to 1.

The Bernoulli distribution models a single trial with two outcomes, coded 1 for success and 0 for failure. In the formula, $p$ is the success probability and $k \in \{0, 1\}$, so the expression returns $p$ when $k = 1$ and $1 - p$ when $k = 0$. The mean is $p$ and the variance is $p(1 - p)$, which peaks at 0.25 when $p = 0.5$, the point of greatest uncertainty. A sigmoid output is exactly a Bernoulli parameter, binary cross-entropy is its negative log-likelihood, and dropout masks units with Bernoulli draws. The model breaks down when outcomes are not binary, or when trials share a hidden cause such as one spam campaign.

Count the successes in $n$ independent Bernoulli trials that share the same success probability $p$, and you get the Binomial distribution. In the formula, $\binom{n}{k}$ counts the orderings with exactly $k$ successes, $p^k$ is the probability of those successes, and $(1-p)^{n-k}$ covers the failures. The mean is $np$ and the variance is $np(1-p)$. The shape is symmetric at $p = 0.5$ and approaches a Gaussian as $n$ grows. The number of correct predictions on a test set of $n$ examples is a Binomial count. The model fails when trials are dependent, for example when near-duplicate test images make errors arrive in clusters.

Many events are individually rare but have many chances to happen: requests hitting a server, typos on a page, photons striking a sensor pixel. The Poisson distribution counts such events in a fixed interval. In the formula, $\lambda$ is the average number of events per interval, $k$ is the observed count, and $k!$ accounts for orderings. Both the mean and the variance equal $\lambda$, and independent Poisson counts add their rates. It is the limit of a Binomial with large $n$ and small $p$, with $\lambda = np$. The telltale failure is overdispersion: if the observed variance far exceeds the mean, events are clustered and a Poisson model underestimates spikes.

How many attempts until the first success? If each independent trial succeeds with probability $p$, the formula says the first success lands on trial $k$ after $k - 1$ failures, each with probability $1-p$. The mean is $1/p$, the variance is $(1-p)/p^2$, and the CDF is $P(X \le k) = 1 - (1-p)^k$. The Geometric is the only discrete distribution that is memoryless: $P(X > s + t \mid X > s) = P(X > t)$, so ten failed trials do not bring success any closer. Watch the convention. Some libraries count the failures before the first success, which shifts the mean to $(1-p)/p$.

The PMF $p(k) = P(X = k)$ answers point questions, and the CDF $F(k)$ in the formula answers "at most" questions by summing the PMF up to $k$. Each can be recovered from the other: $P(X = k) = F(k) - F(k-1)$, and any interval probability is a difference, $P(a < X \le b) = F(b) - F(a)$. A CDF never decreases and climbs from 0 to 1. Percentiles, p-values, and calibration curves are all CDF lookups. The standard off-by-one error is using $F(b) - F(a)$ for $P(a \le X \le b)$, which silently drops the mass at $a$ for a discrete variable.

A density is probability per unit length. The first line of the formula says an interval's probability is the area under $f$ between $a$ and $b$. The other two lines are the density's rules: $f$ is never negative, and its total area is 1. Because a single point has zero width, $P(X = x) = 0$, which is why densities compare relative likelihoods rather than giving probabilities. The CDF $F(x) = \int_{-\infty}^{x} f(t)\,dt$ accumulates area, and differentiating it recovers $f$. The common error is reading $f(x)$ as a probability. A Gaussian with standard deviation 0.1 has a peak density of about 4, which is perfectly valid.

When all you know is that a value lies between $a$ and $b$, the Uniform distribution spreads density evenly. The formula gives height $1/(b-a)$ inside the interval and zero outside, so the area is 1. The mean is the midpoint $(a+b)/2$, the variance is $(b-a)^2/12$, and the CDF is linear, $(x-a)/(b-a)$. Among distributions confined to $[a, b]$ it has the maximum entropy. Random crops, random hyperparameter search, and Glorot initialization all use it. The trap is scale. A uniform range over a quantity spanning several orders of magnitude, such as a learning rate, puts almost all samples in the top decade.

In a Poisson process with rate $\lambda$ events per unit time, the wait until the next event follows an Exponential distribution. The density $\lambda e^{-\lambda x}$ in the formula starts at height $\lambda$ and decays, so short waits are the most common. The CDF is $1 - e^{-\lambda x}$, the mean is $1/\lambda$, and the variance is $1/\lambda^2$. Like the Geometric, the Exponential is memoryless: $P(X > s + t \mid X > s) = P(X > t)$. That fits radioactive decay and request arrivals. It does not fit hardware that wears out, where an old part is more likely to fail soon, so Weibull or Gamma models are used there.

How do you describe uncertainty about a probability itself, such as a click-through rate? The Beta distribution lives on $[0, 1]$. In the formula, $\alpha - 1$ and $\beta - 1$ act like counts of successes and failures, and the beta function $B(\alpha, \beta)$ normalizes the area to 1. The mean is $\alpha/(\alpha + \beta)$. $\text{Beta}(1, 1)$ is uniform, equal parameters above 1 peak at 0.5, and parameters below 1 give a U shape. After $h$ successes and $t$ failures, a Beta prior updates to $\text{Beta}(\alpha + h, \beta + t)$. A trap is reading the mean alone: $\text{Beta}(9, 5)$ and $\text{Beta}(90, 50)$ share it, but their standard deviations differ threefold.

The Gamma distribution models positive quantities built from several waits or small contributions, such as the time until the fifth phone call. In the formula, $\alpha$ is the shape, $\beta$ is the rate, and $\Gamma(\alpha)$ normalizes the area. The mean is $\alpha/\beta$ and the variance is $\alpha/\beta^2$. With $\alpha = 1$ it reduces to the Exponential, and for integer $\alpha$ it is the sum of $\alpha$ independent Exponential waits with rate $\beta$. The chi-squared distribution is a special case, and the Gamma is the conjugate prior for a Poisson rate. The main pitfall is parameterization: NumPy and SciPy take a scale $1/\beta$, not a rate.

Why is the Gaussian everywhere? The Central Limit Theorem says that the average $\bar{X}_n$ of $n$ independent draws with mean $\mu$ and finite variance $\sigma^2$ becomes approximately normal, with mean $\mu$ and variance $\sigma^2/n$, whatever the original distribution looks like. The arrow marked $d$ in the formula means convergence in distribution. The spread shrinks like $\sigma/\sqrt{n}$, so four times the data halves the error bar. This justifies Gaussian confidence intervals for test accuracy and the Gaussian model of mini-batch gradient noise. It fails for heavy-tailed data without a finite variance, such as the Cauchy, and converges slowly for strongly skewed data.

A joint distribution assigns probability to every combination of values. For two discrete variables it is a table of $P(X = x, Y = y)$, and for continuous variables it is a density surface $f(x, y)$. The formula states the two requirements: every entry is non-negative, and the double sum over $x$ and $y$ equals 1. The joint contains everything, since marginals, conditionals, and correlations all follow from it. In supervised learning, the training data are samples from an unknown joint $P(X, Y)$. The cost is size. A joint table over 20 binary features needs about a million entries, which is why later concepts factor it into smaller pieces.

Often you care about one variable, but the model is written for several. Marginalizing sums the joint over the variable you want to ignore. As the formula shows, $P(X = x)$ adds $P(X = x, Y = y)$ across every $y$, and for densities an integral over $y$ does the same job. This is the law of total probability from the first topic, written for random variables. It is the operation behind a mixture model's likelihood, where the hidden cluster is summed out, and behind a VAE's evidence $p(x)$, where the latent code is integrated out. The warning is that marginals lose information, because many different joint tables share the same marginals.

Observing $X = x$ picks out one slice of the joint table, the row where $X = x$. That slice does not sum to 1, so the formula divides it by the marginal $P(X = x)$ to renormalize. The result, $P(Y \mid X = x)$, is a full distribution over $Y$ for that fixed $x$. A classifier learns exactly this object, and regression learns its mean, $E[Y \mid X]$. Comparing conditionals across values of $x$ is a direct test of dependence: if every slice has the same shape, $Y$ is independent of $X$. The division fails when $P(X = x) = 0$, which is why unseen feature combinations need smoothing.

Any joint distribution can be written as a product of conditionals, one variable at a time. In the formula, each $X_i$ is conditioned on every earlier variable $X_1, \ldots, X_{i-1}$. This is an identity, not an assumption. It follows from applying the product rule $P(A, B) = P(A) P(B \mid A)$ repeatedly. Autoregressive language models use it directly, scoring a sentence as the product of next-token probabilities. Assumptions enter only when you drop conditioning variables, as a Bayesian network or a Markov chain does. Any variable ordering is valid, but some make the conditionals much easier to learn. For long products, multiply in log space to avoid underflow.

The multivariate Gaussian extends the bell curve to $d$ dimensions. In the formula, $\boldsymbol{\mu}$ is the mean vector, $\boldsymbol{\Sigma}$ is the $d \times d$ covariance matrix, and $|\boldsymbol{\Sigma}|$ is its determinant, which scales the normalizer. The exponent is a squared Mahalanobis distance, so contours of equal density are ellipses aligned with the eigenvectors of $\boldsymbol{\Sigma}$. Diagonal entries are variances, and off-diagonal entries are covariances. Marginals and conditionals are Gaussian again, which makes Gaussian processes and Kalman filters tractable. The model fails when data are multimodal or heavy-tailed, and a nearly singular $\boldsymbol{\Sigma}$ makes the inverse numerically unstable.

A point estimator turns a sample into a single guess. In the formula, $g$ is any function of the observations $X_1, \ldots, X_n$. The sample mean estimates $\mu$, the sample proportion estimates $p$, and the sample variance $S^2 = \frac{1}{n-1}\sum (X_i - \bar{X})^2$ estimates $\sigma^2$. Three properties judge an estimator. It is unbiased if it is right on average, consistent if it converges to $\theta$ as $n$ grows, and efficient if its variance is low. The $n - 1$ appears because deviations are measured from $\bar{X}$, which is fitted to the same data and so sits closer to the points than $\mu$ does.

An estimator can miss in two ways. Bias is systematic: $E[\hat{\theta}] - \theta$ is the gap between its average and the truth. Variance is noise, meaning how far $\hat{\theta}$ scatters around its own average across samples. The formula shows that mean squared error adds the squared bias and the variance, so the two can be traded. A slightly biased estimator with much smaller variance can win on MSE, and that is the case for regularization: ridge regression shrinks coefficients toward zero, which adds bias but cuts variance. The trap is judging estimators by unbiasedness alone, which can select a far noisier one.

A confidence interval turns an estimate into a range. In the formula, $\bar{X}$ is the sample mean, $z_{\alpha/2}$ is the normal critical value (1.96 for 95%), and $\sigma/\sqrt{n}$ is the standard error. With $\sigma$ unknown and small $n$, you use the sample $S$ and a $t$ critical value with $n - 1$ degrees of freedom. The 95% describes the procedure: across repeated samples, 95% of intervals built this way cover the true mean. It does not mean this particular interval contains $\mu$ with probability 0.95, which is a Bayesian credible-interval statement. The width scales as $1/\sqrt{n}$, so halving it takes four times the data.

For many statistics, such as a median, a ratio, or an F1 score, there is no convenient formula for the standard error. The bootstrap estimates it by simulation. Treat the observed sample as a stand-in for the population, draw $B$ resamples of size $n$ with replacement, and recompute the statistic on each. In the formula, $\hat{\theta}^*_b$ is the value from resample $b$. The spread of these replicates approximates the sampling distribution, and their 2.5th and 97.5th percentiles give a percentile interval. Efron introduced the method in 1979. It fails for extreme statistics such as the maximum, and for dependent data unless you resample whole blocks.

Maximum likelihood picks the parameter under which the observed data are most probable, $\arg\max_\theta P(X \mid \theta)$. MAP multiplies the likelihood by a prior $P(\theta)$ before maximizing, as the formula shows. Taking logs turns the prior into an additive penalty. A Gaussian prior on weights gives an L2 penalty, as in ridge regression, and a Laplace prior gives L1, as in the lasso. A uniform prior makes MAP identical to MLE. With little data the prior matters a great deal, and with lots of data the likelihood dominates and the two agree. The failure mode of pure MLE is overfitting small samples.

A test starts by naming what counts as no effect. The null hypothesis $H_0$ fixes the parameter at $\theta_0$, and the alternative $H_1$ describes the effect you want to detect. The formula shows the two-sided alternative $\theta \neq \theta_0$, which catches differences in either direction. A one-sided alternative such as $\theta < \theta_0$ looks in one direction only and has more power there. The null is what you try to reject, so failing to reject it is not evidence that it is true. Choose the direction before seeing the data, because picking the side after a peek halves the p-value dishonestly.

The p-value in the formula is the probability, computed as if $H_0$ were true, of a test statistic at least as extreme as the observed one. A small value means the data would be surprising under the null, so you reject $H_0$ when the p-value is below $\alpha$, conventionally 0.05. Three misreadings are common. The p-value is not the probability that $H_0$ is true, it is not the chance of an error, and it says nothing about effect size, since a huge sample makes a trivial difference significant. The American Statistical Association's 2016 statement lists these points. Report effect sizes and intervals alongside p-values.

Every test can err in two ways. A Type I error rejects a true null, a false positive, and happens with probability $\alpha$. A Type II error keeps a false null, a missed effect, and happens with probability $\beta$. Power, $1 - \beta$, is the probability of detecting a real effect. Lowering $\alpha$ reduces false positives but raises $\beta$ unless the sample grows. Power increases with sample size, effect size, and $\alpha$, and falls as noise grows. In model evaluation, a Type I error ships a model that is no better, while a Type II error discards a real improvement. Run the power analysis before the experiment, not after.

When $\sigma$ is unknown, you estimate it with the sample standard deviation $S$, which adds noise. The formula's statistic $t$ divides the difference $\bar{X} - \mu_0$ by $S/\sqrt{n}$ and follows a $t$-distribution with $n - 1$ degrees of freedom. That distribution has heavier tails than the normal, so small samples need larger statistics to reach significance, and it approaches the normal as $n$ grows. The one-sample form compares a mean with $\mu_0$, the two-sample form compares two groups, and the paired form tests differences within pairs. It assumes roughly normal, independent observations, and strong skew in tiny samples breaks it.

The chi-squared test compares observed counts $O_i$ with the counts $E_i$ expected under $H_0$. Each term in the formula squares the gap and divides by $E_i$, so a miss of 5 matters more when only 10 were expected than when 1,000 were. The statistic is never negative, and large values are evidence against the null. For goodness of fit across $k$ categories there are $k - 1$ degrees of freedom, and for independence in an $r \times c$ table there are $(r-1)(c-1)$. Use it to check whether class frequencies drifted or whether two categorical features are related. The approximation breaks when expected counts fall below about 5.

A/B testing applies the two-proportion z-test to live traffic. The formula divides the difference in success rates by a standard error built from the pooled rate $\hat{p}$ and the group sizes $n_A$ and $n_B$. The statistics are simple, and most mistakes happen in the experiment design. Checking results daily and stopping at the first significant value inflates the false-positive rate well above $\alpha$. Testing $m$ variants calls for a correction, such as the Bonferroni threshold $\alpha/m$. A significant 0.1% gain may not be worth deploying. Fix the metric, sample size, and stopping rule in advance, and report an interval for the lift.

Bayes' theorem in the formula combines two ingredients. The prior $P(\theta)$ is what you believe about the parameter before seeing data, and the likelihood $P(D \mid \theta)$ says how probable the data are under each parameter value. Their product, normalized by the evidence $P(D)$, is the posterior $P(\theta \mid D)$. The evidence is the law of total probability applied over $\theta$, and it is often the hard part to compute. Updates chain, so today's posterior becomes tomorrow's prior. Thompson sampling uses posteriors to decide which option to explore. The failure mode is a prior that rules out the truth, because zero prior probability stays zero whatever the data say.

A prior is conjugate to a likelihood when the posterior stays in the prior's family, so updating reduces to arithmetic on parameters. The formula shows the main example: a $\text{Beta}(\alpha, \beta)$ prior and $k$ successes in $n$ Binomial trials give $\text{Beta}(\alpha + k, \beta + n - k)$, so the prior parameters act as pseudo-counts. Other standard pairs are a Gamma prior with Poisson counts, which adds the total count to the shape and $n$ to the rate; a Normal prior with a Normal likelihood of known variance; and a Dirichlet prior with Multinomial counts, used in LDA topic models. Conjugacy is a convenience, not a truth, and it can misstate beliefs.

Often you want one number rather than a whole posterior. MAP estimation picks the posterior's mode, the parameter value with the highest posterior density. Because the evidence does not depend on $\theta$, the formula maximizes likelihood times prior, and taking logs turns the product into a sum. A Gaussian prior $\theta \sim N(0, \sigma^2)$ contributes $-\theta^2/(2\sigma^2)$, an L2 penalty, and a Laplace prior contributes an L1 penalty. A regularized network is therefore doing MAP estimation, and a tighter prior means stronger regularization. MAP discards uncertainty, and the mode changes under reparameterization, so it can land on an unrepresentative spike.

To predict a new observation, you could plug in one estimate of $\theta$, but that ignores how unsure you are about $\theta$. The posterior predictive in the formula averages the prediction $P(x_{\text{new}} \mid \theta)$ over every parameter value, weighted by the posterior $P(\theta \mid D)$. Its spread therefore includes both the noise in the data and the uncertainty in the parameter. It is wider than a plug-in prediction when data are scarce and matches it when data are plentiful. Bayesian optimization and active learning use this predictive uncertainty to choose where to sample next. Outside conjugate models the integral rarely has a closed form, so it is estimated by sampling.