03 - Continuous Random Variables: Uniform, Exponential, Normal

You are waiting for a tram that comes every ten minutes, but you do not know when the last one went by. What is the probability of waiting exactly three minutes? Zero. What is the probability of waiting less than three minutes? 30%. The two answers seem to contradict each other until you understand how a variable that can take any value in an interval works: probability is not assigned to points, it is assigned to intervals, computed as an area.

This article builds the tool (density and distribution function), applies it to the three most used distributions, uniform, exponential and normal, and shows how the density changes when a variable is transformed. As in the article on discrete variables, every result is computed twice, with the formula and with a fixed-seed simulation, and the simulation's tolerance is declared before you read its outcome. You need to be able to compute an elementary integral.

Density and distribution function

A random variable is continuous if there is a density f(x) ≥ 0, with total area 1, such that the probability of falling in an interval is the area under the density: P(a < X ≤ b) = ∫_a^b f(x) dx.

The distribution function (CDF) F(x) = P(X ≤ x) accumulates that area from the left, so P(a < X ≤ b) = F(b) − F(a), and where F is differentiable the density is its derivative.

Two consequences need to be understood well. First: the density is not a probability, it can exceed 1 (a uniform density on an interval of length 0.1 is worth 10) and is only the probability per unit of length. Second: P(X = x) = 0 for every single value x, because the area of an interval shrinking to a point tends to zero; this does not make the value impossible, because one of the values has to come up.

Expected value and variance are defined as for discrete variables, with the integral in place of the sum: E[X] = ∫ x f(x) dx and Var(X) = ∫ (x − E[X])² f(x) dx. The same rules hold: the mean of a sum is the sum of the means, the variance adds for independent variables.

The tram is the simplest case, the uniform distribution on [a, b]: constant density 1/(b − a) in that interval and zero outside, mean (a + b)/2 and variance (b − a)²/12. The program numerically integrates the tram's density with scipy.integrate.quad to check area, mean and variance, computes the probabilities with the CDF, shows how the probability of an interval of width h around 3 minutes shrinks with h, and simulates a million arrivals.

from fractions import Fraction as F
import numpy as np
from scipy.integrate import quad

# You wait for the tram: the last one passed at a time you do not know, the next arrives within 0 to 10 minutes.
a, b = 0.0, 10.0
f = lambda x: 1 / (b - a) if a <= x <= b else 0.0      # density
Fc = lambda x: min(max((x - a) / (b - a), 0.0), 1.0)    # distribution function (CDF)

area, _ = quad(f, a, b)
media, _ = quad(lambda x: x * f(x), a, b)
var, _ = quad(lambda x: (x - media) ** 2 * f(x), a, b)
print(f"area of the density = {area:.6f}  mean = {media:.6f}  variance = {var:.6f}  (formulas: {(a+b)/2}, {(b-a)**2/12:.6f})")
print(f"P(wait < 3 min) = {Fc(3):.4f}   P(2 < wait <= 7) = {Fc(7) - Fc(2):.4f}")
print("P(wait = exactly 3 min) = 0, but P(3 <= wait <= 3 + h) is h/10 =", [f"{h}: {Fc(3 + h) - Fc(3):.4f}" for h in (1, 0.1, 0.001)])

N = 10**6
rng = np.random.default_rng(21)
x = rng.uniform(a, b, size=N)
for nome, vera, oss in (("P(<3)", 0.3, np.mean(x < 3)), ("P(2<..<=7)", 0.5, np.mean((x > 2) & (x <= 7)))):
    tol = 3 * (vera * (1 - vera) / N) ** 0.5
    print(f"{nome}: Monte Carlo {oss:.5f} expected {vera} +/- {tol:.5f}")
    assert abs(oss - vera) < tol
print(f"sample mean {x.mean():.4f}, sample variance {x.var():.4f}")

