Skip to content
VibeFormer
Advanced24 min

Sums of Random Variables and Convolution

Distribution of a sum via convolution and MGFs, and the closure properties of common families.

Sums of Random Variables and Convolution

Intuition first

Sums are everywhere: total claim cost, cumulative waiting time, aggregate demand, the numerator of every sample mean. So the distribution of X+YX + Y is one of the most frequently needed objects in probability.

Its mean and variance are easy — linearity and independence handle them. The distribution is harder, because X+Y=sX + Y = s can happen in many ways: XX small and YY large, or the reverse, or anything between. Summing over all those combinations is what convolution does.

There are two routes. Convolve directly, which is an integral (or a sum) and is often unpleasant. Or multiply moment generating functions, which turns the problem into algebra — and for the standard families lets you recognise the answer instantly. The second route is why MGFs earn their keep.

Convolution

Discrete:

P(X+Y=s)=kP(X=k)P(Y=sk)\Prob(X + Y = s) = \sum_k \Prob(X = k)\,\Prob(Y = s - k)

Continuous:

fX+Y(s)=fX(x)fY(sx)dxf_{X+Y}(s) = \int_{-\infty}^{\infty} f_X(x)\,f_Y(s - x)\,dx

Both assume independence.

Solved problem 1 · Two dice, by discrete convolution

Two fair dice. Find the distribution of the sum by convolution.

Step 1 — set up

Each die has P(X=k)=1/6\Prob(X = k) = 1/6 for k=1,,6k = 1, \dots, 6. For the sum ss:

P(S=s)=kP(X=k)P(Y=sk)=k1616\Prob(S = s) = \sum_{k} \Prob(X = k)\Prob(Y = s - k) = \sum_k \frac{1}{6}\cdot\frac{1}{6}

over all kk with both 1k61 \leq k \leq 6 and 1sk61 \leq s - k \leq 6.

Step 2 — determine the valid range of k

The two constraints give

max(1,s6)kmin(6,s1)\max(1,\, s-6) \leq k \leq \min(6,\, s-1)

The number of valid kk is the number of terms, each contributing 1/361/36.

Step 3 — count terms for a few values

s=2:  k{1}1 term,  P=1/36s = 2: \; k \in \{1\} \quad\Longrightarrow\quad 1 \text{ term}, \; \Prob = 1/36s=5:  k{1,2,3,4}4 terms,  P=4/36s = 5: \; k \in \{1,2,3,4\} \quad\Longrightarrow\quad 4 \text{ terms}, \; \Prob = 4/36s=7:  k{1,,6}6 terms,  P=6/36s = 7: \; k \in \{1,\dots,6\} \quad\Longrightarrow\quad 6 \text{ terms}, \; \Prob = 6/36s=10:  k{4,5,6}3 terms,  P=3/36s = 10: \; k \in \{4,5,6\} \quad\Longrightarrow\quad 3 \text{ terms}, \; \Prob = 3/36

Step 4 — the general formula

P(S=s)=6s736,s=2,,12\Prob(S = s) = \frac{6 - \lvert s - 7 \rvert}{36}, \qquad s = 2, \dots, 12

Step 5 — verify

1+2+3+4+5+6+5+4+3+2+136=3636=1  \frac{1 + 2 + 3 + 4 + 5 + 6 + 5 + 4 + 3 + 2 + 1}{36} = \frac{36}{36} = 1 \;\checkmark

The triangular shape is the convolution of two rectangles — a general fact, and the reason sums of uniforms become bell-shaped surprisingly quickly.

Answer

P(S=s)=6s736\Prob(S = s) = \dfrac{6 - \lvert s-7\rvert}{36}, peaking at 6/366/36 for s=7s = 7.

Closure properties

Some families are closed under addition of independent members. These are worth memorising, because recognising one saves the whole convolution.

