57. Lebesgue integration and expectation

🎬 Video · 15 min
💡 Every code box below is live — edit it and hit Run.

§1 gave two formulas for expectation — a sum for discrete variables, an integral for continuous ones — and left mixed cases unaddressed.

Lebesgue integration replaces both with one definition, and it's the reason the limit theorems of §3 can be proved at all.

Riemann vs. Lebesgue

The difference is which axis you slice.

Riemann partitions the domain: chop the x-axis into intervals, approximate f by its value on each, sum f(x_i)\Delta x.

Lebesgue partitions the range: chop the y-axis into levels, and for each level ask how much of the domain maps there. Sum $y_i \cdot \mu({x : f(x) \approx y_i})$.

The standard analogy: to count a pile of coins, Riemann picks them up in the order they lie; Lebesgue sorts them into denominations first and multiplies.

Why this matters: measuring \{x : f(x) \approx y\} requires only that the set be measurable, not that f be continuous or even close to it. That's what makes the Lebesgue integral vastly more robust.

The construction

Build it in three steps.

1. Simple functions. A simple function takes finitely many values: s = \sum_{i=1}^n a_i \mathbb 1_{A_i} with A_i measurable. Define

\int s\,d\mu = \sum_{i=1}^n a_i\,\mu(A_i)

2. Non-negative functions. Approximate from below by simple functions:

\int f\,d\mu = \sup\left\{\int s\,d\mu \;:\; 0 \le s \le f,\; s \text{ simple}\right\}

Every non-negative measurable function is the increasing limit of simple ones, so this is always defined (possibly +\infty).

3. General functions. Split into positive and negative parts, f = f^+ - f^-:

\int f\,d\mu = \int f^+\,d\mu - \int f^-\,d\mu

defined whenever at least one part is finite. f is integrable if \int|f|\,d\mu < \infty.

Expectation, unified

For a random variable X on (\Omega, \mathcal F, P):

\boxed{\;E[X] = \int_\Omega X\,dP\;}

One definition, no cases. The discrete and continuous formulas from §1 are now consequences rather than definitions:

  • Discrete X: P concentrates on atoms, and the integral collapses to \sum_x x\,P(X = x).
  • Continuous X: if P_X has a density f with respect to Lebesgue measure, the integral becomes \int x f(x)\,dx.
  • Mixed X: handled automatically, with no new machinery. An insurance payout that is 0 with probability 0.9 and continuous otherwise has an expectation that neither §1 formula could express.

More generally, E[g(X)] = \int g\,dP_X — the change-of-variables formula that justifies computing expectations on the distribution rather than the underlying space.

Why Lebesgue wins

1. It integrates more functions. The indicator of the rationals, \mathbb 1_{\mathbb Q} on [0,1], is nowhere Riemann integrable: every subinterval contains both rationals and irrationals, so upper sums are 1 and lower sums are 0.

The Lebesgue integral is trivial: \mathbb Q is a null set (last lesson), so

\int_0^1 \mathbb 1_{\mathbb Q}\,d\lambda = 1 \cdot \lambda(\mathbb Q) = 0

2. It handles limits. This is the decisive advantage. The Riemann integral requires uniform convergence to exchange limit and integral — a demanding condition that fails constantly. Lebesgue needs only monotonicity or domination, which is why §3's convergence theorems work.

3. It's general. The same definition works on any measure space — the real line, a discrete set, a function space, a manifold. Probability inherits all of it for free.

When they agree: if f is Riemann integrable on [a,b], it's Lebesgue integrable with the same value. Lebesgue extends rather than contradicts.

Key properties

Linearity: \int(af + bg)\,d\mu = a\int f + b\int g — this is §1's linearity of expectation, now proved rather than asserted.

Monotonicity: f \le g \implies \int f \le \int g.

Null sets don't matter: if f = g almost everywhere, their integrals are equal. Changing a function on a measure-zero set changes nothing — which is why a density is only ever defined up to null sets.

Jensen's inequality: for convex \varphi,

\varphi\big(E[X]\big) \le E\big[\varphi(X)\big]

We've been using this informally since §4 (E[\sqrt Y] \ne \sqrt{E[Y]}, s not unbiased for \sigma, the mean/median gap for lognormals in §8). It's a theorem about integrals against a probability measure.

Worked example

A mixed random variable. An insurance policy pays 0 with probability 0.7, and otherwise pays an amount that is Uniform(0, 1000). What is E[X]?