# Output:
# area of the density = 1.000000  mean = 5.000000  variance = 8.333333  (formulas: 5.0, 8.333333)
# P(wait < 3 min) = 0.3000   P(2 < wait <= 7) = 0.5000
# P(wait = exactly 3 min) = 0, but P(3 <= wait <= 3 + h) is h/10 = ['1: 0.1000', '0.1: 0.0100', '0.001: 0.0001']
# P(<3): Monte Carlo 0.29924 expected 0.3 +/- 0.00137
# P(2<..<=7): Monte Carlo 0.49996 expected 0.5 +/- 0.00150
# sample mean 5.0026, sample variance 8.3285

The area is 1, the mean 5 minutes and the variance 8.33 minutes squared, equal to the formulas. The last series of numbers in the third line is the answer to the paradox at the start: the interval of width h around 3 minutes has probability h/10, which tends to zero with h, so "exactly 3 minutes" weighs zero, while "less than 3" weighs 0.3.

The simulations confirm within tolerance (deviations 0.00076 and 0.00004 against tolerances of 0.0014 and 0.0015). A caveat about the model: the uniform describes the tram only if it arrives on time and you arrive at a random instant independent of the timetable; if you arrive in order to catch that very tram, the distribution is no longer this one.

Exponential: the waiting time without memory

The exponential distribution with parameter λ > 0 has density f(x) = λe^−λx for x ≥ 0, distribution function F(x) = 1 − e^−λx, mean 1/λ and variance 1/λ². It models a waiting time.

The function S(t) = P(X > t) = e^−λt, the probability of surviving beyond t, has a remarkable property: P(X > s + t | X > s) = e^−λ(s+t)/e^−λs = e^−λt = P(X > t).

Whoever has survived up to s faces the same future as whoever starts from zero: this is memorylessness, the continuous analogue of what you saw for the geometric. It is a strong property and rarely true in real systems: a component that wears out is not memoryless. The exponential is the right model only for failures that arrive at random, without ageing.

The example is a component with a mean life of 1000 hours (λ = 0.001 per hour; hypothetical, chosen to make the numbers readable). The probability that it lasts beyond 1500 hours is e^−1.5, and that of lasting beyond 1500 hours knowing that it has already passed 1000 is e^−0.5, identical to the probability of lasting 500 hours when new.

The program samples by inversion: if U is uniform on (0, 1], then X = −ln(U)/λ is exponential, because it inverts the CDF. The instruction 1 - rng.random(N) serves to avoid the logarithm of zero, since random() returns values in [0, 1). The Kolmogorov-Smirnov test measures the maximum distance between the empirical CDF of the sample and the theoretical one.

import numpy as np
from math import exp, log
from scipy.stats import kstest, expon, uniform

# Life of a component: exponential with mean 1000 hours, so lambda = 1/1000
lam = 1 / 1000
S = lambda t: exp(-lam * t)                      # P(X > t), survival function
print(f"P(X > 1500) = e^-1.5 = {S(1500):.5f}")
print(f"P(X > 1500 | X > 1000) = {S(1500) / S(1000):.5f}  and  P(X > 500) = {S(500):.5f}")

# Sampling by inversion: if U is uniform on (0,1], X = -ln(U)/lambda is exponential
N = 10**6
rng = np.random.default_rng(9)
u = 1 - rng.random(N)                            # rng.random() lies in [0,1): 1-u avoids log(0)
x = -np.log(u) / lam
print(f"mean {x.mean():.2f} (expected 1000, standard error {1000 / N**0.5:.2f}); variance {x.var():.0f} (expected 1000000)")

oss = np.mean(x > 1500)
tol = 3 * (S(1500) * (1 - S(1500)) / N) ** 0.5
print(f"frequency X>1500: {oss:.5f} (expected {S(1500):.5f} +/- {tol:.5f})")
assert abs(oss - S(1500)) < tol

dopo = x[x > 1000]                               # the components still alive at 1000 hours
oss_c = np.mean(dopo > 1500)
tol_c = 3 * (S(500) * (1 - S(500)) / len(dopo)) ** 0.5
print(f"among the {len(dopo)} components alive at 1000 hours, share alive at 1500: {oss_c:.5f} (expected {S(500):.5f} +/- {tol_c:.5f})")
assert abs(oss_c - S(500)) < tol_c