Sum of independentResultCondition
Bin(n1,p)+Bin(n2,p)\text{Bin}(n_1,p) + \text{Bin}(n_2,p)Bin(n1+n2,p)\text{Bin}(n_1+n_2,\, p)Same pp
Poisson(λ1)+Poisson(λ2)\text{Poisson}(\lambda_1) + \text{Poisson}(\lambda_2)Poisson(λ1+λ2)\text{Poisson}(\lambda_1+\lambda_2)Always
N(μ1,σ12)+N(μ2,σ22)\mathcal{N}(\mu_1,\sigma_1^2) + \mathcal{N}(\mu_2,\sigma_2^2)N(μ1+μ2,σ12+σ22)\mathcal{N}(\mu_1+\mu_2,\, \sigma_1^2+\sigma_2^2)Always
Gamma(α1,β)+Gamma(α2,β)\text{Gamma}(\alpha_1,\beta) + \text{Gamma}(\alpha_2,\beta)Gamma(α1+α2,β)\text{Gamma}(\alpha_1+\alpha_2,\, \beta)Same β\beta
χk12+χk22\chi^2_{k_1} + \chi^2_{k_2}χk1+k22\chi^2_{k_1+k_2}Always
Exp(λ)×α\text{Exp}(\lambda) \times \alpha termsGamma(α,1/λ)\text{Gamma}(\alpha,\, 1/\lambda)Same λ\lambda
Uniform+Uniform\text{Uniform} + \text{Uniform}Triangular, not uniform
Proving closure with MGFsAdvanced

For independent XX and YY, MX+Y(t)=MX(t)MY(t)M_{X+Y}(t) = M_X(t)M_Y(t). Since an MGF determines a distribution uniquely, recognising the product identifies the sum.

Poisson. With M(t)=exp(λ(et1))M(t) = \exp(\lambda(e^t-1)):

exp ⁣(λ1(et1))exp ⁣(λ2(et1))=exp ⁣((λ1+λ2)(et1))\exp\!\big(\lambda_1(e^t-1)\big)\exp\!\big(\lambda_2(e^t-1)\big) = \exp\!\big((\lambda_1+\lambda_2)(e^t-1)\big)

which is Poisson(λ1+λ2)(\lambda_1+\lambda_2).

Gamma. With M(t)=(1βt)αM(t) = (1-\beta t)^{-\alpha}:

(1βt)α1(1βt)α2=(1βt)(α1+α2)(1-\beta t)^{-\alpha_1}(1-\beta t)^{-\alpha_2} = (1-\beta t)^{-(\alpha_1+\alpha_2)}

Closure holds because the bases match — which is exactly why the same scale β\beta is required. With different scales the product is not of gamma form, and the sum is not gamma.

Why binomial needs the same pp. With M(t)=(1p+pet)nM(t) = (1-p+pe^t)^n:

(1p+pet)n1(1p+pet)n2=(1p+pet)n1+n2(1-p+pe^t)^{n_1}(1-p+pe^t)^{n_2} = (1-p+pe^t)^{n_1+n_2}

Again the bases must be identical. Different pp gives a product of two different bases, which is not a binomial MGF — the sum is then a Poisson-binomial distribution with no closed form.

Solved problem 2 · Sum of two independent uniforms

X,YUniform(0,1)X, Y \sim \text{Uniform}(0,1) independent. Find the density of S=X+YS = X + Y by convolution.

Step 1 — write the convolution integral

fS(s)=fX(x)fY(sx)dxf_S(s) = \int_{-\infty}^{\infty} f_X(x)\,f_Y(s-x)\,dx

Both densities equal 1 on (0,1)(0,1) and 0 elsewhere, so the integrand is 1 exactly when

0<x<1and0<sx<10 < x < 1 \quad\text{and}\quad 0 < s - x < 1

The second condition rearranges to s1<x<ss - 1 < x < s.

Step 2 — intersect the constraints

max(0,s1)<x<min(1,s)\max(0,\, s-1) < x < \min(1,\, s)

The integral of 1 over this interval is simply its length.

Step 3 — case s in (0,1)

Here max(0,s1)=0\max(0, s-1) = 0 and min(1,s)=s\min(1,s) = s, so the interval is (0,s)(0, s) with length ss:

fS(s)=sf_S(s) = s

Step 4 — case s in [1,2)

Here max(0,s1)=s1\max(0,s-1) = s-1 and min(1,s)=1\min(1,s) = 1, so the interval is (s1,1)(s-1, 1) with length 2s2 - s:

fS(s)=2sf_S(s) = 2 - s

Step 5 — assemble and check