Neither §1 formula applies directly: X isn't discrete (it has a continuous part) and isn't continuous (it has an atom at 0). The Lebesgue integral splits over the two pieces automatically:

E[X] = \int_\Omega X\,dP = \underbrace{0 \times 0.7}_{\text{atom}} + \underbrace{0.3 \times \int_0^{1000}\frac{x}{1000}\,dx}_{\text{continuous part}}

= 0 + 0.3 \times 500 = 150

The decomposition works because the measure P_X splits into a discrete part (a Dirac mass at 0) plus a part with a density, and the integral is additive over that decomposition.

Note also that P(X = 0) = 0.7 while P(X = 500) = 0 — the same variable has both an atom and a continuum, which is exactly what the unified definition handles without comment.

Doing it in Python

Riemann and Lebesgue agreeing on a well-behaved function:

import numpy as np
from scipy.integrate import quad

f = lambda x: x**2

# Riemann: partition the DOMAIN
n = 100_000
xs = np.linspace(0, 1, n + 1)
riemann = np.sum(f(xs[:-1]) * np.diff(xs))

# "Lebesgue": partition the RANGE, measuring the preimage of each level
levels = np.linspace(0, 1, n + 1)          # f maps [0,1] onto [0,1]
# For f(x) = x^2, {x : f(x) > y} = (sqrt(y), 1], of measure 1 - sqrt(y)
lebesgue = np.sum((1 - np.sqrt(levels[:-1])) * np.diff(levels))

print(f"Riemann  (domain slices): {riemann:.8f}")
print(f"Lebesgue (range slices) : {lebesgue:.8f}")
print(f"exact                   : {1/3:.8f}")
print(f"scipy quad              : {quad(f, 0, 1)[0]:.8f}")

That "layer cake" identity — \int f = \int_0^\infty \mu(f > y)\,dy for non-negative f — is the Lebesgue construction in one line, and it's a useful computational tool in its own right.

A function Riemann cannot integrate:

import numpy as np
from fractions import Fraction

# 1_Q on [0,1]: 1 at rationals, 0 elsewhere.
def indicator_rational(x, max_denominator=1000):
    return 1.0 if Fraction(x).limit_denominator(max_denominator) == Fraction(x) else 0.0

# Riemann: every subinterval contains both rationals and irrationals,
# so the upper sum is always 1 and the lower sum always 0.
print("Riemann upper sum on any partition: 1.0")
print("Riemann lower sum on any partition: 0.0")
print("-> they never meet, so the Riemann integral does not exist.\n")

# Lebesgue: Q is a null set, so the integral is 0.
print("Lebesgue: integral = 1 * measure(Q) = 1 * 0 = 0")

# Sampling agrees with the Lebesgue answer
rng = np.random.default_rng(0)
draws = rng.uniform(0, 1, 100_000)
hits = sum(indicator_rational(float(x)) for x in draws[:2000])
print(f"\nsampled average over 2000 uniform draws: {hits/2000:.4f}")

Expectation as one integral, across all three kinds of variable:

import numpy as np
from scipy.integrate import quad

rng = np.random.default_rng(1)
n = 2_000_000

# 1. Discrete
die = rng.integers(1, 7, n)
print(f"discrete   : simulated {die.mean():.4f}   theory {3.5}")

# 2. Continuous
expo = rng.exponential(2.0, n)
print(f"continuous : simulated {expo.mean():.4f}   theory {2.0}")

# 3. Mixed -- the insurance policy
pays = rng.random(n) > 0.7
amount = np.where(pays, rng.uniform(0, 1000, n), 0.0)
print(f"mixed      : simulated {amount.mean():.4f}   theory {150.0}")

print(f"\nP(X = 0)   = {(amount == 0).mean():.4f}  (an atom)")
print(f"P(X = 500) = {(amount == 500).mean():.4f}  (no atom -- continuous there)")
print("\nOne definition, E[X] = integral of X dP, covers all three.")

Jensen's inequality, which shows up all over this course:

import numpy as np

rng = np.random.default_rng(2)
X = rng.uniform(1, 10, 500_000)

