4. Probability and Statistics

Probability

Probability spaces and random variables

Probability theory starts from a mathematical model of an experiment whose outcome is uncertain: tossing a coin, sampling users, or initializing a network with random weights. The model has three ingredients: a set of possible outcomes, a collection of events, and an assignment of probabilities to those events. This section motivates them from observed frequencies, states Kolmogorov’s axioms, and defines random variables as maps out of the resulting probability space. Most later calculations use these definitions only implicitly; Appendix A collects the measure-theoretic details behind them.

From relative frequencies to axioms

Repeat an experiment \(n\) times under the same conditions, and let \(N_n(A)\) count the repetitions in which an event \(A\) occurs. The relative frequency \(f_n(A)=N_n(A)/n\) has three properties at every \(n\). Writing \(\Omega\) for the event that some outcome occurs,

\[0\le f_n(A)\le1, \qquad f_n(\Omega)=1, \qquad f_n(A\cup B)=f_n(A)+f_n(B)\quad\text{when $A$ and $B$ cannot occur together.}\]

For many physical and computational experiments, \(f_n(A)\) also settles near a fixed number as \(n\) grows. The frequentist interpretation identifies the probability of \(A\) with this long-run value.

Relative frequencies of three disjoint events in \(5{,}000\) rolls of a fair die. At every \(n\) they add to one, and each settles near its probability: \(1/3\) for a roll of \(1\) or \(2\), \(1/6\) for a three, and \(1/2\) for \(4\), \(5\), or \(6\).

Defining probability as a limit of relative frequencies is hard to make precise. The limit need not exist for an arbitrary sequence of outcomes, and saying which infinite sequences count as random requires much of the theory being defined. Kolmogorov’s 1933 monograph Grundbegriffe der Wahrscheinlichkeitsrechnung took a different route. It keeps the properties above as axioms for an abstract assignment of probabilities, strengthens additivity to countably many events, and leaves open what probabilities mean. The frequency interpretation then returns as a theorem: under the axioms, the law of large numbers proves that relative frequencies over independent repetitions converge, with probability one, to the probability of the event. The same axioms serve Bayesian inference, where probabilities express degrees of belief rather than long-run frequencies.

Sample spaces and events

Definition (sample space and events). The sample space \(\Omega\) is the set of possible outcomes \(\omega\) of an experiment. An event is a subset \(A\subseteq\Omega\), and it occurs when the realized outcome belongs to \(A\). The events that receive probabilities form a collection \(\mathcal F\), called a σ-algebra, that contains \(\Omega\) and is closed under complements and countable unions.

Set operations express logical ones: \(A^c=\Omega\setminus A\) is “not \(A\)”, \(A\cap B\) is “\(A\) and \(B\)”, \(A\cup B\) is “\(A\) or \(B\)”, and \(A\subseteq B\) means that \(A\) implies \(B\). Events \(A\) and \(B\) are disjoint, or mutually exclusive, when \(A\cap B=\varnothing\). Closure under countable unions matters because statements about limits combine countably many events. For example, “the running average eventually stays within \(\varepsilon\) of \(\mu\)” is a countable union of countable intersections.

For a finite or countable sample space, \(\mathcal F\) contains every subset. For outcomes in \(\mathbb R^d\), it is the Borel σ-algebra, the smallest σ-algebra containing all intervals, or all boxes. It contains every open and closed set and every set obtained from these by countably many set operations. It cannot simply contain every subset, because no uniform law on an interval can assign a probability to all of them (Appendix A).

Probability measures

Definition (probability measure). A probability measure on \((\Omega,\mathcal F)\) is a function \(P:\mathcal F\to[0,1]\) that satisfies Kolmogorov’s axioms:

1. \(P(\Omega)=1\);
2. countable additivity: for pairwise disjoint events \(A_1,A_2,\ldots\in\mathcal F\),

\[P\Big(\bigcup_{i=1}^\infty A_i\Big)=\sum_{i=1}^\infty P(A_i).\]

The triple \((\Omega,\mathcal F,P)\) is a probability space.

Every other rule of probability follows from these axioms. Taking every \(A_i=\varnothing\) gives \(P(\varnothing)=0\), so countable additivity includes finite additivity. For events \(A\) and \(B\) it follows that

\[P(A^c)=1-P(A), \qquad A\subseteq B\ \Rightarrow\ P(A)\le P(B), \qquad P(A\cup B)=P(A)+P(B)-P(A\cap B).\]

The last identity splits \(A\cup B\) into the disjoint pieces \(A\setminus B\), \(A\cap B\), and \(B\setminus A\). The union bound states that \(P(\bigcup_i A_i)\le\sum_iP(A_i)\) for any countable collection of events, without an independence assumption. If each of \(m\) possible failures has probability at most \(\alpha/m\), the probability of any failure is at most \(\alpha\). Countable additivity also makes probability continuous along increasing or decreasing sequences of events, which is what lets probabilities pass to limits (Appendix A).

Two constructions cover most models in these notes. On a finite or countable \(\Omega\), nonnegative weights \(p(\omega)\) that sum to one define \(P(A)=\sum_{\omega\in A}p(\omega)\); equal weights on a finite \(\Omega\) give the classical rule of the next subsection. On \(\mathbb R^d\), a nonnegative function \(p\) with integral one, a density, defines \(P(A)=\int_Ap(x)\,dx\). The uniform law on \([0,1]\), which gives each interval its length, is the simplest example.

Under the uniform law, each single outcome has probability zero, although every outcome is possible. An event of probability zero is a null event; it need not be empty. A statement that holds outside a null event holds almost surely (a.s.), or with probability one.

Equally likely outcomes and counting

When \(\Omega\) is finite and every outcome is equally likely, probabilities are ratios of counts:

\[P(A)=\frac{|A|}{|\Omega|}.\]

The assumption of equal weights carries the whole argument. When the ordered outcomes of an experiment are equally likely, its unordered selections with replacement and its occupancy patterns generally are not, so the sample space must be chosen before counting. Choosing \(k\) items from \(n\) distinct items gives four standard counts:

The multinomial coefficient \(n!/(n_1!\cdots n_m!)\) counts the ways to divide \(n\) distinct items into labeled groups of sizes \(n_1,\ldots,n_m\).

Drawing two cards from a standard deck shows why the replacement convention matters. With replacement, \(P(\text{two aces})=(4/52)^2\approx0.0059\). Without replacement, the product rule of the next section gives \((4/52)(3/51)\approx0.0045\): the first draw changes the deck faced by the second. For more than two events, inclusion–exclusion extends the formula for \(P(A\cup B)\) by alternately adding and subtracting intersections; Appendix B states it and uses it to count permutations with no fixed point.

Random variables and their laws

An outcome can be complicated, such as a sequence of tosses, an image, or a whole sampled dataset, and most questions concern numerical summaries of it. A random variable is such a summary, with one requirement: every question about its value must be an event, so that it has a probability.

Definition (random variable). Let \((\Omega,\mathcal F,P)\) be a probability space, and let \(E\) be a set of values with a σ-algebra \(\mathcal E\), such as \(\mathbb R^d\) with its Borel sets. A random variable with values in \(E\) is a function \(X:\Omega\to E\) that is measurable:

\[X^{-1}(B)=\{\omega\in\Omega:X(\omega)\in B\}\in\mathcal F \qquad\text{for every }B\in\mathcal E.\]

It is a real random variable when \(E=\mathbb R\), and a random vector when \(E=\mathbb R^d\).

A random variable is an ordinary function; the randomness lies in which outcome \(\omega\) occurs. The event \(X^{-1}(B)\) is written \(\{X\in B\}\), and its probability \(P(X\in B)\). Measurability guarantees that these probabilities are defined, and it holds for every function met in practice: sums, products, continuous functions, and limits of random variables are again random variables (Appendix A).

Definition (law). The distribution, or law, of \(X\) is the probability measure \(P_X\) on \((E,\mathcal E)\) given by

\[P_X(B)=P(X\in B)=P\big(X^{-1}(B)\big), \qquad B\in\mathcal E.\]

The law is the pushforward of \(P\) through \(X\), and it makes \((E,\mathcal E,P_X)\) a probability space of its own. Most calculations need only the law, so the underlying \(\Omega\) usually stays in the background. Variables with the same law are equal in distribution, written \(X\overset{d}{=}Y\), even if they are defined on different probability spaces or never take equal values. For one toss of a fair coin, the indicators of heads and of tails are equal in distribution, but they never agree.

For two tosses of a fair coin, \(\Omega=\{HH,HT,TH,TT\}\), \(\mathcal F\) contains all subsets, and \(P\) gives each outcome probability \(1/4\). The number of heads \(X\) maps \(\Omega\) to \(\{0,1,2\}\). The preimages of the three values are \(\{TT\}\), \(\{HT,TH\}\), and \(\{HH\}\), so the law of \(X\) assigns them \(1/4\), \(1/2\), and \(1/4\). Several outcomes can give the same value. Counting the possible values of \(X\) and treating them as equally likely would give a different model: probability weights come from assumptions about the experiment, rather than from the names of its outcomes.

The number of heads in two fair tosses, as a map from the sample space to the real line. The event \(\{X=1\}\) is the preimage \(\{HT,TH\}\), and the law of \(X\) collects the probability of each preimage.

Capital letters denote random quantities; lowercase letters denote their realized values. A sample is \(D=(Z_1,\ldots,Z_n)\), and its observed value is \((z_1,\ldots,z_n)\). The underlying data-generating law will be written \(P_\ast\); a statistical model proposes a family of laws \(\{P_\theta:\theta\in\Theta\}\).

An observation may itself be a pair \(Z=(X,Y)\) containing an input vector and a target. Its joint law describes both how inputs occur and how targets vary with those inputs. A conditional model for \(Y\mid X\) describes only the second component. Specifying a distribution is a mathematical assumption; deciding whether it describes the observed data adequately is part of statistical modeling.

Mass, density, and cumulative probability

Definition (mass function and density). A real random variable \(Z\) is discrete if \(P(Z\in S)=1\) for some countable set \(S\). Its probability mass function is \(p_Z(z)=P(Z=z)\), with \(\sum_{z\in S}p_Z(z)=1\). It is absolutely continuous if there is a probability density function \(p_Z\ge0\) such that

\[P(Z\in B)=\int_B p_Z(z)\,dz \qquad\text{for every }B\in\mathcal B(\mathbb R),\]

and in particular \(\int_{\mathbb R}p_Z(z)\,dz=1\).

The notation \(p\) will denote either a mass function or a density; the support and sum or integral identify which one is meant. A density measures probability per unit of the variable. Its height can exceed one: a uniform variable on \((0,1/2)\) has density \(2\), and \(P(0<Z\le0.2)=2(0.2)=0.4\). An individual point has probability zero under an absolutely continuous law, even where its density is positive. Some laws combine point masses and a continuous component, so a single ordinary density need not describe every distribution.

Definition (distribution function). The cumulative distribution function (CDF) of a real random variable \(Z\) is

\[F_Z(z)=P(Z\le z),\qquad z\in\mathbb R.\]

