Probability distributions, moments and expectation values

The cleanest case first: an experiment with countably many outcomes x₁, x₂, … and a number Pₙ attached to each, meaning the fraction of an unbounded run of identical trials that returns xₙ. Two axioms do all the work: every Pₙ ≥ 0, and Σ Pₙ = 1. The second is a completeness statement — the outcome list misses nothing — and it is the line to check first whenever a distribution is handed to you, because a list summing to 0.98 is silently losing one trial in fifty. The mean weights each outcome by its share, ⟨x⟩ = Σ xₙPₙ, and already carries this topic's central warning: a fair die has ⟨n⟩ = (1+2+3+4+5+6)/6 = 3.5, and no face shows 3.5. The mean is a property of the ensemble, never a prediction for one roll. Quantum mechanics will reuse this case verbatim: measuring an observable with a discrete, non-degenerate spectrum returns eigenvalue aₙ with probability Pₙ = |cₙ|², and the requirement Σ|cₙ|² = 1 is exactly this normalisation, resurfacing in Topic 1.4 as Parseval's identity.

When the outcome is continuous — a position, an arrival time — the probability of any exact value is zero: uncountably many candidates cannot each carry a finite share. The workable object is a density: ρ(x) dx is the probability of landing between x and x + dx, and only integrals of ρ are probabilities, P(a < x < b) = ∫ₐᵇ ρ dx. Two consequences trip students every year. First, ρ carries units — m⁻¹ for a position density, s⁻¹ for a lifetime density — because it must integrate against dx to a pure number. Second, the height of ρ is bounded by nothing: ρ(x) = 3x² on 0 ≤ x ≤ 1 m is perfectly normalised, ∫₀¹ 3x² dx = 1, yet reaches ρ = 3 m⁻¹ at the right edge. Only the area obeys the ≤ 1 rule; the height does whatever it must to hold the area. The Dirac delta bridges the two languages: a discrete outcome at x₀ with weight P is the density term P δ(x − x₀) — the form in which mixed spectra, bound levels plus a continuum, will be written later in the course.

Every number a distribution can export is one integral: ⟨f(x)⟩ = ∫ f(x) ρ(x) dx, the value of f at each outcome weighted by that outcome's share. The powers f = xᵏ are the moments. For ρ = 3x² on [0, 1 m]: ⟨x⟩ = ∫₀¹ 3x³ dx = 3/4 m. Averaging does not commute with nonlinear functions: ⟨x²⟩ = ∫₀¹ 3x⁴ dx = 3/5 = 0.600 m², while ⟨x⟩² = 0.5625 m² — and the gap is no accident, since ⟨x²⟩ − ⟨x⟩² = ⟨(x − ⟨x⟩)²⟩ ≥ 0, with equality only for a δ-sharp weight. Keep mean, median and mode distinct: here the mode is 1 m where ρ peaks, the median m solves m³ = 1/2 giving 0.794 m, and the mean is 0.750 m. Skew pulls the three apart — here, piled against the right wall, in the order mean < median < mode — and a symmetric single-peaked density brings them back together. Symmetry is also the working shortcut: when ρ is even about a point, the integrand (x − x₀)ρ is odd and the mean sits at x₀ by cancellation — the one-line argument that will give ⟨x⟩ = L/2 for every infinite-well eigenstate before any integral is attempted.

The definitional form σ² = ⟨(x − ⟨x⟩)²⟩ says what variance is: the mean squared deviation from the mean. Expanding the square gives ⟨x²⟩ − 2⟨x⟩⟨x⟩ + ⟨x⟩² = ⟨x²⟩ − ⟨x⟩², the computational form, quickest by hand: for ρ = 3x² on [0, 1 m], σ² = 3/5 − 9/16 = 3/80 = 0.0375 m², so σ = 0.194 m. Numerically the shortcut can be a trap: a narrow density parked at a large mean — a 1 nm spread on a 10 μm mean — makes ⟨x²⟩ and ⟨x⟩² agree to eight significant figures, and a float64 subtraction hands back only the eight that remain; numpy.var centres the data before squaring for exactly this reason. Quote σ, not σ², as the physical width: it shares x's unit. And keep it distinct from an instrument's error bar — this σ is intrinsic to the distribution and survives perfect apparatus. The Δx in Δx Δp ≥ ħ/2 is precisely this standard deviation, computed later from the weight |ψ|².

The Lorentzian ρ(x) = (Γ/2π)/[(x − x₀)² + (Γ/2)²] integrates to exactly 1, yet its tails fall only as 1/x²: the variance integrand x²ρ tends to a constant, so σ² diverges, and even ⟨x⟩ exists only as a symmetric principal value. This is no invented pathology — it is the natural line shape of an exponentially decaying state, and it returns as the Breit–Wigner resonance in the particle-physics units, which is why resonances are quoted by their full width at half maximum Γ and never by a standard deviation. The practical symptom is that averaging fails: the mean of N independent Cauchy draws is distributed exactly like a single draw, so a simulation that tries to estimate ⟨x⟩ wanders forever instead of settling. Before quoting any moment, check the tail: ⟨xᵏ⟩ exists only if ρ falls faster than |x|⁻ᵏ⁻¹ as |x| → ∞.

Everything above is classical probability, and that is the point: quantum mechanics changes where the weight comes from, not what you do with it. Born's rule — postulated in Unit 2 and applied to position in Unit 3 — will declare ρ(x) = |ψ(x)|² for position and Pₙ = |cₙ|² for a discrete spectrum; every formula on this page then applies verbatim, merely rewritten there in operator dress, ⟨x⟩ = ⟨ψ|x̂|ψ⟩. What the theory predicts is therefore ensemble statistics: prepare N identical systems, measure each once, and the histogram converges on ρ while the sample mean approaches ⟨x⟩ with statistical scatter dying as σ/√N. No formula here addresses a single run — not a week-one gap to be repaired later, but the permanent shape of the theory's output. One NumPy line makes it concrete: rng.choice([1, 2, 3], p=[0.2, 0.5, 0.3], size=10000) has sample means scattering about the exact ⟨a⟩ = 2.1 with spread σ/√N = 0.7/100 = 0.007 — the 1/√N crawl every Monte Carlo lab in this course will exhibit.