Encode a distribution as a polynomial: g(z) = ∑ P(X = k) zᵏ. Adding independent random variables becomes multiplying their polynomials. Convolution is just polynomial multiplication.
Ordinary generating function
For a discrete random variable X taking values 0, 1, 2, …, the ordinary generating function (OGF) is g(z) = E[zᴺ] = ∑ₖ P(X = k) zᵏ. The coefficients are the probabilities. Evaluating g(1) = 1 (probabilities sum to 1). The derivative g′(1) = E[X].
Scheme
; Generating functions as polynomial coefficient lists; g(z) = [p0, p1, p2, ...] means p0 + p1*z + p2*z^2 + ...; Multiply two polynomials (= convolve two distributions)
(define (poly-mul p q)
(let* ((m (length p)) (n (length q))
(result-len (- (+ m n) 1)))
(let loop ((k 0) (acc '()))
(if (= k result-len) (reverse acc)
(let inner ((j 0) (s 0))
(if (= j m) (loop (+ k 1) (cons s acc))
(if (and (>= (- k j) 0) (< (- k j) n))
(inner (+ j 1) (+ s (* (list-ref p j) (list-ref q (- k j)))))
(inner (+ j 1) s))))))))
; Fair coin: P(X=0) = 1/2, P(X=1) = 1/2
(define coin '(0.50.5))
; Sum of 2 coins
(define two-coins (poly-mul coin coin))
(display "2 coins: ") (display two-coins) (newline)
; (0.25 0.5 0.25) => P(0)=1/4, P(1)=1/2, P(2)=1/4; Sum of 3 coins
(define three-coins (poly-mul two-coins coin))
(display "3 coins: ") (display three-coins) (newline)
; (0.125 0.375 0.375 0.125); Fair die: P(X=k) = 1/6 for k in 1..6
(define die '(00.16670.16670.16670.16670.16670.1667))
(define two-dice (poly-mul die die))
(display "2 dice: ") (display two-dice)
; Coefficients give P(sum=2), P(sum=3), ..., P(sum=12)
Python
import numpy as np
def poly_mul(p, q):
"""Multiply polynomials = convolve distributions."""returnlist(np.convolve(p, q))
coin = [0.5, 0.5]
print("2 coins:", [round(x, 4) for x in poly_mul(coin, coin)])
print("3 coins:", [round(x, 4) for x in poly_mul(poly_mul(coin, coin), coin)])
die = [0] + [1/6]*6
two_dice = poly_mul(die, die)
for s inrange(2, 13):
print(f" P(sum={s:>2}) = {two_dice[s]:.4f}")
Moments from derivatives
The generating function encodes all moments. E[X] = g′(1). E[X(X−1)] = g″(1), so Var(X) = g″(1) + g′(1) − (g′(1))². No summation needed: just differentiate the polynomial and evaluate at z = 1.
Scheme
; Extract moments from the generating function; g'(1) = E[X], g''(1) = E[X(X-1)]; Differentiate a polynomial
(define (poly-deriv p)
(let loop ((i 1) (acc '()))
(if (>= i (length p)) (reverse acc)
(loop (+ i 1) (cons (* i (list-ref p i)) acc)))))
; Evaluate polynomial at z
(define (poly-eval p z)
(let loop ((i 0) (acc 0) (zk 1))
(if (>= i (length p)) acc
(loop (+ i 1) (+ acc (* (list-ref p i) zk)) (* zk z)))))
; Fair die generating function
(define die '(00.16670.16670.16670.16670.16670.1667))
(define die-prime (poly-deriv die))
(define die-double (poly-deriv die-prime))
(define mean (poly-eval die-prime 1))
(define factorial-moment (poly-eval die-double 1))
(define variance (+ factorial-moment mean (- (* mean mean))))
(display "E[X] = ") (display mean) (newline)
(display "Var(X) = ") (display variance) (newline)
; E[X] ~ 3.5, Var(X) ~ 2.917
Moment generating function
The moment generating function (MGF) is M(t) = E[eᵗᴺ]. It is the OGF evaluated at z = eᵗ. Its advantage: M′(0) = E[X], M″(0) = E[X²], and in general M⁽ᵏ⁾(0) = E[Xᵏ]. The MGF of a sum of independent variables is the product of their MGFs, exactly like the OGF.
Polynomials are stored as coefficient lists indexed from 0. The polynomial multiplication is exactly discrete convolution. For continuous distributions, the OGF becomes an integral transform, but the algebraic structure is the same: independence of random variables corresponds to multiplication of their generating functions. The textbook uses this machinery to prove the CLT: take the log of the MGF, expand to second order, and recover the normal MGF.
Ready for the real thing? Read Chapter 10 of Grinstead & Snell.
Probability-generating function
The Handshake