How to Calculate Moment Generating Functions

Key Insights

Moment generating functions (MGFs) provide a powerful alternative to directly computing expectations by transforming probability distributions into exponential form, making complex calculations—especially for sums of independent variables—algebraically simple.
The MGF uniquely identifies a distribution and encodes all moments through derivatives at zero: the nth derivative of M(t) at t=0 equals E[X^n], turning moment calculation into a differentiation problem.
While MGFs don’t exist for all distributions (heavy-tailed distributions often lack convergence), they excel at handling common distributions and proving distributional results through their multiplication property for independent sums.

Introduction to Moment Generating Functions

The moment generating function (MGF) of a random variable X is defined as:

M(t) = E[e^(tX)]

This deceptively simple definition transforms probability calculations into a different mathematical space where many operations become easier. The MGF exists when this expectation converges for t in some neighborhood around zero.

Why should you care about MGFs? Three reasons: First, they uniquely determine probability distributions—if two random variables have the same MGF, they have the same distribution. Second, they “generate” moments through differentiation, avoiding tedious integration. Third, they convert convolution operations (sums of independent variables) into simple multiplication.

Characteristic functions (φ(t) = E[e^(itX)]) serve a similar purpose but always exist due to the bounded complex exponential. MGFs are often more convenient when they exist because they avoid complex numbers.

import numpy as np
import matplotlib.pyplot as plt

# Visualize the e^(tX) transformation
X = np.linspace(-3, 3, 100)
t_values = [0, 0.5, 1.0, 1.5]

plt.figure(figsize=(10, 6))
for t in t_values:
    plt.plot(X, np.exp(t * X), label=f't={t}')
plt.xlabel('X')
plt.ylabel('e^(tX)')
plt.title('The Exponential Transformation in MGFs')
plt.legend()
plt.grid(True)
plt.show()

Calculating MGFs from Probability Distributions

For discrete random variables, the MGF is a sum:

M(t) = Σ e^(tx) P(X = x)

For continuous random variables, it’s an integral:

M(t) = ∫ e^(tx) f(x) dx

Let’s derive the MGF for a Bernoulli(p) distribution. X takes value 1 with probability p and 0 with probability 1-p:

M(t) = e^(t·0)(1-p) + e^(t·1)(p) = (1-p) + pe^t

For a Normal(μ, σ²) distribution, the calculation requires completing the square:

M(t) = ∫_{-∞}^{∞} e^(tx) · (1/√(2πσ²)) · e^(-(x-μ)²/(2σ²)) dx

After algebraic manipulation (completing the square in the exponent):

M(t) = exp(μt + σ²t²/2)

import numpy as np
from scipy import integrate

def mgf_discrete(values, probabilities, t):
    """Calculate MGF for discrete distribution."""
    return np.sum(probabilities * np.exp(t * values))

def mgf_continuous(pdf, support, t):
    """Calculate MGF for continuous distribution numerically."""
    def integrand(x):
        return np.exp(t * x) * pdf(x)
    
    result, error = integrate.quad(integrand, support[0], support[1])
    return result

# Example: Bernoulli distribution
p = 0.6
values = np.array([0, 1])
probs = np.array([1-p, p])
t = 0.5
print(f"Bernoulli MGF at t={t}: {mgf_discrete(values, probs, t):.4f}")
print(f"Analytical: {(1-p) + p*np.exp(t):.4f}")

# Example: Exponential distribution
lambda_param = 2.0
pdf_exp = lambda x: lambda_param * np.exp(-lambda_param * x)
t = 0.5
if t < lambda_param:  # MGF only exists for t < lambda
    mgf_val = mgf_continuous(pdf_exp, (0, np.inf), t)
    analytical = lambda_param / (lambda_param - t)
    print(f"Exponential MGF at t={t}: {mgf_val:.4f}")
    print(f"Analytical: {analytical:.4f}")

Deriving Moments from MGFs

The power of MGFs lies in this property: the nth derivative of M(t) evaluated at t=0 gives the nth moment:

E[X^n] = M^(n)(0)

This works because:

M’(t) = d/dt E[e^(tX)] = E[X·e^(tX)]

Setting t=0: M’(0) = E[X]

For variance, we need E[X²] - (E[X])²:

M’’(0) = E[X²]

Var(X) = M’’(0) - (M’(0))²

from sympy import symbols, exp, diff, lambdify

def extract_moments(mgf_expr, order=2):
    """Extract moments from symbolic MGF expression."""
    t = symbols('t')
    moments = {}
    
    for n in range(1, order + 1):
        derivative = diff(mgf_expr, t, n)
        moment = derivative.subs(t, 0)
        moments[f'E[X^{n}]'] = float(moment)
    
    mean = moments['E[X^1]']
    if order >= 2:
        variance = moments['E[X^2]'] - mean**2
        moments['Var(X)'] = variance
    
    return moments

# Example: Exponential(λ) has MGF = λ/(λ-t)
t, lam = symbols('t lambda', positive=True, real=True)
mgf_exp = lam / (lam - t)

# Substitute λ = 2
mgf_exp_2 = mgf_exp.subs(lam, 2)
moments = extract_moments(mgf_exp_2, order=2)

print("Exponential(2) moments:")
for key, value in moments.items():
    print(f"  {key}: {value:.4f}")

MGFs of Common Distributions

Here’s a reference table of MGFs for standard distributions:

