Sum of Random Variables

Learning Objectives

By the end of this section, you will have a deep understanding of how random variables combine when added together. You will be able to:

1. Derive the distribution of $Z = X + Y$ for both discrete and continuous random variables using the convolution formula
2. Apply the convolution integral $f_Z(z) = \int f_X(t) f_Y(z-t)\,dt$ and understand what each part represents
3. Use moment generating functions (MGFs) to find distributions of sums without computing convolutions
4. Recognize which distribution families are closed under convolution (sums stay in the same family)
5. Connect these ideas to the Central Limit Theorem, understanding why sums of random variables become approximately normal
6. Apply sum distributions to AI/ML problems: mini-batch gradient averaging, ensemble predictions, and error aggregation
7. Implement convolution calculations in Python using NumPy and SciPy

The Historical Story: From Gambling to Deep Learning

The study of sums of random variables has a rich history that connects gambling, astronomy, and modern machine learning.

The Gamblers' Problem (17th Century)

When Blaise Pascal and Pierre de Fermat corresponded about gambling problems in 1654, they encountered questions like: "What's the probability that the sum of three dice exceeds 10?" This required understanding how probability distributions combine when random variables are added.

The key insight: to find $P(X + Y = k)$, you must consider all pairs $(x, y)$ where $x + y = k$ and sum their joint probabilities. This is the essence of convolution.

The Astronomical Challenge (18th-19th Century)

Astronomers like Gauss and Laplace needed to combine multiple measurement errors. If each measurement has its own error distribution, what's the distribution of the total error? They discovered that Gaussian (normal) errors have a beautiful property: the sum of Gaussians is Gaussian. This closure property made the normal distribution the centerpiece of statistical theory.

Modern Machine Learning

Today, sums of random variables appear constantly in deep learning:

• Stochastic Gradient Descent: Each mini-batch gradient is a sum of individual gradients
• Ensemble Methods: Averaging predictions from multiple models is summing random variables
• Dropout: The output of a dropout layer is a random sum of activations
• Attention Mechanisms: Weighted sums of value vectors involve random weights

"The whole is more predictable than its parts." — The Central Limit Theorem summarized

Why Sums of Random Variables Matter

Given two independent random variables $X$ and $Y$, we often need to find the distribution of their sum $Z = X + Y$. This is more complex than it might seem.

Simple Example: Two Dice

Let $X$ = first die roll, $Y$ = second die roll. What's $P(X + Y = 7)$?

We need to count all pairs $(x, y)$ where $x + y = 7$:

• (1, 6): $P(X=1) \cdot P(Y=6) = \frac{1}{6} \cdot \frac{1}{6} = \frac{1}{36}$
• (2, 5): $P(X=2) \cdot P(Y=5) = \frac{1}{36}$
• (3, 4): $P(X=3) \cdot P(Y=4) = \frac{1}{36}$
• (4, 3): $P(X=4) \cdot P(Y=3) = \frac{1}{36}$
• (5, 2): $P(X=5) \cdot P(Y=2) = \frac{1}{36}$
• (6, 1): $P(X=6) \cdot P(Y=1) = \frac{1}{36}$

Summing: $P(X + Y = 7) = 6 \times \frac{1}{36} = \frac{6}{36} = \frac{1}{6}$

This is exactly what the convolution formula does systematically.

Discrete Case: The Convolution Formula

Mathematical Definition

For two independent discrete random variables $X$ and $Y$, the PMF of $Z = X + Y$ is:

or equivalently:

The $*$ symbol denotes convolution.

Unpacking the Formula

• $p_Z(k)$: The probability that the sum equals $k$
• $\sum_x$: Sum over all possible values $x$ that $X$ can take
• $p_X(x)$: Probability that $X = x$
• $p_Y(k - x)$: Probability that $Y = k - x$ (so that $x + (k-x) = k$)
• Product: Uses independence: $P(X=x, Y=k-x) = P(X=x) \cdot P(Y=k-x)$
• Sum: Different $(x, k-x)$ pairs are mutually exclusive, so we add their probabilities

Interactive: Discrete Convolution

The visualization below shows how the convolution formula works step by step. Select a distribution type and a sum value $k$ to see which pairs of $(x, y)$ contribute to $P(X + Y = k)$.

Continuous Case: The Convolution Integral

Mathematical Definition

For two independent continuous random variables $X$ and $Y$, the PDF of $Z = X + Y$ is:

Unpacking the Integral

• $f_Z(z)$: The probability density of the sum at value $z$
• $\int_{-\infty}^{\infty} \cdots dt$: Integrate over all possible values $t$ that $X$ can take
• $f_X(t)$: Density of $X$ at value $t$
• $f_Y(z - t)$: Density of $Y$ at value $z - t$

Example: Sum of Two Uniform[0,1] Variables