cases = [
    ("x^2 (convex)",   lambda x: x**2,      "<="),
    ("e^x (convex)",   lambda x: np.exp(x/5), "<="),
    ("sqrt (concave)", lambda x: np.sqrt(x), ">="),
    ("log (concave)",  lambda x: np.log(x),  ">="),
]
print(f"{'g':>16} {'g(E[X])':>12} {'E[g(X)]':>12} {'relation':>10}")
for name, g, rel in cases:
    lhs, rhs = g(X.mean()), g(X).mean()
    actual = "<=" if lhs <= rhs else ">="
    print(f"{name:>16} {lhs:>12.4f} {rhs:>12.4f} {actual:>10}")

print("\nConvex: g(E[X]) <= E[g(X)]. Concave: the reverse.")
print("This is why s is not unbiased for sigma (§4), and why the mean of a")
print("lognormal exceeds its median (§8).")

Your turn

1. What is \int_0^1 \mathbb 1_{\mathbb Q}\,d\lambda, and why can't Riemann compute it?

2. If f = g almost everywhere, what can you say about their integrals?

3. X is 0 with probability 0.5 and Uniform(0,10) otherwise. Find E[X].

Solutions

1. The integral is 0.

\mathbb 1_{\mathbb Q} is 1 on the rationals and 0 elsewhere. Since \lambda(\mathbb Q) = 0 (last lesson),

\int_0^1 \mathbb 1_{\mathbb Q}\,d\lambda = 1 \cdot \lambda(\mathbb Q) + 0 \cdot \lambda(\text{irrationals}) = 0

Riemann fails because both the rationals and the irrationals are dense. Every subinterval, however small, contains points where the function is 1 and points where it is 0. So every upper sum equals 1 and every lower sum equals 0, and they never converge to a common value.

Lebesgue sidesteps this by never subdividing the domain — it asks only "what is the measure of the set where f = 1?", and that set is measurable with measure 0.

2. Their integrals are equal.

f = g \text{ a.e.} \implies \int f\,d\mu = \int g\,d\mu

The set where they differ has measure zero, so it contributes nothing.

Two important consequences:

  • Densities are not unique. Change a PDF at finitely (or countably) many points and it's still a valid density for the same distribution. That's why the value of a density at a point has no meaning — only integrals over sets do.
  • A.e. equal random variables have the same expectation, variance, and distribution. They are the same object for all statistical purposes.

3. Split by the atom and the continuous part:

E[X] = 0.5 \times 0 + 0.5 \times E[\text{Uniform}(0,10)] = 0.5 \times 5 = 2.5

Neither §1 formula alone handles this. It has an atom at 0 (so no density exists there) and a continuum elsewhere (so it isn't discrete). The measure-theoretic definition needs no special case — the integral is simply additive over the decomposition of P_X into its discrete and continuous parts.

Mixed distributions are common in practice: insurance claims with a mass at zero, rainfall with a mass at "no rain", censored survival times. This is the framework that handles them without improvisation.

Check yourself in code

Verify the layer-cake identity — that integrating over range slices gives the same answer as over domain slices — for f(x) = x^2 on [0,1].

Print exactly this:

riemann 0.3333
lebesgue 0.3333
exact 0.3333
agree: True

Use 100000 slices for both. For the range-slicing version use \lambda(\{x : x^2 > y\}) = 1 - \sqrt y. Round every value to 4 decimal places.

import numpy as np

n = 100_000

xs = np.linspace(0, 1, n + 1)
riemann = np.sum(xs[:-1] ** 2 * np.diff(xs))
print("riemann", round(riemann, 4))

# Compute the same integral by slicing the RANGE: sum over levels y of
# lambda({x : f(x) > y}) = 1 - sqrt(y). Then compare both with 1/3.
import numpy as np

n = 100_000

xs = np.linspace(0, 1, n + 1)
riemann = np.sum(xs[:-1] ** 2 * np.diff(xs))
print("riemann", round(riemann, 4))

levels = np.linspace(0, 1, n + 1)
lebesgue = np.sum((1 - np.sqrt(levels[:-1])) * np.diff(levels))
print("lebesgue", round(lebesgue, 4))
print("exact", round(1 / 3, 4))
print("agree:", round(riemann, 4) == round(lebesgue, 4) == round(1 / 3, 4))

Lebesgue integration slices the range instead of the domain, so it only ever needs sets to be measurable rather than functions to be smooth. That gives one definition of expectation, E[X] = \int X\,dP, covering discrete, continuous and mixed variables alike — and it makes limits interchangeable with integrals under conditions weak enough to be useful.

Next: exactly which conditions those are.