Distribution	Parameters	MGF M(t)
Bernoulli	p	1-p + pe^t
Binomial	n, p	(1-p + pe^t)^n
Poisson	λ	exp(λ(e^t - 1))
Exponential	λ	λ/(λ-t), t < λ
Normal	μ, σ²	exp(μt + σ²t²/2)
Gamma	α, β	(1 - βt)^(-α), t < 1/β

Let’s verify these with scipy.stats:

from scipy import stats
import numpy as np

def verify_mgf(dist, mgf_func, t_value, params):
    """Verify MGF against scipy moments."""
    # Calculate moments from MGF derivatives (numerical)
    epsilon = 1e-6
    mgf_mean = (mgf_func(t_value + epsilon, *params) - 
                mgf_func(t_value - epsilon, *params)) / (2 * epsilon)
    
    # Get moments from scipy
    scipy_mean = dist(*params).mean()
    
    return mgf_mean, scipy_mean

# Poisson distribution
def poisson_mgf(t, lam):
    return np.exp(lam * (np.exp(t) - 1))

lam = 3.0
mgf_moment, scipy_moment = verify_mgf(
    stats.poisson, poisson_mgf, 0.0, (lam,)
)
print(f"Poisson(λ={lam}):")
print(f"  MGF-derived mean: {mgf_moment:.4f}")
print(f"  Scipy mean: {scipy_moment:.4f}")

# Normal distribution
def normal_mgf(t, mu, sigma):
    return np.exp(mu * t + 0.5 * sigma**2 * t**2)

mu, sigma = 5.0, 2.0
mgf_moment, scipy_moment = verify_mgf(
    stats.norm, normal_mgf, 0.0, (mu, sigma)
)
print(f"\nNormal(μ={mu}, σ={sigma}):")
print(f"  MGF-derived mean: {mgf_moment:.4f}")
print(f"  Scipy mean: {scipy_moment:.4f}")

Properties and Applications

The most powerful property of MGFs is for independent random variables X and Y:

M_{X+Y}(t) = M_X(t) · M_Y(t)

This turns convolution into multiplication. For example, if X ~ Normal(μ₁, σ₁²) and Y ~ Normal(μ₂, σ₂²) are independent:

M_{X+Y}(t) = exp(μ₁t + σ₁²t²/2) · exp(μ₂t + σ₂²t²/2) = exp((μ₁+μ₂)t + (σ₁²+σ₂²)t²/2)

This is the MGF of Normal(μ₁+μ₂, σ₁²+σ₂²), proving that sums of independent normals are normal.

import numpy as np
from scipy import stats

def demonstrate_mgf_convolution():
    """Show MGF multiplication for sum of independent variables."""
    # Two Poisson random variables
    lambda1, lambda2 = 3.0, 5.0
    
    # Generate samples
    n_samples = 10000
    X = stats.poisson(lambda1).rvs(n_samples)
    Y = stats.poisson(lambda2).rvs(n_samples)
    Z = X + Y
    
    # MGF of sum should equal MGF of Poisson(lambda1 + lambda2)
    def poisson_mgf(t, lam):
        return np.exp(lam * (np.exp(t) - 1))
    
    t = 0.3
    mgf_product = poisson_mgf(t, lambda1) * poisson_mgf(t, lambda2)
    mgf_sum = poisson_mgf(t, lambda1 + lambda2)
    
    print(f"Poisson({lambda1}) + Poisson({lambda2}):")
    print(f"  Product of MGFs: {mgf_product:.6f}")
    print(f"  MGF of Poisson({lambda1 + lambda2}): {mgf_sum:.6f}")
    print(f"  Sample mean of Z: {Z.mean():.2f}")
    print(f"  Theoretical mean: {lambda1 + lambda2:.2f}")

demonstrate_mgf_convolution()

Practical Implementation Tips

MGFs don’t always exist. Heavy-tailed distributions like Cauchy or log-normal may have MGFs that diverge for all t ≠ 0. Always check convergence:

import numpy as np
from scipy import integrate
import warnings

def safe_mgf_calculation(pdf, support, t, max_t=10):
    """Calculate MGF with convergence checking."""
    if abs(t) > max_t:
        return None, "t value too large"
    
    def integrand(x):
        return np.exp(t * x) * pdf(x)
    
    try:
        with warnings.catch_warnings():
            warnings.simplefilter("error")
            result, error = integrate.quad(
                integrand, support[0], support[1],
                limit=100
            )
            
        if error > 1e-3 * abs(result):
            return None, f"High integration error: {error}"
        
        if not np.isfinite(result):
            return None, "MGF diverges"
            
        return result, None
    
    except (integrate.IntegrationWarning, RuntimeWarning):
        return None, "Integration failed to converge"

# Test with exponential (should work)
lambda_param = 2.0
pdf_exp = lambda x: lambda_param * np.exp(-lambda_param * x)

for t in [0.5, 1.5, 1.9, 2.5]:
    result, error = safe_mgf_calculation(pdf_exp, (0, 20), t)
    if error:
        print(f"t={t}: {error}")
    else:
        analytical = lambda_param / (lambda_param - t) if t < lambda_param else np.inf
        print(f"t={t}: MGF={result:.4f}, Analytical={analytical:.4f}")

When MGFs fail, consider cumulant generating functions K(t) = log(M(t)) or stick with characteristic functions. For numerical work, always validate against known distributions and check for numerical stability.

The MGF is an elegant tool that transforms probability theory into calculus. Master it, and you’ll find distribution theory much more tractable.