BrevFeed
Sign in

11. Moments and moment generating functions

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

Expectation and variance are the first two entries in an infinite family called moments. Higher moments describe finer and finer features of a distribution's shape — and remarkably, a single function packages every one of them at once.

Moments

The k-th moment of X is simply

\mu_k' = E[X^k]

  • k = 1: E[X] — the mean, the centre.
  • k = 2: E[X^2] — feeds directly into variance via \operatorname{Var}(X) = E[X^2] - \mu^2.
  • k = 3, 4, \dots: progressively finer shape information.

More useful in practice are the central moments, taken about the mean:

\mu_k = E\big[(X - \mu)^k\big]

Then \mu_1 = 0 always, and \mu_2 = \operatorname{Var}(X).

The third and fourth get standardised into named shape statistics:

\text{skewness} = \frac{\mu_3}{\sigma^3}, \qquad \text{kurtosis} = \frac{\mu_4}{\sigma^4}

Skewness measures lopsidedness — exactly the asymmetry we saw in §0's salary data. Positive means a long right tail (mean pulled above the median); negative means a long left tail; zero for anything symmetric.

Kurtosis measures tail heaviness — how much probability sits far from the centre. A Normal has kurtosis 3, so excess kurtosis (kurtosis − 3) is reported instead, making the Normal the zero point.

Dividing by \sigma^k makes both dimensionless, so they don't change if you switch from rupees to dollars, or from metres to feet. That's the point of standardising.

The moment generating function

Here is the function that ties them all together:

M_X(t) = E\big[e^{tX}\big]

Plug in any real t and this one function has every moment of X hiding inside it.

Why it works

Expand e^{tX} as its Taylor series:

e^{tX} = 1 + tX + \frac{t^2X^2}{2!} + \frac{t^3X^3}{3!} + \cdots

Take expectations term by term (linearity again):

M_X(t) = 1 + t\,E[X] + \frac{t^2}{2!}E[X^2] + \frac{t^3}{3!}E[X^3] + \cdots

Every moment appears, each multiplied by a power of t. It's a power series whose coefficients are the moments.

Now differentiate k times and set t = 0. Every term below the k-th differentiates to zero, every term above still carries a factor of t and vanishes at 0, and the k! from differentiating t^k cancels the k! in the denominator — leaving exactly one survivor:

\boxed{\;M_X^{(k)}(0) = E[X^k]\;}

Differentiate at zero, and you dig the moments out one at a time. Hence the name: it generates moments.

The two properties that make it valuable

1. Uniqueness. If two random variables have the same MGF (on an interval around 0), they have the same distribution. This turns "prove these two things are identically distributed" — normally hard — into "compute two functions and compare", which is often routine.

2. Sums become products. If X and Y are independent:

M_{X+Y}(t) = E\big[e^{t(X+Y)}\big] = E\big[e^{tX}e^{tY}\big] = E\big[e^{tX}\big]E\big[e^{tY}\big] = M_X(t)\,M_Y(t)

The middle step is where independence is used, and it's the whole trick. Adding random variables — genuinely difficult, it's a convolution — turns into multiplying functions. Combined with uniqueness, this is how you prove that a sum of independent Normals is Normal, or that a sum of independent Poissons is Poisson.

Two other handy rules:

M_{aX+b}(t) = e^{bt}M_X(at), \qquad M_X(0) = 1 \;\text{always}

That last one is a free correctness check on any MGF you derive.

Worked example

Derive the MGF of a Bernoulli(p) variable and recover its mean and variance.

X is 1 with probability p and 0 with probability 1-p. Directly from the definition:

M_X(t) = E[e^{tX}] = e^{t \cdot 0}(1-p) + e^{t \cdot 1}p = (1-p) + pe^t

Check: M_X(0) = (1-p) + p = 1. ✓

Differentiate once:

M_X'(t) = pe^t \implies M_X'(0) = p = E[X] \quad\checkmark

Differentiate again:

M_X''(t) = pe^t \implies M_X''(0) = p = E[X^2]

That E[X^2] = p is not a coincidence — since X only takes values 0 and 1, X^2 = X identically.

Now the variance:

\operatorname{Var}(X) = E[X^2] - (E[X])^2 = p - p^2 = p(1-p) \quad\checkmark

Exactly the known Bernoulli variance, recovered mechanically.

Now get the Binomial for free. A Binomial(n, p) is a sum of n independent Bernoullis, so by the product rule:

M(t) = \big[(1-p) + pe^t\big]^n

Differentiating once at 0 gives np; a little more work gives \operatorname{Var} = np(1-p). Deriving those by summing \sum_k k\binom{n}{k}p^k(1-p)^{n-k} is a genuinely unpleasant exercise; via the MGF it's a one-liner.

The catch

The MGF has a real weakness: it may not exist. E[e^{tX}] is an integral that can diverge, and for heavy-tailed distributions it does — for every t \ne 0.

The standard example is the Cauchy distribution, whose MGF is infinite for all t \ne 0. Its mean doesn't exist either, so there are no moments to generate.

When the MGF fails, the characteristic function takes over — that's the next lesson.

Doing it in Python

You can compute an MGF numerically and differentiate it to recover moments:

import numpy as np