Let $X, Y \sim \text{Uniform}(0, 1)$ independently. What is $f_{X+Y}(z)$?

The indicator function $\mathbf{1}_{0 \leq z-t \leq 1}$ is 1 only when $z-1 \leq t \leq z$. Combined with $0 \leq t \leq 1$, this gives:

This is the triangular distribution on [0, 2]. Two flat (uniform) distributions convolve to create a peaked distribution!

Interactive: Continuous Convolution

The visualization below shows the convolution integral in action. Adjust the value of $z$ to see how the integrand $f_X(t) \cdot f_Y(z-t)$ changes. The area under the integrand curve gives $f_Z(z)$.

The MGF Method: Convolution Becomes Multiplication

Convolution integrals can be difficult to compute. Fortunately, there's a powerful shortcut: the moment generating function (MGF).

The Key Theorem

For independent random variables $X$ and $Y$:

The MGF of the sum equals the product of the MGFs.

Why This Works

The third equality uses independence: for independent $X, Y$, $E[g(X) \cdot h(Y)] = E[g(X)] \cdot E[h(Y)]$.

Example: Sum of Normals

Let $X \sim N(\mu_1, \sigma_1^2)$ and $Y \sim N(\mu_2, \sigma_2^2)$ independently.

MGFs:

Product:

This is the MGF of $N(\mu_1 + \mu_2, \sigma_1^2 + \sigma_2^2)$. Therefore:

No convolution integral required!

Closure Properties: Which Sums Stay in the Family?

Some distribution families have the beautiful property that sums of independent members stay in the same family. These are called closed under convolution.

Central Limit Theorem Preview

The most profound result about sums of random variables is the Central Limit Theorem (CLT):

The standardized sum of $n$ independent, identically distributed (i.i.d.) random variables converges to a standard normal distribution as $n \to \infty$, regardless of the original distribution (given finite variance).

Formal Statement

Let $X_1, X_2, \ldots, X_n$ be i.i.d. with mean $\mu$ and variance $\sigma^2 < \infty$. Let $S_n = X_1 + X_2 + \cdots + X_n$. Then:

Why This Happens

The MGF approach gives insight into why CLT works:

1. Sum MGF: $M_{S_n}(t) = [M_X(t)]^n$ (product of n identical MGFs)
2. Taylor expansion: For small $t$, $M_X(t) \approx 1 + \mu t + \frac{\sigma^2 + \mu^2}{2} t^2 + \cdots$
3. Standardization: The standardized variable's MGF converges to $\exp(t^2/2)$, the MGF of $N(0,1)$

The full proof is covered in Chapter 10. For now, let's see CLT in action:

AI/ML Applications

1. Mini-batch Gradient Averaging

In stochastic gradient descent, we compute the gradient on a mini-batch of $n$ samples:

Each $g_i$ is the gradient for sample $i$. The average gradient $\bar{g}$ is a sum of random variables (divided by $n$). By CLT, this average has approximately normal distribution around the true gradient, with variance $\propto 1/n$.

2. Ensemble Model Predictions

When averaging predictions from $n$ models:

If each model's prediction $\hat{y}_i$ has variance $\sigma^2$, and models are independent, then:

Ensemble variance decreases with number of models! This is why ensemble methods (random forests, bagging) are so effective.

3. Dropout Regularization

During dropout, each neuron is randomly zeroed with probability $p$. The output is a random sum of Bernoulli-weighted activations:

By CLT, for many neurons, this sum is approximately normal. This helps explain dropout's regularization effect through the lens of sum distributions.

4. Aggregate Prediction Error

Total prediction error over a test set is a sum:

Understanding the distribution of this sum helps construct confidence intervals for model performance and decide if observed differences between models are statistically significant.

Python Implementation

Computing Convolutions

