Sets up the triple of sample space, sigma-field and probability measure on which every later probabilistic statement is made, derives the elementary consequences of the axioms including continuity along monotone sequences, and defines conditioning and independence, from which the partition theorem and Bayes' theorem follow in two lines each.
Note(What the axioms are for)Original
Probability is a measure-theoretic subject wearing different clothes. A probability space is a measure space of total mass one, an event is a measurable set, and a random variable is a measurable function; the names change because the questions do. Folland (Folland (1999), §10.1, p. 314) sets the two vocabularies side by side in a table, and the translation is exact.
What the axioms buy is modest and indispensable. They do not say what the probability of any particular event is — that is part of the description of the experiment, not of the theory — and they do not define probability as a limiting frequency. They say only that probabilities are numbers in [0,1] that add over disjoint alternatives, including over countably many of them. Everything in the rest of this part is deduced from that, together with the two derived notions of this unit: conditioning, which restricts attention to a sub-population, and independence, which says that a restriction changes nothing.
Definition(Sample space, sigma-field, probability measure)AdaptedAdapted from Grimmett & Stirzaker (2020), Ch. 1, the recapitulation on p. 4 and Definition (1.3.1), p. 5; Folland (1999), §10.1, the dictionary on p. 314
Let Ω be a set, whose elements ω are the possible outcomes of an experiment; Ω is the sample space.
A collection F of subsets of Ω is a σ-field (measure theorists say σ-algebra) if it contains ∅, is closed under complementation, and is closed under countable unions. Its members are the events; to say that A is an event is exactly to say A∈F. Closure under countable intersections follows by De Morgan.
A probability measure on (Ω,F) is a function P:F→[0,1] with
P(∅)=0 and P(Ω)=1;
countable additivity: whenever A1,A2,… are members of F that are pairwise disjoint, meaning Ai∩Aj=∅ for i=j,
P(i=1⋃∞Ai)=i=1∑∞P(Ai).
The triple (Ω,F,P) is a probability space. It is a measure space (see Sigma-algebras and measures) whose total mass is 1.
An event A with P(A)=0 is called null; if P(A)=1 we say A occurs almost surely, abbreviated a.s. A null event is not the same as the impossible event ∅: a dart thrown at a target hits some particular point, and the event that it hits that point had probability zero.
The elementary case is a finite Ω with all outcomes equally likely, F the collection of all subsets, and P(A)=∣A∣/∣Ω∣. Then probability is counting, and the whole difficulty is combinatorial.
Let (Ω,F,P) be a probability space and let A,B,A1,A2,… be events.
Complement.P(Ac)=1−P(A).
Monotonicity. If A⊆B then P(B)=P(A)+P(B∖A)≥P(A).
Addition law.P(A∪B)=P(A)+P(B)−P(A∩B).
Inclusion-exclusion. For any n≥1,
P(i=1⋃nAi)=i∑P(Ai)−i<j∑P(Ai∩Aj)+i<j<k∑P(Ai∩Aj∩Ak)−⋯+(−1)n+1P(A1∩⋯∩An),
each sum running over unordered index sets of the stated size.
Countable subadditivity.P(⋃i≥1Ai)≤∑i≥1P(Ai), with no disjointness assumed.
Continuity along monotone sequences. If A1⊆A2⊆⋯ and A=⋃i≥1Ai, then P(Ai)→P(A) as i→∞. If B1⊇B2⊇⋯ and B=⋂i≥1Bi, then P(Bi)→P(B).
Property 6 is the one that distinguishes countable from merely finite additivity: given finite additivity, countable additivity holds if and only if P is continuous in this sense.
Idea of the proof(Everything is a disjoint union)Original
The only tool the axioms give is additivity over disjoint events, so every proof below consists of rewriting the event of interest as a disjoint union and then adding.
For the addition law, A∪B is A together with the part of B outside A. For a nested increasing sequence, A is A1 together with the successive shells Ai+1∖Ai, which are disjoint, so countable additivity turns P(A) into a telescoping series whose partial sums are exactly P(An). For subadditivity one cannot disjointify without losing mass, so one disjointifies downwards: replace Ai by the part of it not already covered, which is a subset, and use monotonicity to pay for the replacement.
Conditioning is a renormalisation. Restricting attention to the outcomes in B makes B the new sure event, so every probability must be divided by P(B) to restore total mass one; A then occurs exactly when A∩B does. The partition theorem and Bayes' theorem are what this renormalisation looks like read in the two possible directions: forwards, summing causes to get the probability of an effect; backwards, dividing to recover the probability of a cause from an observed effect.
Proof(Consequences 1-6)AdaptedAdapted from Grimmett & Stirzaker (2020), the proofs of Lemma (1.3.4) and Theorem (1.3.5), pp. 6-7
Throughout, disjoint unions are written out explicitly and axiom 2 of the definition is applied; a finite disjoint union is covered by taking all but finitely many of the Ai to be ∅, which is legitimate because P(∅)=0.
1.Ω=A∪Ac is a disjoint union, so 1=P(Ω)=P(A)+P(Ac).
2. If A⊆B then B=A∪(B∖A) is a disjoint union, so P(B)=P(A)+P(B∖A), and the second term is non-negative because P takes values in [0,1].
3. Write A∪B=A∪(B∖A), again disjoint, so P(A∪B)=P(A)+P(B∖A). Now B∖A=B∖(A∩B) and A∩B⊆B, so step 2 gives P(B∖A)=P(B)−P(A∩B). Substituting gives the addition law.
4. Induction on n. The case n=2 is step 3. Suppose the formula holds for n−1 events and put C=A1∪⋯∪An−1. By step 3,
P(C∪An)=P(C)+P(An)−P(C∩An),
and C∩An=⋃i≤n−1(Ai∩An) is a union of n−1 events, so the inductive hypothesis applies to both P(C) and P(C∩An). Collecting terms of each size: a set S⊆{1,…,n} with n∈/S contributes once, from P(C), with sign (−1)∣S∣+1; a set containing n with ∣S∣≥2 contributes once, from −P(C∩An), where it appears as S∖{n} with sign (−1)∣S∣, so after the minus sign its sign is (−1)∣S∣+1 as required; and S={n} contributes the term P(An).
5. Disjointify downwards: set D1=A1 and Di=Ai∖(A1∪⋯∪Ai−1) for i≥2. The Di are pairwise disjoint with the same union as the Ai, and Di⊆Ai. Hence by countable additivity and then monotonicity,
P(i⋃Ai)=i∑P(Di)≤i∑P(Ai).
6. Let A1⊆A2⊆⋯ with union A. Put C1=A1 and Ci+1=Ai+1∖Ai for i≥1. These are pairwise disjoint, since the Ai increase, and ⋃i≤nCi=An while ⋃i≥1Ci=A. Countable additivity gives
the last equality by finite additivity. For a decreasing sequence B1⊇B2⊇⋯, apply what has just been proved to the increasing sequence Bic, whose union is Bc, and use step 1 twice:
1−P(B)=P(Bc)=ilimP(Bic)=1−ilimP(Bi).□Definition(Conditional probability, independence)AdaptedAdapted from Grimmett & Stirzaker (2020), Definition (1.4.1), p. 9, Definition (1.5.1) and the following remark, p. 13, and eq. (1.5.5), p. 14; Riley, Hobson & Bence (2006), eqs. (30.19)-(30.22), pp. 1129-1130; Folland (1999), §10.1, the conditional probability measure and the definition of an independent collection, p. 315
Let (Ω,F,P) be a probability space.
Conditional probability. If B is an event with P(B)>0, the probability of AgivenB is
P(A∣B)=P(B)P(A∩B).
For fixed B the map A↦P(A∣B) is itself a probability measure on (Ω,F): it is non-negative, assigns 1 to Ω, and inherits countable additivity from P because (⋃iAi)∩B=⋃i(Ai∩B) preserves disjointness. So every result proved for probability measures holds verbatim for conditional probabilities.
Rearranged, the definition is the chain ruleP(A∩B)=P(B)P(A∣B), which iterates: provided P(A1∩⋯∩An−1)>0,
When P(B)>0 this says exactly P(A∣B)=P(A), which is the meaning; the product form is taken as the definition because it is symmetric in A and B and remains meaningful when P(B)=0.
A family {Ai:i∈I} of events is independent if
P(i∈J⋂Ai)=i∈J∏P(Ai)for every finite J⊆I,
and pairwise independent if this holds for every J of size 2. Pairwise independence is strictly weaker; see the remark below.
Finally, for an event C with P(C)>0, events A and B are conditionally independent given C if P(A∩B∣C)=P(A∣C)P(B∣C), that is, if they are independent under the conditional measure.
Partition. A finite or countable family {Bi} of events is a partition of Ω if the Bi are pairwise disjoint and their union is Ω. Every outcome then lies in exactly one Bi.
Let {B1,B2,…} be a partition of Ω with P(Bi)>0 for every i, and let A be any event.
Partition theorem (law of total probability):
P(A)=i∑P(A∣Bi)P(Bi).
In particular P(A)=P(A∣B)P(B)+P(A∣Bc)P(Bc) whenever 0<P(B)<1.
Bayes' theorem. If also P(A)>0, then for each kP(Bk∣A)=P(A)P(A∣Bk)P(Bk)=∑iP(A∣Bi)P(Bi)P(A∣Bk)P(Bk).
Odds form. For any two events B,C of positive probability and any A with P(A)>0,
P(C∣A)P(B∣A)=P(C)P(B)⋅P(A∣C)P(A∣B),
in which P(A) does not appear. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), eq. (30.29), p. 1133) note that this form is the useful one when only relative probabilities are wanted.
The names attached to the three factors in part 2 are standard: P(Bk) is the prior, P(A∣Bk) the likelihood of the observation A under hypothesis Bk, and P(Bk∣A) the posterior.
1. The sets A∩Bi are pairwise disjoint, because the Bi are, and their union is A∩Ω=A. Countable additivity gives
P(A)=i∑P(A∩Bi).
Each term is rewritten by the chain rule, which is available because P(Bi)>0:
P(A∩Bi)=P(A∣Bi)P(Bi).
Summing gives the claim. For the two-event version take the partition {B,Bc}, both of positive probability by hypothesis.
2. Since P(A)>0 the quantity P(Bk∣A) is defined, and by the definition of conditional probability applied twice to the same intersection,
P(Bk∣A)P(A)=P(A∩Bk)=P(A∣Bk)P(Bk).
Dividing by P(A) gives the first equality. Replacing P(A) by the sum from part 1 gives the second; this is the only place where the family is required to be a partition rather than just a disjoint family.
3. Apply the first equality of part 2 to B and to C and divide the two identities; the common factor P(A) cancels. □
Remark(Where each hypothesis is used)AdaptedAdapted from Grimmett & Stirzaker (2020), the remark after Definition (1.5.1) and Example (1.5.2), pp. 13-14, and the comment on Theorem (1.3.5), p. 7; Folland (1999), §10.1, p. 315
Countable, not finite, additivity. Only the continuity property (consequence 6) uses the full strength of axiom 2; consequences 1-4 need no more than finite additivity. Continuity is what lets a probability be computed as a limit of approximations, and without it the distribution function of a random variable need not be right-continuous and the Borel-Cantelli lemmas have no proof. Convergence of random variables rests on consequence 6 at almost every step.
Positivity of the conditioning event.P(A∣B) is undefined when P(B)=0, and this is not a removable blemish: conditioning on a null event genuinely has no unique answer, and different limiting procedures give different results. Conditional expectation is the construction that handles the case P(B)=0, and it does so by changing the question rather than by patching the quotient. In the partition theorem the hypothesis P(Bi)>0 is needed for each term to be defined; if some P(Bi)=0 the identity still holds with that term read as P(A∩Bi)=0.
Independence is not disjointness. Grimmett and Stirzaker single out the belief that A∩B=∅ makes A and B independent as a common error. It is close to the opposite of the truth: if A and B are disjoint then the occurrence of B rules A out entirely, and P(A∣B)=0. Disjoint events of positive probability are never independent.
Pairwise is weaker than mutual. Let X and Y be the results of two independent fair coin tosses, coded 0 and 1, and let Z=X+Ymod2. Put A={X=1}, B={Y=1}, C={Z=1}, each of probability 21. Any two of the three determine the third, so each pairwise intersection has probability 41=21⋅21 and the three events are pairwise independent. But A∩B∩C=∅, since X=Y=1 forces Z=0, so
P(A∩B∩C)=0=81=P(A)P(B)P(C),
and the family is not independent. This is why the definition quantifies over every finite subfamily and not merely over pairs; Grimmett and Stirzaker (Grimmett & Stirzaker (2020), Example (1.5.2), p. 14) give a different construction on a nine-point sample space, and Folland (Folland (1999), §10.1, Exercise 1, p. 315) makes the same point.
Conditional independence is a separate condition. It neither implies nor is implied by independence. In the example above, A and B are independent but are not conditionally independent given C: conditionally on C, knowing A determines B.
Method(Conditioning on a partition)Original
Most elementary probability calculations are one of the following two, and the work is in choosing the partition.
Identify the experiment. Write down Ω explicitly only if it is small; otherwise it is enough to be confident that it could be written down. Fix the events of interest.
Choose a partition {Bi} that resolves whatever makes the problem hard: the unknown branch taken, the unobserved cause, the value of a hidden parameter, the result of the first step. A good partition is one for which both P(Bi) and P(A∣Bi) are immediate.
Forwards (prediction). Compute P(A)=∑iP(A∣Bi)P(Bi). Check the partition by verifying ∑iP(Bi)=1.
Backwards (inference). If instead A has been observed and the branch is wanted, compute the unnormalised posteriors P(A∣Bi)P(Bi), then divide each by their sum. The sum is P(A), so step 3 is contained in step 4 and is a useful arithmetic check.
Check the direction.P(A∣B) and P(B∣A) are equal only when P(A)=P(B). Confusing them is the single most common error in applied probability; it has a name, the prosecutor's fallacy.
For a sequence of steps, use the chain rule rather than a partition: multiply the conditional probability of each step given everything before it.
Example(A three-module detector: forwards and backwards)Original
A tracking detector routes each incoming event to exactly one of three reconstruction modules. Module 1 handles a fraction 0.50 of the events, module 2 a fraction 0.30, and module 3 a fraction 0.20. Each module occasionally misreconstructs an event; from calibration runs the misreconstruction probabilities are
P(E∣M1)=0.02,P(E∣M2)=0.05,P(E∣M3)=0.10,
where Mi is the event that module i handled the event and E is the event that it was misreconstructed. The modules partition the sample space: P(M1)+P(M2)+P(M3)=0.50+0.30+0.20=1.
Forwards: the overall error rate. By the partition theorem,
Backwards: which module is to blame? The three unnormalised posteriors are the three products just computed, 0.010, 0.015 and 0.020, and dividing each by their sum 0.045 gives
The three posteriors sum to 2/9+3/9+4/9=1, as they must.
What the numbers say. Module 3 handles only one event in five but is responsible for 4/9 of the bad ones, while module 1 handles half the events and causes only 2/9 of them. The prior ordering 0.50>0.30>0.20 has been reversed by the likelihoods. Note also the gap between the two conditional probabilities that are easy to confuse: P(E∣M3)=0.10, but P(M3∣E)=0.444, larger by a factor of about 4.4.
The odds form, as a check. Comparing modules 3 and 1 without computing P(E) at all:
and indeed (4/9)/(2/9)=2. The prior odds of 0.4 against module 3 are overturned by a likelihood ratio of 5 in its favour.
A second question, by the chain rule. Two events arrive and are routed independently. What is the probability that both are handled by module 1 and that neither is misreconstructed? Independence across events, and the chain rule within each event, give
[P(M1)P(Ec∣M1)]2=[(0.50)(0.98)]2=0.492=0.2401.Application(Where this machinery is actually used)Original
Inference from data. Bayes' theorem is the engine of Bayesian inference: the partition is over hypotheses or parameter values, the likelihood is supplied by the measurement model, and the posterior is the output. The detector example above is the discrete skeleton of every such calculation, including the medical-screening example that Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §30.2.3, p. 1132) work through, in which a test that is 99.99% sensitive still leaves a positive result only one-third likely to be a true positive because the prior prevalence is low. Misreading P(positive∣ill) as P(ill∣positive) is the same error that Grimmett and Stirzaker label the prosecutor's fallacy.
Reliability and fault trees. A system is modelled as a network whose components fail independently; the probability that the system works is computed by inclusion-exclusion over the paths, or by conditioning on the state of one critical component. Consequence 4 and the partition theorem are the whole method, and independence of components is the modelling assumption that makes the numbers finite in count.
Statistical mechanics. The microcanonical ensemble is the equally-likely-outcomes measure P(A)=∣A∣/∣Ω∣ on the set of accessible microstates, and the counting of states that makes up much of the subject is exactly the combinatorics of such a space. Conditioning on the total energy of a system plus reservoir, and applying the partition theorem over the energy of the subsystem, is how the canonical distribution is derived.
Sequential experiments and Markov models. The chain rule is the definition of a stochastic process specified step by step: the probability of a whole trajectory is the product of one-step conditional probabilities. When each step depends only on the previous state, the chain rule collapses to the Markov property and one is in the setting of Discrete-time Markov chains. The product-space construction that Grimmett and Stirzaker sketch is what guarantees that a probability space carrying infinitely many independent trials exists at all, which is what Laws of large numbers silently assumes.
Conditional independence as a modelling language. Measurements are usually modelled as independent given the state of the system, not independent outright; the state is the common cause that induces the correlation. This single conditional-independence assumption underlies the likelihood function of Maximum likelihood, the recursion of the Kalman filter, and the graphical models used throughout data assimilation.
Defines a random variable as a measurable function on a probability space and its distribution function as the object that carries all the probabilistic information about it, characterises which functions are distribution functions, separates the discrete and continuous cases, and catalogues the standard distributions with the mechanism each one models.
Note(Why a distribution function and not a list of probabilities)Original
A random variable is just a numerical readout of an experiment: the sample space Ω may be a set of card hands, molecular configurations or detector hit patterns, but what gets recorded is a number. The obvious way to describe the number is to list P(X=x) for each possible value x, and for a die or a photon count that is exactly right. For a length, a time or a temperature it is useless, because every single value has probability zero.
The device that covers both cases is the distribution function F(x)=P(X≤x), an increasing function from 0 to 1. It exists for every random variable, it determines all the probabilities of interest, and it is a function on R rather than on Ω — so once it is written down the original probability space can be forgotten. Two experiments with entirely different sample spaces can share a distribution function, and for most purposes they are then the same experiment.
Definition(Random variable, distribution function, distribution)AdaptedAdapted from Grimmett & Stirzaker (2020), §2.1, Definitions (2.1.3) and (2.1.4) and the discussion preceding them, pp. 29-30; Folland (1999), §10.1, pp. 314-315
such an X is called F-measurable. The condition is exactly what is needed for the probabilities below to be defined; it is automatic when F consists of all subsets of Ω, as it does for a finite or countable sample space.
Following the usual abuse, the event {ω:X(ω)≤x} is written {X≤x} and its probability P(X≤x).
The distribution function of X is
F:R→[0,1],F(x)=P(X≤x),
written FX when the variable must be named. Measure-theoretically, X pushes P forward to a Borel probability measure on the line,
PX(B)=P(X−1(B))(B⊆R Borel),
called the distribution of X; and F determines PX completely. Two random variables with the same distribution function are identically distributed, whether or not they are defined on the same space, and no purely distributional statement can tell them apart.
The two tails of X are T1(x)=P(X>x)=1−F(x) and T2(x)=P(X≤−x)=F(−x) for large positive x; how fast they decay governs which moments exist.
Definition(Discrete and continuous random variables)AdaptedAdapted from Grimmett & Stirzaker (2020), §2.3, Definitions (2.3.1) and (2.3.2) and the surrounding discussion, p. 35, and §4.1, Definition (4.1.2), p. 99; Riley, Hobson & Bence (2006), §30.4.1-§30.4.2, eqs. (30.38)-(30.44), pp. 1139-1142
X is discrete when its range lies inside a countable subset {x1,x2,…} of the line. Its (probability) mass function is
f:R→[0,1],f(x)=P(X=x),
which vanishes off the countable set and satisfies ∑if(xi)=1. Then F(x)=∑xi≤xf(xi) is a step function, constant between the xi and jumping by f(xi) at xi; such a distribution is called atomic, and the points xi with f(xi)>0 are its atoms. Riley, Hobson and Bence call f the probability function and F the cumulative probability function.
X is continuous when some integrable f:R→[0,∞), called a (probability) density function for X, reproduces F by integration:
F(x)=∫−∞xf(u)duat every x∈R.
Then F is continuous, P(X=x)=0 for every x, and
P(a<X≤b)=∫abf(u)du.
The density is not unique — changing f at finitely many points changes no integral — but wherever F is differentiable one takes f(x)=F′(x). Physicists usually define the density the other way round, by P(x<X≤x+dx)=f(x)dx, which is the same statement in infinitesimal form.
Two cautions. First, the word "continuous" describes F, not the function X on Ω. Second, these two classes are not exhaustive: there exist singular distribution functions, continuous but with no density, and there exist mixtures of all three kinds. A convenient device for handling discrete and continuous variables in one formalism, used throughout physics, is to allow densities containing Dirac deltas, writing f(x)=∑ipiδ(x−xi) for a discrete variable; see The Dirac delta.
Definition(Which functions are distribution functions, and what they tell you)AdaptedAdapted from Grimmett & Stirzaker (2020), §2.1, Lemma (2.1.6) with the remark following it on p. 31, and Lemma (2.1.11), p. 32
Let X be a random variable with distribution function F.
F is non-decreasing: x<y implies F(x)≤F(y).
F(x)→0 as x→−∞ and F(x)→1 as x→+∞.
F is right-continuous: F(x+h)→F(x) as h↓0.
Conversely, any F:R→[0,1] having properties 1-3 arises as FY for some random variable Y on some probability space, constructed in step 5 of the proof. So these three properties characterise distribution functions exactly.
Moreover, writing F(x−)=limy↑xF(y) for the left limit, which exists by monotonicity,
In particular F is continuous at x if and only if P(X=x)=0, so the atoms of X are precisely the jumps of F, and there are at most countably many of them because the jumps at distinct points sum to at most 1.
Idea of the proof(The three properties are the three continuity statements)Original
Each property of F is the continuity theorem for probability measures applied to a different monotone sequence of events, and the side on which F is continuous is decided by the choice of ≤ rather than < in the definition.
Monotonicity needs no limit: {X≤x} sits inside {X≤y} when x<y, and a probability measure is monotone. The two limits at ±∞ come from the sequences {X≤−n}, which decrease to ∅, and {X≤n}, which increase to Ω.
Right-continuity is the same argument at a finite point: the events {X≤x+1/n} decrease, and their intersection is {X≤x} — because a real number that is ≤x+1/n for every n is ≤x. Had we defined F(x)=P(X<x), the corresponding intersection would have been {X≤x} again and F would have come out left-continuous instead. The mass at a point is then the gap between the two conventions, P(X=x)=P(X≤x)−P(X<x)=F(x)−F(x−).
Proof(Properties of the distribution function)AdaptedAdapted from Grimmett & Stirzaker (2020), §2.1, the proofs of Lemma (2.1.6) and Lemma (2.1.11), pp. 30-32
Write A(x)={X≤x}, an event by the measurability requirement in the definition.
1. Monotone. If x<y then A(y)=A(x)∪{x<X≤y}, a disjoint union of events, so by additivity
F(y)=F(x)+P(x<X≤y)≥F(x).
This also proves the second displayed formula of the statement.
2. Limits at ±∞. The events Bn=A(−n) decrease as n grows, and ⋂nBn=∅, since no real number X(ω) satisfies X(ω)≤−n for every n. By continuity of P along decreasing sequences (Probability spaces, consequence 6), F(−n)=P(Bn)→P(∅)=0. As F is monotone, convergence along the integers forces F(x)→0 as x→−∞. Symmetrically, A(n) increases with union Ω, so F(n)→1 and hence F(x)→1 as x→+∞.
3. Right-continuity. The events Cn=A(x+1/n) decrease, and ⋂nCn=A(x): if X(ω)≤x+1/n for all n then X(ω)≤x, and the reverse inclusion is clear. Continuity along decreasing sequences gives F(x+1/n)→F(x), and monotonicity upgrades this to F(x+h)→F(x) as h↓0.
4. Mass at a point. Let Dn={x−1/n<X≤x}. These decrease with ⋂nDn={X=x}. By part 1, P(Dn)=F(x)−F(x−1/n), and continuity along decreasing sequences gives
P(X=x)=n→∞lim[F(x)−F(x−1/n)]=F(x)−F(x−).
The first displayed formula, P(X>x)=1−F(x), is the complement rule.
5. The converse. Given F with properties 1-3, take Ω=(0,1) with Lebesgue measure and set
X(ω)=inf{x∈R:F(x)≥ω},
the generalised inverse of F. Right-continuity of F makes the infimum attained, so F(X(ω))≥ω; and monotonicity makes X(ω)≤x equivalent to ω≤F(x). Hence P(X≤x)=P((0,F(x)])=F(x), and property 2 guarantees that X is finite for every ω∈(0,1). □
Remark(Where each hypothesis is used, and what goes wrong without it)AdaptedAdapted from Grimmett & Stirzaker (2020), §2.1, the discussion of measurability on p. 29, §2.3, the remarks on singular and mixed variables and Example (2.3.5), pp. 35-36
Measurability. Nothing above works if {X≤x} is not an event: F(x) is then simply undefined. On a countable sample space with F the full power set the condition is vacuous, which is why it can be ignored in elementary work. It becomes real on Ω=R with the Borel σ-field, where non-measurable functions exist (see Measurable functions).
Countable additivity. Properties 2 and 3 are the continuity theorem, which as noted in Probability spaces is equivalent to countable rather than finite additivity. A finitely additive set function can produce a monotone "distribution function" that fails to tend to 1, or that jumps without having an atom.
Right- rather than left-continuity is a convention, forced by writing F(x)=P(X≤x). Some applied texts define F(x)=P(X<x) and obtain a left-continuous function. The two differ exactly at the atoms, and every formula involving F at a jump has to be read in the right convention; we use ≤ throughout, as do all four of our sources.
Discrete and continuous do not exhaust the possibilities. The worked example below is a variable that is neither: it has a density on part of the line and an atom at one point. Grimmett and Stirzaker give a similar example in which a coin toss decides between a point mass at −1 and a uniform angle. There are also singular distributions, continuous but with F′=0 almost everywhere; the standard construction is supported on the Cantor set. When a theorem below is stated for continuous variables, the hypothesis is that F is an integral of a density, not merely that F is a continuous function.
Non-uniqueness of the density. Because f enters only through its integrals, altering f on a set of Lebesgue measure zero changes nothing. So "the" density is really an equivalence class, and statements such as "f is largest at the mode" require choosing the continuous representative when one exists.
Note(A catalogue of the standard distributions)AdaptedAdapted from Grimmett & Stirzaker (2020), §3.5, Examples (1)-(6), pp. 64-65, and §4.4, Examples (1)-(8), pp. 106-108
In each row the parameters are as shown and the mass or density function is zero outside the stated range. The mean and variance are included for reference; they are defined, and the method for computing them given, in Expectation, variance and covariance and Generating functions.
Discrete.
Distribution
f(k)=P(X=k)
Range
Mean
Variance
Mechanism
Bernoulli (p)
pk(1−p)1−k
k∈{0,1}
p
p(1−p)
one trial; the indicator of an event
Binomial bin(n,p)
(kn)pk(1−p)n−k
0≤k≤n
np
np(1−p)
successes in n independent trials
Poisson (λ)
e−λλk/k!
k≥0
λ
λ
rare events: n→∞, p→0, np→λ
Geometric (p)
p(1−p)k−1
k≥1
1/p
(1−p)/p2
trial number of the first success
Negative binomial (r,p)
(r−1k−1)pr(1−p)k−r
k≥r
r/p
r(1−p)/p2
trial number of the rth success
Continuous.
Distribution
f(x)
Range
Mean
Variance
Mechanism
Uniform [a,b]
1/(b−a)
a<x<b
(a+b)/2
(b−a)2/12
no preference within an interval
Exponential (λ)
λe−λx
x≥0
1/λ
1/λ2
waiting time with no memory
Normal N(μ,σ2)
2πσ21exp[−2σ2(x−μ)2]
x∈R
μ
σ2
sums of many small independent effects
Gamma Γ(λ,t)
λtxt−1e−λx/Γ(t)
x≥0
t/λ
t/λ2
waiting time for the tth event
Chi-squared χ2(d)
Γ(21,2d) density
x≥0
d
2d
sum of d squared standard normals
Cauchy
1/[π(1+x2)]
x∈R
none
none
ratio of two independent normals
Beta (a,b)
xa−1(1−x)b−1Γ(a+b)/[Γ(a)Γ(b)]
0<x<1
a/(a+b)
(a+b)2(a+b+1)ab
a random proportion
Three conventions to watch. The geometric distribution is here the number of trials up to and including the first success, so it starts at k=1 and has mean 1/p; many texts instead count the failures before the first success, shifting everything by one and giving mean (1−p)/p. Grimmett and Stirzaker parametrise the gamma distribution as Γ(λ,t) with λ a rate, so that Γ(λ,1) is the exponential and χ2(d)=Γ(21,2d); writers who use a scale parameter θ=1/λ invert it. The Cauchy distribution has no mean, not a mean of zero: the defining integral is not absolutely convergent.
For the standard normal the density and distribution function have their own symbols,
ϕ(x)=2π1e−x2/2,Φ(y)=∫−∞yϕ(v)dv,
and if X is N(μ,σ2) then (X−μ)/σ is N(0,1), so that FX(x)=Φ((x−μ)/σ). This standardisation is why one table of Φ serves for all normal distributions.
Example(A clipped decay time: a variable that is neither discrete nor continuous)Original
A counter records the time T at which an unstable state decays, but the acquisition window closes after c=3 hours; if no decay has occurred by then the counter writes down 3 and stops. Model the true decay time as exponential with rate λ=21 per hour, and let X=min(T,3) be the recorded value.
The distribution function. For x<0, F(x)=0. For 0≤x<3 the recorded value is at most x exactly when the true decay time is, so
F(x)=P(T≤x)=1−e−x/2.
For x≥3, F(x)=1, because X never exceeds 3. The function is non-decreasing, tends to 0 and 1 at the ends, and is right-continuous — at x=3 it jumps up to the value 1 and stays there, which is continuity from the right, not from the left.
The atom. By the formula P(X=x)=F(x)−F(x−),
P(X=3)=1−(1−e−3/2)=e−3/2=0.2231301…,
which is just the probability that the state survives the whole window. At every other point F is continuous, so X has exactly one atom and X is not continuous. Nor is it discrete, since the set of values it takes, (0,3], is uncountable and F is not a step function. It is a mixture: a density f(x)=21e−x/2 on (0,3), of total mass 1−e−3/2=0.7768698…, together with a point mass of 0.2231301… at 3. The two add to 1, as they must.
Interval probabilities. Between one and two hours,
Quantiles, and why the generalised inverse is needed. The median solves F(m)=21, that is e−m/2=21, giving
m=2ln2=1.3862944 hours,
comfortably inside the window; the upper quartile solves e−q/2=41, giving q=4ln2=2.7725887, still inside. But the 90th percentile would need e−x/2=0.1, that is x=2ln10=4.6051702, which lies beyond the window. There is no x with F(x)=0.9: the function leaps from 0.7768698 to 1 at x=3. The generalised inverse inf{x:F(x)≥0.9} used in step 5 of the proof returns 3, which is the right answer — nine runs in ten do finish by the time the counter stops, because in fact more than nine in ten do.
What the clipping costs. The recorded distribution agrees with the true one below 3 hours and lumps everything above into a single spike. Any analysis that treats X as exponential will therefore underestimate the mean lifetime; the correct treatment keeps the atom, and this is the elementary case of what survival analysis calls right-censoring.
Application(Distributions as models)Original
Counting. Radioactive decays in a fixed interval, photons on a pixel, cosmic rays through a detector and raindrops in a gauge are all modelled as Poisson, for the reason in the catalogue: a very large number of potential contributors each with a very small probability. The quality of the model is a statement about independence, not about the counts, and when the contributors cluster the data are overdispersed relative to Poisson, which is itself a diagnostic.
Waiting and lifetimes. The exponential distribution is the unique continuous distribution with no memory, and is therefore the universal model for the time to the next event of a process with a constant hazard rate: spontaneous emission, nuclear decay, the time between arrivals in a queue. The gamma distribution is the waiting time for the tth such event, which is why it describes the total dead time of a counter that must record several hits.
Measurement error. Gaussian errors are assumed throughout experimental science, on the strength of The central limit theorem: an error that is the accumulation of many small independent contributions is approximately normal regardless of their individual laws. The assumption fails when one contribution dominates, or when the tails matter, and the Cauchy distribution — the ratio of two normals, with no mean at all — is the standard warning of how badly it can fail. The Lorentzian line shape of a resonance is the same function.
Quantum mechanics. A wavefunction is a probability density in exactly the sense of this unit: Arfken, Weber and Harris (Arfken, Weber & Harris (2013), Example 23.2.2, p. 1135) treat the radial position of the hydrogen 1s electron as a continuous random variable with density 4r2a−3e−2r/a, obtained from ∣ψ∣2 by integrating over angles. Normalisation of the wavefunction is the condition ∫f=1, and the Born rule is the assertion that ∣ψ∣2 is a density and not merely a weight.
Simulation. The generalised inverse constructed in step 5 of the proof is the inverse-transform sampling method: feed a uniform random number on (0,1) into F−1 and out comes a sample from F. It works for any distribution function, mixed ones included, and it is how the clipped variable of the worked example would be simulated. More efficient alternatives for specific distributions are in Joint distributions and transformations.
Defines the expectation of a random variable and the three derived numbers that summarise a distribution and the link between two of them, proves the rules that make expectation linear and variance quadratic, and establishes the Cauchy-Schwarz bound that confines the correlation coefficient to the interval from minus one to one.
A distribution function contains everything, and is therefore hard to compare, hard to communicate and hard to estimate from data. In practice a distribution is summarised by a few numbers: where it sits, how wide it is, and — for two variables at once — how much one moves with the other.
The first is the expectation, the centre of mass of the distribution. The second is the variance, its moment of inertia about that centre. The third, covariance, is the mixed second moment of a pair, and the fourth, correlation, is the covariance rescaled so that it lies in [−1,1] and can be quoted without units. The rescaling is not cosmetic: a covariance of 10−6 may mean total dependence or none at all, depending on the units of the two variables, and a correlation of −0.9 means the same thing whatever they are measured in.
Everything in this unit follows from two facts: expectation is linear, and E(Z2)≥0.
If X is discrete with mass function f, its expectation (mean, expected value) is
E(X)=x∑xf(x),
the sum running over the countably many x with f(x)>0, provided the sum converges absolutely. If X is continuous with density f,
E(X)=∫−∞∞xf(x)dx,
provided ∫∣x∣f(x)dx<∞. When the stated condition fails, X has no mean.
Both are instances of one definition. The distribution PX of X is a measure on R, and
E(X)=∫RtdPX(t)=∫ΩXdP,
a Lebesgue integral (The Lebesgue integral); expectation is integration, and in Folland's dictionary the two words name the same operation. Grimmett and Stirzaker write this uniformly as E(X)=∫xdFX(x), which also covers variables that are neither discrete nor continuous. We write μ or μX for E(X) where it is convenient, and omit the parentheses as EX.
Moments. For a positive integer k, the kth moment of X is mk=E(Xk) and the kth central moment is σk=E((X−m1)k), each defined when the corresponding expectation exists.
Definition(Variance, covariance, correlation)AdaptedAdapted from Grimmett & Stirzaker (2020), Definition (3.3.5) and the computation of $\\sigma_2$, p. 54, Definition (3.3.10), p. 56, and Definition (3.6.7), p. 68; Riley, Hobson & Bence (2006), eq. (30.48) and §30.5.3, p. 1146, and eqs. (30.133)-(30.134), p. 1200; Arfken, Weber & Harris (2013), eqs. (23.26)-(23.27) and (23.36), pp. 1139, 1142; Folland (1999), §10.1, p. 314
Let X and Y be random variables with finite second moments, that is E(X2)<∞ and E(Y2)<∞; equivalently X,Y∈L2, in the language of Lp spaces. Write μX=EX and μY=EY.
The variance of X is its second central moment,
var(X)=E((X−μX)2)≥0,
and the standard deviation is σX=var(X), which carries the same units as X. A variance is never negative.
The covariance of X and Y is their mixed central second moment,
cov(X,Y)=E((X−μX)(Y−μY)),
so that cov(X,X)=var(X): covariance generalises variance. X and Y are uncorrelated when cov(X,Y)=0, equivalently when E(XY)=E(X)E(Y).
When σX and σY are both non-zero, the correlation coefficient is
ρ(X,Y)=σXσYcov(X,Y),
a dimensionless number. If ρ>0 the variables are positively correlated, if ρ<0 negatively correlated, and if ρ=0 uncorrelated.
Definition(The rules)AdaptedAdapted from Grimmett & Stirzaker (2020), Lemma (3.3.3), p. 54, Theorems (3.3.8) and (3.3.11) and Lemma (3.3.9), pp. 55-56, Lemma (3.6.6), p. 68, Lemma (3.6.8) and Theorem (3.6.9), p. 69, Theorem (4.3.3) and Lemma (4.3.4), pp. 103-104; Riley, Hobson & Bence (2006), the properties (i)-(iii) on p. 1144, eq. (30.51), p. 1147, and eq. (30.135), p. 1201; Arfken, Weber & Harris (2013), property 2 on p. 1139 and the theorem on p. 1143
Let X and Y be random variables and a,b,c∈R; each statement is asserted whenever the expectations it names exist.
Law of the unconscious statistician. For g:R→R,
E(g(X))=x∑g(x)f(x)(discrete),E(g(X))=∫−∞∞g(x)f(x)dx(continuous),
provided the sum converges absolutely or the integral exists. The density of g(X) is never needed. For two variables with joint mass function or density fX,Y and g:R2→R, the same formula holds with fX,Y in place of f and a double sum or integral. Grimmett and Stirzaker call this the change of variable formula and record that it is also known by the longer name used here.
Positive linear operator.E(1)=1; if X≥0 then E(X)≥0; and
E(aX+bY)=aE(X)+bE(Y).
No independence is required. In particular E(a)=a for a constant a.
Products under independence. If X and Y are independent then E(XY)=E(X)E(Y), so independent variables are uncorrelated. The converse is false; see the remark.
Behaviour under a linear change.var(a)=0 for a constant, and
var(aX+b)=a2var(X),cov(aX+b,cY+d)=accov(X,Y),ρ(aX+b,cY+d)=sgn(ac)ρ(X,Y)
for a,c=0. Variance is quadratic and blind to shifts; correlation is invariant up to the sign of the scalings.
Variance of a sum.var(X+Y)=var(X)+var(Y)+2cov(X,Y), which reduces to var(X)+var(Y) when X and Y are uncorrelated — in particular when they are independent.
Cauchy-Schwarz, and the range of ρ.(E(XY))2≤E(X2)E(Y2), with equality if and only if P(aX=bY)=1 for some real a,b not both zero. Consequently
∣cov(X,Y)∣≤σXσY,so−1≤ρ(X,Y)≤1,
and ∣ρ(X,Y)∣=1 holds exactly when P(aX+bY=c)=1 for some a,b,c with a,b not both zero: that is, exactly when Y is almost surely an affine function of X. The sign of ρ is then the sign of the slope.
The mean minimises mean squared error.a↦E((X−a)2) is a quadratic in a with its minimum at a=E(X), and the minimum value is var(X).
Tail formula. If X≥0 with distribution function F, then E(X)=∫0∞[1−F(x)]dx.
Idea of the proof(Linearity, then one square)Original
Only two mechanisms are at work.
Linearity comes from the fact that a sum is integrated term by term; in the discrete case one writes both variables as combinations of indicators of the joint events {X=x}∩{Y=y} and reorganises a double sum. Everything about means follows, and so does the computing formula for the variance: expand (X−μ)2 and apply linearity to the three terms, two of which collapse because E(X)=μ is a constant. Rules 4, 5, 6 and 8 are this expansion carried out in slightly different arrangements.
Positivity of a square does the rest. The quantity E((aU+cV)2) is non-negative for every pair of real coefficients, and expanding it gives a quadratic form in (a,c) whose coefficients are E(U2), E(UV) and E(V2). A quadratic that never goes negative has non-positive discriminant, and that inequality is Cauchy-Schwarz. Applying it to the centred variables U=X−μX, V=Y−μY bounds the covariance by the product of the standard deviations, which is what confines ρ to [−1,1]. Equality forces the quadratic to have a real root, that is, to vanish for some (a,c)=(0,0); and a non-negative random variable with zero mean is zero almost surely, so aU+cV=0 with probability one. That is the linear relation in the equality case.
Independence enters only in rule 3, and only through the factorisation of the joint mass function or density, which turns a double sum into a product of single sums.
Proof(The rules)AdaptedAdapted from Grimmett & Stirzaker (2020), the proofs of Theorems (3.3.8) and (3.3.11), Lemma (3.3.9), pp. 55-56, Theorem (3.6.9), p. 69, and Lemma (4.3.4), p. 104; Arfken, Weber & Harris (2013), the quadratic-form argument on p. 1143; Riley, Hobson & Bence (2006), the expansion (30.135), p. 1201, and the proof that independence gives $E[XY]=E[X]E[Y]$, p. 1200
We give the discrete case where the two cases differ only by replacing sums with integrals.
Rule 1 (law of the unconscious statistician). Let Y=g(X) and let y run over the values of Y. Grouping the values of X by the value they are sent to,
the regrouping being legitimate because the double sum converges absolutely. In the continuous case the corresponding statement for g≥0 follows from rule 9 applied to g(X) together with an exchange of the order of integration, and the general case follows by splitting g=g+−g− into positive and negative parts.
Rule 2 (linearity). Put Ax={X=x} and By={Y=y}. The events Ax∩By are disjoint and exhaust Ω, and on Ax∩By the variable aX+bY takes the value ax+by. Hence by rule 1 applied to the pair,
The inner sums are the marginals: ∑yP(Ax∩By)=P(Ax), since the By partition Ω, and similarly in the other order. This is aE(X)+bE(Y). Positivity and E(1)=1 are immediate from the definition.
Rule 3 (products). Independence of X and Y means P(Ax∩By)=P(Ax)P(By) for all x,y, so the double sum factorises:
This is the only step in the whole unit that uses independence.
Rule 4 (computing formulas). By linearity, with μ=EX a constant,
E((X−μ)2)=E(X2)−2μE(X)+μ2=E(X2)−μ2.
The same expansion on (X−μX)(Y−μY) gives E(XY)−μXE(Y)−μYE(X)+μXμY=E(XY)−μXμY.
Rule 5 (linear change).E(aX+b)=aμX+b by linearity, so (aX+b)−E(aX+b)=a(X−μX): the shift cancels. Squaring and taking expectations gives var(aX+b)=a2var(X), and multiplying the two centred variables gives the covariance rule. Since σaX+b=∣a∣σX, dividing yields ρ(aX+b,cY+d)=acρ(X,Y)/(∣a∣∣c∣)=sgn(ac)ρ(X,Y).
Rule 6 (variance of a sum). Centre both variables and expand the square:
Rule 7 (Cauchy-Schwarz). Take any U,V in L2 and any real coefficients a,c. Since (aU+cV)2≥0 pointwise, positivity and linearity of E give
Q(a,c)=a2E(U2)+2acE(UV)+c2E(V2)≥0for all a,c.
Suppose first E(U2)>0. Fixing c=1 makes Q a quadratic in a with positive leading coefficient that is never negative, so its discriminant is non-positive:
4(E(UV))2−4E(U2)E(V2)≤0,
which is the inequality. If instead E(U2)=0 then U=0 almost surely, both sides vanish, and the equality case holds with (a,c)=(1,0).
For the equality case, suppose E(U2)>0 and the discriminant is zero. Then the quadratic has a real root a0, so E((a0U+V)2)=0. A non-negative random variable with zero expectation is zero almost surely, so P(a0U+V=0)=1. Conversely a linear relation makes Q vanish at that point, forcing the discriminant to zero.
The bound on ρ. Apply the inequality to the centred variables U=X−μX and V=Y−μY: then E(UV)=cov(X,Y), E(U2)=var(X) and E(V2)=var(Y), so cov(X,Y)2≤var(X)var(Y) and dividing by σX2σY2>0 gives ρ2≤1. Equality holds exactly when a0(X−μX)+(Y−μY)=0 almost surely, that is when P(aX+bY=c)=1 for suitable constants; solving for Y shows Y=−a0X+(μY+a0μX) almost surely, an affine function of X with slope −a0, and rule 5 gives ρ=sgn(−a0).
Rule 8 (the mean minimises). By rule 4 applied to X−a,
E((X−a)2)=var(X−a)+(E(X−a))2=var(X)+(μX−a)2,
a quadratic in a minimised at a=μX with value var(X).
the exchange of the order of integration being legitimate by Tonelli's theorem (Fubini and Tonelli) because the integrand is non-negative. □
Remark(Where the hypotheses are used)AdaptedAdapted from Grimmett & Stirzaker (2020), the absolute-convergence proviso in Definition (3.3.1) and Example (3.3.12), pp. 53, 57, the remark after Theorem (3.3.8), p. 56, and the Cauchy density, p. 108; Riley, Hobson & Bence (2006), the statement that uncorrelated does not imply independent and the symmetry criterion, p. 1201; Folland (1999), §10.1, p. 314
Absolute convergence, not convergence. The mean of a discrete variable is a sum over an unordered set, so it must not depend on the order of summation; this is why absolute convergence is demanded rather than convergence. Grimmett and Stirzaker's example f(k)=Ak−2 for k=±1,±2,… has ∑k∣k∣f(k)=A∑k=0∣k∣−1=∞, so no mean exists even though the distribution is symmetric and one is tempted to call its mean zero. The Cauchy density 1/[π(1+x2)] is the continuous counterpart: symmetric, and with no mean at all.
Finite second moment, for everything after the definition of the mean. Variance, covariance, correlation and the Cauchy-Schwarz inequality all need X,Y∈L2. Folland states the convention that σ2(X)=∞ when X∈/L2. If var(X)=∞ then rule 6 is vacuous and ρ is undefined; the first place this bites is the central limit theorem, whose finite-variance hypothesis cannot be dropped.
Non-degeneracy, for the correlation coefficient.ρ requires σX,σY>0. A constant random variable is uncorrelated with everything and has no correlation coefficient with anything.
Independence is used only in rule 3, and only in one direction. Uncorrelated does not imply independent. Let X be uniform on {−1,0,1} and set Y=X2. Then EX=0 and E(XY)=E(X3)=31(−1+0+1)=0, so cov(X,Y)=0−0⋅EY=0 and the two are uncorrelated. They are as dependent as they could be — Y is a function of X — and indeed
P(X=0,Y=1)=0=31⋅32=P(X=0)P(Y=1).
What has happened is that correlation measures only the linear part of a dependence, and the dependence here is purely quadratic. Riley, Hobson and Bence give the general reason: any joint distribution symmetric in x about μX has zero correlation, whatever else it does. The one important case in which uncorrelated does imply independent is the joint normal, treated in The multivariate normal distribution.
Rule 6 needs only uncorrelatedness, not independence, which is why variances of sums are so widely usable: pairwise uncorrelated errors suffice, and nothing about higher-order dependence matters.
The exchange of integrals in rule 9 needs the integrand to be non-negative, which is why that rule is stated for X≥0. For a general X the correct version splits into positive and negative parts and both must converge.
Method(Computing a variance or covariance without losing it to rounding)Original
Centre first if you can. Work with U=X−m for any convenient constant m near the mean; by rule 5, var(X)=var(U) and cov(X,Y)=cov(U,Y−m′), so shifting is free and costs no accuracy.
Then use whichever of the two formulas is cheaper.E(U2)−(EU)2 needs one pass over the distribution; E((U−μU)2) needs the mean first, hence two. With exact arithmetic they agree.
Beware the difference of two large numbers. When σX≪∣μX∣ the formula E(X2)−(EX)2 subtracts two nearly equal quantities and loses most of its significant digits; this is precisely what step 1 prevents, and it is the reason numerical libraries subtract a provisional mean before accumulating squares.
For a function of X, do not find its distribution. Apply rule 1 directly: E(g(X))=∑xg(x)f(x) or ∫g(x)f(x)dx. Deriving the density of g(X) first, by the methods of Joint distributions and transformations, is extra work and extra opportunity for error.
For a linear combination, use rules 2, 5 and 6 rather than constructing the distribution of the combination: var(∑iaiXi)=∑iai2var(Xi)+2∑i<jaiajcov(Xi,Xj).
Check the answer. A variance must be non-negative; a correlation must lie in [−1,1]; and var(X+Y) computed from rule 6 must agree with a direct computation from the distribution of X+Y. Each of these catches a different arithmetic slip.
Example(Two detector channels: a complete calculation)Original
Two photon-counting channels A and B view the same weak source through a beam splitter. In a fixed 1ms gate, let X be the number of counts in channel A, taking values 0,1,2, and Y the number in channel B, taking values 0,1. Calibration gives the joint mass function (all entries in twentieths):
20f(x,y)
y=0
y=1
20fX(x)
x=0
4
1
5
x=1
3
5
8
x=2
2
5
7
20fY(y)
9
11
20
The six entries sum to 20, so f is a mass function; the margins are the row and column sums.
Means. By the definition applied to the marginals,
comfortably inside [−1,1], as rule 7 demands. The two channels are moderately, not strongly, coupled: a little over a third of the way to a deterministic linear relation.
Check by rule 6. The predicted variance of the total count S=X+Y is
var(S)=0.59+0.2475+2(0.145)=1.1275.
Computing it directly, S takes the values 0,1,2,3 with probabilities 4/20, 4/20, 7/20, 5/20, so E(S)=(0+4+14+15)/20=33/20=1.65 — which equals 1.1+0.55, as linearity requires — and E(S2)=(0+4+28+45)/20=77/20=3.85, giving var(S)=3.85−1.652=3.85−2.7225=1.1275. The two agree exactly.
What would independence have given? Under independence the covariance would vanish, so var(S) would be 0.8375 rather than 1.1275: the positive coupling inflates the variance of the total by 35%. Treating correlated channels as independent is the standard way to underestimate an error bar.
The best linear predictor. The slope of the least-squares line predicting Y from X is cov(X,Y)/var(X)=0.145/0.59=29/118=0.245763, and the residual variance is var(Y)(1−ρ2)=0.2475(1−0.143982)=0.211864. Knowing X removes about 14.4% of the variance of Y; the exact conditional prediction, which does better, is the subject of Conditional expectation.
Application(Where the second moment is the answer)Original
Error bars and propagation. The standard deviation of an estimate is its error bar, and rules 5 and 6 are the whole of linear error propagation: for Z=∑iaiXi,
σZ2=i∑ai2σi2+2i<j∑aiajcov(Xi,Xj).
The cross terms are the part routinely and wrongly dropped; see Error propagation. For a non-linear Z=g(X1,…,Xn) the same formula is applied to the first-order Taylor expansion, with ai=∂g/∂xi evaluated at the means.
Why averaging helps, and exactly how much. If X1,…,Xn are uncorrelated with common variance σ2, rules 5 and 6 give var(Xˉ)=σ2/n, so the standard error falls like n−1/2. This single calculation is the quantitative content of Laws of large numbers at finite n, the basis of the n−1/2 convergence rate of Monte Carlo integration, and the reason correlated samples — from a Markov chain, say — buy less than their number suggests: the cross terms in rule 6 do not vanish, and the effective sample size is smaller.
Fluctuations in statistical mechanics. In the canonical ensemble the energy E is a random variable with E(E)=−∂lnZ/∂β and var(E)=∂2lnZ/∂β2, where β=1/kBT. Differentiating the mean with respect to temperature gives
CV=∂T∂E(E)=kBT2var(E),
so a heat capacity is literally a variance. The same identity in other variables relates the isothermal compressibility to particle-number fluctuations and the magnetic susceptibility to magnetisation fluctuations; it is the elementary form of the fluctuation-dissipation theorem, and it is why a response function diverges at a critical point, where fluctuations do.
Regression and least squares. The slope of the best linear predictor of Y from X is cov(X,Y)/var(X) and the fraction of variance it explains is ρ2; see Least-squares regression. Arfken's observation that the mean is the minimiser of the sum of squared deviations, which is rule 8, is the one-parameter case of the same principle.
Covariance matrices. For a vector of variables the covariances form a symmetric positive semi-definite matrix, and the Cauchy-Schwarz inequality of rule 7 is exactly the statement that its 2×2 minors are non-negative. Diagonalising it is principal component analysis (EOF and principal component analysis); propagating it through a linear dynamical model is the Kalman filter; and its role as the parameter of a Gaussian is the subject of The multivariate normal distribution.
Tail bounds. The variance is what converts into a probability statement through Chebyshev's inequality, P(∣X−μ∣≥kσ)≤k−2; see Markov, Chebyshev, Jensen and Chernoff inequalities. That one line is why a standard deviation is a useful summary and not merely a conventional one.
Gives the working recipe for turning the joint density of several random variables into the density of functions of them — invert the map, multiply by the modulus of the Jacobian, integrate out what is not wanted — and derives from it the convolution formula for the density of a sum of independent variables.
Almost every calculation in applied probability has the shape: I know the distribution of the inputs; what is the distribution of the output? The inputs are random variables X1,…,Xn with a known joint density, the output is Y=g(X1,…,Xn), and what is wanted is the density of Y.
There is one recipe. Complete g to an invertible map of Rn to itself by adjoining whatever extra coordinates are convenient; transform the density by the rule for changing variables inside a multiple integral, which contributes the modulus of a Jacobian determinant; then integrate out the coordinates that were only there to make the map invertible. Everything else in this unit — the formula for a monotone function of one variable, the convolution formula for a sum, the chi-squared density, the Box-Muller sampler — is this recipe with a particular choice of the extra coordinates.
The Jacobian is not a technicality. A density is a probability per unit volume, so when a map stretches volume the density must shrink in compensation; the Jacobian is exactly the local volume factor.
Definition(Joint distributions, marginals, independence)AdaptedAdapted from Grimmett & Stirzaker (2020), Definition (3.6.2) and the marginals on p. 67, and §4.5, Definitions (4.5.1)-(4.5.2) with items (3)-(7), pp. 110-112; Riley, Hobson & Bence (2006), §30.4.3 p. 1143 and eq. (30.130), p. 1199; Folland (1999), §10.1, pp. 315-316
Let X and Y be random variables on one probability space. Their joint distribution function is
F(x,y)=P(X≤x,Y≤y),(x,y)∈R2.
If X and Y are both discrete, their joint mass function is f(x,y)=P(X=x,Y=y), and the marginal mass functions are recovered by summing out the other variable,
fX(x)=y∑f(x,y),fY(y)=x∑f(x,y),
because the events {Y=y} partition the sample space.
The pair is jointly continuous with joint densityf:R2→[0,∞) if
F(x,y)=∫−∞y∫−∞xf(u,v)dudvfor all x,y,
in which case f=∂2F/∂x∂y wherever F is twice differentiable, and for reasonable B⊆R2,
P((X,Y)∈B)=∬Bf(x,y)dxdy.
Physicists read f(x,y)dxdy directly as the probability that (X,Y) lies in the rectangle of sides dx,dy at (x,y). The marginal densities are
fX(x)=∫−∞∞f(x,y)dy,fY(y)=∫−∞∞f(x,y)dx.
Expectations are computed without any further construction: E(g(X,Y))=∬g(x,y)f(x,y)dxdy, the two-variable law of the unconscious statistician of Expectation, variance and covariance.
Call X and Yindependent when, for every choice of Borel sets B1 and B2, the two events {X∈B1} and {Y∈B2} are independent in the sense of Probability spaces. Folland states this equivalently: the joint distribution on R2 is the product of the two individual distributions. All of this extends verbatim to n variables, with F(x1,…,xn), an n-fold integral, and n−1 integrations to obtain one marginal.
Discrete. Independence of X and Y is equivalent to the identity f(x,y)=fX(x)fY(y) holding at every point (x,y), and also to the weaker-looking requirement that f(x,y)=g(x)h(y) at every point, for some pair of one-variable functions g,h.
Jointly continuous. Independence of X and Y is equivalent to F(x,y)=FX(x)FY(y) at every point, and this in turn to the possibility of choosing the joint density to be fX(x)fY(y).
Factorisation suffices. If a joint density factorises as f(x,y)=g(x)h(y) on the whole plane, then X and Y are independent and fX,fY are g,h normalised. The qualifier on the whole plane is essential: a density supported on a region that is not a product of intervals does not factorise even if its formula appears to.
Functions preserve independence. If X1,…,Xn are independent and the index set is split into disjoint blocks, then Borel functions of separate blocks are independent. In particular, if X and Y are independent, so are u(X) and v(Y) for any Borel u,v.
Method(The transformation and convolution formulas)AdaptedAdapted from Grimmett & Stirzaker (2020), §4.7, Theorem (4.7.3), Corollary (4.7.4), p. 122, and eq. (4.7.16), p. 126; Grimmett & Stirzaker (2020), §4.8, Theorem (4.8.1) and the following specialisation, p. 129; Riley, Hobson & Bence (2006), §30.14, pp. 1206-1207, and eq. (30.79), p. 1160
1. Change of variables in n dimensions. Let X=(X1,…,Xn) have joint density f, and let
T:D→R,y=T(x)
be a one-to-one map of an open set D⊆Rn containing the support of f onto R⊆Rn, with inverse x=T−1(y) whose partial derivatives ∂xi/∂yj are continuous. Let
J(y)=det(∂yj∂xi)i,j=1n
be the Jacobian determinant of the inverse map. Then Y=T(X) has joint density
fY(y)={f(T−1(y))J(y),0,y∈R,otherwise.
2. One variable, piecewise monotone. Let X have density fX and let I1,…,Im be intervals partitioning R on each of which g is strictly monotone and continuously differentiable, with inverse hi defined on g(Ii). Then Y=g(X) has density
fY(y)=i=1∑mfX(hi(y))hi′(y),
a term being omitted wherever its hi is undefined. For a single increasing or decreasing branch this is the n=1 case of part 1; for the affine map g(x)=ax+b with a=0 it reads fY(y)=∣a∣−1fX((y−b)/a).
3. Sums. If X and Y have joint density f, then X+Y has density
fX+Y(z)=∫−∞∞f(x,z−x)dx.
If moreover X and Y are independent, this is the convolution of the two marginal densities,
For independent integer-valued X and Y the same statement reads P(X+Y=n)=∑rP(X=r)P(Y=n−r).
Idea of the proof(Why a Jacobian, and why a sum convolves)Original
A joint density is a probability per unit area. If a map T carries a small region A near x onto a small region T(A) near y=T(x), then (X,Y)∈A exactly when the transformed pair lies in T(A), so the two regions carry the same probability. Probability equals density times area, and the areas differ by the local stretching factor of T — which is ∣detDT∣ going forwards, hence ∣J∣=∣detDT−1∣ going backwards. So the new density is the old one evaluated at the preimage, divided by the forward stretch, which is the same as multiplied by ∣J∣. Nothing probabilistic has happened: this is Multiple integrals and change of variables, applied to the integral that defines a probability.
The sum formula is this with a deliberately wasteful choice of map. To find the law of Z=X+Y one needs an invertible map, so adjoin a second coordinate that will be discarded: take (x,y)↦(z,w)=(x+y,x). This is a shear; it preserves area, so ∣J∣=1, and the joint density of (Z,W) is f(w,z−w). Integrating out the unwanted w gives the marginal density of Z, which is the stated formula. Independence then factorises the integrand, and the integral of a product of two shifted functions is by definition a convolution.
That the density of a sum of independent variables is a convolution is the reason transforms dominate the subject: a convolution becomes a product under the Fourier or Laplace transform, so the sum of independent variables becomes multiplication of characteristic functions. The piecewise formula of part 2 is the only genuine addition to the picture, and it says that when a map folds the line, each sheet of the fold contributes its own term.
Proof(The three formulas)AdaptedAdapted from Grimmett & Stirzaker (2020), the proof of Corollary (4.7.4), pp. 122-123, and of Theorem (4.8.1), p. 129; Riley, Hobson & Bence (2006), the probability-element derivation in §30.14, p. 1206
Part 1. Let B⊆R be a region and A=T−1(B)⊆D its preimage, so that Y∈B holds exactly when X∈A. Then
P(Y∈B)=P(X∈A)=∫Af(x)dx.
By the change-of-variables theorem for multiple integrals, applied to the one-to-one continuously differentiable map T−1:B→A,
∫Af(x)dx=∫Bf(T−1(y))∣J(y)∣dy.
Thus P(Y∈B)=∫Bf(T−1(y))∣J(y)∣dy for every such B. Comparing with the definition of a joint density for Y — namely that P(Y∈B)=∫BfY for all B — identifies the integrand, and outside R the probability is zero so the density may be taken to be zero. Riley, Hobson and Bence obtain the same formula in one line by equating the probability elements ∣fdx1⋯dxn∣ and ∣gdy1⋯dyn∣; that is the infinitesimal version of this argument.
Part 2. Fix i and suppose first that g is increasing on Ii. For y∈g(Ii),
P(g(X)≤y,X∈Ii)=P(X≤hi(y),X∈Ii),
because g increasing makes g(x)≤y equivalent to x≤hi(y) within Ii. Differentiating in y and using the chain rule gives the contribution fX(hi(y))hi′(y), and hi′>0 here so this is fX(hi(y))∣hi′(y)∣. If g is decreasing on Ii the inequality reverses, the contribution is −fX(hi(y))hi′(y), and hi′<0, so again the contribution is fX(hi(y))∣hi′(y)∣. The events {X∈Ii} are disjoint with union of full probability, since the Ii partition R, so
P(g(X)≤y)=i∑P(g(X)≤y,X∈Ii),
and differentiating term by term gives the stated sum.
Part 3. Apply part 1 to the map T(x,y)=(x+y,x)=(z,w) of R2 onto itself. Its inverse is x=w, y=z−w, with
so ∣J∣=1 and the joint density of (Z,W) is fZ,W(z,w)=f(w,z−w). Integrating out w gives the marginal density of Z=X+Y:
fX+Y(z)=∫−∞∞f(w,z−w)dw,
which is the stated formula after renaming w as x. If X and Y are independent then f(x,y)=fX(x)fY(y) by criterion 2, and the integrand becomes fX(x)fY(z−x); substituting y=z−x exchanges the two roles and gives the second form. The exchange of the order of integration implicit in computing a marginal is justified by Tonelli's theorem (Fubini and Tonelli), the integrand being non-negative.
In the discrete case the same argument needs no Jacobian: the events {X=r,Y=n−r}, as r ranges over the integers, are disjoint with union {X+Y=n}, so their probabilities add. □
Remark(Where each hypothesis is used)AdaptedAdapted from Grimmett & Stirzaker (2020), the hypotheses of Theorem (4.7.3) and Corollary (4.7.4), p. 122, the discussion of many-one maps preceding eq. (4.7.16), p. 126, Example (4.7.9), p. 125, and the footnote and Exercise (4.5.1) on marginals, p. 111
One-to-one, on the support only. The transformation rule needs T to be injective where the density lives; what it does elsewhere is irrelevant, and a failure of injectivity on a set of Lebesgue measure zero is harmless because such a set carries no probability. The polar map of the first example below is not injective at the origin, and that costs nothing. When injectivity genuinely fails on a set of positive measure, the map folds, and part 2 — or its n-dimensional analogue, a sum over the sheets — is required. Using part 1 on a folding map loses the contributions of all but one sheet; applying it to Y=X2 gives half the right answer.
A non-vanishing Jacobian, and continuous partial derivatives. These are the hypotheses of the change-of-variables theorem for multiple integrals, and they are used in exactly one place, the second display of the proof of part 1. Where J=0 the map collapses volume and the formula would assign density zero to a point that may carry probability; such points must be excluded, and again a measure-zero set of them is harmless.
Getting the range right is half the work. The formula is fY(y)=f(T−1(y))∣J(y)∣on R=T(D) and zero elsewhere, and the commonest error in practice is to write the formula down correctly and then integrate it over the wrong region. Grimmett and Stirzaker's Example (4.7.8) is the standard trap: for independent exponentials, the pair X=X1, S=X1+X2 has joint density λ2e−λs, which looks like a product of a function of x with a function of s — but only on the region 0≤x≤s, which is not a rectangle, so X and S are not independent. This is why criterion 3 above insists that the factorisation hold on the whole plane.
Independence, used only in the last step of part 3. The formula fX+Y(z)=∫f(x,z−x)dx holds for any jointly continuous pair, dependent or not; it is only the collapse to a convolution of marginals that needs f=fXfY. Applying the convolution formula to dependent variables is a genuine error, not a conservative approximation: it can give a density with the wrong variance in either direction, since var(X+Y) differs from varX+varY by 2cov(X,Y).
This machinery does not condition on null events. It computes the law of a transformed variable, not the law of one variable given that another takes an exact value. Grimmett and Stirzaker's Example (4.7.9) shows what goes wrong if the two are confused: for independent exponentials X1,X2, "the conditional density of X1+X2 given X1=X2" comes out as λ2ye−λy if the condition is approached through the ratio X1/X2→1 and as λe−λy if it is approached through the difference X1−X2→0. The event has probability zero and the answer depends on the approach. Conditioning of this kind belongs to Conditional expectation.
Joint continuity is stronger than continuity of each variable. A jointly continuous pair has continuous marginals in the sense of having densities, but those densities need not be continuous functions; Grimmett and Stirzaker record a density, due to Clarke, whose first marginal is discontinuous. Conversely, two variables each with a density need not be jointly continuous at all: if Y=X then the pair is supported on a line, which has zero area, and no joint density exists.
Method(Finding the density of a function of random variables)Original
Given the joint density f of X=(X1,…,Xn) and a function Y1=g(X) whose density is wanted:
Complete the map. Choose n−1 further functions Y2,…,Yn of X so that T=(g,Y2,…,Yn) is one-to-one on the support of f. Choose them for algebraic convenience: a coordinate of X itself is usually the easiest, and that choice turns step 5 into an integral you can do.
Invert. Solve y=T(x) for x=T−1(y) explicitly. If the solution is multi-valued, the map folds: split the support into pieces on which it is one-to-one, carry steps 2-4 out on each, and add the results.
Jacobian. Form the matrix of partial derivatives ∂xi/∂yj of the inverse map and take the modulus of its determinant, ∣J(y)∣. (Differentiating the forward map instead and taking the reciprocal of the modulus of its determinant gives the same number, and is often less work.)
Substitute. Write fY(y)=f(T−1(y))∣J(y)∣.
Determine the range and integrate out. Express the image R=T(D) as explicit inequalities on y — this is the step that is easy to get wrong — and then integrate fY over y2,…,yn within that range to obtain the marginal density of Y1.
Shortcuts worth knowing.
For Y=g(X) with a single variable and a monotone or piecewise monotone g, use the formula of part 2 directly.
For a sum of independent variables, use the convolution formula rather than re-deriving it; for a sum of many, use characteristic functions instead, since convolution is associative but tedious.
If the problem has a symmetry, exploit it: for a function of ∑iXi2 with a spherically symmetric joint density, integrating over a shell of radius r and thickness dr replaces the whole construction, as Riley, Hobson and Bence do in deriving the chi-squared density.
When only the distribution function of Y is wanted, computing P(g(X)≤y) as an integral of f over the region {g≤y} is often quicker than finding the density.
Check. Verify that fY≥0, that ∫fY=1, and that the mean and variance agree with what linearity and rule 6 of Expectation, variance and covariance predict. These three checks between them catch a missing Jacobian, a wrong range and an algebraic slip.
Example(Two velocity components to a speed: the two-dimensional Maxwell distribution)Original
A gas is confined to a thin film so that its molecules move in a plane. The two velocity components Vx and Vy of a molecule are independent and each N(0,σ2), where σ2=kBT/m. What is the distribution of the speed S=Vx2+Vy2?
Step 1: the joint density and the completed map. By independence,
so ∣J∣=s. The map is one-to-one on {s>0,0≤θ<2π}; it fails only at the origin, a single point of zero area, which carries no probability.
Step 4: substitute. Since vx2+vy2=s2,
fS,Θ(s,θ)=2πσ21e−s2/2σ2⋅s=a function of sσ2se−s2/2σ2⋅a function of θ2π1,s>0,0≤θ<2π.
Step 5: range and marginal. The range is the product region {s>0}×[0,2π), so the factorisation is valid on the whole of it and criterion 3 applies: S and Θ are independent, with Θ uniform on [0,2π) and
fS(s)=σ2se−s2/2σ2,s>0,
the Rayleigh density. Had the Jacobian factor s been forgotten, the "density" would have been largest at s=0, which is visibly wrong: there are more directions available at larger radius.
Step 7: checks and numbers. Substituting u=s2/2σ2 gives ∫0∞fS=∫0∞e−udu=1. The distribution function is FS(s)=1−e−s2/2σ2, so
P(S≤σ)=1−e−1/2=0.393469,P(S≤2σ)=1−e−2=0.864665.
The most probable speed is at the maximum of fS: fS′(s)=(1−s2/σ2)e−s2/2σ2/σ2 vanishes at s=σ. The mean and root-mean-square speeds are
the second obtained from linearity of expectation without any integration, and the ordering σ<1.2533σ<1.4142σ of mode, mean and root-mean-square is the familiar one for a right-skewed distribution. For nitrogen at 300K, with m=28.013u, one finds σ=kBT/m=298.40ms−1, a mean speed E(S)=373.99ms−1, a root-mean-square speed of 422.00ms−1, and P(S≤298.40ms−1)=0.3935.
The same calculation in three dimensions, with ∣J∣=r2sinϕ, produces the factor s2 instead of s and so the full Maxwell-Boltzmann speed distribution.
Example(A decay time plus a readout delay: convolution)Original
An unstable state decays at a time D measured from the start of a gate, exponentially distributed with mean lifetime τ; the electronics then timestamps the event after a delay U that is uniformly distributed between 0 and τ, independently of D. Work in units of τ, so that D has density e−d for d>0 and U has density 1 on (0,1). What is the distribution of the recorded time T=D+U?
Convolution. By part 3, with fU(u)=1 on (0,1) and zero elsewhere,
the upper limit being min(1,t) because fD vanishes unless t−u>0, that is unless u<t, while fU vanishes unless u<1. Carrying out the elementary integral ∫0ae−(t−u)du=e−t(ea−1) with a=min(1,t):
Checks. The two branches agree at t=1, where both give 1−e−1=0.632121, so fT is continuous; it rises from 0 on (0,1) and decays exponentially thereafter, with its maximum at t=1. Its total mass is
neither of which required the density at all — a useful independent check on it.
A checkable number. The probability that the recorded time falls within one mean lifetime of the start of the gate is
P(T≤1)=∫01(1−e−t)dt=1−(1−e−1)=e−1=0.367879.
For comparison, the undelayed decay satisfies P(D≤1)=1−e−1=0.632121. The delay has moved a large part of the probability past the one-lifetime mark: the uniform smearing, though small in mean, is comparable to the lifetime and badly distorts the early part of the spectrum. Fitting a pure exponential to T would overestimate τ, and the standard remedy — deconvolving the known instrument response from the measured spectrum — is exactly the inverse of the calculation just performed.
Application(Where transformations of densities are the work)Original
Instrument response. Every measured spectrum is the convolution of the true spectrum with the response function of the apparatus: a pulse-height distribution smeared by electronic noise, a spectral line broadened by a finite slit, a time-of-flight peak widened by a trigger jitter. Part 3 is the statement of this, and the Gaussian-on-Gaussian case explains why resolutions add in quadrature, σobs2=σtrue2+σinst2. Unfolding or deconvolution reverses it, and is hard precisely because convolution destroys the high-frequency content that would be needed to invert it.
Derived quantities and error propagation. A measured quantity is almost never the one of interest: one measures a period and wants a frequency, measures two lengths and wants an aspect ratio, measures counts and wants a rate. Part 2 gives the exact density of the derived quantity, and the first-order version of it — the affine case fY(y)=∣a∣−1fX((y−b)/a) with a=g′ — is the linearised rule of Error propagation. The exact rule matters when the function is strongly curved over the width of the input distribution, which is when the linear rule silently fails.
Sampling distributions. The chi-squared, Student and Fisher distributions of Sampling distributions are each produced by one application of part 1 to independent normals: χn2=∑iZi2 by the spherical-shell version of the method, as Riley, Hobson and Bence carry it out; Student's t as a ratio of a normal to the square root of an independent chi-squared; and F as a ratio of two chi-squareds. The Rayleigh density of the first example is the n=2 case, in disguise.
Random number generation. Inverse-transform sampling is part 2 run backwards: if U is uniform on (0,1) then F−1(U) has distribution function F. The Box-Muller algorithm is the first worked example read backwards — generate the Rayleigh radius from a uniform by inverting 1−e−s2/2, generate an independent uniform angle, and the Cartesian components are two independent standard normals. The whole of Monte Carlo rests on such transformations, and the Jacobian is what makes the output density correct.
Statistical mechanics and kinetic theory. The step from a distribution over velocity vectors to one over speeds, over energies, or over a relative velocity between two molecules is in every case an application of this method; the familiar density-of-states factors s, s2, E are Jacobians and nothing else. The same is true of the change from Cartesian to action-angle variables in phase space, where Liouville's theorem is the statement that the relevant Jacobian is 1 along the flow.
Sums and limit theorems. Repeated convolution is what the central limit theorem is about, and the reason the proof is conducted with characteristic functions is that they turn the convolution of part 3 into multiplication. The stable distributions — normal, Cauchy, Lévy — are exactly the fixed points of repeated convolution followed by rescaling.
Defines the conditional distribution and the conditional expectation E(Y∣X) as a random variable, proves the tower, pull-through and conditional-variance rules, shows that E(Y∣X) is the function of X closest to Y in mean square, and works a mixture example whose checks are log2 and 21.
Two quantities are measured together and only one of them is observed. A detector records an energy deposit but not the particle that caused it; a survey records a household's income but not its wealth; a simulation stores a coarse field but not the microstate behind it. The question in every case is the same: given what was seen, what is the average of what was not?
The answer is a number for each observed value, and therefore a function of the observed variable, and therefore itself a random variable. That one shift of viewpoint - from E(Y∣X=x), a number, to E(Y∣X), a random variable - is what makes conditional expectation a tool rather than a notation. It gives the tower property, which computes hard expectations in two easy stages; it gives the conditional variance formula, which splits a variance into a part explained by X and a part not; and it gives the best mean-square prediction of Y from X, which is the whole content of regression.
This unit covers the case where X and Y are both discrete or both continuous. The general construction, which conditions on a σ-field rather than on a variable, is stated at the end and used in Martingales.
Definition(Conditional distributions, conditional expectation and conditional variance)AdaptedAdapted from Grimmett & Stirzaker (2020), §3.7, Definitions (3.7.1), (3.7.2) and (3.7.3), pp. 71-72; Grimmett & Stirzaker (2020), §4.6, Definitions (4.6.1) and (4.6.2) with the preceding limit argument, and the definition of $\psi$ following Example (4.6.3), pp. 117-118; Riley, Hobson & Bence (2006), §30.11.3, p. 1199
Let X and Y be random variables on a probability space (Ω,F,P), with joint mass function fX,Y or joint density fX,Y in the sense of Joint distributions and transformations, and write fX for the marginal of X.
Discrete case. For any x with P(X=x)>0, the conditional mass function of Y given X=x is
fY∣X(y∣x)=P(Y=y∣X=x)=fX(x)fX,Y(x,y),
and the conditional distribution function is FY∣X(y∣x)=P(Y≤y∣X=x). Both are undefined at any x with P(X=x)=0, because the elementary conditional probability of Probability spaces conditions only on events of positive probability.
Continuous case. Here {X=x} has probability zero and the elementary definition is unavailable. The way round it is to condition on the thickened event {x≤X≤x+h}, which has probability ≈fX(x)h, and let h↓0:
One formula, fY∣X=fX,Y/fX, serves both cases. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §30.11.3, p. 1199) describe it as taking the slice of the joint density at the observed value and dividing by whatever is needed to make it integrate to 1; since ∫fX,Y(x,v)dv=fX(x), that divisor is exactly fX(x).
Conditional expectation. Suppose E∣Y∣<∞. For each x with fX(x)>0 put
Then the conditional expectation of Y given X is the random variable
E(Y∣X)=ψ(X),
the composition of the function ψ with the random variable X. Its value at a sample point ω is ψ(X(ω)). Grimmett and Stirzaker emphasise that although the phrase sounds like a number, the object is a random variable; the number is E(Y∣X=x) for a fixed x.
Conditional variance. With ψ as above and E(Y2)<∞, set v(x)=E((Y−ψ(x))2∣X=x), the variance of the conditional distribution of Y given X=x, and define the random variable var(Y∣X)=v(X). Equivalently var(Y∣X)=E(Y2∣X)−(E(Y∣X))2.
Everywhere below, an identity between conditional expectations is an identity between random variables and holds with probability 1; we follow Grimmett and Stirzaker in not repeating the qualification.
Definition(The rules of conditional expectation)AdaptedAdapted from Grimmett & Stirzaker (2020), §3.7, Theorem (3.7.4), Theorem (3.7.6) with eq. (3.7.7), and Exercises (3.7.1)(a)-(f), (3.7.2) and (3.7.4), pp. 72-74; Grimmett & Stirzaker (2020), §7.7, Lemmas (7.7.10) and (7.7.11), p. 381
Let X, Y, Z be random variables with E∣Y∣,E∣Z∣<∞, let a,b∈R, and let g be a function for which the expectations written below exist.
Partial averaging.E(ψ(X)g(X))=E(Yg(X)), where ψ(X)=E(Y∣X).
Tower property (first form). Taking g≡1 in rule 1,
E(E(Y∣X))=E(Y),that isE(Y)=x∑E(Y∣X=x)P(X=x)
in the discrete case and E(Y)=∫E(Y∣X=x)fX(x)dx in the continuous case.
Linearity.E(aY+bZ∣X)=aE(Y∣X)+bE(Z∣X), and E(1∣X)=1.
Positivity. If Y≥0 then E(Y∣X)≥0.
Pull-through.E(Yg(X)∣X)=g(X)E(Y∣X): a function of the conditioning variable behaves as a constant.
Independence. If X and Y are independent then E(Y∣X)=E(Y), a constant.
Relabelling.E(Y∣h(X))=E(Y∣X) whenever h is one-to-one.
Tower property (general form).E(E(Y∣X,Z)X)=E(Y∣X)=E(E(Y∣X)X,Z). Conditioning twice, coarsely after finely or finely after coarsely, gives the coarser answer.
Almost sure uniqueness. If φ(X) and ψ(X) both satisfy rule 1 for every admissible g, then P(φ(X)=ψ(X))=1. Consequently rule 1 characterisesE(Y∣X) and may be taken as its definition; Grimmett and Stirzaker say so explicitly and use that route for the general theory.
Conditional variance formula. If E(Y2)<∞ then
var(Y)=E(var(Y∣X))+var(E(Y∣X)).
Definition(Conditional expectation is the best mean-square prediction)AdaptedAdapted from Grimmett & Stirzaker (2020), §7.9, Definition (7.9.5), Theorem (7.9.9), Projection theorem (7.9.14) and Theorem (7.9.17), pp. 388-391
Write ∥U∥2=E(U2), and let
H={h(X):h:R→R,E(h(X)2)<∞}
be the set of square-integrable functions of X. An element Y∈H carries the name best predictor (in full, minimum mean-squared-error predictor, which is Grimmett and Stirzaker's term) when ∥Y−Y∥2≤∥Y−Z∥2 for every Z∈H. Suppose E(Y2)<∞. Then:
(a) Existence and uniqueness. A best predictor exists, and it is unique up to a null event: if Y and Y both attain the minimum then P(Y=Y)=1.
(b) Orthogonality. For M∈H, the two statements
E((Y−M)Z)=0for all Z∈Hand∥Y−M∥2≤∥Y−Z∥2for all Z∈H
are equivalent. This is the projection theorem: minimising the distance and being perpendicular to H are the same condition.
(c) Identification. The best predictor of Y given X is E(Y∣X), and the minimum value of the mean squared error is
E((Y−E(Y∣X))2)=E(var(Y∣X)).
Because (c) holds for every pair X,Y with E(Y2)<∞, whatever their types, Grimmett and Stirzaker point out that it can be adopted as the definition of E(Y∣X): the conditional expectation is the orthogonal projection of Y onto the functions of X. That is the formulation which generalises.
Definition(The general definition, by conditioning on a sub-σ-field)AdaptedAdapted from Grimmett & Stirzaker (2020), §7.9, Definition (7.9.20) with eqs. (7.9.21)-(7.9.23), pp. 391-392; Folland (1999), §3.2, Exercise 17, p. 93
Let G⊆F be a sub-σ-field. For Y with E(Y2)<∞, Grimmett and Stirzaker define E(Y∣G) to be the G-measurable random variable with
E((Y−E(Y∣G))Z)=0for all G-measurable Z with E(Z2)<∞,
which they then show may equivalently be required only for the indicators Z=1G, G∈G; taking G=Ω recovers E(E(Y∣G))=E(Y). Taking G to be the smallest σ-field making a family {Xi} measurable returns the earlier notion, E(Y∣G)=E(Y∣Xi,i∈I).
Folland (Folland (1999), §3.2, Exercise 17, p. 93) gives the construction that removes the second-moment hypothesis. On a σ-finite measure space (Ω,F,μ) with sub-σ-algebra G and ν=μ∣G, every f∈L1(μ) admits a G-measurable g∈L1(ν) with
∫Gfdμ=∫Ggdνfor every G∈G,
and g is determined up to a ν-null set. Folland records that g is what probability theory calls the conditional expectation of f on G. Specialised to a probability measure this is E(Y∣G), defined for every Y∈L1 by exactly the partial-averaging property of rule 1, and existence comes from the Radon-Nikodym theorem (The Radon-Nikodym theorem) rather than from a projection.
Idea of the proofOriginal
Everything in this unit comes from one identity. Summing or integrating fY∣X(y∣x)fX(x) over x rebuilds the joint law, because by construction that product isfX,Y(x,y). So averaging the conditional average over the conditioning variable undoes the conditioning: that is the tower property, and the pull-through and partial-averaging rules are the same computation with an extra factor g(x) carried along, which can leave the inner average because it does not depend on y.
The prediction statement needs a second idea, geometric rather than algebraic. Square-integrable random variables form an inner-product space with (U,V)=E(UV), the functions of X form a closed subspace H of it, and the best predictor is the foot of the perpendicular from Y to H. Existence of that foot is the standard closest-point argument, in which the parallelogram law turns a minimising sequence into a Cauchy sequence. Identifying the foot as E(Y∣X) is then one line: partial averaging says precisely that Y−E(Y∣X) is orthogonal to every g(X).
Proof(Partial averaging, the tower property and the algebraic rules)AdaptedAdapted from Grimmett & Stirzaker (2020), §3.7, proofs of Theorems (3.7.4) and (3.7.6), pp. 72-73
We give the discrete argument; the continuous one is the same with sums replaced by integrals, which is legitimate because the conditional density was defined so that fY∣X(y∣x)fX(x)=fX,Y(x,y) at every x with fX(x)>0, and the set {x:fX(x)=0} contributes nothing to any of the integrals.
using fY∣XfX=fX,Y. The last double sum is E(Yg(X)), again by the law of the unconscious statistician, now for the function (x,y)↦yg(x) of the pair. The interchange of the two sums is justified by absolute convergence, which is the hypothesis that the expectations exist.
Rule 2. Put g≡1 in rule 1. The right-hand side is E(Y) and the left-hand side is E(ψ(X)), which written out is ∑xE(Y∣X=x)P(X=x).
Rules 3, 4, 6. Each holds at every fixed x and therefore holds after substituting X for x. Linearity holds because E(⋅∣X=x) is an expectation with respect to the probability mass function fY∣X(⋅∣x), and expectations are linear; the same remark gives E(1∣X=x)=1 and E(Y∣X=x)≥0 for Y≥0. For rule 6, independence makes fX,Y=fXfY, hence fY∣X(y∣x)=fY(y) for every x with fX(x)>0, so ψ(x)=E(Y) does not depend on x.
Rule 5. Fix x with fX(x)>0. Conditional on {X=x} the factor g(X) equals the constant g(x), so
E(Yg(X)∣X=x)=y∑yg(x)fY∣X(y∣x)=g(x)E(Y∣X=x).
Substituting X for x gives the claim.
Rule 7. If h is one-to-one then the events {X=x} and {h(X)=h(x)} are the same event, so the two conditional distributions, and hence the two conditional expectations, coincide. Grimmett and Stirzaker make exactly this observation.
Rule 8. Write χ(x,z)=E(Y∣X=x,Z=z), so that E(Y∣X,Z)=χ(X,Z). Conditioning on {X=x} and using the definition twice,
The factor fX,Z(x,z) cancels - this is the only step in the proof, and it is why the property holds - leaving ∑yy∑zfX,Y,Z(x,y,z)/fX(x)=∑yyfY∣X(y∣x)=ψ(x). For the second equality in rule 8, note that E(Y∣X) is already a function of X alone, hence a function of the pair (X,Z), so rule 5 with Y replaced by 1 returns it unchanged.
Rule 9. If E(φ(X)g(X))=E(ψ(X)g(X)) for all admissible g, write D=φ(X)−ψ(X) and choose g(x)=sign(φ(x)−ψ(x)), which is bounded. Then E∣D∣=E(Dg(X))=0, and a non-negative random variable with zero mean vanishes with probability 1.
Rule 10. Abbreviate ψ=E(Y∣X) and v=var(Y∣X)=E(Y2∣X)−ψ2. Take expectations and use rule 2 twice:
E(v)=E(Y2)−E(ψ2),var(ψ)=E(ψ2)−(Eψ)2=E(ψ2)−(EY)2.
Adding the two displays cancels E(ψ2) and leaves E(Y2)−(EY)2=var(Y). □
Proof(Best mean-square prediction)AdaptedAdapted from Grimmett & Stirzaker (2020), §7.9, proofs of Theorem (7.9.9), the Projection theorem (7.9.14) and Theorem (7.9.17), pp. 389-391
Step 1: H is a closed linear space. A linear combination of square-integrable functions of X is a square-integrable function of X, by the Minkowski inequality, so H is a linear space. It is closed in the norm ∥⋅∥2: if hn(X) is a Cauchy sequence in H then it converges in mean square to some square-integrable limit, and a mean-square limit of functions of X is again (almost surely equal to) a function of X. This is the completeness of L2 (L^p spaces).
Step 2: a closest point exists. Let d=inf{∥Y−Z∥2:Z∈H} and pick Z1,Z2,⋯∈H with ∥Y−Zn∥2→d. For any A,B∈H the parallelogram law, applied to the vectors Y−A and Y−B, reads
∥A−B∥22=2∥Y−A∥22+2∥Y−B∥22−4Y−21(A+B)22;
expanding both sides as expectations verifies it. Take A=Zm, B=Zn. Since 21(Zm+Zn)∈H (here linearity of H is used), the last term is at least 4d2, so
∥Zm−Zn∥22≤2(∥Y−Zm∥22−d2)+2(∥Y−Zn∥22−d2)⟶0.
Thus (Zn) is Cauchy, and by closedness it converges in mean square to some Y∈H. The triangle inequality gives ∥Y−Y∥2≤∥Y−Zn∥2+∥Zn−Y∥2→d, so ∥Y−Y∥2=d and Y attains the minimum.
Step 3: the minimiser is unique up to a null event. If Y∈H also attains d, apply the parallelogram law with A=Y, B=Y: the first two terms give 4d2 and the last is at least 4d2, so ∥Y−Y∥22≤0. Hence E((Y−Y)2)=0 and the two agree with probability 1.
Step 4: minimising is the same as being orthogonal. Suppose M∈H satisfies E((Y−M)Z)=0 for all Z∈H. For M′∈H the difference M−M′ lies in H, so the cross term in
E((Y−M′)2)=E((Y−M)2)+2E((Y−M)(M−M′))+E((M−M′)2)
vanishes, leaving ∥Y−M′∥22=∥Y−M∥22+∥M−M′∥22≥∥Y−M∥22: M minimises. Conversely suppose M minimises but E((Y−M)Z)=c>0 for some Z∈H; rescaling Z we may take E(Z2)=1. Then M′=M+cZ∈H has
E((Y−M′)2)=E((Y−M)2)−2cE((Y−M)Z)+c2=E((Y−M)2)−c2,
contradicting minimality. (If c<0 replace Z by −Z.) So a minimiser is orthogonal to H.
Step 5: the projection is E(Y∣X). First, ψ(X)=E(Y∣X) lies in H: by the conditional Jensen inequality in the form ψ(x)2=E(Y∣X=x)2≤E(Y2∣X=x), which is the Cauchy-Schwarz inequality applied inside the conditional distribution, taking expectations gives E(ψ(X)2)≤E(Y2)<∞. Second, for any Z=h(X)∈H,
E((Y−ψ(X))h(X))=E(Yh(X))−E(ψ(X)h(X))=0
by rule 1 with g=h. So ψ(X) is orthogonal to H, and by step 4 it minimises. Finally, by rules 2 and 5,
E((Y−ψ(X))2)=E(E((Y−ψ(X))2X))=E(var(Y∣X)),
the inner conditional expectation being the variance of the conditional distribution because ψ(X) is its mean. □
Remark(Where the hypotheses are used, and what fails without them)Original
fX(x)>0 is needed to define anything at x. Conditioning is division by fX(x). At values the conditioning variable never takes, there is no conditional distribution, and E(Y∣X) is simply not specified there. This costs nothing, because the set {fX=0} has probability zero, but it is why every identity above is an identity "with probability 1" and no more.
The limiting definition in the continuous case is a choice, not a theorem.P(Y≤y∣X=x) has no elementary meaning; we defined it as the limit of P(Y≤y∣x≤X≤x+h). Conditioning the same event on a different family of shrinking sets can give a different answer. The standard warning case is a point on a sphere: conditioning on a shrinking band of latitude and conditioning on a shrinking wedge of longitude produce different "conditional distributions" of the other coordinate, because the two families of thick events shrink in different ways. Nothing below is affected, since we fix the family once and for all, but the notation P(⋅∣X=x) hides a convention.
E∣Y∣<∞ is what rules 1-9 need; E(Y2)<∞ is what the prediction theorem needs. Rule 10 and the whole of the projection argument use the second moment, in three places: to put ψ(X) in H, to make the parallelogram law meaningful, and to make d finite. For Y with E∣Y∣<∞ but E(Y2)=∞ - a Pareto variable with index 3/2, say - the conditional expectation still exists by the Folland construction above, and rules 1-9 still hold, but "closest in mean square" is vacuous because every candidate is at infinite distance.
Linearity of H is used, and it is not free. Step 2 needs 21(Zm+Zn)∈H. If H were replaced by a non-convex set of predictors - say, the functions of X taking only the values 0 and 1 - the minimising sequence argument collapses and the closest point can fail to be unique.
Rule 5 is why the proof of rule 8 is a cancellation and not a computation. The only property of g(X) used is that it is already determined by the conditioning variable. Replace g(X) by g(Y) and the rule is false: for X,Y independent standard normals, E(Y⋅Y∣X)=1=YE(Y∣X)=0.
Uncorrelated is weaker than "conditional mean is constant", which is weaker than independence. Rule 6 goes one way only. Let X be uniform on {−1,0,1} and Y=X2. Then E(Y∣X)=X2, not constant, yet cov(X,Y)=E(X3)=0. Conversely E(Y∣X)=E(Y) forces cov(X,Y)=0 by rules 1 and 2, but does not force independence, since it says nothing about E(Y2∣X).
Method(Computing a conditional expectation)Original
Decide which variable is observed. Write down the joint law fX,Y of the observed X and the unobserved Y. If the model is specified the other way round - X has some law and, given X=x, Y has some law - then fX,Y(x,y)=fX(x)fY∣X(y∣x) is already the joint law and step 2 is unnecessary.
Marginalise. Compute fX(x)=∑yfX,Y(x,y) or ∫fX,Y(x,y)dy, and note where it vanishes.
Divide. Form fY∣X(y∣x)=fX,Y(x,y)/fX(x) and check that it integrates (or sums) to 1 in y for each fixed x. This catches most algebra slips at once.
Average in y. Compute ψ(x)=∫yfY∣X(y∣x)dy. Recognising the conditional law as a standard distribution usually saves the integral: if fY∣X(⋅∣x) is exponential with rate λ(x), then ψ(x)=1/λ(x) with no work.
Substitute the random variable.E(Y∣X)=ψ(X). If what is wanted is E(Y), stop here and apply the tower property, E(Y)=∫ψ(x)fX(x)dx, which is usually a one-dimensional integral in place of a two-dimensional one.
Check. The tower property is a free consistency test: compute E(Y) both ways. If the second moment is wanted as well, the conditional variance formula gives a second independent check, since the two terms E(var(Y∣X)) and var(E(Y∣X)) are separately non-negative and must sum to var(Y).
Example(A mixed exponential: conditioning, the tower property and the variance split)Original
A radioactive source is placed at a random distance from a counter, and the waiting time to the next recorded count is exponential with a rate that depends on the distance. Model this by letting the rate X be uniform on [1,2] and, given X=x, letting the waiting time Y be exponential with parameter x. The joint density is then
and zero elsewhere. It is a density: ∫0∞xe−xydy=1 for each x, and ∫121dx=1.
Step 1: the conditional law of Y given X. It was given: fY∣X(y∣x)=xe−xy for y>0, exponential with parameter x. Hence, reading the mean and variance of an exponential off the standard list,
The direct route would be the double integral ∫12∫0∞yxe−xydydx, whose inner integral is 1/x: the same calculation, with the tower property making the first stage invisible.
Step 3: the variance of Y, both ways. For the second moment, E(Y2∣X=x)=2/x2, so
Their sum is 1−(log2)2, matching var(Y) exactly, as rule 10 requires. Only 3.8% of the variance of the waiting time is attributable to not knowing the rate; the rest is the intrinsic randomness of the exponential.
Step 4: the best prediction and its error. By the theorem, the best mean-square prediction of Y from X is Y=1/X, with mean squared error E(var(Y∣X))=21. Compare the best linear prediction a+bX. Here E(XY)=E(XE(Y∣X))=E(X⋅X−1)=1 by rules 5 and 2, so
cov(X,Y)=1−23log2=−0.0397208…,var(X)=121,
and the best linear predictor has mean squared error var(Y)−cov(X,Y)2/var(X)=0.500614… The nonlinear predictor wins, as it must, but only just: over 1≤x≤2 the curve 1/x is nearly straight.
Step 5: a check by a different route. Rule 2 applied to the indicator of an event gives P(Y>1)=E(P(Y>1∣X))=∫12e−xdx=e−1−e−2=0.2325442…. Integrating the joint density directly, ∫12∫1∞xe−xydydx=∫12e−xdx, the same number. (All values here were confirmed symbolically.)
Application(Regression, filtering and coarse-graining)Original
Regression is conditional expectation. The regression function of Y on X is by definition x↦E(Y∣X=x), and the theorem above is the reason least squares is the default fitting criterion: minimising mean squared error over all functions of the predictors returns the conditional mean, and minimising it over affine functions returns the line that statisticians call the regression line. The familiar decomposition of a total sum of squares into "explained" and "residual" parts is rule 10 estimated from data, and the coefficient of determination R2 estimates var(E(Y∣X))/var(Y) - in the worked example, 0.0195/0.5195=3.8%.
Filtering and data assimilation. A Kalman filter propagates exactly one object: the conditional distribution of the unobserved state given all observations so far. The update step is Bayes' theorem on conditional densities, the forecast step is the tower property, and the reported state estimate is the conditional mean because that is the estimate with least mean squared error. Weather and ocean models assimilate observations by the same logic on a very large state space, with ensembles standing in for the conditional distribution.
Coarse-graining in statistical physics. Integrating a Gibbs distribution over fast degrees of freedom at fixed slow ones produces a conditional distribution, and the resulting effective Hamiltonian is built from conditional expectations of the microscopic one. The conditional variance formula is the statement that the fluctuation of an observable splits into a within-cell part and a between-cell part - the decomposition behind block-spin renormalisation and behind the separation of thermal noise from sampling noise in a Monte Carlo estimate.
Variance reduction in simulation. If a quantity Y can be averaged analytically given some auxiliary variable X, replacing the sampled Y by the computed E(Y∣X) leaves the mean unchanged (rule 2) and reduces the variance by exactly E(var(Y∣X)) (rule 10). This is conditional Monte Carlo, and it is free whenever the inner average is available in closed form.
Martingales. The defining property of a martingale, E(Sn+1∣X1,…,Xn)=Sn, is a statement about conditional expectation, and the whole theory in Martingales runs on the tower property. The same is true of the Markov property in Discrete-time Markov chains, which says that conditioning on the whole past is the same as conditioning on the present.
Defines the N(μ,V) density on Rn, proves its normalisation, mean and covariance matrix by orthogonal diagonalisation, shows that linear images, marginals and conditionals are all normal with explicit parameters, and works a three-dimensional example whose conditional law is N(9/5,13/5).
Three facts make the multivariate normal the working distribution of the physical sciences. It is the limit law of sums, by the central limit theorem, so it describes any quantity assembled from many small independent contributions. It is closed under every linear operation - sums, projections, rotations, marginalisation and conditioning all return normal laws - so a linear model with normal inputs has normal outputs and the bookkeeping never leaves the family. And it is determined entirely by a mean vector and a covariance matrix, so quoting a best-fit parameter set with its error matrix is a complete statement of the result, not a summary of one.
The price is one hypothesis that is easy to forget: jointly normal is much stronger than normal one coordinate at a time. The remark below gives a pair of standard normal variables that is not bivariate normal at all, and for which zero correlation does not give independence.
Definition(The multivariate normal density)AdaptedAdapted from Grimmett & Stirzaker (2020), §4.9, the discussion of (4.9.1)-(4.9.3) and Definition (4.9.4), pp. 131-132; Riley, Hobson & Bence (2006), §30.15.2, the form preceding (30.148) and eq. (30.148), pp. 1209-1210
Throughout, x, μ and t are column vectors in Rn and MT is the transpose of M. A real symmetric n×n matrix V is positive definite if uTVu>0 for every u=0; by the spectral theorem (Spectral theorem for Hermitian matrices) this holds exactly when all eigenvalues of V are strictly positive, and then V is invertible with detV>0.
A random vector X=(X1,…,Xn)T has the multivariate normal (or multinormal) distribution N(μ,V), with μ∈Rn and V symmetric positive definite, if its joint density is
The case μ=0, V=I is the standard multivariate normal, a vector of n independent N(0,1) variables. For n=1 the formula is the familiar N(μ,σ2) density with V=(σ2). For n=2, writing V=(σ12ρσ1σ2ρσ1σ2σ22) with ∣ρ∣<1, it is the bivariate normal density
which is the form Grimmett and Stirzaker (Grimmett & Stirzaker (2020), §4.5, Example (4.5.9), p. 113) write out; their standard case is μ=0, σ1=σ2=1.
A convention to watch. Grimmett and Stirzaker use row vectors and write the exponent as −21(x−μ)V−1(x−μ)′, with x′ for the transpose; Riley, Hobson and Bence use column vectors and write −21(x−μ)TV−1(x−μ), which is what we use. The two expressions are the same scalar, but every matrix identity below changes shape between them: a linear map is Y=XD in the row convention and Y=DTX in the column convention, and the induced covariance is D′VD either way. We follow the column convention of Riley, Hobson and Bence throughout and translate the statements we take from Grimmett and Stirzaker; no result depends on the choice.
The singular case detV=0 is excluded by the definition, because V−1 does not exist. It is covered by the alternative definition below.
Definition(The parameters, and the equivalent definition by linear combinations)AdaptedAdapted from Grimmett & Stirzaker (2020), §4.9, Theorem (4.9.5), Definition (4.9.7) and the discussion of the singular case, pp. 132-133; Grimmett & Stirzaker (2020), §5.8, Example (5.8.6), pp. 209-210; Riley, Hobson & Bence (2006), §30.15.2, the moment generating function and the identification $A=V^{-1}$, pp. 1209-1210
Let X∼N(μ,V) with V positive definite. Then
E(X)=μ, that is E(Xi)=μi for each i;
V is the covariance matrix: vij=cov(Xi,Xj), so the diagonal entries are the variances and V=E((X−μ)(X−μ)T);
the moment generating function and the characteristic function are
MX(t)=E(etTX)=exp(tTμ+21tTVt),φX(t)=E(eitTX)=exp(itTμ−21tTVt),
both valid for every t∈Rn (see Generating functions and Characteristic functions);
the quadratic form in the exponent is a chi-squared variable: (X−μ)TV−1(X−μ)∼χ2(n). Riley, Hobson and Bence add that if the Xi are required to satisfy m independent linear constraints ∑iciXi=0 then the same quadratic form is χ2(n−m).
An equivalent definition. Grimmett and Stirzaker also define the multivariate normal by the property that
aTX=a1X1+⋯+anXn has a univariate normal distribution for every a∈Rn.
Every X with the density above satisfies this. The converse holds with one exception: if some non-trivial linear combination is degenerate, P(aTX=b)=1 with a=0, then aTVa=var(aTX)=0, so V is singular and no density exists. Such an X is said to have a singular multivariate normal distribution; it lives on a proper affine subspace of Rn. When the distribution is non-singular the two definitions agree. The linear-combination definition is the more useful one in practice, because it needs no inverse and survives degeneracy.
Definition(Linear images, marginals and conditionals are normal)AdaptedAdapted from Grimmett & Stirzaker (2020), §4.9, Theorem (4.9.6) with the remark that the case $m<n$ is omitted, p. 133, and Exercises (4.9.3), (4.9.7) and (4.9.9), pp. 134-135; Grimmett & Stirzaker (2020), §4.6, Example (4.6.7), p. 119; Riley, Hobson & Bence (2006), §30.12, eqs. (30.140)-(30.141), p. 1205
Let X∼N(μ,V) on Rn.
(a) Affine images. Let D be an n×m matrix of rank m≤n, let c∈Rm, and set Y=DTX+c. Then
Y∼N(DTμ+c,DTVD)on Rm.
In particular every aTX is univariate normal with mean aTμ and variance aTVa, and the covariance law VY=SVXST that Riley, Hobson and Bence derive for general linearly related variables is the special case of (a) in their notation. Grimmett and Stirzaker state (a) for m≤n and prove it for m=n, remarking that the proof for m<n is harder and referring elsewhere for it; we follow them and prove the m=n case, then obtain the cases we need from the characteristic function.
(b) Marginals. Every sub-vector of X is multivariate normal, with the corresponding sub-vector of μ and the corresponding sub-matrix of V. In particular each Xi∼N(μi,vii).
(c) Conditionals. Split X, μ and V into blocks of sizes p and q=n−p,
X=(X1X2),μ=(μ1μ2),V=(V11V21V12V22),V21=V12T.
Then, conditional on X2=x2, the vector X1 is multivariate normal:
X1∣{X2=x2}∼N(μ1+V12V22−1(x2−μ2),V11−V12V22−1V21).
The conditional mean is affine in the observed x2; the conditional covariance S=V11−V12V22−1V21, the Schur complement of V22 in V, does not depend on x2 at all. Hence, in the language of Conditional expectation,
E(X1∣X2)=μ1+V12V22−1(X2−μ2),var(X1∣X2)=S(a constant matrix),
so for jointly normal variables the best mean-square predictor is already the best linear predictor.
(d) Independence equals zero correlation.X1 and X2 are independent if and only if V12=0. Componentwise: jointly normal variables are independent exactly when they are uncorrelated.
(e) Sampling. If C is any n×n matrix with CCT=V and W is standard multivariate normal, then μ+CW∼N(μ,V). Such a C is automatically invertible, since (detC)2=detV>0. Two standard choices of C are the symmetric square root of V and the Cholesky factor. Grimmett and Stirzaker state the Cholesky version in their row convention: with L lower triangular and V=L′L, the row vector Z=μ+XL is N(μ,V); transposing gives C=LT here.
Idea of the proofOriginal
One change of variables does nearly everything. Because V is real symmetric, the spectral theorem supplies an orthogonal Q with V=QΛQT and Λ=diag(λ1,…,λn). Setting y=QT(x−μ) rotates the quadratic form into ∑iyi2/λi and leaves the volume element alone, because ∣detQ∣=1. The density then factorises into n one-dimensional Gaussians with variances λi, and every statement about the general case becomes a statement about independent one-dimensional Gaussians: the normalising constant is a product of n Gaussian integrals, the mean is 0 by symmetry, the covariance matrix is Λ rotated back, and positive definiteness is exactly the condition that makes all n integrals converge.
Conditioning needs a second, purely algebraic device. The matrix V factorises as a product of a lower-triangular, a block-diagonal and an upper-triangular matrix - block Gaussian elimination - and reading the factorisation off in the exponent completes the square in x1 with x2 held fixed. What is left over is a quadratic form in x2 alone, which cancels against the marginal density, and the determinant splits in the same way. The conditional law drops out with no integration at all.
Proof(Normalisation, mean and covariance)AdaptedAdapted from Grimmett & Stirzaker (2020), §4.9, the evaluation of $\int K\exp(-\tfrac12 Q)$ and the symmetry argument, p. 132; Riley, Hobson & Bence (2006), §30.15.2, the worked integral (30.149), p. 1210
Let V be symmetric positive definite with eigenvalues λ1,…,λn>0. By the spectral theorem (Spectral theorem for Hermitian matrices) there is an orthogonal Q with V=QΛQT, Λ=diag(λi), whence V−1=QΛ−1QT and detV=λ1⋯λn.
Step 1: the integral. Substitute y=QT(x−μ), so x=μ+Qy. The Jacobian determinant is detQ=±1, so dx=dy. The exponent becomes
So the stated constant makes f integrate to 1. This is the step that needs positive definiteness: if some λi≤0 the corresponding factor diverges (for λi<0) or is undefined (for λi=0), and no constant can repair it.
Step 2: the mean. With Y=QT(X−μ), step 1 shows that the density of Y is ∏i(2πλi)−1/2e−yi2/2λi, a product of N(0,λi) densities. So the Yi are independent, each with mean 0 and variance λi; each is integrable because a Gaussian has all moments. Hence E(Y)=0 and, by linearity, E(X)=μ+QE(Y)=μ. (Grimmett and Stirzaker reach the same conclusion from the symmetry f(x)=f(−x) of the centred density.)
Step 3: the covariance matrix.cov(Y)=Λ, since the Yi are independent with variances λi. Then X−μ=QY, and for any matrix A and vector Z with mean 0, E(AZ(AZ)T)=AE(ZZT)AT. Therefore
E((X−μ)(X−μ)T)=QΛQT=V,
that is cov(Xi,Xj)=vij.
Step 4: the transforms. For the moment generating function, complete the square. With u=x−μ,
as expanding the last bracket confirms, using the symmetry of V−1. Integrating f(x)etTx over x therefore leaves exp(tTμ+21tTVt) times the integral of a density centred at μ+Vt, which is 1. This is the form Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §30.15.2, p. 1209) obtain, with their A=V−1. The characteristic function follows by analytic continuation to imaginary argument, exactly as in the one-dimensional case (Characteristic functions); it is the result Grimmett and Stirzaker (Grimmett & Stirzaker (2020), §5.8, Example (5.8.6), pp. 209-210) obtain by this same rotation.
Step 5: the chi-squared statement. In the coordinates of step 1, (X−μ)TV−1(X−μ)=∑iYi2/λi, and Yi/λi are independent N(0,1). A sum of n squared independent standard normals is χ2(n) by definition. □
Proof(Affine images, marginals, conditionals and independence)AdaptedAdapted from Grimmett & Stirzaker (2020), §4.9, proof of Theorem (4.9.6) for $m=n$, p. 133
(a) Affine images, invertible case. Let D be n×n and invertible, Y=DTX+c. The map x↦y=DTx+c is a bijection of Rn with inverse x=(DT)−1(y−c) and Jacobian determinant ∣detD∣−1, so by the change of variables formula (Multiple integrals and change of variables) the density of Y is f((DT)−1(y−c))∣detD∣−1. Put W=DTVD and ν=DTμ+c. Then W−1=D−1V−1(DT)−1 and, with z=y−ν,
((DT)−1z)TV−1((DT)−1z)=zTD−1V−1(DT)−1z=zTW−1z,
while detW=(detD)2detV, so ∣detD∣−1(detV)−1/2=(detW)−1/2. The density of Y is therefore the N(ν,W) density. Positive definiteness of W is inherited: uTWu=(Du)TV(Du)>0 for u=0 since D is invertible.
(a) Affine images, general rank. For D of rank m<n, compute the characteristic function of Y=DTX+c directly from step 4 above: for s∈Rm,
which is the characteristic function of N(DTμ+c,DTVD); the matrix DTVD is positive definite because Ds=0 for s=0 when D has rank m. By the uniqueness half of the inversion theorem (Characteristic functions) the distributions agree. Grimmett and Stirzaker state this case without proof, referring to another text; the characteristic-function argument is the standard way to supply it.
(b) Marginals. Take D in (a) to be the n×p matrix whose columns are the standard basis vectors ei1,…,eip of the chosen coordinates and c=0. It has rank p, DTμ is the corresponding sub-vector of μ, and DTVD is the corresponding p×p sub-matrix of V.
(c) Conditionals. Write u=x1−μ1, v=x2−μ2 and S=V11−V12V22−1V21. Block Gaussian elimination factorises V as
V=(I0V12V22−1I)(S00V22)(IV22−1V210I),
as multiplying out confirms, the (1,1) block being S+V12V22−1V21=V11. Two consequences follow at once, because the outer factors are triangular with unit diagonal and so have determinant 1 and easily inverted:
Applying the second display to (u,v) and writing w=u−V12V22−1v gives the completed square
(x−μ)TV−1(x−μ)=wTS−1w+vTV22−1v.
(S is positive definite: it is the inverse of the top-left block of V−1, which is positive definite because V is.) The marginal density of X2 is the N(μ2,V22) density, by (b). Dividing, and using detV=detSdetV22 to split the constant,
which is the N(μ1+V12V22−1(x2−μ2),S) density, since w=x1−(μ1+V12V22−1(x2−μ2)). No integration was performed: the whole content is the factorisation of V.
(d) Independence. If V12=0 then S=V11 and the conditional density of X1 given X2=x2 does not depend on x2, so f=fX1fX2 and the blocks are independent. Conversely independent blocks have cov(X1i,X2j)=0, which is V12=0.
(e) Sampling. If W∼N(0,I) and CCT=V then μ+CW has mean μ and covariance CICT=V. Taking determinants in CCT=V gives (detC)2=detV, which is strictly positive because V is positive definite, so C is invertible and (a) applies with DT=C: the law is N(μ,V). (If V is only non-negative definite, C may be singular and (a) no longer applies, but μ+CW is still multivariate normal in the linear-combination sense, since aT(μ+CW)=aTμ+(CTa)TW is univariate normal for every a.) □
Remark(Where the hypotheses are used, and the trap of marginal normality)Original
Positive definiteness is needed for the density and nowhere else. It enters exactly once, in step 1, to make the n Gaussian integrals converge; everything downstream inherits it. Merely non-negative definite V still defines a perfectly good distribution through the characteristic function exp(itTμ−21tTVt), but that distribution is concentrated on an affine subspace of dimension rankV and has no density on Rn. The simplest instance is X2=X1 with X1∼N(0,1): the law of (X1,X2) sits on a line, V=(1111) is singular, and detV=0 kills the normalising constant.
Invertibility of V22, not of V, is what conditioning needs. The factorisation in (c) uses V22−1 and nothing else. If the conditioning block is itself degenerate, the observed coordinates carry a deterministic relation and one of them must be dropped before conditioning.
Jointly normal is strictly stronger than marginally normal.Grimmett and Stirzaker (Grimmett & Stirzaker (2020), §4.9, Exercise (4.9.5), pp. 134-135) set this exercise: let X∼N(0,1), fix a>0, and let Y=X when ∣X∣<a and Y=−X when ∣X∣≥a. Then Y is again N(0,1), since the transformation reflects a symmetric density about the origin on part of its range; yet ∣Y∣=∣X∣ always, so the pair is confined to two lines and has no bivariate normal density. Their ρ(a)=cov(X,Y) runs continuously from −1 (at a=0, where Y=−X) to +1 (as a→∞, where Y=X), so some intermediate a makes X and Y uncorrelated while they remain completely dependent. Statement (d) is therefore about jointly normal variables and collapses without that hypothesis.
Normality of every linear combination is the right hypothesis, and pairwise checks do not suffice. Verifying that X1, X2 and X1+X2 are each normal does not make (X1,X2) bivariate normal; the definition quantifies over alla, and the example above fails it at a=(1,1)T or a=(1,−1)T depending on the cut-off.
The conditional covariance being constant is special to the normal family. For a general joint law var(Y∣X) depends on X, and the best predictor is a curve rather than a line. Here both degeneracies happen at once, which is precisely why least-squares linear regression is exactly optimal under a normal model and only approximately so otherwise.
A check on the conditional formula. For the standard trivariate normal with unit variances and correlations ρ1=ρ(X,Y), ρ2=ρ(Y,Z), ρ3=ρ(Z,X), statement (c) with X1=Z and X2=(X,Y)T gives V12=(ρ3,ρ2) and V22−1=(1−ρ12)−1(1−ρ1−ρ11), hence
E(Z∣X,Y)=1−ρ12(ρ3−ρ1ρ2)X+(ρ2−ρ1ρ3)Y,var(Z∣X,Y)=1−ρ121−ρ12−ρ22−ρ32+2ρ1ρ2ρ3,
which is exactly what Grimmett and Stirzaker (Grimmett & Stirzaker (2020), §4.9, Exercise (4.9.9), p. 135) ask the reader to show. Their ρ2 and ρ3 are labelled cyclically, as that agreement fixes. The bivariate case p=q=1 with unit variances reduces to E(Y∣X)=ρX and var(Y∣X)=1−ρ2, which is their Example (4.6.7) (Grimmett & Stirzaker (2020), §4.6, Example (4.6.7), p. 119).
Method(Working with a multivariate normal)Original
Write down μ and V, and check V. It must be symmetric; it must be positive definite if a density is wanted. Cheap checks: every diagonal entry positive, every 2×2 principal minor positive, detV>0. A full check is Sylvester's criterion or a Cholesky factorisation, which fails exactly when V is not positive definite (Quadratic forms and the Rayleigh quotient).
Reduce to the standard case when you need an integral. Either rotate, Y=QT(X−μ) with Q from the spectral theorem, or whiten, W=C−1(X−μ) with CCT=V. Both turn the problem into one about independent N(0,1) variables.
For any linear function, do not integrate. Read the answer off (a): DTX+c is normal with mean DTμ+c and covariance DTVD. Sums, differences, averages and projections are all of this form.
For a marginal, delete rows and columns. Keep the entries of μ and V belonging to the coordinates retained; discard the rest. No integration.
For a conditional, use the Schur complement. Order the coordinates so that the observed ones come last, form S=V11−V12V22−1V21, and read off the conditional mean μ1+V12V22−1(x2−μ2) and covariance S.
Check. Three identities cost nothing and catch most errors: detV=detS⋅detV22; the conditional variance of any coordinate must be no larger than its unconditional variance; and setting V12=0 must return the marginal.
To simulate, factor once. Compute the Cholesky factor C with CCT=V once, then each sample is μ+CW with W a fresh vector of independent standard normals.
Example(A three-dimensional normal: marginal, conditional and a tail probability)Original
Take n=3 with
μ=1−12,V=421231112.
Step 1: is V a legitimate covariance matrix? It is symmetric. Its leading principal minors are 4>0, det(4223)=8>0 and
so V is positive definite by Sylvester's criterion. (Numerically its eigenvalues are 1.30798, 1.64310, 6.04892; they sum to 9=trV and multiply to 13=detV, two checks at once.) Hence X∼N(μ,V) has the density of the definition, with normalising constant (2π)−3/213−1/2=0.0176100.
Step 2: a marginal. By (b), (X1,X2)T∼N((1−1),(4223)), so X1∼N(1,4), X2∼N(−1,3) and ρ(X1,X2)=2/4⋅3=1/3=0.57735.
Step 3: condition X1 on the other two. Take p=1, q=2, so
Observe x2=(0,3)T, so x2−μ2=(1,1)T and the conditional mean is
μ1+V12V22−1(x2−μ2)=1+51(3+1)=59=1.8.
Therefore
X1∣{X2=0,X3=3}∼N(59,513).
Step 4: the checks. The determinant identity holds: detS⋅detV22=513⋅5=13=detV. The conditional variance 2.6 is smaller than the unconditional variance 4, as observing correlated coordinates must make it. And the conditional mean has moved up from 1 to 1.8, consistent with X2 and X3 both having been observed above their means while both covariances with X1 are positive. A Monte Carlo run of 4×106 draws, keeping those with X2 and X3 within 0.05 of the conditioning values, gave a sample mean 1.819 and sample variance 2.685 from 2106 retained draws, agreeing with 1.8 and 2.6 to within the sampling error.
Step 5: a probability. With the conditional law in hand,
The unconditional answer is Φ((3−1)/2)=Φ(1)=0.841345. Conditioning has therefore made the value 3less surprising, not more: the upper tail P(X1>3) rises from 0.158655 to 0.228375, and the standardised distance from the centre falls from z=1 to z=0.744208. Two effects compete, and the one that wins is the shift of the centre. Observing X2 and X3 above their means pulls the mean of X1 up from 1 to 1.8, which moves 3 closer; it also shrinks the standard deviation from 2 to 13/5=1.612452, which on its own would move 3 further away in units of the spread. The break-even point is the conditional mean that would leave z at 1 with the new spread, namely 3−1.612452=1.387548; the actual conditional mean 1.8 is larger than that, so the mean shift dominates. All numerical values here were computed and cross-checked symbolically.
Application(Error ellipses, least squares, Gaussian fields and filters)Original
Error matrices and confidence regions. A fit reports parameters θ^ with a covariance matrix V, and the standard interpretation is that θ^ is approximately N(θ,V). Statement 4 above then says that (θ^−θ)TV−1(θ^−θ) is χ2(n), which is where contours of constant Δχ2 as confidence regions come from: the 1σ, 2σ ellipsoids are level sets of the quadratic form, with axes along the eigenvectors of V and semi-axes proportional to λi. Marginalising a parameter away is statement (b) - delete its row and column - while fixing it is statement (c), and the two give different error bars on the remaining parameters, the conditional one always smaller. Confusing the two is a common source of understated uncertainties.
Linear least squares. With normal errors, the maximum-likelihood estimate minimises a quadratic form and is itself a linear function of the data, hence normal by (a); its covariance matrix is the usual (ATΣ−1A)−1. Statement (c) is the reason linear regression is optimal and not merely convenient in this setting.
Gaussian random fields. A field ϕ(x) whose finite-dimensional distributions are multivariate normal is specified completely by a mean function and a two-point correlation function; this is the standard model for the cosmic microwave background, for turbulence at fixed scale, and for the Gaussian free field. Conditioning on measurements at some points to predict the field elsewhere is statement (c) applied to the covariance kernel - in geostatistics and machine learning the same formula is called kriging, or Gaussian process regression.
Kalman filtering and data assimilation. Because a normal prior through a linear model with normal noise gives a normal posterior, the entire filtering recursion reduces to updating a mean vector and a covariance matrix by statements (a) and (c). This is why the Kalman filter is finite-dimensional while general nonlinear filtering is not.
Equilibrium statistical mechanics. A system of harmonic degrees of freedom at temperature T has Boltzmann weight exp(−21qTKq/kBT), which is exactly the N(0,kBTK−1) density. Equipartition, the fluctuation-dissipation relation between the response matrix K−1 and the covariance of the displacements, and normal-mode decomposition are statements 2, 3 and the rotation of step 1 respectively. The Ornstein-Uhlenbeck process of Stochastic differential equations has multivariate normal finite-dimensional laws for the same reason.
Encodes a distribution in its probability, moment and cumulant generating functions, turns convolution into multiplication and random sums into composition, and applies the composition rule to branching processes, extracting the extinction probability 1/2 for a worked offspring law.
A distribution on {0,1,2,…} is an infinite list of numbers. Packing that list into the coefficients of a power series replaces three awkward operations by three easy ones:
operation on distributions
operation on generating functions
convolution (adding independent variables)
multiplication
summing a random number of copies
composition
extracting moments
differentiating at one point
The third is a convenience. The first two are the reason the method exists: a convolution is a sum over all the ways a total can be split, and a power series multiplies by doing exactly that bookkeeping. Everything in this unit, including the entire theory of branching processes, follows from those two lines of the table.
Three transforms are in use, and the difference between them is only a change of variable. The probability generating functionG(s)=E(sX) handles variables on the non-negative integers. The moment generating functionM(t)=E(etX), obtained by putting s=et, handles any variable but may be infinite. The characteristic functionφ(t)=E(eitX), which exists always, is the subject of Characteristic functions.
Definition(The three generating functions)AdaptedAdapted from Grimmett & Stirzaker (2020), §5.1, Definition (5.1.10) and the change of variable at eq. (5.1.21), pp. 170, 172, and §5.7, Definition (5.7.1), p. 203; Riley, Hobson & Bence (2006), §30.7.1, eq. (30.71), p. 1158, §30.7.2, eq. (30.85), p. 1162, and §30.7.4, eq. (30.92), p. 1166; Arfken, Weber & Harris (2013), §23.2, eq. (23.32), p. 1141
Let X be a random variable with mass function f or density f.
Probability generating function. If X takes values in {0,1,2,…}, its probability generating function is
GX(s)=E(sX)=n=0∑∞snP(X=n).
The series has non-negative coefficients summing to 1, so it converges absolutely for ∣s∣≤1 at least; its radius of convergence R satisfies R≥1. Riley, Hobson and Bence write it ΦX(t) and Grimmett and Stirzaker GX(s); the object is the same.
Moment generating function. For any real random variable,
defined at those real t for which the sum or integral is finite. Always MX(0)=1. Where both exist, MX(t)=GX(et). Physicists write MX(t)=⟨etX⟩, as Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §23.2, eq. (23.32), p. 1141) do. Grimmett and Stirzaker note that MX is a Laplace transform of the distribution up to the sign of the exponent: the usual Laplace transform of a density carries e−tx, this one carries e+tx.
Cumulant generating function. Wherever MX is finite and positive on a neighbourhood of the origin,
KX(t)=logMX(t)=κ1t+κ22!t2+κ33!t3+⋯,
and the coefficients κk are the cumulants of X. There is no constant term because MX(0)=1.
Joint generating functions. For two variables on the non-negative integers, GX1,X2(s1,s2)=E(s1X1s2X2), and likewise MX1,X2(t1,t2)=E(et1X1+t2X2) for general variables.
Method(The rules)AdaptedAdapted from Grimmett & Stirzaker (2020), §5.1, eq. (5.1.16), Theorem (5.1.18) with eq. (5.1.19), the uniqueness statement (5.1.13)-(5.1.14), Theorems (5.1.23), (5.1.25) and (5.1.29), pp. 170-174, and §5.7, properties (a)-(c) and Remark (5.7.10), pp. 203, 206; Riley, Hobson & Bence (2006), §30.7.1, eqs. (30.74)-(30.76) and (30.80)-(30.84), pp. 1158-1162, and §30.7.2, eqs. (30.86)-(30.90), pp. 1163-1164
Let X, Y be random variables and let a,b be real constants.
For probability generating functions (X,Y on the non-negative integers):
Boundary values.GX(0)=P(X=0) and GX(1)=1.
Recovery.P(X=n)=GX(n)(0)/n!. In particular GX determines the distribution: two variables with the same probability generating function on some interval (−R′,R′) with R′>0 have the same mass function.
Factorial moments.E(X(X−1)⋯(X−k+1))=GX(k)(1), where G(k)(1) means lims↑1G(k)(s) when the radius of convergence is exactly 1. Hence
E(X)=GX′(1),var(X)=GX′′(1)+GX′(1)−GX′(1)2.
Sums. If X and Y are independent then GX+Y(s)=GX(s)GY(s); for n independent summands the generating functions multiply, and for n independent copies of one variable GSn=(GX)n.
Random sums. Let X1,X2,… be independent, identically distributed with common generating function GX, and let N≥0 be an integer-valued random variable independent of all the Xi, with generating function GN. Then S=X1+⋯+XN (an empty sum being 0) has
GS(s)=GN(GX(s)).
Grimmett and Stirzaker record the mnemonic that fixes the order of composition: if N≡n is constant then GN(s)=sn and the formula returns GX(s)n, as rule 4 demands.
Independence. The factorisation GX1,X2(s1,s2)=GX1(s1)GX2(s2), holding at every pair (s1,s2), is equivalent to independence of X1 and X2 - not merely implied by it.
Defective variables. If P(X=∞)>0 then GX(s)=E(sX) still converges for ∣s∣<1 and lims↑1GX(s)=1−P(X=∞)<1. A probability generating function with G(1−)<1 is the signature of a defective distribution.
For moment generating functions (any real X, wherever the functions are finite):
Moments. If MX is finite on an open interval containing 0, then all moments exist, MX has the convergent expansion MX(t)=∑k≥0E(Xk)tk/k! there, and E(Xk)=MX(k)(0); in particular var(X)=MX′′(0)−MX′(0)2.
Affine changes.MaX+b(t)=ebtMX(at). Taking b=−μ with μ=EX produces the generator of the moments about the mean.
Sums. If X1,…,Xm are independent, with m a fixed integer, then M∑i≤mciXi(t)=∏i≤mMXi(cit). For a random number of summands, with N and the identically distributed Xi as in rule 5, MS(t)=GN(MX(t)).
Cumulants.KX+Y=KX+KY for independent X,Y, and κ1=E(X), κ2=var(X). For X∼N(μ,σ2), KX(t)=μt+21σ2t2, so κ1=μ, κ2=σ2, and every cumulant from the third on vanishes.
Uniqueness, with a caveat. If MX=MY and both are finite on a neighbourhood of 0, then X and Y have the same distribution. Riley, Hobson and Bence state this and say a rigorous proof lies outside their scope. Grimmett and Stirzaker add the warning that the moments alone do not always determine a distribution: there exist distinct distribution functions all of whose moments agree, the standard example being built from the log-normal. Finiteness of M near the origin is their simplest sufficient condition for a distribution to be determined by its moments.
Method(Branching processes)AdaptedAdapted from Grimmett & Stirzaker (2020), §5.4, the modelling assumptions (a)-(b), Theorem (5.4.1), Lemma (5.4.2) and Theorem (5.4.5), pp. 192-195
A branching process is a population in which each individual, independently of all others, produces a random number of offspring with a common mass function f and probability generating function G. Write Zn for how many individuals generation n contains, take Z0=1, and put Gn(s)=E(sZn), μ=E(Z1)=G′(1) and σ2=var(Z1).
Iteration.Gm+n(s)=Gm(Gn(s)), so Gn is the n-fold composite of G with itself and G1=G.
Extinction. As n→∞, P(Zn=0) increases to the probability η of ultimate extinction, and η is the smallest non-negative root of s=G(s). Moreover η=1 if μ<1, η<1 if μ>1, and in the critical case μ=1 extinction is certain, η=1, so long as the offspring law is not concentrated at a single value, that is, so long as σ2>0.
Idea of the proofOriginal
Why a product? Expanding GX(s)GY(s) and collecting the coefficient of sn gives ∑kP(X=k)P(Y=n−k), which is precisely the number of ways X+Y can equal n weighted by their probabilities - the convolution. The power series does the sum over splittings automatically, because that is what multiplying polynomials means. Equivalently and more quickly: sX and sY are independent when X and Y are, so the expectation of their product factorises.
Why a composition? For a random sum, condition on N. Given N=n the sum has generating function GX(s)n, so averaging over n gives ∑nGX(s)nP(N=n) - and that sum isGN evaluated at the number GX(s). Substituting a generating function into a generating function is the same operation as randomising the number of terms.
Branching processes are then immediate: generation n+m is a random sum over the individuals of generation m, each contributing an independent copy of an n-generation descendant count, so generation sizes compose. Extinction is the fixed point of that composition, and the geometry of a convex increasing curve meeting the diagonal decides whether the fixed point is 1 or something smaller - which is why the mean family size μ=G′(1), the slope at the diagonal point s=1, is what separates the cases.
Proof(The rules)AdaptedAdapted from Grimmett & Stirzaker (2020), §5.1, the convolution computation (5.1.6), the proof of Theorem (5.1.18) by Abel's theorem, eq. (5.1.19), and the proofs of Theorems (5.1.23) and (5.1.25), pp. 169-174; Riley, Hobson & Bence (2006), §30.7.1, eqs. (30.74)-(30.76), p. 1158
Rules 1 and 2. Setting s=0 in the series leaves the constant term, P(X=0); setting s=1 sums the mass function, giving 1. Rule 2 is the standard uniqueness statement for power series: inside the disc of convergence a power series may be differentiated term by term any number of times, and evaluating the nth derivative at 0 kills every term but one (Power series).
Rule 3. For ∣s∣<1, differentiating k times term by term,
GX(k)(s)=i≥k∑i(i−1)⋯(i−k+1)si−kP(X=i).
Every coefficient here is non-negative, so Abel's theorem applies: as s↑1 the sum tends to ∑ii(i−1)⋯(i−k+1)P(X=i), whether that is finite or +∞. That sum is E(X(X−1)⋯(X−k+1)) by the law of the unconscious statistician (Expectation, variance and covariance). Abel's theorem is exactly what is needed here, because when R=1 one has no right to evaluate the series at the endpoint without it. For the variance, take k=1,2:
Rule 4. If X and Y are independent then so are the random variables sX and sY, being functions of independent variables, and the expectation of a product of independent variables is the product of the expectations:
GX+Y(s)=E(sX+Y)=E(sXsY)=E(sX)E(sY).
The same conclusion comes from multiplying the two series and recognising the coefficient of sn as the convolution ∑k=0nP(X=k)P(Y=n−k), which is P(X+Y=n); Riley, Hobson and Bence present the argument in that order, justifying the interchange of the two summations by the geometry of the index set. Independence is used once and is indispensable: see the remark.
where the inner conditional expectation was replaced by an unconditional one because N is independent of the Xi. By rule 4 applied n−1 times the nth term is GX(s)nP(N=n), and summing over n gives ∑nGX(s)nP(N=n)=GN(GX(s)), since 0≤GX(s)≤1 for 0≤s≤1 puts the argument inside the disc of convergence of GN.
Rule 6. If X1,X2 are independent the argument of rule 4 applies to s1X1 and s2X2. Conversely, if the joint generating function factorises then comparing the coefficient of s1is2j on the two sides gives P(X1=i,X2=j)=P(X1=i)P(X2=j) for every i,j, which is independence.
Rule 7. The series ∑nsnP(X=n) omits the value +∞ entirely, since s∞=0 for ∣s∣<1; letting s↑1 and applying Abel's theorem to the non-negative coefficients gives ∑n<∞P(X=n)=1−P(X=∞).
Rules 8-11. Rule 8 is the substitution s=et together with the expansion
MX(t)=E(k≥0∑k!(tX)k)=k≥0∑k!tkE(Xk),
the interchange of expectation and sum being justified when MX is finite on a neighbourhood of 0; this is the derivation Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §23.2, eq. (23.32), p. 1141) give. Rule 9 is E(et(aX+b))=ebtE(e(at)X). Rule 10 repeats the argument of rules 4 and 5 with etX in place of sX. Rule 11 follows from rule 10 by taking logarithms, and the Gaussian values from K(t)=logexp(μt+21σ2t2), which is already a polynomial, so the coefficients beyond t2 are zero. □
Proof(The branching process theorems)AdaptedAdapted from Grimmett & Stirzaker (2020), §5.4, proofs of Theorem (5.4.1), Lemma (5.4.2) and Theorem (5.4.5), pp. 193-195
Iteration. Each member of generation m+n descends from exactly one member of generation m. Writing Wi for the number of generation-(m+n) descendants of the ith member of generation m, we have Zm+n=W1+⋯+WZm. The Wi are independent, by the independence assumption on family sizes; they are identically distributed, each with the law of Zn started from a single ancestor; and Zm is a function of the family sizes in generations 0 to m−1 only, whereas each Wi is a function of the family sizes of generation m onwards, so Zm is independent of the whole sequence W1,W2,…. Rule 5 therefore applies and gives Gm+n=Gm∘Gn. Taking m=1 and iterating, Gn=G∘G∘⋯∘G (n factors).
Moments. Differentiate Gn(s)=G(Gn−1(s)) by the chain rule and let s↑1. Since Gn−1(1)=1,
Gn′(1)=G′(1)Gn−1′(1)=μGn−1′(1),
so E(Zn)=μn by induction from E(Z0)=1. Differentiating once more,
Write vn=var(Zn)=Gn′′(1)+μn−μ2n, substitute Gn′′(1)=vn−μn+μ2n and G′′(1)=σ2−μ+μ2 into the last display, and the terms in μn, μ2n−1 and μ2n all cancel, leaving the one-term recursion
vn=μvn−1+σ2μ2n−2,v0=0.
Unrolling it gives vn=σ2(μ2n−2+μ2n−3+⋯+μn−1)=σ2μn−1(1+μ+⋯+μn−1), which is nσ2 when μ=1 and σ2μn−1(μn−1)/(μ−1) otherwise.
Extinction. The events {Zn=0} increase with n, because an empty generation has empty descendants; so ηn=P(Zn=0) is non-decreasing and, by continuity of probability along increasing sequences, ηn↑η=P(some Zn=0)=P(extinction). Now ηn=Gn(0)=G(Gn−1(0))=G(ηn−1), and G is continuous on [0,1], so letting n→∞ gives
η=G(η):
the extinction probability is a fixed point. It is the smallest non-negative one: if ψ≥0 satisfies ψ=G(ψ), then, because G has non-negative coefficients and so is non-decreasing on [0,1], η1=G(0)≤G(ψ)=ψ, and inductively ηn=G(ηn−1)≤G(ψ)=ψ; letting n→∞ gives η≤ψ.
For the three cases, note that G′′(s)=E(Z1(Z1−1)sZ1−2)≥0 on [0,1], so G is convex there, and G(1)=1. A convex function on [0,1] that equals 1 at s=1 meets the diagonal y=s in at most one further point.
If μ=G′(1)>1 then G(s)<s for s slightly below 1, while G(0)=P(Z1=0)≥0; so G(s)−s changes sign on [0,1) and there is a root η<1.
If μ<1 then G(s)>s for all s∈[0,1) - by convexity the chord from (s,G(s)) to (1,1) has slope at most G′(1)=μ<1, which forces G(s)≥1−μ(1−s)>s - so the only fixed point in [0,1] is 1 and η=1.
If μ=1 and σ2>0 then G is strictly convex (strict convexity is exactly σ2>0 here, since σ2=0 with μ=1 means Z1≡1 and G(s)=s), so G(s)>s on [0,1) and again η=1. In the degenerate case σ2=0, μ=1, every individual has exactly one child, G(s)=s, every point is a fixed point and the population never dies: the smallest non-negative root is 0, which is the correct extinction probability. □
Remark(Where the hypotheses are used, and where the method breaks)Original
Rule 4 needs independence, not merely zero correlation. Let X be Bernoulli with parameter 21 and Y=1−X. Then GX(s)=GY(s)=21(1+s), so the product is 41(1+s)2, the generating function of a bin(2,21) variable. But X+Y≡1, with generating function s. The two disagree everywhere except at s=1.
Rule 5 needs N independent of the Xi, and the proof shows exactly why. The step that replaces E(sX1+⋯+Xn∣N=n) by E(sX1+⋯+Xn) is the only use of that hypothesis, and it fails if the number of terms is chosen after looking at them. If X1,X2,… are independent Bernoulli(21) and N=1+X1, then S=X1+⋯+XN can be enumerated in four equally likely cases. If X1=0 then N=1 and S=X1=0; if X1=1 then N=2 and S=1+X2, which is 1 or 2 according to X2. Hence P(S=0)=21, P(S=1)=41, P(S=2)=41 and GS(s)=21+41s+41s2, whereas GN(u)=21(u+u2) gives GN(GX(s))=21(21+s+4(1+s)2)=83+21s+81s2. The two agree at s=1, as both must, and nowhere else: the formula is wrong here because N was chosen after looking at X1.
The moment generating function can fail to exist.MX(t)=∞ for every t=0 when X is Cauchy, and for every t>0 when X is log-normal or has any heavy tail; rules 8-12 are then empty, and the characteristic function of Characteristic functions is the right tool. Riley, Hobson and Bence note that boundedness of X guarantees existence for all t; the honest general statement is that MX is finite on an interval containing 0, possibly the single point {0}.
Moments do not always determine a distribution. Grimmett and Stirzaker state this explicitly and point to the log-normal: there are distinct distribution functions with identical moments of every order. So "matching all the moments" is not by itself a proof that two laws agree; what is needed is rule 12's hypothesis that M is finite near the origin, or a characteristic-function argument.
A terminology clash to be aware of.Arfken, Weber and Harris (Arfken, Weber & Harris (2013), §23.2, p. 1141) define the higher moments about the mean in their eq. (23.31), then introduce ⟨etX⟩ and remark that the moments it generates are taken relative to x=0 rather than to the expectation value, calling those quantities central moments. Everywhere else in this handbook, and in Expectation, variance and covariance, central moment means a moment about the mean, E((X−EX)k), and MX as defined above generates the moments about the origin; rule 9 with b=−μ is what converts between them. We follow the handbook's usage and record Arfken, Weber and Harris's wording so that a reader moving between the two is not caught out.
The critical branching case is genuinely delicate. With μ=1 and σ2>0 the population dies out with probability 1 even though E(Zn)=1 for every n; the mean is held up by ever rarer, ever larger survivals, and var(Zn)=nσ2 grows without bound. This is why the theorem must say "provided the variance is strictly positive": the one excluded case, Z1≡1, has μ=1 and never dies.
Abel's theorem is doing real work in rule 3. For a variable with E(X)=∞ - say P(X=n)∝n−2 - the radius of convergence is exactly 1, G is finite on [0,1], and G′(s)→+∞ as s↑1. Rule 3 reports this correctly as E(X)=+∞ only because Abel's theorem permits the limit to be taken.
Method(Using a generating function)Original
Choose the transform. Non-negative integer values and questions about probabilities: use G(s)=E(sX). Real values, or questions about moments: use M(t)=E(etX), after checking it is finite near 0. Heavy tails or a limit theorem: use the characteristic function instead.
Build it from the pieces. Write the quantity of interest as a sum of independent variables (multiply), as a sum of a random number of them (compose), or as an affine function (rule 9). Almost never sum the defining series directly.
Recognise the answer. Keep a short table: Bernoulli(p) gives q+ps; bin(n,p) gives (q+ps)n; Poisson(λ) gives eλ(s−1); geometric on {1,2,…} gives ps/(1−qs); N(μ,σ2) has M(t)=exp(μt+21σ2t2); exponential(λ) has M(t)=λ/(λ−t) for t<λ. By uniqueness, matching one of these identifies the distribution.
Extract what you need. Probabilities: P(X=n)=G(n)(0)/n!, or read off the coefficient of sn in a closed form expanded as a series. Moments: differentiate at s=1 (factorial moments) or at t=0 (ordinary moments). Cumulants: expand logM.
Check.G(1)=1 and M(0)=1 always; G(0) must equal P(X=0); G′(1) must be the mean you expect; and in a composite model the answer must reduce to the known one when the randomness is switched off (for example, N constant in rule 5).
For a branching process, write down G, compute μ=G′(1), and solve s=G(s) on [0,1]; the smallest non-negative root is the extinction probability, and μ≶1 tells you in advance whether it is 1.
Example(A branching process with a quadratic offspring generating function)Original
A neutron entering a block of fissile material is absorbed, and the absorption produces 0, 1 or 2 further neutrons with probabilities
p0=41,p1=41,p2=21.
Each neutron behaves independently and in the same way, so the generations of neutrons form a branching process with Z0=1 and
The five coefficients sum to 3211+324+329+324+324=1, as they must, and they are the mass function of Z2: for instance P(Z2=3)=81. The degree is 4, correct because at most two neutrons each produce at most two.
Step 3: moments of later generations, two ways. From G2 directly, G2′(1)=81+2⋅329+3⋅81+4⋅81=1625, and G2′′(1)+G2′(1)−G2′(1)2=256495. From the general formulas,
with roots s=21 and s=1. The smallest non-negative root is
η=21.
So a single initial neutron sets off a chain reaction that dies out with probability 21 and grows without bound with probability 21. This is consistent with μ>1, which guarantees η<1.
Step 5: a check on the convergence. The theorem says ηn=Gn(0) increases to η. Iterating ηn=G(ηn−1) from η0=0:
increasing towards 0.5, as predicted. Note how slowly: the convergence rate is G′(η)=41+η=43 per generation. (The value η2=3211 agrees with the constant term of G2 computed in step 2, an independent check. All the numbers in this example were verified symbolically.)
Step 6: what happens if the material is made slightly less reactive. Replace p2 by 41 and p1 by 21: then G(s)=41+21s+41s2, μ=21+21=1, and σ2=21+1−1=21>0. The fixed-point equation becomes 4s=1+2s+s2, that is (s−1)2=0, with the single root s=1: extinction is now certain, although the expected population size is exactly 1 in every generation. The transition between η=21 and η=1 is entirely controlled by μ crossing 1 - the criticality condition of a chain reaction.
Application(Chain reactions, queues, polymers and epidemics)Original
Nuclear chain reactions. The example above is the standard neutron-multiplication model: μ is the effective multiplication factor k of reactor physics, and the three regimes μ<1, μ=1, μ>1 are subcritical, critical and supercritical. The extinction probability η is the chance that a single stray neutron fails to start a divergent chain, which is why reactor startup is a statistical rather than a deterministic event and why source neutrons are deliberately supplied. Electromagnetic and hadronic cascades in detectors, and photon multiplication in a photomultiplier dynode chain, are modelled the same way, with the generating function of the per-stage gain composed once per stage.
Compound Poisson totals. Rule 5 with N Poisson is the compound Poisson model: total charge deposited by a Poisson number of independent scattering events, total claim size in insurance, total shot noise in a photodiode. The random sum formula gives the transform immediately, and expanding logMS gives the cumulants of the total in one line - the basis of Campbell's theorem for shot noise.
Queues and recurrence times. The equilibrium queue length of an M/G/1 queue is found by writing a balance relation between generating functions and solving it; the first-return time of a random walk is obtained in the same way, and Random walks uses the convolution rule to turn the renewal relation between return probabilities into an algebraic one.
Polymers and percolation. The size distribution of a branched polymer cluster is a branching process in the mean-field approximation, and the gel point - the onset of an infinite connected cluster - is precisely the critical point μ=1. The same calculation gives the percolation threshold on a tree and the giant-component threshold in a random network, where μ is the mean number of further edges reached from an edge.
Epidemics. In the early stage of an outbreak, infections form a branching process whose mean offspring number is the basic reproduction number R0. The extinction theorem says that an outbreak started by one case fizzles out with probability 1 if R0≤1, and with probability η<1 given by the smallest root of s=G(s) if R0>1 - a statement about the whole offspring distribution, not just its mean, which is why superspreading (large σ2 at fixed μ) raises the chance of early extinction while leaving R0 unchanged.
Defines φX(t)=E(eitX), proves that it is bounded, uniformly continuous and non-negative definite and that it turns independent sums into products and affine maps into phases, states the inversion and continuity theorems and the moment expansion, and inverts sin2t/t2 to the triangular density as a worked check.
The moment generating function M(t)=E(etX) of Generating functions does everything one wants - until the integral diverges, which it does for the Cauchy distribution at every t=0 and for any heavy-tailed law on half the line. Replacing the real exponential by eitX cures this at a stroke: ∣eitX∣=1, so the expectation exists for every random variable and every real t, with no hypothesis whatever. Riley, Hobson and Bence (Riley, Hobson & Bence (2006), §30.7.3, p. 1166) give exactly this as the reason for preferring the characteristic function.
What is gained is not only safety. The characteristic function is the Fourier transform of the distribution - Folland (Folland (1999), §10.1, p. 314) lists the two names as translations of one another in his dictionary between analysts' and probabilists' usage - so the whole apparatus of Fourier inversion becomes available: the distribution can be recovered from φ, convolution becomes multiplication, and convergence of distributions becomes pointwise convergence of transforms. That last equivalence, the continuity theorem, is what makes the proof of the central limit theorem three lines long.
One warning about words. In analysis, "characteristic function" means the indicator function χE of a set; Folland's dictionary records the clash. In this chapter the term always means E(eitX).
Definition(The characteristic function)AdaptedAdapted from Grimmett & Stirzaker (2020), §5.7, Definition (5.7.2) with the surrounding discussion and Definition (5.7.7), pp. 204-205; Riley, Hobson & Bence (2006), §30.7.3, eq. (30.91) and the relation $C_X(t)=M_X(it)$, pp. 1165-1166
The characteristic function of a real random variable X is the function φX:R→C,
φX(t)=E(eitX)=E(costX)+iE(sintX),
the second expression being the definition of the expectation of a complex-valued random variable and showing that φX(t) is defined for every t, since costX and sintX are bounded. Equivalently φX(t)=∫ReitxdF(x), which reads ∑jeitxjP(X=xj) in the discrete case and ∫eitxf(x)dx in the continuous case. Riley, Hobson and Bence write it CX(t) and record the relation to the moment generating function,
φX(t)=MX(it)wherever MX is defined,
so that any identity for M transfers by substituting it for t - with the caveat of the remark below that this substitution needs justification, since M may be finite nowhere but the origin.
Relation to the Fourier transform. If X has density f then, in the convention of The Fourier transform where f^(k)=∫f(x)e−ikxdx,
The probabilists' sign convention is the opposite of the handbook's forward transform, and this is the only difference. Reed and Simon carry an extra (2π)−n/2: their transform of a measure is μ^(λ)=(2π)−n/2∫e−iλxdμ(x), so for the distribution of X in one dimension μ^(λ)=(2π)−1/2φX(−λ).
Joint characteristic function. For random variables X and Y on the same probability space,
φX,Y(s,t)=E(eisX+itY),
and similarly for a random vector X on Rn, φX(t)=E(eitTX) for t∈Rn.
Non-negative definiteness. A function ϕ:R→C is non-negative definite if
j,k=1∑nϕ(tj−tk)zjzk≥0
for every n, all real t1,…,tn and all complex z1,…,zn. Reed and Simon call a bounded continuous ϕ with this property a function of positive type, and record three immediate consequences of the definition: ϕ(0)≥0, ϕ(−x)=ϕ(x) and ∣ϕ(x)∣≤ϕ(0).
Theorem(Basic properties)AdaptedAdapted from Grimmett & Stirzaker (2020), §5.7, Theorem (5.7.3)(a)-(c) with the attribution of the converse to Bochner, Theorems (5.7.5) and (5.7.6), Theorem (5.7.8) with eq. (5.7.9), pp. 204-206, and Problem (5.12.22), p. 234; Reed & Simon (1980), §IX.2, Theorem IX.9, p. 330
Let X, Y be random variables and a,b∈R.
Normalisation and bound.φX(0)=1 and ∣φX(t)∣≤1 for every t.
Uniform continuity.φX is uniformly continuous on R - for everyX, with no moment hypothesis.
Non-negative definiteness.φX is non-negative definite in the sense above.
Bochner's theorem (the converse). Properties 1-3 are not merely necessary but sufficient: whenever ϕ is continuous, non-negative definite and normalised by ϕ(0)=1, some random variable has ϕ for its characteristic function - and by the uniqueness statement below that variable's distribution is the only one with this ϕ. Grimmett and Stirzaker state this and give no proof. Reed and Simon (Reed & Simon (1980), §IX.2, Theorem IX.9, p. 330) prove the measure-theoretic version: on Rn, the cone they call the functions of positive type is precisely the image of the finite positive measures under their Fourier transform. Translating by μ^(λ)=(2π)−1/2φ(−λ) and normalising at the origin turns their theorem into the statement above.
Affine images.φaX+b(t)=eitbφX(at).
Independent sums. If X and Y are independent then φX+Y(t)=φX(t)φY(t).
Independence criterion. Independence of X and Y is equivalent to the factorisation φX,Y(s,t)=φX(s)φY(t) holding at every pair of real arguments s,t. It is not enough that φX+Y(t)=φX(t)φY(t) for all t, which is the statement with s=t only; Grimmett and Stirzaker flag this explicitly and set the construction of a counterexample as an exercise.
Symmetry.φX(−t)=φX(t) always; and X has the same distribution as −X if and only if φX is real-valued, that is ImφX≡0.
Stability of the class. If φ is a characteristic function then so are φ, φ2, ∣φ∣2 and Reφ; but ∣φ∣ need not be one.
Theorem(Moments)AdaptedAdapted from Grimmett & Stirzaker (2020), §5.7, Theorem (5.7.4)(a)-(b), p. 205, Remark (5.7.10) and Theorem (5.7.12), p. 206, and Problem (5.12.30), p. 234; Riley, Hobson & Bence (2006), §30.7.3, the moment formula, p. 1166
Let φ=φX and let k≥1 be an integer.
(a) From derivatives to moments. If φ(k)(0) exists then E∣Xk∣<∞ when k is even, and E∣Xk−1∣<∞ when k is odd. Note what this does and does not promise: an odd derivative buys only the moment one order below it.
(b) From moments to the expansion. If E∣Xk∣<∞ then
φ(t)=j=0∑kj!E(Xj)(it)j+o(tk)(t→0),andφ(k)(0)=ikE(Xk).
Riley, Hobson and Bence write the second relation as E(Xn)=(−i)nφ(n)(0), which is the same statement.
(c) Joint moments.im+nE(XmYn)=∂sm∂tn∂m+nφX,Ys=t=0, subject to the appropriate differentiability.
(d) The moment problem. The moments do not always determine the distribution: Grimmett and Stirzaker record that distinct distribution functions can have identical moments of every order, the usual example being built from the log-normal, and that no simple necessary and sufficient condition is known. The simplest condition they offer as sufficient: M(t)=E(etX) should be finite throughout some neighbourhood of the origin. They also state, as a theorem whose proof they omit, that for a>0 the three conditions "∣M(t)∣<∞ for ∣t∣<a", "φ is analytic on the strip ∣Imz∣<a" and "limsupk(∣mk∣/k!)1/k≤1/a" are equivalent, and that under any of them M extends analytically with φ(t)=M(it).
Grimmett and Stirzaker state (a) and (b) without proof, remarking that it is Taylor's theorem for a function of a complex variable and referring the reader to other texts; we follow them.
Theorem(Inversion and continuity)AdaptedAdapted from Grimmett & Stirzaker (2020), §5.9, Theorem (5.9.1) with the integrability condition, Inversion theorem (5.9.2), Corollary (5.9.3) with its footnote, the two-dimensional formula, and Continuity theorem (5.9.5), pp. 212-214; Riley, Hobson & Bence (2006), §30.7.3, the inverse transform, p. 1166; Folland (1999), §8.3, Theorem 8.26 and Corollary 8.27, pp. 252-253
Inversion, density form. If X is continuous with density f and characteristic function φ, then
f(x)=2π1∫−∞∞e−itxφ(t)dt
at every x at which f is differentiable. This is the Fourier inversion theorem of Fourier inversion written in probabilists' notation; where the integral does not converge absolutely it is read as a principal value. Absolute integrability, ∫−∞∞∣φ(t)∣dt<∞, suffices to place φ among the transforms of variables possessing a density; it is not required. Folland (Folland (1999), §8.3, Theorem 8.26, p. 252) states the underlying L1 theorem with the hypothesis that both the function and its transform are integrable, the conclusion holding almost everywhere.
Inversion, general form. Let X have distribution function F and characteristic function φ, and define F(x)=21(F(x)+limy↑xF(y)), the midpoint value at a jump. Then for all a<b,
F(b)−F(a)=N→∞lim∫−NN2πite−iat−e−ibtφ(t)dt.
No assumption is made about the type of X: this holds for discrete, continuous and mixed distributions alike.
Uniqueness. Agreement of two characteristic functions at every real argument is equivalent to agreement of the two distribution functions everywhere. Grimmett and Stirzaker add a footnote worth keeping: agreement of φX and φY on a bounded interval does not force FX=FY, and they give an example.
Joint inversion. If X,Y have joint density f and joint characteristic function φ, then at every point of differentiability of f,
f(x,y)=4π21∬R2e−isxe−ityφ(s,t)dsdt,
from which the independence criterion 7 follows.
Continuity theorem. Let F1,F2,… be distribution functions with characteristic functions φ1,φ2,….
(a) Suppose Fn→F, where F is a distribution function whose characteristic function is φ. Then φn(t)→φ(t) at every real t.
(b) Conversely, suppose the pointwise limit φ(t)=limn→∞φn(t) exists at every real t, and that this limit function is continuous at the origin. Then there is a distribution function F whose characteristic function is φ, and Fn→F.
Here Fn→F means convergence at every continuity point of F, the mode of Convergence of random variables. Grimmett and Stirzaker refer the proofs of the inversion and continuity theorems to other texts; we do the same.
Theorem(The standard characteristic functions)AdaptedAdapted from Grimmett & Stirzaker (2020), §5.8, Examples (5.8.1)-(5.8.5) and (5.8.8), pp. 208-211, for the Bernoulli, binomial, exponential, Cauchy, normal, gamma and chi-squared rows; Grimmett & Stirzaker (2020), §5.8, Example (5.8.6), pp. 209-210, for $N(0,V)$, extended to $N(\mu,V)$ by Theorem (5.7.6); Grimmett & Stirzaker (2020), §5.1, Examples (5.1.17)(a) and (d), p. 171, for the point-mass and Poisson probability generating functions $s^{c}$ and $e^{\lambda(s-1)}$, converted by $\varphi_X(t)=G_X(e^{it})$
distribution
density or mass function
φ(t)
point mass at c
P(X=c)=1
eitc
Bernoulli(p)
P(X=1)=p, P(X=0)=q
q+peit
bin(n,p)
(kn)pkqn−k
(q+peit)n
Poisson(λ)
e−λλk/k!
exp(λ(eit−1))
exponential(λ)
λe−λx, x≥0
λ/(λ−it)
gamma Γ(λ,s)
λsxs−1e−λx/Γ(s), x≥0
(λ/(λ−it))s
χ2(d)=Γ(21,2d)
(1−2it)−d/2
N(0,1)
(2π)−1/2e−x2/2
e−t2/2
N(μ,σ2)
exp(iμt−21σ2t2)
Cauchy
1/(π(1+x2))
e−∣t∣
N(μ,V) on Rn
exp(itTμ−21tTVt)
Two entries deserve comment. The Cauchy characteristic function e−∣t∣ has a corner at the origin: its one-sided derivatives there are +1 from the left and −1 from the right, so φ′(0) does not exist. Statement (b) read backwards then confirms that the Cauchy law has no mean, since E∣X∣<∞ would have forced φ′(0)=iE(X) to exist. The implication runs in that direction only: statement (a) at k=1 is vacuous (see the remark below), so the absence of a mean does not by itself prevent a characteristic function from being differentiable at 0, and the corner here has to be read off the formula rather than deduced from the missing moment. And the normal entry is obtained, as Grimmett and Stirzaker emphasise, not by substituting s=it into M(s)=es2/2 as if i were real, but by analytic continuation, which their Remark (5.7.11) licenses; they warn twice against treating i as a real number, noting that for the exponential distribution the illegitimate substitution happens to give the right answer and for the Cauchy distribution it gives nothing at all.
Idea of the proofOriginal
Each of the three basic properties has a one-line reason. Boundedness: the integrand has modulus 1, so the average does too at worst. Uniform continuity: the increment φ(t+h)−φ(t) is E(eitX(eihX−1)), whose modulus is bounded by E∣eihX−1∣, a quantity depending on h alone and tending to 0 because the integrand is bounded by 2 and tends to 0 pointwise. Non-negative definiteness: the double sum is literally the expected squared modulus of ∑jzjeitjX, and a squared modulus cannot be negative.
The rules for sums and affine maps are the exponential law ea+b=eaeb plus the factorisation of expectations over independent variables - the same two facts that made generating functions multiply in Generating functions. The deep results, inversion and continuity, are not proved here: inversion is the Fourier inversion theorem, and continuity rests on it. What this unit supplies is the dictionary that lets those analytic theorems be read as statements about distributions.
Proof(The basic properties, and the rules for sums and affine maps)AdaptedAdapted from Grimmett & Stirzaker (2020), §5.7, proofs of Theorem (5.7.3)(a)-(c), Theorem (5.7.5) and Theorem (5.7.6), pp. 204-205; Reed & Simon (1980), §IX.2, the positivity computation and its three consequences, p. 329
Property 1.φX(0)=E(e0)=E(1)=1. For the bound, ∣E(Z)∣≤E∣Z∣ for a complex random variable Z, and ∣eitx∣=1 for real t,x, so ∣φX(t)∣≤E∣eitX∣=1.
using ∣eitX∣=1 again. The right-hand side does not involve t, so uniform continuity will follow once ϵ(h)→0. Set Yh=∣eihX−1∣. Then 0≤Yh≤2 for every h, and Yh→0 pointwise as h→0 because u↦eiu is continuous at 0. The bounded convergence theorem (Dominated convergence) therefore gives E(Yh)→0. This is where the boundedness of the integrand does the work that no moment hypothesis could: nothing was assumed about X.
Property 3. Let t1,…,tn be real and z1,…,zn complex, and put W=∑jzjeitjX. Then
∣W∣2=WW=j,k∑zjzkei(tj−tk)X,
since eitkX=e−itkX. Taking expectations term by term - a finite sum, so no interchange issue arises - gives
0≤E(∣W∣2)=j,k∑zjzkφX(tj−tk),
which is the claim. Reed and Simon run the same computation for a finite measure and then read off the three consequences of positivity by specialising: n=1 gives φ(0)≥0, and n=2 with t1=x, t2=0 forces the 2×2 matrix (φ(0)φ(−x)φ(x)φ(0)) to be self-adjoint with non-negative determinant, giving φ(−x)=φ(x) and ∣φ(x)∣≤φ(0). With φ(0)=1 this reproves property 1 and the first half of property 8.
Property 5.φaX+b(t)=E(eit(aX+b))=eitbE(ei(at)X)=eitbφX(at), the constant eitb coming out of the expectation because it is not random.
Property 6. If X and Y are independent then so are the complex random variables eitX and eitY, being functions of X and of Y separately. Expanding both into real and imaginary parts, multiplying out, applying E(UV)=E(U)E(V) to each of the four real products, and reassembling gives
φX+Y(t)=E(eitXeitY)=E(eitX)E(eitY)=φX(t)φY(t).
This is the route Grimmett and Stirzaker indicate; the detour through real and imaginary parts is needed because the factorisation theorem for expectations is stated for real variables.
Property 7, forward direction. The argument of property 6 with eisX and eitY gives the factorisation of φX,Y. For the converse one inverts the joint characteristic function and finds that the joint distribution function factorises; this uses the two-dimensional inversion formula and we do not reproduce it.
Property 8.φX(t)=E(eitX)=E(e−itX)=φX(−t). By property 5 with a=−1, b=0, the characteristic function of −X is φX(−t)=φX(t). So X and −X have the same characteristic function exactly when φX=φX, that is when φX is real; by uniqueness, that happens exactly when X and −X have the same distribution.
Property 9.φ=φ−X by property 8. If X′ is an independent copy of X then φ2=φX+X′ and ∣φ∣2=φφ=φX−X′, both by property 6. Finally Reφ=21(φ+φ) is the characteristic function of the variable equal to X or to −X with probability 21 each, independently of X, since a mixture of distributions has the corresponding mixture of characteristic functions. That ∣φ∣ need not be a characteristic function is Grimmett and Stirzaker's exercise, and is why the square in ∣φ∣2 matters. □
Remark(Where the hypotheses are used, and four traps)Original
Property 6 uses independence and property 7 shows that the converse needs more. The single-variable identity φX+Y=φXφY is the factorisation of the joint characteristic function evaluated only on the diagonal s=t, and a function of two variables is not determined by its values on a line. Grimmett and Stirzaker make this point at their eq. (5.7.9) and ask the reader to construct dependent X, Y satisfying it. So "the transform of the sum is the product of the transforms" is evidence of independence, never a proof of it.
Continuity of the limit at t=0 is indispensable in the continuity theorem, and failure is not a technicality. Let Z∼N(0,1) and Xn=nZ. By property 5, φn(t)=e−n2t2/2, which converges pointwise to the function equal to 1 at t=0 and 0 elsewhere. That limit is not continuous at 0, and correspondingly there is no limiting distribution: Fn(x)→21 for every fixed x, and 21 is not a distribution function. The probability mass has escaped to ±∞, and the discontinuity of the limit at the origin is exactly the signal that it did.
An odd derivative buys one moment less than you expect. Statement (a) says that the existence of φ′(0) gives only E∣X0∣<∞, which is no information at all; it does not give E∣X∣<∞. Even derivatives are different: φ′′(0) existing does force E(X2)<∞. The asymmetry is real and is why the hypothesis of the central limit theorem is a finite variance and not a finite mean plus a differentiable transform.
Characteristic functions that agree on an interval need not agree. Uniqueness requires agreement on all of R. Grimmett and Stirzaker's footnote to Corollary (5.9.3) records that agreement on a bounded interval is not enough and points to their Problem (5.12.27)(b) (Grimmett & Stirzaker (2020), §5.12, Problem 27(b), p. 234), which supplies the pair: one continuous distribution and one discrete distribution whose characteristic functions coincide on [−π,π] and then part company, the discrete one being periodic. So matching transforms over a finite band of frequencies identifies nothing.
The modulus is not a characteristic function. Property 9 gives ∣φ∣2 but stops there; Grimmett and Stirzaker's exercise asks for an example where ∣φ∣ fails. The class of characteristic functions is a convex set, closed under products and under conjugation, but not under square roots. It is not a cone: every characteristic function satisfies φ(0)=1, so λφ is again one only for λ=1. The cone is Reed and Simon's larger class of functions of positive type, of which the characteristic functions are exactly the members normalised at the origin.
φ(t)=M(it) is a theorem, not a definition. Riley, Hobson and Bence state the relation flatly, and it is correct wherever M is finite. But M may be finite only at t=0 - the Cauchy case - and then the formula says nothing, while φ(t)=e−∣t∣ exists and is perfectly well behaved. Grimmett and Stirzaker's Theorem (5.7.12) states the precise conditions under which the analytic continuation is legitimate: ∣M(t)∣<∞ for ∣t∣<a is equivalent to φ being analytic on the strip ∣Imz∣<a. Treating i as a real number inside an integral is not a proof even when the answer is right, a point they make twice.
Two notational hazards. The analysts' "characteristic function" is the indicator function, as Folland's dictionary records; and the probabilists' sign convention ∫e+itxdF is opposite to this handbook's forward Fourier transform, while Reed and Simon add a factor (2π)−n/2. Any formula imported from a Fourier text must be converted before use, by the dictionary given in the definition block.
Method(Using characteristic functions)Original
Write down φ for each ingredient, from the table or by direct integration. For a density supported on a half-line, the integral is usually a contour integral or a Laplace transform in disguise; for a bounded discrete variable it is a finite sum.
Combine. Independent sums multiply (property 6); affine maps contribute a phase and rescale the argument (property 5); mixtures mix. A standardised sum (Sn−nμ)/(σn) therefore has characteristic function e−itnμ/σ[φX(t/(σn))]n, which is the opening line of the proof of the central limit theorem.
Identify or invert. If the result is in the table, uniqueness finishes the job. Otherwise, if ∫∣φ∣<∞, recover the density from f(x)=2π1∫e−itxφ(t)dt; for a lattice variable, use the general inversion formula or extract the coefficient directly.
For a limit, take the pointwise limit of φn and check continuity at 0. If the limit is continuous at the origin, the continuity theorem converts it into convergence in distribution. If it is not, there is no limit law and mass is escaping.
For moments, expand near 0. Match coefficients against ∑jE(Xj)(it)j/j!; the coefficient of tj is ijE(Xj)/j!.
Check.φ(0)=1 and ∣φ(t)∣≤1 must hold; if X is real and symmetric, φ must be real and even; the small-t expansion must reproduce the mean and variance you already know; and a corner or other non-smoothness at t=0 signals an infinite moment of the corresponding order.
Example(The sum of two uniforms, by inversion)Original
Let X and Y be independent, each uniform on [−1,1], and let S=X+Y.
Step 1: the characteristic function of one summand. The density is 21 on [−1,1], so for t=0
φX(t)=∫−1121eitxdx=21⋅iteit−e−it=tsint,
and φX(0)=1, which is the limit of sint/t, so φX(t)=sinct is continuous everywhere. Two checks: φX is real and even, as property 8 requires for a symmetric variable, and ∣φX(t)∣≤1.
Step 2: the characteristic function of the sum. By property 6,
φS(t)=φX(t)φY(t)=t2sin2t.
Step 3: moments from the expansion. Using sint/t=1−t2/6+t4/120−⋯ and squaring,
φS(t)=1−3t2+452t4+O(t6).
Matching against ∑jE(Sj)(it)j/j!: the coefficient of t2 is −E(S2)/2, so E(S2)=2/3; the coefficient of t4 is E(S4)/24, so
E(S4)=24⋅452=1516=1.06667.
Both agree with a direct computation: var(X)=31 so E(S2)=31+31=32, and E(S4)=E(X4)+6E(X2)E(Y2)+E(Y4)=51+6⋅91+51=1516, the odd cross terms vanishing by symmetry.
Step 4: recover the density by inversion. Here ∣φS(t)∣≤min(1,t−2) is integrable, so the density form of the inversion theorem applies:
the sine part dropping out because the rest of the integrand is even. Evaluating,
fS(z)=42−∣z∣for ∣z∣≤2,fS(z)=0 otherwise.
This is the triangular density on [−2,2], which is also what the convolution of two uniform densities gives directly - the independent confirmation that the inversion was done correctly.
Step 5: numerical checks. At z=0 the formula predicts fS(0)=21, and indeed 2π1∫sin2t/t2dt=2π1⋅π=21. At z=21 it predicts (2−0.5)/4=0.375, and numerical quadrature of the inversion integral gives 0.37493, agreeing to the accuracy of the quadrature. Finally
(Every value in this example was verified symbolically and, where an integral was involved, numerically as well.)
Step 6: a contrast. Repeat the calculation with X and Y independent Cauchy variables. Then φX(t)=e−∣t∣, so φS(t)=e−2∣t∣ and, by property 5, the averageZ=21S has φZ(t)=φS(t/2)=e−∣t∣: the average of two independent Cauchy variables is again Cauchy, with no narrowing whatever. Grimmett and Stirzaker set this as a problem and ask why it does not violate the law of large numbers; the answer is that the Cauchy law has no mean, so Laws of large numbers does not apply. The corner of e−∣t∣ at t=0 is where that failure is visible.
Application(Limit theorems, lineshapes and scattering)Original
The central limit theorem and its relatives. The standard proof - expand φ to second order, raise to the nth power, take the limit, apply the continuity theorem - is carried out in The central limit theorem. The same three steps give the Poisson limit of the binomial, the convergence of a standardised Poisson or gamma variable to the normal, and the stable laws that arise when the variance is infinite; the Cauchy calculation in step 6 above is the simplest stable law other than the normal.
Infinitely divisible laws and Lévy processes. A law is infinitely divisible when φ=(φn)n for a characteristic function φn, at every n; the Lévy-Khinchin formula classifies all such φ. These are exactly the laws that can be the increment distributions of a process with stationary independent increments, which is how Brownian motion and the jump processes of mathematical finance and of anomalous diffusion are built.
Spectral lineshapes. The intensity profile of a spectral line is the Fourier transform of the correlation function of the emitting dipole, which is a characteristic function in disguise. Doppler broadening, which convolves the emission with a Maxwellian velocity distribution, multiplies transforms and gives a Gaussian profile; collisional broadening gives a Lorentzian, which is the Cauchy density, with the non-differentiability of e−∣t∣ at the origin corresponding to the infinite second moment of a Lorentzian line. A Voigt profile is the convolution of the two, and is handled by multiplying e−σ2t2/2 by e−γ∣t∣.
Scattering and structure factors. In diffraction from a disordered medium, the scattered amplitude at wavevector q is E(eiq⋅R) for the distribution of scatterer positions R - a characteristic function evaluated at the momentum transfer. The Debye-Waller factor e−q2⟨u2⟩/2 is exactly the characteristic function of a Gaussian displacement, and the intermediate scattering function measured in neutron spin-echo experiments is φΔR(τ)(q) for the displacement over time τ.
Random walks and Fourier methods. The position after n steps of a random walk has characteristic function φstep(t)n, so lattice Green functions, recurrence criteria and first-passage results in Random walks are all obtained by contour integration of a power of one transform. The same device, in the form of the structure function, underlies the numerical solution of master equations by Fourier transform.