It is nondecreasing and right-continuous, with limits zero and one at the two ends of the real line, and it determines the law \(P_Z\) completely (Appendix A). Conversely, every nondecreasing, right-continuous function with these limits is the CDF of some random variable; the inverse-CDF construction later in the chapter builds one from a uniform variable. For a density, \(F_Z(z)=\int_{-\infty}^z p_Z(u)\,du\) and \(F_Z'(z)=p_Z(z)\) almost everywhere. CDF differences give interval probabilities:

\[P(a<Z\le b)=F_Z(b)-F_Z(a),\qquad a<b.\]

Define the left limit \(F_Z(z^-)=\lim_{u\uparrow z}F_Z(u)=P(Z<z)\). The jump at \(z\) is therefore

\[P(Z=z)=F_Z(z)-F_Z(z^-).\]

These identities apply to continuous, discrete, and mixed laws. For the uniform example, integrating the density gives

\[F_Z(z)= \begin{cases} 0,&z\le0,\\ 2z,&0<z<1/2,\\ 1,&z\ge1/2. \end{cases}\]

The shaded area under the Unif\((0,1/2)\) density, \(2\times0.2=0.4\), is the height of the CDF at \(0.2\).

The generalized quantile is

\[F_Z^{-1}(q)=\inf\{z:F_Z(z)\ge q\},\qquad 0<q<1.\]

The median is a \(1/2\) quantile. The survival function \(1-F_Z(z)\) gives \(P(Z>z)\) and is useful for tail probabilities.

Blitzstein and Hwang’s Introduction to Probability gives a systematic treatment of these objects, conditioning, and distribution families.

Conditioning and dependence

Definition (conditional probability). For an event \(B\) with \(P(B)>0\), the conditional probability of \(A\) given \(B\) is

\[P(A\mid B)=\frac{P(A\cap B)}{P(B)}.\]

Conditioning restricts attention to outcomes consistent with the observation and renormalizes their weights. For fixed \(B\), the map \(A\mapsto P(A\mid B)\) is again a probability measure on \((\Omega,\mathcal F)\), so every rule of the previous section applies to conditional probabilities. Rearranging gives the product rule, \(P(A\cap B)=P(A\mid B)P(B)\). If positive-probability events \(B_1,\ldots,B_k\) partition the sample space up to a null set, summing over the disjoint cases gives total probability. When \(P(A)>0\), it also gives Bayes’ rule. Zero-probability partition cells can be omitted:

\[P(A)=\sum_j P(A\mid B_j)P(B_j),\qquad P(B_j\mid A)=\frac{P(A\mid B_j)P(B_j)}{\sum_iP(A\mid B_i)P(B_i)}.\]

For example, suppose one percent of messages are spam. A filter flags \(90\%\) of spam and \(5\%\) of ordinary messages. Among \(10{,}000\) messages, the expected flagged counts are \(90\) spam and \(495\) ordinary messages. Therefore

\[P(\text{spam}\mid\text{flagged}) =\frac{0.90(0.01)}{0.90(0.01)+0.05(0.99)} \approx0.154.\]

The prevalence, or base rate, matters alongside the filter’s behavior within each class. Reversing the conditioning event changes the question.

Of 10,000 messages, 585 are flagged, and only 90 of those are spam.

Bayes’ rule is often clearest in odds form. For two hypotheses \(H_1,H_2\) and data \(D\),

\[\frac{P(H_1\mid D)}{P(H_2\mid D)} =\frac{P(H_1)}{P(H_2)}\cdot\frac{P(D\mid H_1)}{P(D\mid H_2)}.\]

Posterior odds are prior odds multiplied by the likelihood ratio. In the filter example, the prior odds of spam are \(1/99\), and a flag multiplies them by \(0.90/0.05=18\). The posterior odds \(18/99\) correspond to the probability \(18/117\approx0.154\) found above. The likelihood ratio measures what the observation says; the prior odds say how rare the hypothesis was to begin with.

Joint laws, marginalization, and the chain rule

A joint distribution describes variables together. If \((X,Y)\) has a joint mass function, then \(p_X(x)=\sum_y p_{X,Y}(x,y)\); for a joint density, replace the sum with an integral. This operation is marginalization. A conditional density is

\[p_{Y\mid X}(y\mid x)=\frac{p_{X,Y}(x,y)}{p_X(x)},\qquad p_X(x)>0.\]

It remains meaningful for continuous \(X\), even though \(P(X=x)=0\): it describes the conditional law through densities, rather than dividing two point probabilities. Conditional laws are only determined up to changes on sets of conditioning values having probability zero.

Repeated factorization gives the probability chain rule:

\[p(z_1,\ldots,z_n)=p(z_1)\prod_{i=2}^n p(z_i\mid z_1,\ldots,z_{i-1}).\]

This identity requires no independence. An autoregressive language model uses this factorization to describe a sequence by distributions for successive tokens.

For a small discrete example, suppose two binary variables have joint probabilities

Summing rows gives \(P(X=0)=P(X=1)=0.5\); summing the second column gives \(P(Y=1)=0.4\). Within the \(X=1\) row, renormalization gives \(P(Y=1\mid X=1)=0.3/0.5=0.6\). The variables are dependent because knowing \(X=1\) changes the probability of \(Y=1\). Conversely, total probability recovers \(0.4\) as \(0.2(0.5)+0.6(0.5)\). Marginalization and conditioning are distinct operations on the same joint law.

Independence and exchangeability

Definition (independence). Events \(A\) and \(B\) are independent if \(P(A\cap B)=P(A)P(B)\). Events \(A_1,\ldots,A_n\) are mutually independent if

\[P\Big(\bigcap_{i\in S}A_i\Big)=\prod_{i\in S}P(A_i) \qquad\text{for every subset }S\subseteq\{1,\ldots,n\}.\]

Random variables \(X_1,\ldots,X_n\) are independent if \(P(X_1\in B_1,\ldots,X_n\in B_n)=\prod_iP(X_i\in B_i)\) for all measurable sets \(B_1,\ldots,B_n\): any events determined by different variables are independent, and the joint law is the product of the marginal laws. An infinite family is independent when every finite subfamily is. Independence of two variables is written \(X\perp Y\).

When \(P(B)>0\), independence of \(A\) and \(B\) says that \(P(A\mid B)=P(A)\): learning that \(B\) occurred leaves the probability of \(A\) unchanged. For variables with a joint density or mass function, independence is the factorization \(p_{X,Y}(x,y)=p_X(x)p_Y(y)\) almost everywhere, and for real variables it is equivalent to factorization of the joint CDF. Independence is different from disjointness: two disjoint events of positive probability are dependent, since observing one rules out the other.

The variables \(Z_1,\ldots,Z_n\) are independent and identically distributed, abbreviated iid, when their joint law is the product of one common marginal law. Pairwise independence alone does not imply joint independence. For independent fair bits \(U,V\), the triple \((U,V,U\mathbin{\oplus}V)\) is pairwise independent, but the third bit is determined by the first two.

The two parts of iid make separate claims. Variables can be independent with different distributions, as when measurements have different known noise variances. They can also have identical marginal distributions while being dependent, as with repeated measurements of the same individual. Sampling without replacement from a finite population creates dependence even when every sampled position has the same marginal law. Independence is often a useful approximation; whether it is appropriate depends on how data were collected.

A sequence is exchangeable when its joint law is unchanged by any permutation of the indices. Iid sequences are exchangeable, and so are draws without replacement from a finite population, even though those draws are dependent. Permutation tests later in the chapter rest on exactly this symmetry.

Conditional independence, \(X\perp Y\mid C\), means that their conditional joint law factors given \(C\). For example, repeated observations may be independent after conditioning on a shared parameter but dependent when that parameter is integrated out. A \(\operatorname{Bernoulli}(q)\) variable equals one with probability \(q\) and zero with probability \(1-q\). If \(Q\) is random and \(X,Y\mid Q\) are independent \(\operatorname{Bernoulli}(Q)\) variables, observing \(X=1\) favors larger values of \(Q\), which in turn favor \(Y=1\). A common random success probability can therefore create marginal dependence. Conditional independence is a property of the specified joint model; it does not, by itself, establish a causal relationship.

Simpson’s paradox

Conditioning on a third variable can reverse a comparison. The table gives success counts for two kidney-stone treatments, split by stone size (Charig et al., 1986):

Treatment A has the higher success rate for small stones and for large stones, yet the lower rate overall. There is no arithmetic error. Writing \(S\) for success, \(T\) for treatment, and \(G\) for stone size, total probability expresses each overall rate as a weighted average of the group rates:

\[P(S\mid T=t)=\sum_g P(S\mid T=t,G=g)\,P(G=g\mid T=t).\]

The weights differ between treatments: \(263\) of the \(350\) patients given A had large stones, the harder cases, against \(80\) of the \(350\) given B. This reversal is Simpson’s paradox.

Treatment A succeeds more often within each stone size, but three quarters of its patients had large stones.

Which comparison answers a question depends on how treatments were assigned and on what the question is. Here stone size influenced the choice of treatment and also the chance of success, so the within-size comparisons are the relevant ones for choosing a treatment. The probabilities alone do not settle whether a variable should be conditioned on; that judgment requires assumptions about how the data arose.

Expectations, variation, and conditional averages

Definition (expectation). The expectation of a real random variable \(Z\) is its probability-weighted average, the integral of \(Z\) over the sample space:

\[\mathbb EZ=\int_\Omega Z(\omega)\,P(d\omega).\]

The integral is the Lebesgue integral of \(Z\) with respect to \(P\). It is always defined for \(Z\ge0\), possibly as \(+\infty\), and a general \(Z\) is integrable when \(\mathbb E|Z|<\infty\), in which case its expectation is finite. Appendix A defines this integral, compares it with the Riemann integral, and proves that it equals an integral over the values of \(Z\), as in the next paragraph.

The integral over \(\Omega\) can be moved to the law of \(Z\): for measurable \(g\) with \(\mathbb E|g(Z)|<\infty\), \(\mathbb E[g(Z)]=\int_{\mathbb R}g(z)\,P_Z(dz)\), an integral over the real line (Appendix A). This law of the unconscious statistician becomes a sum or an ordinary integral for a mass function or a density:

\[\mathbb E[g(Z)]=\sum_z g(z)p_Z(z) \quad\text{or}\quad \mathbb E[g(Z)]=\int g(z)p_Z(z)\,dz.\]

For a law with neither a mass function nor a density, the same integral is taken against the CDF and written \(\int g(z)\,dF_Z(z)\). In every case there is no need to first find the distribution of \(g(Z)\). Expectations are linear: \(\mathbb E[\sum_i a_iZ_i]=\sum_i a_i\mathbb E Z_i\) whenever each \(Z_i\) is integrable. Independence is unnecessary. In particular, an indicator \(\mathbf1_A\) equals one on \(A\) and zero elsewhere, so \(\mathbb E\mathbf1_A=P(A)\). Expected counts are therefore sums of event probabilities, including when the counted events are dependent.

For a uniformly random ordering of \(n\) items, let \(I_i\) indicate that item \(i\) stays in its original position. The number of fixed points \(\sum_iI_i\) has expectation \(n\cdot(1/n)=1\) for every \(n\), although the indicators are dependent. Indicators also give the tail-sum formula. A variable with values in \(\{0,1,2,\ldots\}\) satisfies \(X=\sum_{k\ge1}\mathbf1\{X\ge k\}\), so

\[\mathbb EX=\sum_{k\ge1}P(X\ge k), \qquad\text{and}\qquad \mathbb EX=\int_0^\infty P(X>t)\,dt\]

for any nonnegative variable.

For \(\mu=\mathbb E Z\) and finite second moment,

\[\operatorname{Var}(Z)=\mathbb E[(Z-\mu)^2] =\mathbb E[Z^2]-\mu^2, \qquad \operatorname{SD}(Z)=\sqrt{\operatorname{Var}(Z)}.\]

Variance has squared units; standard deviation has the same units as \(Z\). For finite second moments,

\[\operatorname{Cov}(X,Y)=\mathbb E[(X-\mathbb EX)(Y-\mathbb EY)],\]

and correlation divides this by the two standard deviations, when both are positive. Covariance measures linear association. Independence implies zero covariance, but the converse fails: for \(X\sim\operatorname{Unif}(-1,1)\) and \(Y=X^2\), symmetry gives \(\operatorname{Cov}(X,Y)=0\) although \(Y\) is determined by \(X\).

Expanding a squared sum gives

\[\operatorname{Var}\!\left(\sum_i a_iZ_i\right) =\sum_i a_i^2\operatorname{Var}(Z_i) +2\sum_{i<j}a_i a_j\operatorname{Cov}(Z_i,Z_j).\]

Thus averaging iid observations with variance \(\sigma^2\) gives variance \(\sigma^2/n\). Correlation changes this reduction; duplicating the same observation \(n\) times does not supply \(n\) independent pieces of information.

For instance, if each observation has variance \(\sigma^2\) and every distinct pair has correlation \(\rho\), then

\[\operatorname{Var}(\bar Z_n)=\frac{\sigma^2}{n}\bigl(1+(n-1)\rho\bigr).\]

For fixed positive \(\rho\), this variance approaches \(\rho\sigma^2\) rather than zero. This calculation is one reason to distinguish the number of recorded rows from the amount of independent information. It also shows why an iid standard-error formula can be misleading for clustered or temporally correlated observations.

For a column vector \(Z\in\mathbb R^d\), define

\[\mu=\mathbb EZ,\qquad \Sigma=\mathbb E[(Z-\mu)(Z-\mu)^\top].\]

The covariance matrix is symmetric positive semidefinite because \(a^\top\Sigma a=\operatorname{Var}(a^\top Z)\ge0\). An affine transformation satisfies \(\mathbb E[AZ+b]=A\mu+b\) and \(\operatorname{Cov}(AZ+b)=A\Sigma A^\top\).

Conditional expectation and total variance

Definition (conditional expectation). For an integrable \(Y\), let \(m(x)=\mathbb E[Y\mid X=x]\) be the mean of the conditional law of \(Y\) given \(X=x\). With a conditional density,

\[m(x)=\int y\,p_{Y\mid X}(y\mid x)\,dy,\]

and with a conditional mass function the integral becomes a sum. The conditional expectation \(\mathbb E[Y\mid X]\) is the random variable \(m(X)\). It is a function of \(X\) rather than a number: its value depends on the observed \(X\).

The tower property, or law of total expectation, says

\[\mathbb E[\mathbb E(Y\mid X)]=\mathbb EY.\]

It follows by integrating first over \(y\) and then over \(x\). Appendix A gives the general definition, which needs no density and can condition on any body of information, such as a whole past. It characterizes \(m(X)\) as the function of \(X\) whose average over every event determined by \(X\) equals that of \(Y\).

Assume now that \(\mathbb E[Y^2]<\infty\). Write \(Y= m(X)+R\), where \(R=Y-m(X)\) satisfies \(\mathbb E[R\mid X]=0\). Expanding \((Y-\mathbb EY)^2=(m(X)-\mathbb EY+R)^2\) and taking expectations eliminates the cross term, yielding the law of total variance:

\[\boxed{\operatorname{Var}(Y) =\mathbb E[\operatorname{Var}(Y\mid X)] +\operatorname{Var}(\mathbb E[Y\mid X]).}\]

The first term measures average variation within groups defined by \(X\); the second measures variation between their means. The same residual argument shows that, for any square-integrable prediction rule \(h(X)\),

\[\mathbb E[(Y-h(X))^2] =\mathbb E[\operatorname{Var}(Y\mid X)] +\mathbb E[(m(X)-h(X))^2].\]

More explicitly, \(Y-h(X)=R+[m(X)-h(X)]\). After squaring, the cross term has expectation

\[\mathbb E\{R[m(X)-h(X)]\} =\mathbb E\{[m(X)-h(X)]\mathbb E[R\mid X]\}=0.\]

The factor depending only on \(X\) can be taken outside the inner conditional expectation. This orthogonality of the residual to functions of \(X\) is stronger than having zero ordinary covariance with \(X\). A population linear predictor that minimizes expected squared error has residuals orthogonal to its chosen features in expectation. Ordinary least squares on one sample instead gives empirical orthogonality to the columns of its design matrix. The true conditional mean removes every square-integrable predictable component.

Consequently, the conditional mean minimizes expected squared prediction error. This is a population statement: estimating that mean from finite data remains a statistical problem.

There is also a law of total covariance:

\[\operatorname{Cov}(U,V) =\mathbb E[\operatorname{Cov}(U,V\mid C)] +\operatorname{Cov}(\mathbb E[U\mid C],\mathbb E[V\mid C]).\]

It follows by subtracting the conditional means and expanding, just as for total variance. Returning to the shared random probability \(Q\), conditionally independent Bernoulli draws \(U,V\mid Q\) have conditional covariance zero and conditional means both equal to \(Q\). Consequently, \(\operatorname{Cov}(U,V)=\operatorname{Var}(Q)\). Variation in a shared hidden variable explains their positive marginal association.

For a concrete mixture, let \(C\) indicate group membership, with \(P(C=1)=0.3\) and \(P(C=0)=0.7\). Let \(G\) be an independent standard Normal variable, which has mean zero and variance one, and set \(Y=2C+G\). Conditional on \(C\), the mean of \(Y\) is \(2C\) and its variance is one. Thus \(\mathbb EY=0.6\) and \(\operatorname{Var}(Y)=1+4(0.3)(0.7)=1.84\). Simulation reproduces both sources of variation:

import numpy as np

rng = np.random.default_rng(7)
c = rng.binomial(1, 0.3, size=200_000)
conditional_mean = 2.0 * c
y = conditional_mean + rng.normal(size=c.size)
print("mean:", y.mean(), "theory:", 0.6)
print("variance:", y.var(), "theory:", 1.84)
print("within + between:", 1.0 + conditional_mean.var())

The mixture \(Y=2C+G\): its variance \(1.84\) is the within-group variance \(1\) plus the variance \(0.84\) of the group means \(0\) and \(2\).

Two moment inequalities

Two useful moment inequalities are Cauchy–Schwarz, \(|\mathbb E[UV]|\le\sqrt{\mathbb E[U^2]\mathbb E[V^2]}\), and Jensen’s inequality, \(\phi(\mathbb EZ)\le\mathbb E[\phi(Z)]\) for convex \(\phi\) when the relevant expectations exist. Jensen explains why applying a nonlinear function before averaging generally changes the result; for example, \(\mathbb E[Z^2]\ge(\mathbb EZ)^2\).

For differentiable convex \(\phi\), Jensen follows directly from the supporting-line inequality \(\phi(z)\ge\phi(\mu)+\phi'(\mu)(z-\mu)\) with \(\mu=\mathbb EZ\): expectation removes the linear term. The direction reverses for concave functions. In particular, averaging positive quantities and then taking a logarithm is at least as large as averaging their logarithms, provided the relevant expectations exist. Expectation preserves linear operations, while convexity controls certain nonlinear ones.

Distribution families

A distribution family specifies its support, parameters, and probability rule. Write \(\operatorname{Gamma}(\alpha,\lambda)\) with shape \(\alpha>0\) first and rate \(\lambda>0\) second, and \(\operatorname{Beta}(\alpha,\beta)\) with its two positive shape parameters. The gamma convention is shape–rate throughout these notes; another common convention uses shape and scale, where scale is \(1/\lambda\). For a normal law, the second argument of \(\mathcal N(\mu,\sigma^2)\) is the variance, not the standard deviation. Geometric and negative binomial variables count the failures before the first or the \(r\)th success; some sources count trials instead, which adds \(1\) or \(r\).

Discrete observations

A binomial variable counts successes in \(n\) independent Bernoulli trials. The categorical distribution describes a class or token, while its one-hot representation is a random vector. Counts of the \(K\) outcomes across \(n\) iid categorical draws have a multinomial distribution:

\[P(N_1=n_1,\ldots,N_K=n_K) =\frac{n!}{\prod_k n_k!}\prod_kq_k^{n_k}, \qquad \sum_k n_k=n.\]

The binomial formula follows by counting which \(k\) of the \(n\) trials succeed. Each such sequence has probability \(q^k(1-q)^{n-k}\), and there are \(\binom nk\) such sequences. Its moments follow more simply by writing the count as a sum of indicators: linearity gives mean \(nq\), and independence removes covariance terms to give variance \(nq(1-q)\). The same random variable can therefore be analyzed through its probability mass function or through its construction, whichever makes a calculation clearer.

A geometric variable counts the failures before the first success in independent trials with success probability \(q\). Since \(P(X\ge k)=(1-q)^k\), it is memoryless:

\[P(X\ge a+b\mid X\ge a)=P(X\ge b).\]

Failures already observed do not change the distribution of the remaining wait. The number of failures before the \(r\)th success is a sum of \(r\) independent geometric counts and has the negative binomial distribution.

Sampling without replacement changes the binomial. Draw \(n\) items from a population of \(N\) items, \(K\) of which are successes. The number of successes drawn is hypergeometric. Its mean \(nK/N\) matches the binomial mean with \(q=K/N\), but its variance carries the finite-population correction \((N-n)/(N-1)\). The draws are negatively dependent: each success drawn leaves fewer behind. The correction is negligible when \(n\) is a small fraction of \(N\).

The Poisson family is the law of rare events. If \(X_n\sim\operatorname{Binomial}(n,\lambda/n)\) with \(\lambda\) fixed, then for each \(k\),

\[\binom nk\Big(\frac\lambda n\Big)^k\Big(1-\frac\lambda n\Big)^{n-k}\longrightarrow e^{-\lambda}\frac{\lambda^k}{k!}.\]

Many independent opportunities, each unlikely, produce a moderate total count. Independent Poisson counts add their rates. Conversely, if \(X\sim\operatorname{Poisson}(\lambda)\) and \(Y\sim\operatorname{Poisson}(\mu)\) are independent, then \(X\mid X+Y=m\sim\operatorname{Binomial}(m,\lambda/(\lambda+\mu))\). Equality of mean and variance is a modeling restriction, not a property of counts in general; counts whose variance exceeds their mean are called overdispersed.

Binomial counts with mean \(4\) approach the Poisson\((4)\) law as the number of trials grows (left). Ten draws without replacement from twenty items, eight of them successes, vary less than ten draws with replacement (right).

Continuous observations and uncertain probabilities

The normalizing functions are

\[\Gamma(\alpha)=\int_0^\infty t^{\alpha-1}e^{-t}\,dt, \qquad B(\alpha,\beta)=\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}.\]

The gamma function \(\Gamma(\alpha)\) is distinct from the distribution name \(\operatorname{Gamma}(\alpha,\lambda)\). Setting the shape to one gives \(\operatorname{Gamma}(1,\lambda)=\operatorname{Exp}(\lambda)\).

The SciPy gamma distribution uses shape and scale: the law written here as \(\operatorname{Gamma}(\alpha,\lambda)\) is gamma(a=alpha, scale=1/rate) in SciPy. NumPy draws it with rng.gamma(shape=alpha, scale=1/rate, size=...). For example, shape \(3\) and rate \(2\) give mean \(3/2\) and variance \(3/4\):

import numpy as np
from scipy.stats import gamma

alpha, rate = 3.0, 2.0
law = gamma(a=alpha, scale=1.0 / rate)
print("mean and variance:", law.mean(), law.var())  # 1.5, 0.75
assert np.allclose([law.mean(), law.var()], [alpha/rate, alpha/rate**2])

The exponential distribution is the continuous memoryless law:

\[P(Z>s+t\mid Z>s)=\frac{e^{-\lambda(s+t)}}{e^{-\lambda s}}=P(Z>t),\]

and among continuous distributions on \([0,\infty)\) it is the only one with this property. It is the limit of a geometric waiting time when time is divided into short intervals, each containing an event with probability about \(\lambda\) times its length. For independent \(Z_i\sim\operatorname{Exp}(\lambda_i)\), the minimum satisfies \(P(\min_iZ_i>t)=\prod_ie^{-\lambda_it}\), so it is \(\operatorname{Exp}(\sum_i\lambda_i)\), and \(Z_j\) is the smallest with probability \(\lambda_j/\sum_i\lambda_i\). Sums of independent gamma variables with a common rate add their shapes: \(\operatorname{Gamma}(\alpha_1,\lambda)+\operatorname{Gamma}(\alpha_2,\lambda)\) is \(\operatorname{Gamma}(\alpha_1+\alpha_2,\lambda)\). In particular, a sum of \(k\) independent \(\operatorname{Exp}(\lambda)\) variables is \(\operatorname{Gamma}(k,\lambda)\).

These facts fit together in the Poisson process, a model for events occurring at rate \(\lambda\) in continuous time. The gaps between successive events are iid \(\operatorname{Exp}(\lambda)\), so the time of the \(k\)th event is \(\operatorname{Gamma}(k,\lambda)\). The number of events in any interval of length \(t\) is \(\operatorname{Poisson}(\lambda t)\), and counts in disjoint intervals are independent. The waiting-time and counting descriptions are two views of one model: the first event comes after time \(t\) exactly when the interval \([0,t]\) contains no event, and both probabilities equal \(e^{-\lambda t}\). Appendix D adds renewal processes, Markov chains, and martingales, the other basic models of randomness that unfolds over time.

One run of a rate-\(1\) Poisson process. The gaps between events are exponential, the count \(N(t)\) rises by one at each event, and the counts in the four windows of length \(2\) are independent Poisson\((2)\) draws.

A beta distribution models a probability or fraction on \((0,1)\). It also arises from ratios: for independent \(X\sim\operatorname{Gamma}(\alpha,\lambda)\) and \(Y\sim\operatorname{Gamma}(\beta,\lambda)\), the fraction \(X/(X+Y)\) is \(\operatorname{Beta}(\alpha,\beta)\) and is independent of the total \(X+Y\) (Appendix B).

The Dirichlet distribution generalizes beta to a random probability vector \(Q\sim\operatorname{Dirichlet}(\alpha_1,\ldots,\alpha_K)\). For \(q_k>0\), \(\sum_kq_k=1\), and \(\alpha_k>0\), its density in the first \(K-1\) coordinates is

\[p_Q(q_1,\ldots,q_{K-1}) =\frac{\Gamma(\alpha_0)}{\prod_k\Gamma(\alpha_k)}\prod_kq_k^{\alpha_k-1}, \quad q_K=1-\sum_{k<K}q_k,\quad \alpha_0=\sum_k\alpha_k.\]

Its mean is \(\mathbb E Q_k=\alpha_k/\alpha_0\). The total \(\alpha_0\) controls concentration around the mean when the ratios \(\alpha_k/\alpha_0\) are held fixed. Normalizing independent \(\operatorname{Gamma}(\alpha_k,1)\) variables by their sum gives a \(\operatorname{Dirichlet}(\alpha_1,\ldots,\alpha_K)\) vector; with every \(\alpha_k=1\), it is uniform on the simplex.

Unlike a Normal law, whose shape is fixed and only shifts and stretches with its parameters, gamma and beta distributions change shape. A gamma density with shape \(\alpha\le1\) is largest at zero, while larger shapes move its mode away from zero and make it more symmetric. For the beta family, \(\alpha=\beta=1\) gives a uniform density, \(\alpha=\beta>1\) concentrates toward the middle, \(\alpha=\beta<1\) places more density near the endpoints, and unequal parameters tilt it toward one end.

Gamma densities with rate \(1\) (left) and beta densities (right). The shape parameters change the form of each density, not only its location and scale.

Distributions derived from Normal samples

Several distributions arise from Normal samples and reappear throughout inference. If \(G_1,\ldots,G_\nu\) are iid standard normals, then \(V=\sum_iG_i^2\sim\chi^2_\nu=\operatorname{Gamma}(\nu/2,1/2)\), with mean \(\nu\) and variance \(2\nu\). If \(G\sim\mathcal N(0,1)\) independently of \(V\), then \(G/\sqrt{V/\nu}\) has Student’s \(t\) distribution with \(\nu\) degrees of freedom. It has heavier tails than a normal; its mean is zero for \(\nu>1\) and its variance is \(\nu/(\nu-2)\) for \(\nu>2\). As \(\nu\to\infty\), it approaches \(\mathcal N(0,1)\). The \(t_1\) law is the Cauchy distribution, with density \(1/[\pi(1+x^2)]\); it is also the law of a ratio of two independent standard normals, and it has neither a finite mean nor a finite variance. Moment assumptions in limit theorems therefore matter.

For independent \(U\sim\chi^2_m\) and \(V\sim\chi^2_k\), the ratio \((U/m)/(V/k)\) has the \(F\) distribution \(F_{m,k}\), used in variance-ratio and regression tests. If \(G\sim\mathcal N(\mu,\sigma^2)\), then \(L=e^G\) is lognormal, with mean \(e^{\mu+\sigma^2/2}\). A lognormal variable has finite moments of every order, yet its right tail is heavy enough that \(\mathbb E[e^{sL}]\) is infinite for every \(s>0\).

Student \(t\) densities approach the standard Normal as the degrees of freedom grow (left). On a logarithmic scale (right), their tail probabilities stay far above the Normal tail.

For \(Z\sim\mathcal N(\mu,\sigma^2)\), standardizing gives \((Z-\mu)/\sigma\sim\mathcal N(0,1)\). Thus \(P(Z\le z)=\Phi((z-\mu)/\sigma)\), where \(\Phi\) is the standard normal CDF: all one-dimensional normal probabilities reduce to one reference CDF. Write \(\Phi^{-1}(q)\) for its \(q\) quantile.

Moment generating functions

The moment generating function (MGF) of \(Z\) is \(M_Z(s)=\mathbb E[e^{sZ}]\), where this expectation is finite. When it is finite on an open interval around zero, it determines the distribution, and its derivatives at zero are the moments: \(M_Z^{(k)}(0)=\mathbb E[Z^k]\). Its most useful property concerns sums. For independent \(X\) and \(Y\),

\[M_{X+Y}(s)=\mathbb E[e^{sX}e^{sY}]=M_X(s)\,M_Y(s).\]

Multiplying entries of this table proves several facts stated above: independent binomials with a common \(q\) add their trial counts, Poisson counts add their rates, common-rate gammas add their shapes, and independent Normals add their means and variances. Not every distribution has an MGF that is finite near zero; the Cauchy and lognormal laws do not. The characteristic function \(\mathbb E[e^{isZ}]\), with \(i=\sqrt{-1}\), always exists and serves the same purposes; Appendix C uses it to prove the central limit theorem. Exponential moments also give the sharp tail bounds later in this chapter.

Functions of random variables

Many quantities are computed from other random variables: a total, a transformed measurement, a sorted value. Their distributions follow from the joint law by the probability bookkeeping above, and three cases recur: sums, smooth invertible maps, and order statistics.

Sums and convolution

If \(X\) and \(Y\) are independent with densities \(p_X\) and \(p_Y\), then \(S=X+Y\) has density

\[p_S(s)=\int p_X(x)\,p_Y(s-x)\,dx.\]

This is the convolution of the two densities: to reach the total \(s\), the first term takes some value \(x\) and the second must equal \(s-x\). For integer-valued variables, the integral becomes a sum. Two independent \(\operatorname{Unif}(0,1)\) variables have the triangular density \(\min(s,2-s)\) on \([0,2]\), and each further term smooths the density again. Repeated convolution is why sums of many independent terms become approximately Normal, the subject of the central limit theorem below. Moment generating functions turn convolution into multiplication, which is often the easier calculation.

Adding independent uniform variables smooths the density: two give a triangle, and three already lie close to the Normal density with the same mean and variance (dashed).

Change of variables

Let \(Y=g(X)\), where \(g\) is strictly increasing and differentiable with a differentiable inverse. The events \(\{X\le x\}\) and \(\{Y\le g(x)\}\) coincide, so \(F_Y(y)=F_X(g^{-1}(y))\). Differentiating gives

\[p_Y(y)=p_X(g^{-1}(y))\left|\frac{d}{dy}g^{-1}(y)\right| =\frac{p_X(x)}{|g'(x)|},\qquad x=g^{-1}(y).\]

The absolute value makes the same formula hold for decreasing maps. The formula conserves probability: a short interval of length \(dx\) near \(x\) carries probability \(p_X(x)\,dx\), and \(g\) maps it to an interval of length \(|g'(x)|\,dx\) near \(y\). Where \(g\) stretches intervals, density falls; where it compresses them, density rises.

Passing \(X\sim\mathcal N(0,1.8^2)\) through the logistic function. The two shaded bands carry equal probability; the flat end of the curve squeezes its band of \(x\) values into a short range of \(y\), so the density of \(Y\) piles up near \(0\) and \(1\).

The same argument works in \(d\) dimensions. Suppose \(g:\mathbb R^d\to\mathbb R^d\) is invertible and continuously differentiable with a continuously differentiable inverse on the relevant open domains. With the numerator-layout Jacobian \(J_g(x)=(\partial g_i/\partial x_j)_{ij}\), the transformed density is

\[p_Y(y)=p_X(g^{-1}(y))\left|\det J_{g^{-1}}(y)\right| =\frac{p_X(x)}{|\det J_g(x)|},\qquad x=g^{-1}(y).\]

The determinant corrects for local volume change: a region expanded by a factor \(c\) must have its density divided by \(c\) to preserve probability. For a map with several inverse branches, add the contributions of the branches. For \(Y=X^2\) and \(y>0\),

\[p_Y(y)=\frac{p_X(\sqrt y)+p_X(-\sqrt y)}{2\sqrt y}.\]

These formulas assume an absolutely continuous \(X\); point masses require separate accounting. Maps between spaces of different dimensions need an extra step, given in Appendix B. Change of variables reappears with MAP estimates in Bayesian inference, and generative models use it to turn simple noise into complicated data.

The CDF transform and inverse-CDF sampling

The CDF is itself a useful transformation. If \(F\) is continuous and \(X\sim F\), then \(U=F(X)\) is uniform on \((0,1)\): \(P(F(X)\le u)=u\) for \(0<u<1\). Conversely, let \(U\sim\operatorname{Unif}(0,1)\) and let \(F^{-1}\) be the generalized quantile function. Since \(F^{-1}(u)\le x\) exactly when \(u\le F(x)\),

\[P\{F^{-1}(U)\le x\}=P\{U\le F(x)\}=F(x),\]

so \(F^{-1}(U)\) has CDF \(F\). This converse holds for discrete and mixed laws as well. It is inverse-CDF sampling: uniform random numbers become draws from any distribution whose quantile function can be evaluated. For \(\operatorname{Exp}(\lambda)\), solving \(u=1-e^{-\lambda x}\) gives \(F^{-1}(u)=-\log(1-u)/\lambda\).

Evenly spaced levels \(u\) map through the Exp\((1)\) quantile function to values that crowd where the density is high (left). Applied to uniform random numbers, the same map produces exponential draws (right).

Order statistics

Sorting iid observations \(X_1,\ldots,X_n\) with CDF \(F\) gives the order statistics \(X_{(1)}\le\cdots\le X_{(n)}\). The extremes have simple laws, because the maximum is at most \(x\) exactly when every observation is:

\[P(X_{(n)}\le x)=F(x)^n,\qquad P(X_{(1)}>x)=[1-F(x)]^n.\]

When \(F\) has density \(p\), the \(k\)th smallest value has density

\[p_{X_{(k)}}(x)=\frac{n!}{(k-1)!\,(n-k)!}\,F(x)^{k-1}[1-F(x)]^{n-k}\,p(x):\]

\(k-1\) observations fall below \(x\), one falls at \(x\), and \(n-k\) fall above it, and the factorials count the ways to choose which. For uniform observations on \((0,1)\), \(X_{(k)}\sim\operatorname{Beta}(k,n+1-k)\), with mean \(k/(n+1)\). Sample medians and other sample quantiles are order statistics, and so is the sample maximum that estimates an unknown endpoint in Appendix E. Gaps between order statistics, extremes, and records appear in Appendix B.

The sorted values of five uniform observations have Beta\((k,6-k)\) densities, with means \(k/6\).

Gaussian vectors and conditioning

A vector \(Z\in\mathbb R^d\) is jointly Gaussian if every linear combination \(a^\top Z\) is Gaussian, allowing constant combinations. With mean \(\mu\) and positive definite covariance \(\Sigma\), its density is

\[p(z)=\frac{\exp[-\tfrac12(z-\mu)^\top\Sigma^{-1}(z-\mu)]} {(2\pi)^{d/2}\det(\Sigma)^{1/2}}.\]

The quadratic form measures distance relative to the covariance: directions with greater variance are penalized less. If \(Z\sim\mathcal N(\mu,\Sigma)\), then \(AZ+b\sim\mathcal N(A\mu+b,A\Sigma A^\top)\). A singular covariance describes a Gaussian supported on a lower-dimensional affine subspace; the full-dimensional density above does not apply.

An explicit construction links this distribution to linear algebra. If \(G\) has independent standard normal coordinates and \(LL^\top=\Sigma\), then \(Z=\mu+LG\) has the desired Gaussian law. For positive definite \(\Sigma\), a Cholesky factor supplies such an \(L\). The inverse map \(L^{-1}(Z-\mu)\) produces independent standard normal coordinates. For general non-Gaussian vectors, the same covariance transformation can make the coordinates uncorrelated, but it need not make them independent.

Partition a jointly Gaussian vector and its parameters as

\[\begin{pmatrix}X\\Y\end{pmatrix} \sim\mathcal N\!\left( \begin{pmatrix}\mu_X\\\mu_Y\end{pmatrix}, \begin{pmatrix}\Sigma_{XX}&\Sigma_{XY}\\\Sigma_{YX}&\Sigma_{YY}\end{pmatrix} \right).\]

If \(\Sigma_{XX}\) is invertible, then

\[Y\mid X=x\sim\mathcal N\!\left( \mu_Y+\Sigma_{YX}\Sigma_{XX}^{-1}(x-\mu_X), \Sigma_{YY}-\Sigma_{YX}\Sigma_{XX}^{-1}\Sigma_{XY} \right).\]

To see the structure, define \(A=\Sigma_{YX}\Sigma_{XX}^{-1}\) and \(R=Y-\mu_Y-A(X-\mu_X)\). Direct calculation gives \(\operatorname{Cov}(R,X)=0\). Jointly Gaussian blocks with zero cross-covariance are independent, so conditioning on \(X\) changes the affine term but leaves \(R\) unchanged in distribution. Calculating \(\operatorname{Cov}(R)\) gives the displayed conditional covariance. Thus Gaussian conditional means are affine and conditional variances do not depend on the observed value \(x\).

A standard bivariate Normal pair with correlation \(0.7\). Conditioning on \(X=1.2\) takes a vertical slice: the conditional law of \(Y\) is centered on the line \(\mathbb E[Y\mid X=x]=0.7x\), and its standard deviation \(\sqrt{1-0.7^2}\approx0.71\) is the same for every \(x\) and smaller than the marginal standard deviation \(1\).

The joint-Gaussian assumption is essential: normal marginal distributions alone do not imply a joint Gaussian law. Murphy’s Probabilistic Machine Learning: An Introduction, Chapters 2–3, develops these probability models and their multivariate structure.

Averages, limit theorems, and concentration

For iid \(Z_1,\ldots,Z_n\) with mean \(\mu\) and finite variance \(\sigma^2\), let \(\bar Z_n=n^{-1}\sum_iZ_i\). Then \(\mathbb E\bar Z_n=\mu\) and \(\operatorname{Var}(\bar Z_n)=\sigma^2/n\). These exact identities are the starting point for both asymptotic approximations and finite-sample bounds.

Modes of convergence

Limit theorems describe a sequence of random variables \(T_1,T_2,\ldots\) that settles down, and there are several senses of settling:

• almost surely: the numbers \(T_n(\omega)\) converge to \(T(\omega)\) for every outcome \(\omega\) in an event of probability one;
• in probability, \(T_n\xrightarrow{p}T\): \(P(|T_n-T|>\varepsilon)\to0\) for every \(\varepsilon>0\);
• in mean square: \(\mathbb E[(T_n-T)^2]\to0\), and more generally in \(L^r\) when \(\mathbb E|T_n-T|^r\to0\);
• in distribution, \(T_n\Rightarrow T\): the CDFs converge, \(F_{T_n}(t)\to F_T(t)\), at every point \(t\) where \(F_T\) is continuous.

Almost-sure convergence and \(L^r\) convergence each imply convergence in probability, which in turn implies convergence in distribution. None of the reverse implications holds in general, with one exception: convergence in distribution to a constant implies convergence in probability to it. Convergence in distribution concerns only the laws of the \(T_n\), so the variables need not be defined on one probability space. Appendix C gives examples that separate these notions.

Laws of large numbers and the central limit theorem

The iid strong law of large numbers states that

\[\mathbb E|Z_1|<\infty \quad\Longrightarrow\quad \bar Z_n\longrightarrow\mu\quad\text{almost surely}.\]

It implies convergence in probability, the weak law of large numbers. Finite variance is sufficient but not necessary for either law.

The iid central limit theorem (CLT), assuming \(0<\sigma^2<\infty\), states

\[\frac{\sqrt n(\bar Z_n-\mu)}{\sigma}\Rightarrow\mathcal N(0,1).\]

The law of large numbers describes where the average settles. The CLT describes the distribution of its fluctuations after rescaling. It does not say that individual observations become Gaussian. Its usual practical approximation is \(\bar Z_n\approx\mathcal N(\mu,\sigma^2/n)\), and there is no universal sample size at which this approximation becomes accurate. Strong skewness and rare large observations can require much larger samples.

The assumptions cannot be inferred from a large dataset alone. If \(Z_i=Z_1\) for every \(i\), the mean never changes, even if every marginal distribution has finite variance. At the other extreme, iid Cauchy observations remain Cauchy when averaged, so no finite population mean exists for the average to approach. These examples fail for different reasons: dependence in the first and lack of integrability in the second. Many extensions of the limit theorems allow controlled dependence or weaker tail assumptions, but their conclusions require corresponding hypotheses.

Running averages of independent Exp\((1)\) draws settle at the mean \(1\), inside a band that narrows like \(1/\sqrt n\) (left). Running averages of Cauchy draws keep jumping, because the Cauchy law has no mean (right).

The standardized mean \(T_n=\sqrt n(\bar Z_n-1)\) of \(n\) Exp\((1)\) observations loses its skew as \(n\) grows. The histograms summarize 50,000 simulated datasets per panel; the dashed curves are the exact densities, and the solid curve is the standard Normal density.

For \(n=1\), no averaging has occurred: \(T_1=Z_1-1\) is a shifted exponential variable, whose density jumps from zero to one at \(t=-1\). For general \(n\), the support begins at \(-\sqrt n\), since every observation is nonnegative. The exact curves follow from \(\sum_iZ_i\sim\operatorname{Gamma}(n,1)\); as \(n\) grows, the standardized distribution becomes less skewed and its central region approaches the Normal curve.

The following simulation compares the central probability of standardized means with the standard normal value \(P(|G|\le1)\approx0.6827\):

import numpy as np

rng = np.random.default_rng(11)
repetitions = 30_000
for n in (1, 5, 30, 100):
    samples = rng.exponential(scale=1.0, size=(repetitions, n))
    standardized = np.sqrt(n) * (samples.mean(axis=1) - 1.0)
    print(n, "central probability:", np.mean(np.abs(standardized) <= 1))

Here the exponential mean and standard deviation both equal one. Each row is a new sample; variation across rows approximates a sampling distribution.

Bounds with explicit assumptions

The limit theorems describe what happens as the sample grows. The inequalities below hold at every sample size under stated assumptions, and each follows from Markov’s inequality.

Markov’s inequality. For \(W\ge0\) and \(a>0\),

\[P(W\ge a)\le\frac{\mathbb EW}{a}.\]

The proof is \(W\ge a\mathbf1_{\{W\ge a\}}\), followed by expectation. Applying this to \(W=(Z-\mu)^2\) gives Chebyshev’s inequality:

\[P(|Z-\mu|\ge\varepsilon)\le\frac{\sigma^2}{\varepsilon^2}.\]

For iid averages, it gives \(P(|\bar Z_n-\mu|\ge\varepsilon)\le\sigma^2/(n\varepsilon^2)\), proving the weak law when the variance is finite. These bounds hold for every sample size but can be loose. Appendix C gives a one-sided refinement, Cantelli’s inequality.

Chernoff bounds. Markov’s inequality applied to \(e^{sW}\) and to \(e^{-sW}\) with \(s>0\) gives

\[P(W\ge t)\le\inf_{s>0}e^{-st}M_W(s), \qquad P(W\le t)\le\inf_{s>0}e^{st}M_W(-s).\]

For a sum of independent terms, the MGF is a product, so these bounds decay exponentially in the number of terms. For \(W\sim\operatorname{Binomial}(n,q)\) with mean \(\mu=nq\), bounding each factor \(1-q+qe^s\) of the MGF by \(e^{q(e^s-1)}\) and then minimizing gives

\[\begin{aligned} P\{W\ge(1+\delta)\mu\}&\le\left(\frac{e^\delta}{(1+\delta)^{1+\delta}}\right)^{\mu}, &&\delta>0,\\ P\{W\le(1-\delta)\mu\}&\le\left(\frac{e^{-\delta}}{(1-\delta)^{1-\delta}}\right)^{\mu}\le e^{-\mu\delta^2/2}, &&0<\delta<1. \end{aligned}\]

With \(n=1000\) and \(q=0.01\), the probability of at least \(20\) events, twice the mean of \(10\), is at most \((e/4)^{10}\approx0.021\). Chebyshev’s inequality gives only \(9.9/10^2\approx0.099\), and the exact probability is about \(0.0033\). Near the mean, the exponent of a Chernoff bound is close to the Gaussian exponent \(t^2/(2\sigma^2)\) of the Normal approximation. Far from the mean it can grow much more slowly, and the Normal approximation then understates the tail by orders of magnitude.

Hoeffding’s inequality. If \(Z_i\) are independent and \(a_i\le Z_i\le b_i\) almost surely, then for \(t>0\) and \(\sum_i(b_i-a_i)^2>0\),

\[P\!\left(\left|\sum_i(Z_i-\mathbb EZ_i)\right|\ge t\right) \le2\exp\!\left(-\frac{2t^2}{\sum_i(b_i-a_i)^2}\right).\]

For iid observations in \([0,1]\), this becomes \(P(|\bar Z_n-\mu|\ge\varepsilon)\le2e^{-2n\varepsilon^2}\). Thus \(n\ge\log(2/\alpha)/(2\varepsilon^2)\) guarantees an error below \(\varepsilon\) with probability at least \(1-\alpha\). Unlike the CLT approximation, this conclusion is a finite-sample guarantee under the stated boundedness and independence assumptions. Evaluating a fixed classifier’s zero-one loss on iid test examples fits this setting; choosing the classifier using the same examples changes the dependence structure and needs additional analysis.

Hoeffding’s inequality is a Chernoff bound that uses only the ranges of the summands, so it needs no other knowledge of their distributions. The price is that it is far from tight when the variances are much smaller than the ranges allow, as for indicators of rare events. Appendix C proves it, compares it with the Chernoff bound computed from the exact distribution, describes the Gaussian and exponential regimes of Chernoff bounds, and gives bounds that also use the variances of the summands, such as Bernstein’s inequality, together with their extensions to random matrices.

Monte Carlo integration

For an expectation \(I=\mathbb E[h(Z)]\) that is difficult to integrate analytically, Monte Carlo draws iid samples and computes \(\widehat I_n=n^{-1}\sum_i h(Z_i)\). If \(h(Z)\) has finite variance \(v\), then

\[\mathbb E\widehat I_n=I,\qquad \operatorname{Var}(\widehat I_n)=v/n.\]

The usual error scale is therefore \(\sqrt{v/n}\). Multiplying the sample count by four halves this scale. The exponent \(n^{-1/2}\) does not depend on the dimension of \(Z\), although the variance and cost of generating samples can depend strongly on dimension.

For example, \(I=\int_0^1e^{-u^2}\,du=\mathbb E[e^{-U^2}]\) for \(U\sim\operatorname{Unif}(0,1)\):

import numpy as np

rng = np.random.default_rng(19)
u = rng.uniform(size=100_000)
values = np.exp(-(u ** 2))
estimate = values.mean()
estimated_se = values.std(ddof=1) / np.sqrt(values.size)
print("integral estimate:", estimate)
print("estimated standard error:", estimated_se)

The standard error describes variation of the estimate across repeated simulations. It differs from the standard deviation of the integrand values themselves. This distinction carries directly into statistics, where the observations come from data collection and their law is unknown.

Relative error in integrating \(\prod_je^{-u_j^2}\) over the unit cube in \(d=1\) and \(d=10\) dimensions. A midpoint grid with \(m\) points per coordinate uses \(N=m^d\) evaluations, and its error falls like \(N^{-2/d}\); the Monte Carlo root-mean-square error falls like \(N^{-1/2}\) in every dimension.

When the integrand is concentrated where \(Z\) rarely falls, plain Monte Carlo wastes most of its draws. Importance sampling draws from a proposal density \(q\) instead and reweights:

\[\mathbb E_p[h(Z)]=\mathbb E_q\!\left[h(Z)\frac{p(Z)}{q(Z)}\right] \approx\frac1n\sum_{i=1}^n h(Z_i)\frac{p(Z_i)}{q(Z_i)}, \qquad Z_i\overset{\mathrm{iid}}{\sim}q.\]

The estimator is unbiased when \(q>0\) wherever \(h\,p\ne0\). To estimate the Normal tail probability \(P(G>4)\approx3.2\times10^{-5}\), plain sampling needs about \(32{,}000\) draws per tail event, while drawing from \(\mathcal N(4,1)\) puts half the draws in the tail. The weighted terms \(h\,p/q\) can have a large or infinite variance when \(q\) has lighter tails than \(|h|\,p\), and the estimate then becomes unstable. The same reweighting reappears in approximate inference and in evaluating one policy with data collected under another.

Simulation also distinguishes uncertainty about data from computational error. If a fully specified model permits sampling, Monte Carlo error can be reduced by drawing more samples from that model. More simulated data cannot correct a model that poorly represents the real process. In statistical applications, the sampling distribution of an estimator may itself be approximated by simulation; there are then two layers of repetition, one corresponding to hypothetical datasets and one to the computer’s finite approximation of that distribution.

Statistics

Data, models, and inferential targets

Probability describes observations generated from a specified law. Statistical inference uses observations to learn about a law that is incompletely known. Estimation, confidence intervals, tests, and posterior distributions answer different versions of this question. Their interpretation depends on the population of interest and on how the data were obtained.

A population is the collection of units or possible observations to which a conclusion refers. It can be finite, such as the messages received during one month, or represented by a distribution for possible observations. A sample is the collection actually observed. Write

\[D=(Z_1,\ldots,Z_n),\qquad Z_i\sim P_*,\]

where \(P_\ast\) is the true population distribution. Chapter 1 writes this distribution as \(\mathcal D\) when it defines population risk; this chapter uses \(P_\ast\) to keep it visibly distinct from the dataset \(D\). Capital letters denote random observations; \(z_i\) denotes a realized value. The dataset \(D\) is random before collection. After collection, its realized values are fixed, although the same symbol is often used when the distinction is clear. An observation may be a scalar, a vector, or an input–target pair \(Z_i=(X_i,Y_i)\).

Sampling assumptions and study design

The assumption \(Z_1,\ldots,Z_n\overset{\mathrm{iid}}{\sim}P_\ast\) says both that observations are independent and that they have the same distribution. Neither follows from storing observations in separate rows. Two measurements on one person can be dependent. Observations collected before and after a change may have different distributions. A random sample from one population need not represent another population.

Sampling design determines what the dataset can reveal. Uniform random sampling from a finite population produces a different joint law from sampling independently with replacement. Unequal inclusion probabilities can require weighting. Repeated measurements, clustered samples, and time series require uncertainty calculations that retain their dependence. Increasing the number of rows does not eliminate selection bias or make correlated observations independent.

Random sampling supports inference from the sample to a population. Random assignment of treatments supports causal comparisons within an experiment. These are different sources of randomness. Observational associations can be estimated precisely without identifying a causal effect. For example, the difference between the mean outcomes of two self-selected groups is a descriptive contrast; interpreting it as the effect of group membership requires additional assumptions.

Models and identifiability

A statistical model is a family of candidate probability laws,

\[\mathcal P=\{P_\theta:\theta\in\Theta\}.\]

The parameter \(\theta\) indexes a law, and \(p_\theta\) denotes its mass function or density. A model with a fixed, finite-dimensional parameter space is parametric. A nonparametric model allows an unknown distribution or function without a fixed finite-dimensional restriction. A semiparametric model contains a finite-dimensional parameter of interest together with an infinite-dimensional component. All three kinds impose assumptions.

A model is correctly specified if \(P_\ast=P_{\theta_\ast}\) for some \(\theta_\ast\in\Theta\). Otherwise it is misspecified. A model is identifiable when \(P_{\theta_1}=P_{\theta_2}\) implies \(\theta_1=\theta_2\). If two parameter values generate exactly the same observable distribution, no amount of data can distinguish them. A particular target can nevertheless be identifiable even when the full parameter is not: only parameter values that give different target values must be distinguishable.

Suppose the model is \(Z\sim\mathcal N(a+b,1)\). The distribution identifies \(a+b\), but it cannot identify \(a\) and \(b\) separately. This is a property of the model, rather than a failure of an estimation algorithm.

Estimands, statistics, and estimators

An estimand is the quantity the analysis aims to learn. It can be a parameter or a functional of a distribution:

\[\tau=T(P_*).\]

Examples are a population mean, a quantile, the difference between two means, or the expected loss of a fixed predictor. Defining the population and the estimand comes before selecting a formula for estimating it.

The same statistic can serve different purposes. A sample average estimates a probability for Bernoulli data, a rate for Poisson counts, or a mean for general numerical observations. The algebra can be identical while the units, sampling distribution, and interpretation change. A nuisance parameter is an unknown component needed to describe the model but not itself the target; the unknown variance is a nuisance parameter when estimating a Normal mean.

Point estimation returns one value. Interval estimation returns a range with an uncertainty interpretation. A statistical procedure can also return a test decision or a probability distribution. Frequentist inference evaluates procedures through repeated sampling with the population law held fixed. Bayesian inference introduces a probability distribution for unknown parameters and conditions it on the observed data. The likelihood connects both approaches.

Empirical distributions and the plug-in principle

The empirical distribution places probability \(1/n\) on each observed value, counting repetitions:

\[\widehat P_n=\frac1n\sum_{i=1}^n\delta_{Z_i},\]

where \(\delta_z\) is the point mass at \(z\). It replaces an unknown distribution by a distribution constructed entirely from the sample. For scalar observations, its cumulative distribution function is

\[\widehat F_n(x)=\frac1n\sum_{i=1}^n\mathbf1\{Z_i\le x\}.\]

Under iid sampling, each indicator is Bernoulli with probability \(F_\ast(x)\), so, at each fixed \(x\),

\[\mathbb E[\widehat F_n(x)]=F_*(x),\qquad \operatorname{Var}(\widehat F_n(x))=\frac{F_*(x)(1-F_*(x))}{n}.\]

The law of large numbers gives pointwise convergence. The Glivenko–Cantelli theorem strengthens this to simultaneous convergence over all thresholds: \(\sup_x|\widehat F_n(x)-F_\ast(x)|\to0\) almost surely. Estimating a distribution therefore has a concrete meaning even without fitting a named family.

The Dvoretzky–Kiefer–Wolfowitz inequality makes the uniform convergence quantitative at every sample size:

\[P\Big(\sup_x|\widehat F_n(x)-F_*(x)|>\varepsilon\Big)\le2e^{-2n\varepsilon^2}.\]

Setting the right side equal to \(\alpha\) gives a band of half-width \(\sqrt{\log(2/\alpha)/(2n)}\) around \(\widehat F_n\) that contains the entire population CDF with probability at least \(1-\alpha\).

The plug-in principle estimates \(T(P_\ast)\) by \(T(\widehat P_n)\). For the population mean and variance, it gives

\[\widehat\mu=\overline Z=\frac1n\sum_i Z_i, \qquad \widehat v_n=\frac1n\sum_i(Z_i-\overline Z)^2.\]

The plug-in estimate of a quantile \(F_\ast^{-1}(q)\) is \(\widehat F_n^{-1}(q)\), using the generalized inverse defined above. For observations \(2,4,4,7,10\), this gives median \(4\) and \(0.8\)-quantile \(7\). Software sometimes interpolates between order statistics, so finite-sample quantile conventions can differ.

The empirical CDF of \(2,4,4,7,10\) jumps by \(1/5\) at each value, by \(2/5\) at the repeated value \(4\), and first reaches \(0.8\) at \(7\) (left). For 100 standard Normal draws, the 95% DKW band of half-width \(0.136\) contains the whole population CDF (right).

The target determines which summary is meaningful. The data \(0,0,1,1,100\) have mean \(20.4\) and median \(1\). The median is less sensitive to the largest observation, but it estimates a different population quantity. Robustness describes stability under departures from a model or contamination; it does not make different estimands interchangeable. The empirical-distribution viewpoint, estimation, and testing are developed together in Wasserman’s All of Statistics, especially Chapters 6–10.

Sampling distributions and estimator accuracy

The sampling distribution of \(\widehat\tau=t(D)\) is its distribution over repeated datasets generated by the sampling model. It differs from the distribution of an individual observation and from the empirical distribution within one dataset. These three distributions answer different questions: variation among units, variation among estimates, and the observed frequencies in one sample.

Bias, variance, and mean squared error

For a scalar estimand \(\tau\) and an estimator with finite second moment, define

\[\operatorname{Bias}(\widehat\tau)=\mathbb E[\widehat\tau]-\tau, \qquad \operatorname{SE}(\widehat\tau)=\sqrt{\operatorname{Var}(\widehat\tau)}.\]

The expectation holds the population law fixed and averages over datasets. An estimator is unbiased if this bias is zero throughout the model. Its standard error is the standard deviation of its sampling distribution. An estimated standard error replaces unknown features of that distribution by estimates.

The mean squared error is

\[\operatorname{MSE}(\widehat\tau) =\mathbb E[(\widehat\tau-\tau)^2] =\operatorname{Var}(\widehat\tau) +\operatorname{Bias}(\widehat\tau)^2.\]

To obtain the decomposition, write \(\widehat\tau-\tau=(\widehat\tau-\mathbb E\widehat\tau)+(\mathbb E\widehat\tau-\tau)\), expand the square, and observe that the cross term has expectation zero. Bias and variance contribute differently. A modest increase in bias can be worthwhile if it produces a larger decrease in variance.

For iid observations with mean \(\mu\) and finite variance \(\sigma^2\),

\[\mathbb E[\overline Z]=\mu, \qquad \operatorname{Var}(\overline Z)=\frac{\sigma^2}{n}.\]

The observation standard deviation \(\sigma\) measures variation among individual values; the standard error \(\sigma/\sqrt n\) measures variation among sample means. Multiplying the sample size by four halves this standard error. For dependent observations, however,

\[\operatorname{Var}(\overline Z) =\frac1{n^2}\sum_{i,j=1}^n\operatorname{Cov}(Z_i,Z_j),\]

so positive correlations can substantially increase uncertainty.

The plug-in variance \(\widehat v_n\) is biased. The identity

\[\sum_i(Z_i-\overline Z)^2 =\sum_i(Z_i-\mu)^2-n(\overline Z-\mu)^2\]

gives \(\mathbb E[\widehat v_n]=(n-1)\sigma^2/n\). Thus

\[S^2=\frac1{n-1}\sum_i(Z_i-\overline Z)^2\]

is unbiased when \(n\ge2\). Both denominators yield consistent variance estimators. Unbiasedness is one criterion, rather than a universal reason to prefer one estimator.

For iid vector observations \(Z_i\in\mathbb R^d\) with finite second moments, the corresponding sample covariance matrix is

\[\widehat\Sigma=\frac1{n-1}\sum_{i=1}^n(Z_i-\overline Z)(Z_i-\overline Z)^\top.\]

It is unbiased for the population covariance matrix. The empirical-distribution plug-in version uses denominator \(n\). If the centered observations form the rows of a matrix \(X_c\in\mathbb R^{n\times d}\), then \(\widehat\Sigma=X_c^\top X_c/(n-1)\). The sample correlation between coordinates \(j\) and \(k\) is \(\widehat\Sigma_{jk}/\sqrt{\widehat\Sigma_{jj}\widehat\Sigma_{kk}}\) when both sample variances are positive. These estimates summarize variation and linear association within the observed sample; their accuracy depends on the sampling assumptions and sample size.

A scalar example illustrates the bias–variance tradeoff. If \(\widehat p=K/n\) with \(K\sim\operatorname{Binomial}(n,p)\), the smoothed estimator \((K+1)/(n+2)\) has bias \((1-2p)/(n+2)\) and variance \(np(1-p)/(n+2)^2\). Comparing the two mean squared errors shows that smoothing helps exactly when \(p(1-p)>n/[4(2n+1)]\): for a central range of \(p\) the reduction in variance outweighs the squared bias, while near \(p=0\) or \(1\) the sample proportion is better.

import numpy as np

rng = np.random.default_rng(7)
n, repetitions = 20, 100_000
for p in (0.02, 0.5):
    counts = rng.binomial(n, p, size=repetitions)
    estimators = {"sample proportion": counts / n,
                  "smoothed": (counts + 1) / (n + 2)}
    for name, estimates in estimators.items():
        bias = estimates.mean() - p
        variance = estimates.var()  # empirical distribution denominator
        mse = np.mean((estimates - p) ** 2)
        assert np.isclose(mse, variance + bias**2)
        print(p, name, "bias", round(bias, 4), "MSE", round(mse, 5))

Each entry of estimates comes from a fresh dataset. The simulation therefore approximates a sampling distribution, rather than the spread of observations in one sample.

With \(n=20\), smoothing toward \(1/2\) lowers the mean squared error for \(0.14<p<0.86\) and raises it near the endpoints.

Exact distributions for Normal samples

Most sampling distributions are known only approximately, but Normal samples have exact results at every sample size. For iid \(\mathcal N(\mu,\sigma^2)\) observations and \(n\ge2\),

\[\overline Z\sim\mathcal N\Big(\mu,\frac{\sigma^2}{n}\Big), \qquad \frac{(n-1)S^2}{\sigma^2}\sim\chi^2_{n-1}, \qquad \overline Z\ \text{and}\ S^2\ \text{are independent}.\]

There is a geometric reason for these identities. Standardize the sample to a vector \(G=(Z-\mu\mathbf1)/\sigma\) with independent standard Normal coordinates, where \(Z=(Z_1,\ldots,Z_n)^\top\) and \(\mathbf1\) is the all-ones vector. Its projection onto the direction \(\mathbf1/\sqrt n\) is \(\sqrt n(\overline Z-\mu)/\sigma\). Its orthogonal projection onto the \((n-1)\)-dimensional residual subspace has squared norm \(\sum_i(Z_i-\overline Z)^2/\sigma^2\). Orthogonal Gaussian components are independent, and the residual squared norm is a sum of \(n-1\) independent squared standard Normals. This gives both the chi-square law and the independence.

Consequently, the studentized mean has an exact Student \(t\) distribution:

\[\frac{\overline Z-\mu}{S/\sqrt n}\sim t_{n-1},\]

since it is a standard Normal divided by the square root of an independent \(\chi^2_{n-1}/(n-1)\). The independence of \(\overline Z\) and \(S^2\) is special to Normal samples, and the same projection argument gives the exact theory of Normal linear regression later in the chapter. Write \(z_q=\Phi^{-1}(q)\), \(t_{\nu,q}\), and \(\chi^2_{\nu,q}\) for the \(q\) quantiles of the standard Normal, \(t_\nu\), and \(\chi^2_\nu\) distributions.

Consistency, asymptotic normality, and the delta method

An estimator sequence is consistent if \(\widehat\tau_n\xrightarrow{p}\tau\). This is a large-sample property. Unbiasedness is a statement about an expectation at a fixed sample size. An unbiased estimator need not become more accurate: using \(Z_1\) to estimate \(\mu\) is unbiased but ignores every additional observation. Conversely, \(\widehat v_n\) is biased for each finite \(n\) and consistent.

An estimator is asymptotically Normal at rate \(\sqrt n\) if

\[\sqrt n(\widehat\tau_n-\tau)\Rightarrow\mathcal N(0,V),\]

with \(0<V<\infty\). The asymptotic variance \(V\) is the variance of the scaled limiting distribution. The resulting Normal approximation has variance \(V/n\); the stronger conclusion \(n\operatorname{Var}(\widehat\tau_n)\to V\) requires additional control of moments. Slutsky’s theorem allows a quantity converging in probability to a constant to replace that constant in sums, products, and ratios of distributional limits, with a nonzero limiting denominator for ratios. No independence between the converging quantities is required. If \(\widehat V\xrightarrow{p}V\), it gives

\[\frac{\widehat\tau_n-\tau}{\sqrt{\widehat V/n}} \Rightarrow\mathcal N(0,1).\]

This studentization replaces an unknown scale by a consistent estimate. It is the basis of many intervals and tests. A large sample alone does not guarantee the approximation: finite variance, suitable dependence assumptions, and regularity of the estimator matter.

The delta method applies a local linear approximation to a random estimator. If \(\widehat\theta_n\in\mathbb R^d\) satisfies

\[\sqrt n(\widehat\theta_n-\theta)\Rightarrow\mathcal N(0,V)\]

and \(g:\mathbb R^d\to\mathbb R^k\) is differentiable at \(\theta\), then

\[\sqrt n\{g(\widehat\theta_n)-g(\theta)\} \Rightarrow \mathcal N\!\left(0,J_g(\theta)VJ_g(\theta)^\top\right).\]

Here the numerator-layout Jacobian has shape \(k\times d\). The result follows by expanding \(g(\widehat\theta_n)-g(\theta)\) as \(J_g(\theta)(\widehat\theta_n-\theta)\) plus a remainder negligible on the \(n^{-1/2}\) scale. For scalar \(g\), the variance is \(\nabla g(\theta)^\top V\nabla g(\theta)\), with a column gradient.

For a Bernoulli proportion with fixed \(0<p<1\) and \(g(p)=\log[p/(1-p)]\), \(g'(p)=1/[p(1-p)]\). The delta method gives the asymptotic standard-error approximation

\[\operatorname{SE}_{\mathrm{asymp}}\{g(\widehat p)\} \approx\frac1{\sqrt{np(1-p)}}.\]

Replacing \(p\) by \(\widehat p\) gives an estimated standard error when \(0<\widehat p<1\); at \(\widehat p=0\) or \(1\) the log-odds estimate is undefined. If \(J_g(\theta)=0\), the first-order limit is degenerate and a higher-order expansion is needed.

The delta method for the log-odds of a proportion near \(p=0.2\). Over the spread of \(\widehat p\) the logit is nearly linear, so the log-odds are approximately Normal, with standard deviation multiplied by the slope \(1/[p(1-p)]=6.25\).

The sample median is asymptotically Normal as well. If the population has a unique median \(m\) at which its density \(p_Z\) is continuous and positive, then \(\sqrt n(\widehat m_n-m)\Rightarrow\mathcal N(0,1/[4p_Z(m)^2])\). The variance depends on the density at the median: a flat density there means that small changes in the data move the median far.

Constructing point estimators

Maximum likelihood is the main general-purpose way to construct estimators, and the one most directly connected to training with log loss. The method of moments, which matches sample averages to model expectations, is a simpler alternative described in Appendix E.

Likelihood and maximum likelihood

The likelihood evaluates a model on the observed data, treating the parameter as variable. Under iid sampling,

\[L_n(\theta)=\prod_{i=1}^n p_\theta(z_i), \qquad \log L_n(\theta)=\sum_{i=1}^n\log p_\theta(z_i).\]

For dependent observations, the appropriate joint density or conditional factorization replaces the iid product. \(L_n\) is not a probability distribution over \(\theta\) and need not integrate to one. Likelihood values are compared across parameters for the same observed data. Multiplying every likelihood value by the same positive, parameter-independent constant changes none of those comparisons.

The maximum likelihood estimator satisfies

\[\widehat\theta_{\mathrm{MLE}} \in\operatorname*{arg\,max}_{\theta\in\Theta}\log L_n(\theta).\]

A maximum can be nonunique or fail to exist. Solving a derivative equation finds an interior stationary point, which still must be checked against boundaries and other candidates. Minimizing average negative log-likelihood gives the same maximizer; this is the connection between likelihood estimation and empirical risk minimization with log loss. Here \(L_n\) denotes likelihood, whereas \(\ell(y,a)\) elsewhere denotes a prediction loss.

For Bernoulli observations, let \(k=\sum_i z_i\). Then

\[\log L_n(p)=k\log p+(n-k)\log(1-p), \qquad \frac{d}{dp}\log L_n(p)=\frac{k}{p}-\frac{n-k}{1-p}.\]

For \(0<k<n\), setting the derivative to zero gives \(\widehat p=k/n\). The log-likelihood is strictly concave on \((0,1)\), so this is its unique maximum. If all observations are failures or all are successes, the maximum is respectively \(0\) or \(1\) when those endpoints belong to the parameter space.

Eight successes in ten trials and eighty in a hundred give the same estimate, \(0.8\), but the larger sample’s likelihood is more sharply peaked; each curve is divided by its maximum. On the log scale (right), dashed parabolas share each curve’s peak and curvature, and that curvature is the observed information of the next subsection.

For iid \(\mathcal N(\mu,v)\) observations, writing \(v=\sigma^2>0\) for the variance,

\[\log L_n(\mu,v) =-\frac n2\log(2\pi v)-\frac1{2v}\sum_i(z_i-\mu)^2.\]

Differentiating first in \(\mu\) gives \(\widehat\mu=\overline z\). Substituting that value and differentiating in \(v\) gives \(\widehat v=n^{-1}\sum_i(z_i-\overline z)^2\), provided the residual sum of squares is positive. The MLE uses denominator \(n\) because likelihood maximization and unbiasedness are different requirements. If all observations coincide, the likelihood increases without bound as \(v\downarrow0\), so this model has no positive-variance MLE for that dataset.

MLE is equivariant under one-to-one reparameterizations: if \(\phi=g(\theta)\), its MLE is \(g(\widehat\theta)\). Likelihood maximization identifies a fitted probability law, and changing its coordinates does not change that law.

For iid \(\operatorname{Exp}(\lambda)\) observations, \(\log L_n(\lambda)=n\log\lambda-\lambda\sum_iz_i\), which is maximized at \(\widehat\lambda=1/\overline z\). By equivariance, the MLE of the mean \(1/\lambda\) is \(\overline z\).

The Bernoulli, Normal, and exponential models above are exponential families, as are the Poisson and several other families of this chapter. Their likelihoods depend on the data only through a low-dimensional sufficient statistic, such as \(\sum_iz_i\), and their maximum likelihood estimates set the model’s expected value of that statistic equal to its sample average. Appendix E develops sufficiency and exponential families.

Why should maximizing likelihood recover the parameter? By the law of large numbers, \(n^{-1}\log L_n(\theta)\to\mathbb E_{P_\ast}[\log p_\theta(Z)]\) for each fixed \(\theta\). When the model is correct, \(P_\ast=P_{\theta_\ast}\), and

\[\mathbb E_{\theta_*}[\log p_{\theta_*}(Z)]-\mathbb E_{\theta_*}[\log p_\theta(Z)] =D_{\mathrm{KL}}(P_{\theta_*}\,\|\,P_\theta)\ge0,\]

the Kullback–Leibler divergence of chapter 5. For an identifiable model it is zero only at \(\theta=\theta_\ast\), so the limiting average log-likelihood is maximized at the true parameter. Turning this pointwise limit into consistency of the maximizer requires control that is uniform over \(\theta\), which standard regularity conditions supply.

Score, information, and asymptotic uncertainty

The score is the column gradient \(U_n(\theta)=\nabla_\theta\log L_n(\theta)\). The observed information is its negative derivative,

\[J_n(\theta)=-\nabla_\theta^2\log L_n(\theta).\]

For one observation, write \(s_\theta(Z)=\nabla_\theta\log p_\theta(Z)\). The Fisher information is

\[I_1(\theta)=\mathbb E_\theta[s_\theta(Z)s_\theta(Z)^\top].\]

Under conditions permitting differentiation under the integral, with parameter-independent support,

\[\mathbb E_\theta[s_\theta(Z)]=0, \qquad I_1(\theta)=-\mathbb E_\theta[\nabla_\theta^2\log p_\theta(Z)].\]

The first identity follows by differentiating \(\int p_\theta(z)\,dz=1\). Differentiating again gives the second. Independence makes information add: \(I_n(\theta)=nI_1(\theta)\). Curvature measures how rapidly the log-likelihood changes locally; greater information corresponds to smaller uncertainty in regular estimation problems.

In a correctly specified regular model with an interior true parameter, a consistent MLE typically satisfies

\[\sqrt n(\widehat\theta-\theta_*) \Rightarrow\mathcal N(0,I_1(\theta_*)^{-1}).\]

The conditions include identifiability, smoothness, finite nonsingular information, and sufficient control of derivatives and sample averages. A Taylor expansion of the score equation explains the result:

\[0=U_n(\widehat\theta) \approx U_n(\theta_*)-nI_1(\theta_*)(\widehat\theta-\theta_*).\]

The score is a sum of iid mean-zero vectors. Its CLT produces the Normal limit, and the inverse curvature translates score fluctuations into parameter fluctuations. The limiting Normal approximation has covariance \(I_1(\theta_\ast)^{-1}/n\), estimated by \([nI_1(\widehat\theta)]^{-1}\) or, when justified, \(J_n(\widehat\theta)^{-1}\). Boundary and parameter-dependent-support problems can have different limits and rates. The Cramér–Rao bound, which turns information into a lower bound on the variance of unbiased estimators, is in Appendix E.

When the model is only an approximation

Even a misspecified family can provide a useful fitted description. Under suitable consistency conditions, maximum likelihood then targets

\[\theta^\dagger\in\operatorname*{arg\,max}_{\theta\in\Theta} \mathbb E_{P_*}[\log p_\theta(Z)],\]

rather than a true parameter whose model distribution equals \(P_\ast\). This population optimum is sometimes called the pseudo-true parameter. Equivalently, \(\theta^\dagger\) minimizes the Kullback–Leibler divergence \(D_{\mathrm{KL}}(P_\ast\,\|\,P_\theta)\) over the model: it indexes the member of the family closest to the population in that sense. For example, fitting a Normal distribution to a non-Normal population still targets its mean and variance when these are finite. A useful estimator can therefore outlive an inaccurate distributional assumption, but its uncertainty calculation must be reconsidered.

Under misspecification, the information identity need not equate log-likelihood curvature with score variability. Inverse curvature alone can therefore give incorrect standard errors. A sandwich covariance estimate accounts for both quantities under suitable iid regularity assumptions; its formula and limitations appear in Appendix E.

Confidence intervals and repeated-sampling coverage

A \(1-\alpha\) confidence procedure returns a random set \(C(D)\) such that

\[P_\theta\{\tau(\theta)\in C(D)\}\ge1-\alpha\]

for every parameter value covered by its guarantee. Exact procedures satisfy this at the stated sample size; asymptotic procedures approach the nominal coverage under their assumptions. In repeated sampling, the target is fixed and the interval changes. After observation, a particular interval either contains the target or does not. The confidence level describes the generating procedure rather than a posterior probability for that realized interval.

Exact 95% Student \(t\) intervals from forty samples of size \(20\) drawn from \(\mathcal N(0,1)\). Thirty-eight of them cover the true mean \(0\); the 95% describes the procedure over many repetitions, not any single interval.

A pivot is a function of the sample and unknown parameters whose distribution does not depend on those unknown parameters. Inverting a probability statement about a pivot often yields a confidence interval.

Means: Normal and Student intervals

If observations are iid \(\mathcal N(\mu,\sigma^2)\) with known \(\sigma\), then \((\overline Z-\mu)/(\sigma/\sqrt n)\) is exactly standard Normal. Rearranging its central probability statement gives

\[\overline Z\pm z_{1-\alpha/2}\frac{\sigma}{\sqrt n}.\]

For unknown variance, the exact results for Normal samples make the studentized mean a pivot with the \(t_{n-1}\) distribution. The same rearrangement gives the exact interval

\[\overline Z\pm t_{n-1,1-\alpha/2}\frac S{\sqrt n}.\]

The heavier tails of \(t\) account for estimating the scale. For iid non-Normal observations with finite nonzero variance, studentized means have a large-sample Normal limit; the same \(t\) formula is then an approximation rather than an exact finite-sample statement. The chi-square pivot \((n-1)S^2/\sigma^2\) gives a similarly exact interval for a Normal variance, but that interval relies heavily on Normal tails and does not become approximately valid for other distributions (Appendix F). NIST gives the mean-interval derivation and interpretation.

For two measurements \(A_i,B_i\) on each independently sampled unit, analyze \(W_i=B_i-A_i\). The paired mean difference has estimated standard error \(S_W/\sqrt n\). This incorporates within-unit covariance: \(\operatorname{Var}(B-A)=\operatorname{Var}(B)+\operatorname{Var}(A)-2\operatorname{Cov}(A,B)\). Treating paired measurements as independent discards that information.

For genuinely independent groups, the estimated standard error of \(\overline B-\overline A\) is \(\sqrt{S_B^2/n_B+S_A^2/n_A}\). A large-sample Normal interval follows from the CLT. The Welch \(t\) approximation improves the small-sample treatment without imposing equal population variances; its degrees of freedom depend on both estimated variances and sample sizes.

More generally, an asymptotically Normal estimator with an estimated standard error gives the Wald interval \(\widehat\tau\pm z_{1-\alpha/2}\widehat{\operatorname{SE}}(\widehat\tau)\). Its accuracy is only as good as the Normal approximation and the standard error.

Proportions and score inversion

The Wald approximation \(\widehat p\pm z_{1-\alpha/2}\sqrt{\widehat p(1-\widehat p)/n}\) can perform poorly for small samples or probabilities near \(0\) and \(1\). If no successes occur, it collapses to \([0,0]\), despite substantial uncertainty.

The Wilson interval asks, for each candidate \(p\), whether the observed proportion is unusually far from that value relative to the sampling variation implied by \(p\). Keeping the candidates that are not rejected gives an interval. For \(0<p<1\), its Normal approximation retains values satisfying

\[\frac{|\widehat p-p|}{\sqrt{p(1-p)/n}}\le z, \qquad z=z_{1-\alpha/2}.\]

Equivalently, solve \(n(\widehat p-p)^2\le z^2p(1-p)\) over \(0\le p\le1\), including the degenerate boundary cases by continuity. This avoids division by zero at the endpoints. Solving the quadratic inequality yields

\[\frac{\widehat p+z^2/(2n) \ \pm\ z\sqrt{\widehat p(1-\widehat p)/n+z^2/(4n^2)}} {1+z^2/n}.\]

The endpoints stay in \([0,1]\). This is an approximate interval obtained by inverting a score test (Appendix F), and its finite-sample coverage is not uniformly at least the nominal level. Clopper–Pearson intervals invert exact binomial tails and provide at least nominal coverage, generally with some conservatism because counts are discrete. NIST develops both proportion-interval constructions.

import numpy as np
from scipy.stats import norm, binomtest

k, n, alpha = 18, 100, 0.05
phat = k / n
z = norm.ppf(1 - alpha / 2)
denom = 1 + z*z / n
center = (phat + z*z / (2*n)) / denom
halfwidth = z * np.sqrt(phat*(1-phat)/n + z*z/(4*n*n)) / denom
print("estimate:", phat)
print("estimated SE:", np.sqrt(phat*(1-phat)/n))
print("Wilson:", (center-halfwidth, center+halfwidth))
exact = binomtest(k, n).proportion_ci(confidence_level=1-alpha, method="exact")
print("Clopper–Pearson:", (exact.low, exact.high))

Exact coverage of nominal 95% intervals for a binomial proportion with \(n=30\). The Wald interval’s coverage oscillates and collapses near the endpoints; the Wilson interval stays close to \(0.95\) except very near \(0\) or \(1\).

Bootstrap approximation of sampling uncertainty

Analytic standard errors can be difficult for medians, correlations, and other statistics. The nonparametric bootstrap applies the plug-in principle to the sampling distribution itself. Replace \(P_\ast\) by \(\widehat P_n\), sample \(n\) observations independently from that empirical distribution, and recompute the estimator. Sampling from \(\widehat P_n\) is exactly sampling the observed rows with replacement.

Let \(\widehat\tau^\ast\) denote an estimate from such a bootstrap sample. Conditional on the observed dataset, the bootstrap attempts to approximate

\[\operatorname{Law}_{P_*}(\widehat\tau-\tau) \quad\text{by}\quad \operatorname{Law}_{\widehat P_n}(\widehat\tau^*-\widehat\tau).\]

There are two approximations: replacing the population by the empirical distribution, and estimating the resulting law using finitely many simulations. Increasing the number of bootstrap replicates improves the second approximation but does not provide additional original observations. The foundational construction is due to Efron (1979).

The sampling distribution of the mean of \(30\) Exp\((1)\) observations, from fresh datasets (left), and the bootstrap distribution from resampling one of those datasets (right). The bootstrap distribution is centered at that dataset’s mean rather than at the population mean \(1\), but its spread and skew are close.

For \(B\) replicates \(\widehat\tau^{\ast(1)},\ldots,\widehat\tau^{\ast(B)}\), their sample standard deviation \(\widehat{\operatorname{SE}}_{\mathrm{boot}}\) estimates the standard error of \(\widehat\tau\), and \(\widehat\tau\pm z_{1-\alpha/2}\widehat{\operatorname{SE}}_{\mathrm{boot}}\) is the normal bootstrap interval. The replicates’ \(\alpha/2\) and \(1-\alpha/2\) quantiles form a percentile interval. If those quantiles are \(q^\ast_{\alpha/2}\) and \(q^\ast_{1-\alpha/2}\), a basic interval is

\[[2\widehat\tau-q^*_{1-\alpha/2},\ 2\widehat\tau-q^*_{\alpha/2}].\]

It inverts the bootstrap distribution of the centered error. Studentized and bias-corrected-and-accelerated methods refine these constructions; none is automatically reliable for every statistic and sample size.

The bootstrap is most useful when no convenient formula exists. The sample median’s approximate variance \(1/[4np_Z(m)^2]\) requires the population density at the median, which is hard to estimate well. The standard deviation of resampled medians avoids that step. Resampled medians take only a few distinct values determined by the observed sample, however (the observations themselves when \(n\) is odd), so for small \(n\) their distribution is lumpy and the resulting intervals are coarse.

For paired measurements, the observational unit is a pair. Resampling the two columns independently destroys their association. The following example keeps each pair intact and estimates the mean within-unit change.

import numpy as np

rng = np.random.default_rng(18)
n, replicates = 120, 4000
baseline = rng.normal(10, 3, size=n)
followup = baseline + 0.7 + rng.normal(0, 1.5, size=n)
pairs = np.column_stack([baseline, followup])
change = pairs[:, 1] - pairs[:, 0]
estimate = change.mean()

indices = rng.integers(0, n, size=(replicates, n))
resampled = pairs[indices]  # shape: (replicates, observations, 2)
boot = (resampled[:, :, 1] - resampled[:, :, 0]).mean(axis=1)
lower, upper = np.quantile(boot, [0.025, 0.975])
print("mean change:", estimate)
print("analytic SE:", change.std(ddof=1) / np.sqrt(n))
print("bootstrap SE:", boot.std(ddof=1))
print("percentile interval:", (lower, upper))
print("basic interval:", (2*estimate-upper, 2*estimate-lower))

The synthetic measurements share a baseline effect. Keeping the rows paired preserves that effect’s cancellation in the difference. For a correlation or another statistic involving both columns, the same row-resampling principle preserves their joint empirical distribution.

Ordinary iid bootstrap theory works well for many smooth statistics under adequate moment and regularity assumptions. A sample median additionally requires suitable behavior of the density near the population median. Maxima and support endpoints can fail because empirical resamples never contain a value beyond the observed maximum. Infinite-variance means can fail because their limiting behavior is not captured by ordinary Normal and iid-bootstrap approximations. Time series and clusters require appropriate block or cluster resampling. A bootstrap calculation reproduces the assumptions encoded in its resampling scheme; it cannot correct selection bias.

A parametric bootstrap instead simulates data from a fitted model \(P_{\widehat\theta}\) and refits each simulated dataset. It can use structure unavailable to empirical resampling, but its validity depends more directly on the fitted family. Both forms approximate repeated data collection, rather than drawing a Bayesian posterior for the parameter. The jackknife, a deterministic leave-one-out relative of the bootstrap, appears in Appendix E.

Hypothesis tests, evidence, and power

A hypothesis is a restriction on the population distribution. In a parametric model, write

\[H_0:\theta\in\Theta_0, \qquad H_1:\theta\in\Theta_1,\]

with disjoint sets. A simple hypothesis specifies one distribution; a composite hypothesis allows several. A test chooses a rejection region, often through a scalar statistic that measures departure from the null.

A Type I error rejects a true null. A Type II error fails to reject at a specified alternative. A level-\(\alpha\) test satisfies

\[\sup_{\theta\in\Theta_0}P_\theta(\text{reject }H_0)\le\alpha.\]

Its power function is \(\operatorname{Power}(\theta)=P_\theta(\text{reject }H_0)\), and its Type II error probability at an alternative is \(1-\operatorname{Power}(\theta)\). Power is a function of the alternative, not a single property of the test independent of effect size.

Test–interval duality says that keeping every target value whose level-\(\alpha\) test is not rejected produces a \(1-\alpha\) confidence set when the tests are valid. The Wilson interval above is one example, and any family of valid tests yields intervals in this way. The interpretation is still coverage; inversion does not turn a p-value into a posterior probability. Appendix F shows that likelihood-ratio tests are most powerful for simple hypotheses and describes the Wald, likelihood-ratio, and score tests of regular likelihood theory.

What a p-value states

For a simple null and a statistic \(T\) whose larger values indicate stronger disagreement, the tail probability

\[p_{\mathrm{val}}=P_{H_0}\{T(D')\ge T(D_{\mathrm{obs}})\}\]

uses a hypothetical fresh dataset \(D'\) drawn under the null. Two-sided tests choose an appropriate measure of extremeness, often an absolute standardized difference. A valid p-value obeys

\[P_\theta(p_{\mathrm{val}}\le u)\le u \qquad(0\le u\le1,\ \theta\in\Theta_0).\]

This property justifies rejecting when \(p_{\mathrm{val}}\le\alpha\). Continuous exact-null p-values are often uniform; discreteness can make them conservative. Composite hypotheses require validity across nuisance parameters, possibly through conditioning, a pivot, or a worst-case calculation.

The p-value is a probability calculated under the null, not the probability that the null is true. A large p-value can reflect a small effect, large noise, or insufficient data. Failing to reject equality does not establish practical equivalence. Equivalence is a different hypothesis, requiring a prespecified tolerance for negligible differences.

For example, suppose the observations are Normal with known \(\sigma\), so that \(T=\sqrt n(\overline Z-\mu_0)/\sigma\) is standard Normal under \(H_0:\mu=\mu_0\), and the observed value of \(T\) is \(2.10\). The two-sided p-value is

\[P_{H_0}(|T|\ge2.10)=2[1-\Phi(2.10)]\approx0.036.\]

A level-\(0.05\) test rejects, and a level-\(0.01\) test does not. The value \(0.036\) is not the probability that the null hypothesis is true.

The two-sided p-value for an observed statistic of \(2.10\) is the null probability beyond \(\pm2.10\), which lies inside the level-\(0.05\) rejection region beyond \(\pm1.96\).

For \(K\sim\operatorname{Binomial}(n,p)\), testing \(p=p_0\) against \(p>p_0\) gives the exact p-value

\[\sum_{j=k}^n {n\choose j}p_0^j(1-p_0)^{n-j}.\]

No large-sample Normal approximation is needed. For a two-sided discrete test, “at least as extreme” must be specified; doubling the smaller tail and summing outcomes with probability no greater than the observed outcome can produce different valid conventions. SciPy documents its exact binomial implementation in binomtest.

Effect sizes and sample size

In a Normal mean problem with known \(\sigma\), a two-sided level-\(\alpha\) test rejects when \(|\overline Z-\mu_0|>z_{1-\alpha/2}\sigma/\sqrt n\). At a true mean \(\mu=\mu_0+\Delta\), define \(d=\sqrt n\Delta/\sigma\) and \(z=z_{1-\alpha/2}\). Its power is

\[\operatorname{Power}(\mu)=1-\Phi(z-d)+\Phi(-z-d).\]

The standardized shift grows as \(\sqrt n\). A tiny nonzero effect can therefore become statistically significant with enough data, while an important effect can be missed by a small noisy sample. Reporting an effect estimate and uncertainty interval conveys information that a thresholded test result discards.

For a one-sided Normal test against a prespecified positive difference \(\Delta\), with \(0<\alpha,\beta<1/2\), achieving power \(1-\beta\) requires

\[n\ge\frac{\sigma^2(z_{1-\alpha}+z_{1-\beta})^2}{\Delta^2},\]

rounded up to an integer. The calculation specifies an alternative and assumes known variance; planning with an estimated scale adds uncertainty. Choosing the direction of a one-sided test after looking at the sign of the data invalidates its stated error rate.

The standardized statistic under the null and under an alternative with standardized shift \(d=2.5\). The two-sided rejection region has probability \(0.05\) under the null and about \(0.71\) under this alternative; the rest of the alternative’s probability is the Type II error.

Comparing two proportions

A common experiment compares success rates in two independent groups, such as the conversion rates of users randomly shown two versions of a web page. With \(K_A\sim\operatorname{Binomial}(n_A,p_A)\) and \(K_B\sim\operatorname{Binomial}(n_B,p_B)\), the estimated difference \(\widehat p_B-\widehat p_A\) has variance \(p_A(1-p_A)/n_A+p_B(1-p_B)/n_B\). Under \(H_0:p_A=p_B\), both groups share one rate, estimated by the pooled proportion \(\widehat p=(K_A+K_B)/(n_A+n_B)\), and

\[Z=\frac{\widehat p_B-\widehat p_A}{\sqrt{\widehat p(1-\widehat p)\,(1/n_A+1/n_B)}}\]

is approximately standard Normal for large samples. With \(200\) conversions among \(2{,}000\) users in group A and \(250\) among \(2{,}000\) in group B, \(\widehat p=0.1125\), \(Z\approx2.50\), and the two-sided p-value is about \(0.012\). An interval for the difference is not computed under the null, so it uses the unpooled variance: here \(0.025\pm0.020\). The calculation assumes independent users, random assignment, and a sample size fixed in advance. Checking the p-value repeatedly and stopping at the first significant result inflates the Type I error rate. The same test is Pearson’s chi-square test of independence for the \(2\times2\) table of counts; Appendix F gives chi-square tests for general tables.

Permutation and randomization tests

Some null hypotheses imply an exact symmetry. For two independent iid samples from the same distribution, pooled observations are exchangeable under the null: relabeling which observations belong to which group leaves their joint law unchanged. Comparing the observed statistic with all allowed relabelings gives an exact permutation test, with conservatism possible because the reference distribution is discrete.

Equality of two population means alone does not imply exchangeability. If two groups have different variances or shapes, an unstudentized label-permutation test need not test only equality of means at its nominal level. The symmetry being used is an assumption that needs its own justification.

In a paired randomized experiment, swapping treatment labels within pairs gives the randomization distribution under a sharp null of no effect on any unit. For observational paired differences, sign flipping is justified by independent differences whose null distribution is symmetric about zero. Mean zero alone is insufficient. With \(n\) pairs there are \(2^n\) possible sign patterns, making exact enumeration practical for small samples.

import itertools
import numpy as np

# Illustrative paired differences. The null assumes independent,
# symmetric distributions about zero (or a valid paired randomization).
differences = np.array([1.2, 0.7, -0.1, 0.8, 1.0, 0.4, -0.2, 0.9])
observed = abs(differences.mean())
signs = np.array(list(itertools.product([-1, 1], repeat=len(differences))))
null_statistics = np.abs((signs * differences).mean(axis=1))
p_value = np.mean(null_statistics >= observed - 1e-12)
print("mean difference:", differences.mean())
print("exact two-sided sign-flip p-value:", p_value)

The calculation conditions on the observed magnitudes and enumerates every sign pattern. The tolerance includes computationally tied values. For a random subset of permutations, including the observed arrangement through the conventional \((b+1)/(B+1)\) calculation avoids zero p-values and yields a valid Monte Carlo test under the appropriate random-permutation construction. The official SciPy permutation-test documentation distinguishes independent-sample, paired, and pairing permutations.

Bootstrap and permutation procedures generate different reference worlds. The ordinary bootstrap estimates sampling fluctuations around the fitted empirical distribution. A permutation test generates a reference distribution by a symmetry that holds under its null. Resampling code can look similar even though the probability arguments differ.

Multiple comparisons and selection

If \(m\) true nulls each have rejection probability exactly \(\alpha\) and their test decisions are independent, the probability of at least one false rejection is \(1-(1-\alpha)^m\). At \(m=20\) and \(\alpha=0.05\), it is about \(0.64\). If each rejection probability is only bounded above by \(\alpha\), the same expression is an upper bound under independence. Repeatedly trying analyses and reporting only the smallest p-value creates a related selection problem.

The familywise error rate is \(P(V\ge1)\), where \(V\) counts false rejections. Bonferroni rejects an individual hypothesis only when its valid p-value is at most \(\alpha/m\). The union bound gives familywise control under arbitrary dependence. Holm’s procedure sorts p-values \(p_{(1)}\le\cdots\le p_{(m)}\), compares \(p_{(j)}\) sequentially to \(\alpha/(m-j+1)\), and stops at the first failure. It also controls familywise error and is at least as powerful as Bonferroni.

The false discovery rate is instead

\[\operatorname{FDR}=\mathbb E\!\left[\frac{V}{\max(R,1)}\right],\]

where \(R\) is the total number of rejections. The Benjamini–Hochberg procedure finds the largest \(k\) with \(p_{(k)}\le qk/m\) and rejects the first \(k\) hypotheses; it rejects none if no such \(k\) exists. For valid independent p-values it controls FDR at most \(q\). Extensions permit specified positive-dependence conditions, but arbitrary dependence is not covered by the ordinary guarantee. FDR control bounds an expected proportion and does not promise that each realized rejection set contains at most a fraction \(q\) of errors.

Forty sorted p-values, thirty from true null hypotheses and ten from real effects. Benjamini–Hochberg at \(q=0.1\) rejects every hypothesis up to the last p-value below the line \(jq/m\), eight here, all of them real effects; Bonferroni at level \(0.1\) rejects only the five below \(0.1/m\).

Selection can occur across features, outcomes, subgroups, model versions, stopping times, and analysis choices. A reported standard error or p-value generally describes a fixed procedure under its assumptions. Reusing the same data to choose a promising result and then evaluating it as if it had been prespecified changes that procedure. Independent evaluation data, appropriate multiplicity adjustments, and methods designed for sequential or selective inference address different versions of this problem.

Bayesian inference

Frequentist analysis studies how a procedure behaves when the data are repeatedly sampled at a fixed parameter. Bayesian analysis also specifies a probability distribution for uncertainty about that parameter. Both begin with a sampling model; they answer different probability questions.

Let \(D=(Z_1,\ldots,Z_n)\) denote the random sample, \(L_n(\theta)=p_\theta(D)\) its likelihood, and \(\pi(\theta)\) a prior density. Conditioning on the observed value of \(D\), Bayes’ rule gives the posterior density

\[\pi(\theta\mid D) =\frac{L_n(\theta)\pi(\theta)}{m(D)}, \qquad m(D)=\int L_n(u)\pi(u)\,du.\]

The denominator is the marginal likelihood, or evidence. It normalizes the posterior; the integral averages likelihood values using the prior. For discrete parameters, replace the integral by a sum. A proper prior integrates to one, and the posterior requires \(0<m(D)<\infty\).

The prior and likelihood have distinct roles: the prior describes parameter uncertainty before observing \(D\); the likelihood describes how different parameter values explain the same observed data. Neither alone is the posterior. Priors can encode substantive knowledge or regularize poorly identified parameters. Conclusions depend on these choices, especially with limited data. Bayesian Data Analysis, Chapters 1–3 develops this framework and its elementary models.

Bernoulli observations and a Beta prior

Suppose \(Z_i\mid\theta\sim\operatorname{Bernoulli}(\theta)\) independently, with \(s=\sum_i Z_i\) successes. Choose \(\theta\sim\operatorname{Beta}(\alpha_0,\beta_0)\), where \(\alpha_0,\beta_0>0\) are the prior shape parameters. Subscript \(0\) denotes the prior; subscript \(n\) will denote the posterior after \(n\) observations. Multiplication gives

\[\pi(\theta\mid D) \propto \theta^s(1-\theta)^{n-s} \theta^{\alpha_0-1}(1-\theta)^{\beta_0-1} =\theta^{\alpha_0+s-1}(1-\theta)^{\beta_0+n-s-1}.\]

Therefore

\[\theta\mid D\sim\operatorname{Beta}(\alpha_0+s,\beta_0+n-s).\]

The prior is conjugate because updating stays within its distribution family; Gamma priors for Poisson counts and Dirichlet priors for categories work the same way (Appendix H). The posterior mean is

\[\mathbb E[\theta\mid D] =\frac{\alpha_0+s}{\alpha_0+\beta_0+n} =\frac{\alpha_0+\beta_0}{\alpha_0+\beta_0+n}\frac{\alpha_0}{\alpha_0+\beta_0} +\frac{n}{\alpha_0+\beta_0+n}\frac{s}{n}.\]

This is a weighted average of the prior mean and the sample proportion. For this formula, \(\alpha_0+\beta_0\) measures the prior’s weight relative to the sample size; the pseudo-count interpretation depends on whether one is discussing the mean or the mode.

Writing \(\alpha_n=\alpha_0+s\) and \(\beta_n=\beta_0+n-s\), the posterior variance is

\[\operatorname{Var}(\theta\mid D) =\frac{\alpha_n \beta_n}{(\alpha_n+\beta_n)^2(\alpha_n+\beta_n+1)}.\]

The estimate and its uncertainty are separate outputs. Observing eight successes in ten trials and eighty successes in a hundred trials gives similar evidence about the location of \(\theta\), but substantially different posterior concentration. A point estimate alone conceals that difference.

For example, a \(\operatorname{Beta}(2,2)\) prior and eight successes in ten trials yield a \(\operatorname{Beta}(10,4)\) posterior. The likelihood is maximized at \(0.8\), while the posterior mean is \(10/14\approx0.714\): the prior pulls the estimate toward \(0.5\).

A credible set \(C(D)\) satisfies

\[P\{\theta\in C(D)\mid D\}\ge1-\alpha\]

under the specified prior and model. An equal-tailed interval uses posterior quantiles \(\alpha/2\) and \(1-\alpha/2\); for a continuous posterior its probability is exactly \(1-\alpha\). Discrete posteriors can prevent exact equality. Its probability concerns the uncertain parameter after conditioning on the observed sample. A confidence interval instead guarantees repeated-sampling coverage for a fixed parameter. A credible interval need not have its nominal frequentist coverage at every parameter value.

A Beta\((2,2)\) prior updated by eight successes in ten trials becomes a Beta\((10,4)\) posterior; the shaded 95% equal-tailed credible interval runs from \(0.462\) to \(0.909\). Updating multiplies the prior by the likelihood (right) and renormalizes.

The following computes a posterior interval with the Beta distribution implementation in SciPy. Here ppf means inverse CDF; these are posterior quantiles.

import numpy as np
from scipy.stats import beta as beta_dist

alpha_0, beta_0 = 2.0, 2.0
successes, n = 8, 10
alpha_n = alpha_0 + successes
beta_n = beta_0 + n - successes
posterior = beta_dist(alpha_n, beta_n)
credible_interval = posterior.ppf([0.025, 0.975])
predictive_success = posterior.mean()
print("Posterior mean:", predictive_success)
print("95% credible interval:", credible_interval)
assert np.isclose(predictive_success, 10 / 14)
assert np.isclose(np.diff(posterior.cdf(credible_interval))[0], 0.95)

Prediction and shrinkage

Inference about \(\theta\) is different from prediction of a new observation. If \(Z_{\mathrm{new}}\) and \(D\) are conditionally independent given \(\theta\), the posterior predictive distribution is

\[p(z_{\mathrm{new}}\mid D) =\int p_\theta(z_{\mathrm{new}})\pi(\theta\mid D)\,d\theta.\]

For the Bernoulli example, its success probability is \(\mathbb E[\theta\mid D]\). Future trials are independent conditional on \(\theta\), but generally dependent after averaging over their shared uncertain parameter. A posterior predictive interval \(A(D)\) satisfies \(P\{Z_{\mathrm{new}}\in A(D)\mid D\}\ge1-\alpha\). It concerns the next observation, including its own variability, rather than only an unknown parameter.

For a continuous example, let \(Z_i\mid\mu\sim\mathcal N(\mu,\sigma^2)\) independently, with known \(\sigma^2>0\), and let \(\mu\sim\mathcal N(m_0,v_0)\) with \(v_0>0\). Completing the square in the posterior exponent gives

\[\mu\mid D\sim\mathcal N(m_n,v_n), \qquad v_n^{-1}=v_0^{-1}+n\sigma^{-2}, \qquad m_n=v_n\left(\frac{m_0}{v_0}+\frac{n\bar Z}{\sigma^2}\right).\]

Here inverse variance, called precision, adds across the prior and observations. The posterior mean shrinks \(\bar Z\) toward \(m_0\) according to their relative uncertainty. More data increase the weight on \(\bar Z\). The predictive distribution is

\[Z_{\mathrm{new}}\mid D\sim\mathcal N(m_n,\sigma^2+v_n).\]

The extra \(\sigma^2\) is future observation noise. More generally, total variance gives

\[\operatorname{Var}(Z_{\mathrm{new}}\mid D) =\mathbb E_{\theta\mid D}[\operatorname{Var}(Z_{\mathrm{new}}\mid\theta)] +\operatorname{Var}_{\theta\mid D}[\mathbb E(Z_{\mathrm{new}}\mid\theta)].\]

The terms describe variation within a specified parameter value and uncertainty across parameter values. They are often called aleatoric and epistemic uncertainty. This decomposition is model-dependent; it does not automatically include uncertainty about an incorrectly specified model.

A \(\mathcal N(0,1)\) prior for \(\mu\) and four observations with \(\sigma=2\) and mean \(1.5\) give the posterior \(\mathcal N(0.75,0.5)\), halfway between prior and data because their precisions are equal (left). The posterior predictive distribution of a new observation is much wider than the posterior for \(\mu\) (right).

Posterior predictive checks and prior sensitivity analysis examine whether the model and prior are adequate for the data at hand; Appendix H describes both.

Posterior means and MAP estimates

A posterior distribution can be summarized in several ways. The posterior mean minimizes posterior squared-error loss. A maximum a posteriori estimate is a maximizer of the posterior density:

\[\widehat\theta_{\mathrm{MAP}} \in\operatorname*{arg\,max}_\theta \{\log L_n(\theta)+\log\pi(\theta)\}.\]

Thus a negative log prior becomes a penalty. For example, a Gaussian prior on coefficients yields a quadratic penalty. MAP retains one point and does not propagate parameter uncertainty into prediction. Unlike posterior probabilities, a density mode is not preserved under a change of parameterization, so a MAP estimate depends on the coordinates chosen (Appendix H).

Decisions, loss, and risk

An estimate, a prediction, and a test result are all actions chosen from data. Let \(a\) be an action, \(d(D)\) a decision rule, and \(\ell(\tau(\theta),a)\) the loss when the estimand is \(\tau(\theta)\). Its frequentist risk is

\[R(\theta,d) =\mathbb E_\theta[\ell(\tau(\theta),d(D))].\]

The expectation varies the sample while holding \(\theta\) fixed. Under squared-error loss, risk is MSE; the sample mean in the known-variance Normal model has risk \(\sigma^2/n\). Unbiasedness alone cannot rank all estimators, because a biased estimator can have smaller MSE.

Given a prior, the Bayes risk averages this sampling risk:

\[r(\pi,d)=\int R(\theta,d)\pi(\theta)\,d\theta.\]

Interchanging expectations shows that it is minimized by minimizing posterior expected loss for each observed sample:

\[d^*(D)\in\operatorname*{arg\,min}_a \mathbb E[\ell(\tau(\theta),a)\mid D].\]

For a scalar target \(T=\tau(\theta)\), the principal cases are:

For squared loss, write

\[\mathbb E[(T-a)^2\mid D] =\operatorname{Var}(T\mid D) +\big(\mathbb E[T\mid D]-a\big)^2.\]

Only the second term depends on \(a\), proving the first row. For a continuous target, exact zero-one loss gives expected loss one for every individual action under a continuous posterior; it does not justify MAP without changing the loss or taking a suitable small-neighborhood limit.

An action also depends on its costs. Suppose \(Y\in\{0,1\}\) and \(p=P(Y=1\mid\text{available information})\). With false-positive cost \(c_{\mathrm{FP}}>0\), false-negative cost \(c_{\mathrm{FN}}>0\), and zero cost for correct decisions, predicting one costs \(c_{\mathrm{FP}}(1-p)\) in expectation. Predicting zero costs \(c_{\mathrm{FN}}p\). Therefore the optimal threshold is

\[p\ge\frac{c_{\mathrm{FP}}}{c_{\mathrm{FP}}+c_{\mathrm{FN}}}.\]

The familiar threshold \(1/2\) is a consequence of equal costs.

Expected cost of each decision when a false negative costs four times as much as a false positive. Predicting \(1\) becomes the cheaper action once \(p\) exceeds \(1/5\).

A rule chosen without a prior can be judged by its worst case, or by whether any other rule does at least as well at every parameter value. Appendix H defines these minimax and admissibility criteria and describes Stein’s paradox, in which shrinkage dominates the natural unbiased estimator of a mean vector.

Regression as statistical inference

Regression concerns the conditional distribution of a response \(Y\) given a covariate vector \(X\). Its regression function is \(m(x)=\mathbb E[Y\mid X=x]\). A linear conditional-mean model assumes \(m(x)=x^\top\beta\). This is an assumption about a conditional expectation, not a claim that the covariates or responses must themselves have Normal marginal distributions.

For observed covariate column vectors \(x_i\in\mathbb R^d\), define the design matrix \(\mathbf X\in\mathbb R^{n\times d}\) with rows \(x_i^\top\). Stack the random responses into \(\mathbf Y\in\mathbb R^n\), with observed realization \(y\). The model is

\[\mathbf Y=\mathbf X\beta+\varepsilon, \qquad \mathbb E[\varepsilon\mid \mathbf X]=0.\]

An intercept is represented by a column of ones and counts among the \(d\) coefficients. Conditioning on \(\mathbf X\) separates variation in the responses from variation in the observed design. Each coefficient describes a conditional comparison for this particular set of covariates; changing the covariate set can change its meaning.

Least squares and coefficient uncertainty

Assume \(\mathbf X\) has full column rank. Ordinary least squares minimizes \(\|\mathbf Y-\mathbf X\beta\|^2\) and satisfies

\[\mathbf X^\top \mathbf X\widehat\beta=\mathbf X^\top \mathbf Y, \qquad \widehat\beta=(\mathbf X^\top \mathbf X)^{-1}\mathbf X^\top \mathbf Y.\]

Geometrically, \(\mathbf X\widehat\beta\) is the orthogonal projection of \(\mathbf Y\) onto the column space of \(\mathbf X\). Statistically,

\[\widehat\beta-\beta=(\mathbf X^\top \mathbf X)^{-1}\mathbf X^\top\varepsilon.\]

This identity gives conditional unbiasedness. Writing \(\Omega=\operatorname{Cov}(\varepsilon\mid \mathbf X)\) gives

\[\operatorname{Cov}(\widehat\beta\mid \mathbf X) =(\mathbf X^\top \mathbf X)^{-1}\mathbf X^\top\Omega \mathbf X(\mathbf X^\top \mathbf X)^{-1}.\]

If \(\Omega=\sigma^2I_n\), the covariance simplifies to \(\sigma^2(\mathbf X^\top \mathbf X)^{-1}\). Under these mean and covariance assumptions, the Gauss–Markov theorem says that OLS has smallest covariance among estimators that are linear in \(\mathbf Y\) and unbiased for every \(\beta\). It does not claim superiority to every biased estimator, and it does not require Normal errors. Appendix G proves it and develops the projection geometry of the fit.

Exact small-sample inference adds the assumption \(\varepsilon\mid \mathbf X\sim\mathcal N(0,\sigma^2I_n)\) and requires \(n>d\). Set

\[\widehat\sigma^2=\frac{\|\mathbf Y-\mathbf X\widehat\beta\|^2}{n-d}, \qquad \widehat{\operatorname{SE}}(\widehat\beta_j) =\widehat\sigma\sqrt{[(\mathbf X^\top \mathbf X)^{-1}]_{jj}}.\]

Then \((\widehat\beta_j-\beta_j)/\widehat{\operatorname{SE}}(\widehat\beta_j)\) has a Student \(t\) distribution with \(n-d\) degrees of freedom. Consequently,

\[\widehat\beta_j\pm t_{n-d,1-\alpha/2}\widehat{\operatorname{SE}}(\widehat\beta_j)\]

is an exact marginal confidence interval. These intervals are not automatically simultaneous statements about all coefficients.

The degrees of freedom have a geometric explanation. Fitting \(d\) linearly independent columns leaves residuals in an \((n-d)\)-dimensional orthogonal subspace. Under spherical Gaussian noise, the fitted component and residual component are independent, and the squared residual length divided by \(\sigma^2\) has a \(\chi^2_{n-d}\) distribution. A standardized Normal coefficient error divided by the independent estimated noise scale therefore has the stated \(t\) distribution. Without Gaussian errors this exact argument no longer applies; large-sample inference needs its own assumptions.

With \(\mathrm{RSS}=\|\mathbf Y-\mathbf X\widehat\beta\|^2\) and \(\mathrm{TSS}=\sum_i(Y_i-\overline Y)^2\), the coefficient of determination \(R^2=1-\mathrm{RSS}/\mathrm{TSS}\) is the fraction of the sample variation around the mean explained by a model with an intercept. It never decreases when a column is added, so a larger \(R^2\) alone does not justify a larger model. The \(F\) test of Appendix G asks whether added columns improve a nested model by more than chance would.

Estimating a mean versus predicting an observation

At a fixed covariate vector \(x_0\), define

\[\widehat m(x_0)=x_0^\top\widehat\beta, \qquad h_0=x_0^\top(\mathbf X^\top \mathbf X)^{-1}x_0.\]

Under the Gaussian model, a confidence interval for the conditional mean and a prediction interval for an independent future response are, respectively,

\[\underbrace{\widehat m(x_0)\pm t_{n-d,1-\alpha/2}\widehat\sigma\sqrt{h_0}}_{\text{mean response}}, \qquad \underbrace{\widehat m(x_0)\pm t_{n-d,1-\alpha/2}\widehat\sigma\sqrt{1+h_0}}_{\text{future observation}}.\]

The prediction interval is wider because a new observation contains its own random error. The extra one represents this error variance in units of \(\sigma^2\). These are pointwise intervals at a specified \(x_0\), not simultaneous bands for every covariate value. The NIST discussion of mean-response uncertainty illustrates their repeated-sampling interpretation. With only an intercept, \(h_0=1/n\) and \(\widehat\sigma\) is the sample standard deviation \(S_Y\) of the responses, so a new response from a Normal population has the prediction interval \(\overline Y\pm t_{n-1,1-\alpha/2}S_Y\sqrt{1+1/n}\), with exact coverage jointly over the sample and the new observation.

Pointwise 95% intervals for the mean response (inner band) and for a new observation (outer band), computed from the data and fitted line of the NumPy example below. Both bands widen away from the center of the observed covariates.

import numpy as np
from scipy.stats import t

rng = np.random.default_rng(7)
u = np.linspace(-2, 2, 60)
X = np.column_stack([np.ones(u.size), u])
y = X @ np.array([1.0, 2.0]) + rng.normal(size=u.size)
beta_hat = np.linalg.lstsq(X, y, rcond=None)[0]
n, d = X.shape
s2 = np.sum((y - X @ beta_hat) ** 2) / (n - d)
x0 = np.array([1.0, 1.5])
mean_hat = x0 @ beta_hat
h0 = x0 @ np.linalg.solve(X.T @ X, x0)
critical = t.ppf(0.975, df=n - d)
mean_halfwidth = critical * np.sqrt(s2 * h0)
prediction_halfwidth = critical * np.sqrt(s2 * (1 + h0))
print("Mean interval:", mean_hat + np.array([-1, 1]) * mean_halfwidth)
print("Prediction interval:", mean_hat + np.array([-1, 1]) * prediction_halfwidth)
assert prediction_halfwidth > mean_halfwidth

The formulas express assumptions as well as calculations. Heteroskedastic errors can leave OLS unbiased while invalidating the constant-variance standard errors. Correlated observations likewise require an appropriate covariance estimate. Near-collinear covariates inflate coefficient uncertainty. Residual patterns can reveal omitted nonlinear structure or changing noise variance; a small training residual sum of squares does not establish that these assumptions hold. Appendix G collects diagnostics for these problems, including leverage, variance inflation, and robust standard errors.

Omitted variables and the best linear approximation

A coefficient’s meaning depends on the other covariates in the model. Suppose \(Y=\beta_0+\beta_1X+\beta_2W+\varepsilon\) with \(\mathbb E[\varepsilon\mid X,W]=0\), but \(Y\) is regressed on \(X\) alone, with an intercept. The population slope of that shorter regression is

\[\frac{\operatorname{Cov}(X,Y)}{\operatorname{Var}(X)} =\beta_1+\beta_2\frac{\operatorname{Cov}(X,W)}{\operatorname{Var}(X)}.\]

The second term is omitted-variable bias. It vanishes only if the omitted variable has no effect on \(Y\) or is uncorrelated with \(X\). The short regression still estimates a well-defined quantity, the best linear predictor of \(Y\) from \(X\) alone, but not the coefficient \(\beta_1\) of the longer model. Appendix G gives the general form of this bias and the attenuation caused by measurement error in a covariate, and it shows that each coefficient of a multiple regression is a partial association.

If the conditional mean is nonlinear, a best linear approximation can still be a useful estimand. With finite second moments and invertible \(\mathbb E[XX^\top]\), its population coefficient is

\[\beta^*=\operatorname*{arg\,min}_\beta\mathbb E[(Y-X^\top\beta)^2] =\mathbb E[XX^\top]^{-1}\mathbb E[XY].\]

Its residual is orthogonal to the covariates in expectation. This weaker property does not imply a correctly specified linear conditional mean, so the preceding exact Gaussian intervals cannot simply be carried over unchanged.

Finally, an observational coefficient measures association under the chosen conditioning set. Interpreting it as the effect of an intervention requires additional causal assumptions. Statistical significance alone supplies no such argument. Appendix G describes two designs that state such assumptions explicitly, instrumental variables and difference-in-differences.

From inference to learning

Let a fresh pair \((X,Y)\) have the true distribution \(P_\ast\), and let \(f(X)\) be a prediction. Its population risk and the corresponding empirical risk are

\[R_*(f)=\mathbb E_{P_*}[\ell(Y,f(X))], \qquad \widehat R_n(f)=\frac1n\sum_{i=1}^n\ell(Y_i,f(X_i)).\]

Conditioning on \(X=x\) reduces prediction to the decision problem above. Squared loss is minimized by \(\mathbb E[Y\mid X=x]\); absolute loss by a conditional median; classification under zero-one loss by a class of greatest conditional probability. The phrase Bayes predictor refers to these population-optimal rules even when no prior over parameters is used. The Elements of Statistical Learning, Section 2.4 develops this decision-theoretic view.

For a fixed \(f\) chosen independently of an iid evaluation sample, the empirical risk estimates population risk without selection bias. When the same data select \(f\), its training loss generally understates its future prediction loss. Held-out evaluation estimates predictive performance; a coefficient standard error quantifies a different uncertainty. Neither automatically answers the other’s question.

Appendices

import numpy as np
from scipy.stats import norm, chi2_contingency

conversions = np.array([200, 250])
users = np.array([2000, 2000])
rates = conversions / users
pooled = conversions.sum() / users.sum()
z = (rates[1] - rates[0]) / np.sqrt(pooled * (1 - pooled) * (1 / users).sum())
print("z:", z, "two-sided p:", 2 * norm.sf(abs(z)))
se_diff = np.sqrt((rates * (1 - rates) / users).sum())
print("95% interval for the difference:",
      rates[1] - rates[0] + np.array([-1, 1]) * norm.ppf(0.975) * se_diff)

table = np.column_stack([conversions, users - conversions])
chi2, p, dof, expected = chi2_contingency(table, correction=False)
print("Pearson chi-square:", chi2, "z squared:", z**2, "p:", p)