1import numpy as np
2from scipy import signal
3from scipy.stats import norm, expon, gamma
4import matplotlib.pyplot as plt
5
6# Discrete Convolution: Sum of two dice
7def discrete_convolution_dice():
8    """
9    Compute P(X + Y = k) where X, Y are fair dice.
10    Uses numpy convolution which implements:
11    (f * g)[k] = sum_j f[j] * g[k-j]
12    """
13    # PMF of a fair die: P(X=1) = P(X=2) = ... = P(X=6) = 1/6
14    pmf_die = np.array([1/6] * 6)  # Values 1-6
15
16    # Convolve to get PMF of X + Y
17    # This computes P(X+Y=k) for k = 2, 3, ..., 12
18    pmf_sum = np.convolve(pmf_die, pmf_die)
19
20    # Values range from 2 (1+1) to 12 (6+6)
21    values = np.arange(2, 13)
22
23    print("Distribution of X + Y (two dice):")
24    for k, p in zip(values, pmf_sum):
25        print(f"  P(X+Y = {k:2d}) = {p:.4f} = {int(p*36)}/36")
26
27    return values, pmf_sum
28
29# Continuous Convolution: Sum of two exponentials
30def continuous_convolution_exponential():
31    """
32    The sum of two independent Exp(λ) random variables
33    follows a Gamma(2, λ) distribution.
34
35    f_{X+Y}(z) = ∫ f_X(t) * f_Y(z-t) dt
36    """
37    lambda_rate = 1.0
38
39    # Create grid for visualization
40    z = np.linspace(0, 10, 1000)
41
42    # Analytical result: Gamma(2, λ) = λ² * z * exp(-λz)
43    pdf_sum_analytical = gamma.pdf(z, a=2, scale=1/lambda_rate)
44
45    # Numerical convolution for comparison
46    t = np.linspace(0, 10, 1000)
47    dt = t[1] - t[0]
48    f_x = expon.pdf(t, scale=1/lambda_rate)
49
50    pdf_sum_numerical = np.zeros_like(z)
51    for i, z_val in enumerate(z):
52        # f_{X+Y}(z) = ∫ f_X(t) * f_Y(z-t) dt
53        integrand = f_x * expon.pdf(z_val - t, scale=1/lambda_rate)
54        pdf_sum_numerical[i] = np.trapz(integrand, t)
55
56    print(f"\\nExp(1) + Exp(1) ~ Gamma(2, 1)")
57    print(f"Mean: E[X+Y] = 2/λ = 2.0")
58    print(f"Variance: Var(X+Y) = 2/λ² = 2.0")
59
60    return z, pdf_sum_analytical, pdf_sum_numerical
61
62# Sum of Normals: Stays Normal
63def sum_of_normals():
64    """
65    X ~ N(μ₁, σ₁²), Y ~ N(μ₂, σ₂²) independent
66    => X + Y ~ N(μ₁ + μ₂, σ₁² + σ₂²)
67    """
68    mu1, sigma1 = 2, 1
69    mu2, sigma2 = 3, 1.5
70
71    # Sum parameters
72    mu_sum = mu1 + mu2
73    sigma_sum = np.sqrt(sigma1**2 + sigma2**2)
74
75    x = np.linspace(-5, 15, 1000)
76
77    pdf_sum = norm.pdf(x, mu_sum, sigma_sum)
78
79    print(f"\\nN({mu1}, {sigma1}²) + N({mu2}, {sigma2}²)")
80    print(f"= N({mu_sum}, {sigma_sum**2:.2f})")
81    print(f"= N({mu_sum}, {sigma_sum:.3f}²)")
82
83    return x, pdf_sum
84
85# Run examples
86if __name__ == "__main__":
87    discrete_convolution_dice()
88    continuous_convolution_exponential()
89    sum_of_normals()

Using MGFs to Find Sum Distributions

1import numpy as np
2from scipy.stats import norm, expon, poisson
3
4def mgf_normal(t, mu, sigma):
5    """MGF of Normal(μ, σ²): M(t) = exp(μt + σ²t²/2)"""
6    return np.exp(mu * t + 0.5 * sigma**2 * t**2)
7
8def mgf_exponential(t, lam):
9    """MGF of Exponential(λ): M(t) = λ/(λ-t) for t < λ"""
10    if np.any(t >= lam):
11        return np.inf
12    return lam / (lam - t)
13
14def mgf_poisson(t, lam):
15    """MGF of Poisson(λ): M(t) = exp(λ(e^t - 1))"""
16    return np.exp(lam * (np.exp(t) - 1))
17
18# Demonstrate: MGF of sum = product of MGFs
19def verify_sum_property():
20    """
21    For independent X, Y:
22    M_{X+Y}(t) = E[e^{t(X+Y)}] = E[e^{tX}]E[e^{tY}] = M_X(t) * M_Y(t)
23    """
24    t = 0.5  # Evaluate at t = 0.5
25
26    # Example 1: Two normals
27    mu1, sigma1 = 1, 2
28    mu2, sigma2 = 3, 1
29
30    mgf_x = mgf_normal(t, mu1, sigma1)
31    mgf_y = mgf_normal(t, mu2, sigma2)
32    mgf_product = mgf_x * mgf_y
33
34    # Sum should be N(μ1+μ2, σ1²+σ2²)
35    mu_sum = mu1 + mu2
36    sigma_sum = np.sqrt(sigma1**2 + sigma2**2)
37    mgf_sum_direct = mgf_normal(t, mu_sum, sigma_sum)
38
39    print("Normal + Normal:")
40    print(f"  M_X({t}) × M_Y({t}) = {mgf_x:.4f} × {mgf_y:.4f} = {mgf_product:.4f}")
41    print(f"  M_{{X+Y}}({t}) directly = {mgf_sum_direct:.4f}")
42    print(f"  Match: {np.isclose(mgf_product, mgf_sum_direct)}")
43
44    # Example 2: Two Poissons
45    lam1, lam2 = 3, 5
46
47    mgf_p1 = mgf_poisson(t, lam1)
48    mgf_p2 = mgf_poisson(t, lam2)
49    mgf_p_product = mgf_p1 * mgf_p2
50
51    # Sum should be Poisson(λ1 + λ2)
52    mgf_p_sum = mgf_poisson(t, lam1 + lam2)
53
54    print("\\nPoisson + Poisson:")
55    print(f"  M_X({t}) × M_Y({t}) = {mgf_p1:.4f} × {mgf_p2:.4f} = {mgf_p_product:.4f}")
56    print(f"  M_{{X+Y}}({t}) directly = {mgf_p_sum:.4f}")
57    print(f"  Match: {np.isclose(mgf_p_product, mgf_p_sum)}")
58
59# Identify distribution from MGF
60def identify_from_mgf():
61    """
62    Once we have M_{X+Y}(t) = M_X(t) × M_Y(t),
63    we can often identify the sum's distribution by
64    recognizing the functional form.
65    """
66    print("\\nIdentifying distributions from MGF:")
67    print("  M(t) = exp(μt + σ²t²/2) → Normal(μ, σ²)")
68    print("  M(t) = (λ/(λ-t))^n → Gamma(n, λ)")
69    print("  M(t) = exp(λ(e^t - 1)) → Poisson(λ)")
70    print("  M(t) = (p·e^t / (1-q·e^t))^r → NegBinomial(r, p)")
71
72verify_sum_property()
73identify_from_mgf()