# The Kolmogorov-Smirnov test compares the empirical distribution function with the theoretical one
print(f"KS against Exp(1/1000): D = {kstest(x, expon(scale=1000).cdf).statistic:.5f}")
# Probability integral transform: F(X) is uniform on (0,1) when F is the distribution function of X
print(f"KS of F(X) against the uniform: D = {kstest(1 - np.exp(-lam * x), uniform.cdf).statistic:.5f}")

# Output:
# P(X > 1500) = e^-1.5 = 0.22313
# P(X > 1500 | X > 1000) = 0.60653  and  P(X > 500) = 0.60653
# mean 999.88 (expected 1000, standard error 1.00); variance 1003337 (expected 1000000)
# frequency X>1500: 0.22330 (expected 0.22313 +/- 0.00125)
# among the 366897 components alive at 1000 hours, share alive at 1500: 0.60861 (expected 0.60653 +/- 0.00242)
# KS against Exp(1/1000): D = 0.00108
# KS of F(X) against the uniform: D = 0.00108

The two theoretical values, 0.60653, coincide as expected, and the simulated frequencies are 0.22330 against 0.22313 and 0.60861 against 0.60653. This last deviation is 2.6 standard errors, the highest among the forty seeds I tried to check that there was no distortion (mean of the deviations −0.07 in units of standard error, deviation 0.96): it is therefore a normal fluctuation, and an example of why the tolerance is set at three errors and not at one.

The sample variance, 1,003,337, differs by 0.33% from 10⁶, within its error (the estimated variance of such an asymmetric variable converges slowly, with a standard error of about 2,800). The two Kolmogorov-Smirnov distances, 0.00108, coincide because transforming the sample with its CDF is an increasing transformation that does not change the order of the values: if F is the CDF of X, F(X) is uniform on (0, 1), and this is the foundation of sampling by inversion.

Normal: the distribution of measurement error

The normal distribution (or Gaussian) with mean μ and variance σ² has density f(x) = (1/(σ√(2π))) exp(−(x − μ)²/(2σ²)): a symmetric bell centred at μ, with width set by σ. The standard normal is the one with μ = 0 and σ = 1; its CDF is denoted Φ and has no elementary form, so it is computed numerically (with the error function, erf, or with a library).

Two properties hold that you use constantly: if X is normal, so is aX + b, with mean aμ + b and variance a²σ²; and the sum of independent normals is normal, with means and variances that add. This is why it is enough to know the standard one: every normal reduces to it by subtracting the mean and dividing by the standard deviation.

The probability of lying within k standard deviations of the mean depends neither on μ nor on σ: it is erf(k/√2). The program computes it this way, compares it with SciPy and with a million draws, and checks that the area of the density is 1.

Then it applies everything to a hypothetical quality control: parts with a normal diameter of mean 10.00 mm and standard deviation 0.02 mm, tolerance ±0.03 mm. By symmetry the share of parts out of tolerance is 2(1 − Φ(1.5)). Finally it computes the central 99% interval and checks the rule on the sum of two independent normals.

import numpy as np
from math import erf, sqrt
from scipy.integrate import quad
from scipy.stats import norm

phi = lambda t: np.exp(-t * t / 2) / sqrt(2 * np.pi)       # density of the standard normal
area, _ = quad(phi, -np.inf, np.inf)
print(f"integral of the standard density = {area:.10f}")

# Probability within k standard deviations: P(|Z|<k) = erf(k/sqrt(2)), no tables
N = 10**6
rng = np.random.default_rng(13)
z = rng.standard_normal(N)
for k in (1, 2, 3):
    esatta = erf(k / sqrt(2))
    oss = np.mean(np.abs(z) < k)
    tol = 3 * (esatta * (1 - esatta) / N) ** 0.5
    print(f"|Z|<{k}: erf {esatta:.6f}  scipy {norm.cdf(k) - norm.cdf(-k):.6f}  Monte Carlo {oss:.6f} (+/- {tol:.6f})")
    assert abs(oss - esatta) < tol