p = 0.3
M = lambda t: (1 - p) + p * np.exp(t)          # Bernoulli MGF

def deriv(f, x, k, h=1e-4):
    """Numerical k-th derivative by central differences."""
    if k == 0:
        return f(x)
    return (deriv(f, x + h, k - 1) - deriv(f, x - h, k - 1)) / (2 * h)

print("M(0)   =", round(M(0), 6), "(must be 1)")
print("M'(0)  =", round(deriv(M, 0.0, 1), 4), " -> E[X]   (true", p, ")")
print("M''(0) =", round(deriv(M, 0.0, 2), 4), " -> E[X^2] (true", p, ")")

mean = deriv(M, 0.0, 1)
second = deriv(M, 0.0, 2)
print("Var    =", round(second - mean**2, 4), " (true", round(p * (1 - p), 4), ")")

Symbolically it's exact, and SymPy will do the differentiation for you:

import sympy as sp

t, p = sp.symbols("t p", positive=True)
M = (1 - p) + p * sp.exp(t)

m1 = sp.diff(M, t).subs(t, 0)
m2 = sp.diff(M, t, 2).subs(t, 0)
print("E[X]   =", sp.simplify(m1))
print("E[X^2] =", sp.simplify(m2))
print("Var(X) =", sp.simplify(m2 - m1**2))

And SciPy reports skewness and kurtosis directly:

from scipy.stats import norm, expon

for name, dist in [("normal", norm()), ("exponential", expon())]:
    mean, var, skew, kurt = dist.stats(moments="mvsk")
    print(f"{name:12} mean={float(mean):6.3f} var={float(var):6.3f} "
          f"skew={float(skew):6.3f} excess kurt={float(kurt):6.3f}")

The Normal shows skew 0 and excess kurtosis 0 — it's the reference point. The Exponential shows skew 2, reflecting its long right tail.

Your turn

1. Find the MGF of a constant X = c, and use it to get E[X] and \operatorname{Var}(X).

2. X is uniform on \{1, 2, 3\}. Write its MGF and find E[X].

3. If X and Y are independent with MGFs M_X and M_Y, what's the MGF of 2X + 3?

Solutions

1. X = c with probability 1, so

M(t) = E[e^{tX}] = e^{tc}

Then M'(t) = ce^{tc}, so M'(0) = c = E[X]. ✓

And M''(t) = c^2e^{tc}, so M''(0) = c^2 = E[X^2], giving

\operatorname{Var}(X) = c^2 - c^2 = 0

A constant has zero variance — as it must, since it never strays from its mean.

2. Each value has probability 1/3:

M(t) = \tfrac{1}{3}\left(e^{t} + e^{2t} + e^{3t}\right)

Differentiate:

M'(t) = \tfrac{1}{3}\left(e^{t} + 2e^{2t} + 3e^{3t}\right) \implies M'(0) = \tfrac{1 + 2 + 3}{3} = 2

Which is just the average of 1, 2, 3 — as expected for a uniform distribution.

3. Use the affine rule M_{aX+b}(t) = e^{bt}M_X(at) with a = 2, b = 3:

M_{2X+3}(t) = e^{3t}M_X(2t)

Note Y is irrelevant — it doesn't appear in 2X + 3. (The trap here is reaching for the product rule out of habit; it applies to sums of independent variables, not to any expression that happens to mention two of them.)

Check yourself in code

Recover the mean and variance of a Bernoulli(0.3) from its MGF by numerical differentiation at t = 0, and check them against the known formulas p and p(1-p).

Print exactly this:

M(0) = 1.0
E[X] = 0.3
Var(X) = 0.21
matches formulas: True

Round every value to 4 decimal places before printing (so 0.3 and 0.21 appear without trailing noise), and compare to 4 decimal places.

import numpy as np

p = 0.3
M = lambda t: (1 - p) + p * np.exp(t)

def deriv(f, x, k, h=1e-4):
    if k == 0:
        return f(x)
    return (deriv(f, x + h, k - 1, h) - deriv(f, x - h, k - 1, h)) / (2 * h)

print("M(0) =", round(M(0.0), 4))

# Recover E[X] = M'(0) and E[X^2] = M''(0), form the variance,
# and compare both against p and p*(1-p).

import numpy as np

p = 0.3
M = lambda t: (1 - p) + p * np.exp(t)

def deriv(f, x, k, h=1e-4):
    if k == 0:
        return f(x)
    return (deriv(f, x + h, k - 1, h) - deriv(f, x - h, k - 1, h)) / (2 * h)

print("M(0) =", round(M(0.0), 4))

mean = deriv(M, 0.0, 1)
second = deriv(M, 0.0, 2)
var = second - mean**2

print("E[X] =", round(mean, 4))
print("Var(X) =", round(var, 4))
print("matches formulas:", round(mean, 4) == round(p, 4)
      and round(var, 4) == round(p * (1 - p), 4))

Every moment of X is buried inside M(t). Differentiate at zero and you dig them out one at a time. Match the whole function and you've matched the whole distribution — and because sums of independent variables become products, the MGF turns the hardest routine operation in probability into multiplication.

Next: the fix for the one thing the MGF can't do — exist for every distribution.

← Previous Next →