fS(s)={s0s12s1<s20otherwisef_S(s) = \begin{cases} s & 0 \leq s \leq 1 \\ 2 - s & 1 < s \leq 2 \\ 0 & \text{otherwise}\end{cases}

A triangle peaking at s=1s = 1. Normalisation is the area of a triangle with base 2 and height 1:

12×2×1=1  \tfrac12 \times 2 \times 1 = 1 \;\checkmark

Step 6 — moments, as an independent check

E[S]=E[X]+E[Y]=0.5+0.5=1  (the peak, by symmetry)\E[S] = \E[X] + \E[Y] = 0.5 + 0.5 = 1 \;\checkmark \text{(the peak, by symmetry)}Var(S)=112+112=160.1667\Var(S) = \tfrac{1}{12} + \tfrac{1}{12} = \tfrac{1}{6} \approx 0.1667

Answer

Triangular on [0,2][0,2]: fS(s)=sf_S(s) = s for s1s \leq 1 and 2s2-s for s>1s > 1. Two flat distributions convolve into a peaked one — the first step of the CLT visible in one calculation.

Sums of dependent variables

Convolution requires independence. Without it, the mean still adds by linearity, but the variance needs the covariance term and the distribution requires the full joint:

E[X+Y]=E[X]+E[Y](always)\E[X+Y] = \E[X] + \E[Y] \qquad \text{(always)} Var(X+Y)=Var(X)+Var(Y)+2Cov(X,Y)\Var(X+Y) = \Var(X) + \Var(Y) + 2\Cov(X,Y)

A random number of terms

For S=i=1NXiS = \sum_{i=1}^{N} X_i with NN random and independent of the XiX_i — a compound distribution — condition on NN:

E[S]=E[N]E[X],Var(S)=E[N]Var(X)+Var(N)(E[X])2\E[S] = \E[N]\,\E[X], \qquad \Var(S) = \E[N]\Var(X) + \Var(N)\big(\E[X]\big)^2
Deriving the compound varianceAdvanced

Apply the variance decomposition from the conditional expectation lesson. Given N=nN = n, the sum of nn i.i.d. terms has

E[SN]=NE[X],Var(SN)=NVar(X)\E[S \mid N] = N\,\E[X], \qquad \Var(S \mid N) = N\Var(X)

Now:

E[Var(SN)]=E[N]Var(X)\E\big[\Var(S\mid N)\big] = \E[N]\Var(X)Var(E[SN])=Var(NE[X])=(E[X])2Var(N)\Var\big(\E[S\mid N]\big) = \Var\big(N\E[X]\big) = \big(\E[X]\big)^2\Var(N)

Adding gives the stated result. The first term is variability in the individual claim sizes, the second variability in how many claims arrive.

This is the foundation of insurance risk modelling: total annual claims is a compound distribution, usually with NN Poisson (giving Var(N)=E[N]\Var(N) = \E[N]) and XX gamma or lognormal.

python
import numpy as np
from scipy import stats

rng = np.random.default_rng(0)

# Two dice: convolution against simulation.
pmf = np.array([(6 - abs(s - 7)) / 36 for s in range(2, 13)])
rolls = rng.integers(1, 7, (500_000, 2)).sum(axis=1)
print("s   theory   empirical")
for i, s in enumerate(range(2, 13)):
    print(f"{s:2d}  {pmf[i]:.5f}  {np.mean(rolls == s):.5f}")

# Sum of two uniforms is triangular.
S = rng.random(500_000) + rng.random(500_000)
print(f"\nsum of 2 uniforms: mean {S.mean():.4f} (1.0)  var {S.var():.4f} ({1/6:.4f})")
print(f"density near s=1 is ~2x density near s=0.5:")
for s in (0.25, 0.5, 1.0):
    band = np.mean(np.abs(S - s) < 0.01) / 0.02
    print(f"  s={s}: empirical {band:.3f}   theory {s if s <= 1 else 2-s:.3f}")

# Closure: same p works, different p does not.
a = rng.binomial(10, 0.5, 400_000) + rng.binomial(10, 0.5, 400_000)
b = rng.binomial(10, 0.3, 400_000) + rng.binomial(10, 0.7, 400_000)
print(f"\nBin(10,.5)+Bin(10,.5): var {a.var():.3f}  Bin(20,.5) var {20*0.25:.3f}  <- matches")
print(f"Bin(10,.3)+Bin(10,.7): var {b.var():.3f}  Bin(20,.5) var {20*0.25:.3f}  <- does NOT")

# Compound: total claims with Poisson count and gamma severity.
N = rng.poisson(12, 400_000)
total = np.array([rng.gamma(2, 500, n).sum() if n else 0.0 for n in N[:20_000]])
EX, VarX, EN = 2*500, 2*500**2, 12
print(f"\ncompound: mean {total.mean():,.0f}  theory {EN*EX:,.0f}")
print(f"          var  {total.var():,.0f}  theory {EN*VarX + EN*EX**2:,.0f}")

Exercise 1

XPoisson(3)X \sim \text{Poisson}(3) and YPoisson(5)Y \sim \text{Poisson}(5), independent. Find P(X+Y=2)\Prob(X + Y = 2).

Show solution

By the closure property, X+YPoisson(3+5)=Poisson(8)X + Y \sim \text{Poisson}(3 + 5) = \text{Poisson}(8).

P(X+Y=2)=e8822!=e8×642=32e8\Prob(X+Y = 2) = \frac{e^{-8}\,8^2}{2!} = \frac{e^{-8}\times 64}{2} = 32\,e^{-8}e8=0.00033546e^{-8} = 0.00033546=32×0.000335460.010735= 32 \times 0.00033546 \approx 0.010735

About 1.07%.

Verifying by direct convolution — three terms, since XX and YY must sum to 2:

P(X=0)P(Y=2)+P(X=1)P(Y=1)+P(X=2)P(Y=0)\Prob(X=0)\Prob(Y=2) + \Prob(X=1)\Prob(Y=1) + \Prob(X=2)\Prob(Y=0)=e325e52+3e35e5+9e32e5= e^{-3}\cdot\frac{25e^{-5}}{2} + 3e^{-3}\cdot 5e^{-5} + \frac{9e^{-3}}{2}\cdot e^{-5}=e8(12.5+15+4.5)=32e8  = e^{-8}\left(12.5 + 15 + 4.5\right) = 32e^{-8} \;\checkmark

The closure property replaced three terms with one — and would replace 101 terms with one if we wanted P(X+Y=100)\Prob(X+Y = 100).

Exercise 2

An insurer expects 20 claims per year (Poisson) with claim sizes averaging £4,000 and standard deviation £6,000. Find the mean and standard deviation of annual total claims.

Show solution

A compound Poisson. With E[N]=Var(N)=20\E[N] = \Var(N) = 20 (Poisson), E[X]=4000\E[X] = 4000, Var(X)=60002=36,000,000\Var(X) = 6000^2 = 36{,}000{,}000:

Mean:

E[S]=E[N]E[X]=20×4000=£80,000\E[S] = \E[N]\E[X] = 20 \times 4000 = £80{,}000

Variance:

Var(S)=E[N]Var(X)+Var(N)(E[X])2\Var(S) = \E[N]\Var(X) + \Var(N)\big(\E[X]\big)^2=20(36,000,000)+20(16,000,000)= 20(36{,}000{,}000) + 20(16{,}000{,}000)=720,000,000+320,000,000=1,040,000,000= 720{,}000{,}000 + 320{,}000{,}000 = 1{,}040{,}000{,}000SD(S)=1.04×109£32,249\text{SD}(S) = \sqrt{1.04\times10^9} \approx £32{,}249

Two things worth noting. The coefficient of variation is 32,249/80,0000.4032{,}249/80{,}000 \approx 0.40, so year-to-year totals swing substantially — a bad year of £145,000 is under a 2-standard-deviation event.

And the decomposition is informative: 69% of the variance comes from uncertainty in claim sizes and 31% from uncertainty in the number of claims. Reducing severity variability — through policy limits or reinsurance — would cut total risk more than stabilising claim frequency.

For a Poisson count the formula simplifies to Var(S)=E[N]E[X2]\Var(S) = \E[N]\,\E[X^2], which is a useful check: 20×(36,000,000+16,000,000)=20×52,000,000=1.04×10920 \times (36{,}000{,}000 + 16{,}000{,}000) = 20 \times 52{,}000{,}000 = 1.04\times10^9 ✓.


Next: The Multivariate Normal Distribution, the last distribution in the module.