# Quality control (hypothetical case): diameter N(10.00, 0.02^2) mm, tolerance 10.00 +/- 0.03 mm
mu, sigma, tol_mm = 10.00, 0.02, 0.03
fuori = 2 * (1 - norm.cdf(tol_mm / sigma))                  # symmetry: two equal tails
x = rng.normal(mu, sigma, size=N)
oss = np.mean(np.abs(x - mu) > tol_mm)
tol = 3 * (fuori * (1 - fuori) / N) ** 0.5
print(f"parts out of tolerance: exact {fuori:.5f}, Monte Carlo {oss:.5f} (+/- {tol:.5f})")
assert abs(oss - fuori) < tol
print(f"central 99% interval: mean +/- {norm.ppf(0.995):.4f} sigma = [{mu - norm.ppf(0.995)*sigma:.4f}, {mu + norm.ppf(0.995)*sigma:.4f}] mm")

# Sum of independent normals: N(mu1+mu2, s1^2+s2^2)
y = rng.normal(5, 0.015, size=N)
s = x + y
print(f"sum: mean {s.mean():.4f} (15), deviation {s.std():.5f} (expected {sqrt(0.02**2 + 0.015**2):.5f})")

# Output:
# integral of the standard density = 1.0000000000
# |Z|<1: erf 0.682689  scipy 0.682689  Monte Carlo 0.682402 (+/- 0.001396)
# |Z|<2: erf 0.954500  scipy 0.954500  Monte Carlo 0.954419 (+/- 0.000625)
# |Z|<3: erf 0.997300  scipy 0.997300  Monte Carlo 0.997275 (+/- 0.000156)
# parts out of tolerance: exact 0.13361, Monte Carlo 0.13355 (+/- 0.00102)
# central 99% interval: mean +/- 2.5758 sigma = [9.9485, 10.0515] mm
# sum: mean 15.0000 (15), deviation 0.02500 (expected 0.02500)

The probabilities within one, two and three sigma are 0.682689, 0.954500 and 0.997300, and the simulation confirms them within tolerance. In the case of the parts, 13.36% fall outside tolerance: the designer reads in this number that a tolerance of 1.5 sigma is too tight, and that to scrap less than 1% a tolerance of at least 2.58 sigma is needed (the 99% interval goes from 9.9485 to 10.0515 mm).

Adding another independent normal variable of mean 5 and deviation 0.015, the standard deviation of the total is √(0.02² + 0.015²) = 0.025, not 0.035: standard deviations do not add, variances do. The normal is a good model when the error is the sum of many small independent contributions, as the article on the central limit theorem explains; where extreme events are much more frequent than this (heavy tails), its estimate of rare probabilities is too optimistic.

Transformations: how the density changes

If Y = g(X) and g is monotone and differentiable, the density of Y is f_Y(y) = f_X(g⁻¹(y)) · |dg⁻¹/dy|: the factor with the derivative is the "Jacobian" in one variable, and compensates for the fact that g stretches or shrinks the axis. If g is not monotone, the contributions of the branches are added.

The classic example is Y = X² with X standard normal: every positive value y comes from two values of X, +√y and −√y, each with weight φ(√y)/(2√y), whence f_Y(y) = e^−y/2/√(2πy).

This is the chi-square distribution with one degree of freedom. The program checks the formula in three independent ways, which must give the same number: the numerical integral of the density, the CDF obtained directly as 2Φ(√y) − 1 (because Y ≤ y is equivalent to |X| ≤ √y), and SciPy's chi-square CDF; and it compares them with the simulation.

import numpy as np
from math import sqrt, pi, exp
from scipy.integrate import quad
from scipy.stats import norm, chi2

# Y = X^2 with X standard normal. X^2 is not monotone: sum over the two branches X = +sqrt(y), X = -sqrt(y).
phi = lambda t: exp(-t * t / 2) / sqrt(2 * pi)
f_Y = lambda y: 2 * phi(sqrt(y)) / (2 * sqrt(y))        # f_X(sqrt y) |d sqrt(y)/dy| for each branch
print("f_Y(y) = exp(-y/2) / sqrt(2 pi y)?", all(abs(f_Y(y) - exp(-y / 2) / sqrt(2 * pi * y)) < 1e-15 for y in (0.1, 1, 4)))