Demonstrating the Central Limit Theorem

1import numpy as np
2from scipy.stats import norm
3import matplotlib.pyplot as plt
4
5def demonstrate_clt(dist_sampler, n_samples, n_vars, true_mean, true_var):
6    """
7    Demonstrate CLT: (S_n - nμ) / √(nσ²) → N(0, 1)
8
9    Parameters:
10    - dist_sampler: Function that returns random samples
11    - n_samples: Number of sample means to generate
12    - n_vars: Number of variables in each sum
13    - true_mean, true_var: μ and σ² of the base distribution
14    """
15    sums = []
16    for _ in range(n_samples):
17        # Sum of n_vars independent random variables
18        sample_sum = np.sum(dist_sampler(n_vars))
19        sums.append(sample_sum)
20
21    sums = np.array(sums)
22
23    # Standardize: Z = (S_n - nμ) / √(nσ²)
24    standardized = (sums - n_vars * true_mean) / np.sqrt(n_vars * true_var)
25
26    # Compare to standard normal
27    x = np.linspace(-4, 4, 100)
28
29    print(f"CLT Check (n = {n_vars}):")
30    print(f"  Sample mean of standardized sums: {np.mean(standardized):.4f} (should be ≈ 0)")
31    print(f"  Sample std of standardized sums: {np.std(standardized):.4f} (should be ≈ 1)")
32
33    return standardized
34
35# Example: Uniform(0, 1) - definitely NOT normal!
36uniform_sampler = lambda n: np.random.uniform(0, 1, n)
37uniform_mean = 0.5
38uniform_var = 1/12
39
40print("Base distribution: Uniform(0, 1)")
41print(f"  Mean: {uniform_mean}, Variance: {uniform_var:.4f}")
42
43for n in [1, 2, 5, 10, 30]:
44    print(f"\\n--- n = {n} variables summed ---")
45    demonstrate_clt(uniform_sampler, 10000, n, uniform_mean, uniform_var)
46
47# Example: Exponential(1) - skewed!
48print("\\n" + "="*50)
49print("Base distribution: Exponential(1)")
50exp_sampler = lambda n: np.random.exponential(1, n)
51exp_mean = 1.0
52exp_var = 1.0
53
54for n in [1, 5, 30]:
55    print(f"\\n--- n = {n} variables summed ---")
56    demonstrate_clt(exp_sampler, 10000, n, exp_mean, exp_var)

Common Pitfalls

Summary

You've mastered one of the most important topics in probability theory: understanding how random variables combine when added. Here's what you've learned:

Core Concepts

• Discrete Convolution: $P(X+Y=k) = \sum_x P(X=x) \cdot P(Y=k-x)$
• Continuous Convolution: $f_{X+Y}(z) = \int f_X(t) \cdot f_Y(z-t)\,dt$
• MGF Method: $M_{X+Y}(t) = M_X(t) \cdot M_Y(t)$ for independent RVs
• Closure Properties: Some families (normal, Poisson, gamma) are closed under convolution
• CLT Preview: Sums of many independent RVs tend toward normal distribution

Key Formulas

Next Steps

In the next section, we'll study Order Statistics — the distribution of the minimum, maximum, and other ranked values from a sample. This is essential for understanding extreme values, percentiles, and robust statistics.