N = 10**6
rng = np.random.default_rng(17)
y = rng.standard_normal(N) ** 2
print(f"{'interval':>12} {'integral of f_Y':>17} {'2Phi(sqrt)-1':>14} {'chi2(1)':>9} {'Monte Carlo':>12}")
for a, b in ((0.0, 1.0), (1.0, 4.0), (4.0, 9.0)):
    integr, _ = quad(f_Y, a, b)
    dalla_cdf = (2 * norm.cdf(sqrt(b)) - 1) - (2 * norm.cdf(sqrt(a)) - 1)
    chi = chi2.cdf(b, 1) - chi2.cdf(a, 1)
    oss = np.mean((y > a) & (y <= b))
    tol = 3 * (dalla_cdf * (1 - dalla_cdf) / N) ** 0.5
    assert abs(integr - dalla_cdf) < 1e-9 and abs(chi - dalla_cdf) < 1e-12 and abs(oss - dalla_cdf) < tol
    print(f"({a:.0f}, {b:.0f}] {integr:17.6f} {dalla_cdf:14.6f} {chi:9.6f} {oss:12.6f}")
print(f"E[Y] = {y.mean():.4f} (exact 1), Var(Y) = {y.var():.3f} (exact 2)")

# Output:
# f_Y(y) = exp(-y/2) / sqrt(2 pi y)? True
#     interval   integral of f_Y   2Phi(sqrt)-1   chi2(1)  Monte Carlo
# (0, 1]          0.682689       0.682689  0.682689     0.682203
# (1, 4]          0.271810       0.271810  0.271810     0.272241
# (4, 9]          0.042800       0.042800  0.042800     0.042816
# E[Y] = 1.0011 (exact 1), Var(Y) = 2.000 (exact 2)

The three calculations agree to the sixth decimal place in every interval (0.682689 in (0, 1], 0.271810 in (1, 4], 0.042800 in (4, 9]) and the assertions check that they agree to 10⁻⁹ as well; the simulation is within tolerance. The first value coincides with the probability of lying within one standard deviation of the normal, and not by chance: X² ≤ 1 is equivalent to |X| ≤ 1.

The mean of Y is 1, the variance 2, values that follow from the moments of the normal (E[X²] = 1, E[X⁴] = 3, hence Var(X²) = 3 − 1 = 2). The typical mistake is to forget the negative branch and obtain a density that integrates to 1/2: the area check, as in the first program, reveals it at once.

What is left out

Missing here are the moment generating function and the characteristic function, the Poisson process (whose times between successive arrivals the exponential describes), the families derived from the normal such as Student's t, the sampling of the normal with the Box-Muller method and transformations in several dimensions with the Jacobian determinant. For transformations of vectors go to article 04, where the multivariate normal appears. To return to the discrete case and its links with these models, article 01 contains the axioms from which these integrals also descend.

Frequently asked questions

Why is the probability of a single value zero?

Because probability is the area under the density, and the area of an interval of width zero is zero. This does not mean that the value cannot come up: it means that, among infinitely many possible values, none has a positive probability on its own.

Can a density be greater than 1?

Yes. The density is a probability per unit of length: on a short interval it can exceed 1, as long as the total area stays 1. A uniform on an interval of length 0.1 has density 10.

What does it mean that the exponential has no memory?

That the probability of waiting a further time t does not depend on how long you have already waited. It holds for purely random failures, not for components that wear out.

What you can do now

You can go from density to CDF and back, compute the probability of an interval as an area and understand why a single point weighs zero. You can use the uniform, exponential and normal, sample by inversion, reduce any normal to the standard one and find the density of a function of a variable, checking every result with a second method. In the next article you move from one variable to several together: the joint distribution, covariance and the multivariate normal.

The same series also includes Exploratory Data Analysis (EDA) with Python. The full list is on the series page, Probability and Statistics.