What is a random variable?
Let's start with two problems.
- You're an insurance company that provides car and house insurance. You have a forecast of a range of claims for cars and a forecast of a range of claims for houses. How do you work out the overall range of claims' forecast?
- You are a manufacturing company producing widgets in a multi-stage process. For each stage, you know the range of defect rates. How do you combine them to get an overall range of defect rates?
As with probability distributions, there are two types of random variables: continuous and discrete. I'm going to focus on continuous random variables here.
(Gemini)
Random variables in pictures
Let's imagine we're modeling car insurance claims. What's the next claim value going to be? We don't know exactly what it is (remember, the random part), but we do know it follows a probability distribution and very often, we know what the probability distribution is. Here's a normal probability distribution I'm going to use for illustrative purposes.
If our insurance claims follows this normal probability distribution, we don't know what the next claim will be, but we can assign it probabilities. For example, there's a 25% probability it will be between $3,026 and $4,371.
Some notation and definitions
A continuous random variable \(X\) is made up of a set of measurements \(\{x_0, x_1, \ldots\}\) that follow a probability distribution \(f(x)\) and a cumulative distribution function \(F\). Here are some things that must be true:
- It must follow a probability distribution function (PDF) defined as \(P(a \le X \le b) = \int_a^b f(x), dx\). The PDF must have these properties:
- \(f(x) \ge 0\) for all \(x\)
- The total area "under the curve" for the PDF is 1, or: \(\int_{-\infty}^{\infty} f(x), dx = 1\)
- There must be a cumulative distribution function (CDF) defined as \(F(x) = P(X \le x)\)
The expectation value is the mean value and is denoted by \(E(X)\). The variance is denoted by \(\text{Var}(X)\).
Summing random variables
Let's go back to the insurance company selling car and house insurance. Let's assume for now that the chart on the left is the probability distribution for the value of car claims in one year and the chart on the right is the value of house claims in a year. We'll also assume car and house claims are independent. How do we get the total claims distribution?
The problem is, any car claim value can combine with any house claim value. The depth of the problem's a bit easier to understand if we switch to the discrete domain. Let's say there are only three possible claim values for the house (\(f\)) and claims insurance (\(g\)) If every house insurance possibility can occur with any car insurance claim, then we can write the sum of the random variables as:
\[(f*g)[z] = \sum_{m=0}^{2} f[m]\,g[z-m]\]
That's fine in the discrete case, but there are an infinite number of values in the continuous case; we would have to multiply an infinite number of numbers by an infinite number of numbers, Or would we? It turns out, there's a short-cut that transforms convolution into multiplication.
Let's start with some basics. If we represent the sum as \(S = X_0 + X_1 + ...\) then:
\[ S_n = \sum_{i=1}^{n} X_i \]
by definition. By some math I won't go into, the mean (expectation value) is always the given by the sum of the expectation values:
\[E[S_n] = \sum_{i=1}^{n} E[X_i]\]
the formula for the sum of the variances is a little more complicated. If all the \(X\) variables are independent, then:
\[\text{Var}(S_n) = \sum_{i=1}^{n} \text{Var}(X_i)\]
if they aren't, then:
\[\text{Var}(S_n) = \sum_{i=1}^n \text{Var}(X_i) + 2\sum_{i<j} \text{Cov}(X_i, X_j)\]
So far, this is standard. Let's say each \(X\) has a minimum value of \(min_i\) and a maximum of \(max_i\) then:
\[\sum_{i=1}^{n} \min_i \ \le\ S_n \ \le\ \sum_{i=1}^{n} \max_i\]
which is a fancy way of saying the minimum is the sum of the minimum values and the maximum is the sum of the maximum values.
We also know that:
\[\int_{-\infty}^{\infty} f_{S_n}(z)\, dz = 1\]
The bottom line is, we have some information about the sum probability distribution, but we don't have the actual probability distribution of the sum. To get the probability distribution of the sum, we're going to turn to convolution.
Convolution to get the combined probability distribution
In probability theory, there are a couple of quantities defined for probability distributions. You can use these quantities to reason about the behavior of probability distributions and do math with them. There are two quantities that concern us: the moment generating function and the characteristic function.
- The moment generating function encodes the mean, the variance, the skewness etc. (which are called moments). It's a convenient packaging of a distribution that you can use to get the moments and do math with distributions. The problem is, it doesn't always exist for all probability functions.
- The characteristic function completely characterizes the probability distribution; you can get moments from it and do math with probability distributions, but it is harder to work with than the moment generating functions. By contrast with the moment generating function, it exists for all probability functions.
unless you're deeply into math, you probably won't recognize the equation on the right, so let me tell you what it is; it's the Fourier Transform of the probability distribution. To complete the definitions, you can get \(f\) from \(\varphi\) via the inverse Fourier Transform:
\[f_X(x) = \frac{1}{2\pi} \int_{-\infty}^{\infty} e^{-itx}\, \varphi_X(t)\, dt\]
I came across Fourier Transforms in Physics. They were originally developed to model heat flow and have expanded into signal processing. They're all about decomposing a signal into a series of overlapping waves. So why have they popped up in probability theory? The answer is subtle, but it relates to this term \(e^{itX}\) in the definition of the characteristic function. We can expand it to \(e^{itX} = cos(tX) + i·sin(tX)\) which looks like waves - which in turn provides the links to Fourier analysis. This argument is a bit hand-wavy, but the it would take pages to get into the detail.
Let's turn to adding two random variables. To make it easier, I'll go back to the insurance example of a company offering car and house insurance. How do we add these two distributions to get overall losses?
The process of summing two random variables is called convolution. Here's how it's defined. The symbol for convolution is \(*\).
\[f_Z(z) = f_X(z) * f_Y(z) = (f_X * f_Y)(z) = \int_{-\infty}^{\infty} f_X(x)\, f_Y(z - x)\, dx\]
Once again, we're back in the world of Fourier Transforms. If we use the symbol \(\mathcal{F}\) to represent a Fourier Transform and \(\mathcal{F}^{-1}\) to represent an inverse Fourier Transform, then:
\[f_Z = f_X * f_Y = \mathcal{F}^{-1}\Big\{\mathcal{F}\{f_X\} \cdot \mathcal{F}\{f_Y\}\Big\}\]
What does going into the world of Fourier Transforms achieve? Two things:
- It changes the calculation from convolution into multiplication of two functions; a much easier operation.
- It means we have access to the huge amount of Fourier Transform literature and methods. When it comes to doing this on a computer, we can use the vast amount of work others have done to create Fourier Transform libraries.
Of course, we can sum more than two random variables. Here are the relevant equations:
\[\varphi_{X_i}(t) = \int_{-\infty}^{\infty} e^{itx}\, f_{X_i}(x)\, dx\]
\[\varphi_{S_n}(t) = \prod_{i=1}^{n} \varphi_{X_i}(t)\]
\[f_{S_n}(z) = \frac{1}{2\pi} \int_{-\infty}^{\infty} e^{-itz}\, \varphi_{S_n}(t)\, dt\]
or more succinctly:
\[f_{S_n} = \mathcal{F}^{-1}\left\{ \prod_{i=1}^{n} \mathcal{F}\{f_{X_i}\} \right\}\]
Here's what you do in words and not symbols:
- Get the probability distributions for each random variable.
- Calculate the Fourier Transform for each random variable.
- Multiply the Fourier Transforms.
- Do an inverse Fourier Transform on the result.
Problems in the real world
This sounds easy, but it isn't. Although you can find the Fourier Transform for many probability distributions (search for characteristic function), multiplying them can be tricky. Even worse, transforming the product using an inverse Fourier Transform can be very, very hard.
I'm going to take one example, let's imagine we want to sum two random variables, both with a beta distribution. The Fourier Transforms of a beta function looks like this:
\[\varphi_X(t) = {}_1F_1(\alpha;\ \alpha+\beta;\ it)\]
where:
\[{}_1F_1(a; b; z) = \sum_{k=0}^{\infty} \frac{a^{(k)}}{b^{(k)}} \cdot \frac{z^k}{k!}\]
Now, imagine multiplying two of these beasts and transforming them with an inverse Fourier Transform!
In fact, for many of these types of calculations, no closed-form solution exists.
How you solve this problem is a little complex, but first, let's take a detour and talk about some easier calculations.
Some short cuts exist
Normal distributions are easy. The convolution of a normal with another normal is a normal! The mean and variance of the summed normal distribution is:
\[\mu_{S_n} = E[S_n] = \sum_{i=1}^{n} \mu_i = \sum_{i=1}^{n} E[X_i]\]
\[\sigma_n = \sum_{i=1}^{n} \sigma_i^2\]
Uniform distributions are a bit harder. There are three cases:
- All the distributions are exactly the same (i.e. iid or independent identically distributed) and are Uniform(0,1), then the sum follow an Irwin-Hall distribution.
\[f_{S_n}(z) = \frac{1}{(n-1)!} \sum_{k=0}^{\lfloor z \rfloor} (-1)^k \binom{n}{k} (z-k)^{n-1}, \qquad 0 \le z \le n\] - All the random variables follow the same distribution, but it isn't Unform(0,1), instead is it's Uniform(a, b).
\[ f_{S_n}(z) = \frac{1}{(b-a)(n-1)!} \sum_{k=0}^{\left\lfloor \frac{z-na}{b-a} \right\rfloor} (-1)^k \binom{n}{k} \left(\frac{z-na}{b-a} - k\right)^{n-1}, \qquad na \le z \le nb\] - The random variables follow different Uniform distributions. In this case, there is no convenient closed-form expression for the sum.
Here are some of the other shortcuts.
\begin{array}{|l|l|l|} \hline \textbf{Distribution} & \textbf{Condition} & \textbf{Sum is} \\ \hline \text{Normal} & X_i \sim N(\mu_i, \sigma_i^2) & N\left(\sum \mu_i,\ \sum \sigma_i^2\right) \\ \hline \text{Poisson} & X_i \sim \text{Pois}(\lambda_i) & \text{Pois}\left(\sum \lambda_i\right) \\ \hline \text{Binomial} & X_i \sim \text{Bin}(n_i, p),\ \text{same } p & \text{Bin}\left(\sum n_i,\ p\right) \\ \hline \text{Exponential (same rate)} & X_i \sim \text{Exp}(\lambda),\ \text{i.i.d.} & \text{Gamma}(n, \lambda) \\ \hline \text{Gamma (same rate)} & X_i \sim \text{Gamma}(k_i, \lambda) & \text{Gamma}\left(\sum k_i,\ \lambda\right) \\ \hline \text{Chi-squared} & X_i \sim \chi^2_{k_i} & \chi^2_{\sum k_i} \\ \hline \text{Negative Binomial (same } p\text{)} & X_i \sim \text{NB}(r_i, p) & \text{NB}\left(\sum r_i,\ p\right) \\ \hline \text{Cauchy} & X_i \sim \text{Cauchy}(x_{0,i}, \gamma_i) & \text{Cauchy}\left(\sum x_{0,i},\ \sum \gamma_i\right) \\ \hline \end{array}
What happens if you're trying to sum random variables with different distributions or your sum isn't one of the ones listed here? You have two choices: slog through the math yourself or step into the digital world.
Summing random variables in practice
Often, there's no closed form expression for the sum, so you have to do thing a little differently. Let's assume you're trying to sum several random variables with very different probability distributions. What can you do?
The first step is to go as far as you can with the standard formula I gave above. For example, if you have multiple normally distributed random variables, sum them.
In many cases, you can figure out the characteristic function for the distributions; these are the Fourier Transforms. You can now write out the Fourier Transform of your sum as a series of characteristic function multiplications. Now, evaluate that expression on a grid. For example, if every characteristic function is a function of x, choose x values from a minimum to a maximum with a step size. You may need to use hundreds of thousands or even millions of x values. Finally, use a digital inverse FFT (Fast Fourier Transform) to transform back and get the distribution of the sum.
This process sometimes needs a great amount of computing power. If the power you need is greater than what's available there is another approach. This involves Monte Carlo simulation, or better still, Quasi Monte Carlo simulation. I'm not going to go into that here, it's way too involved.
The overall process of summing random variables is obviously more involved than what I've described, but you get the gist.
Zero-inflation is a problem
So far, I've only handled the case where the probabilities are 1, i.e., there will be car insurance claims or house claims. What happens if that isn't the case? What happens if there's a probability of no losses whatsoever? This is a more complicated situation I'll deal with in another blog post.
Summary
To sum two or more random variables, we use convolution, and Fourier Transforms are very much part of that process. In some cases, there are nice closed-form expressions for the sum, but most of the time, there are not. In practice, we need to digitize the data and step into the world of the discrete Fourier Transform.
AI statement
All the writing was don by me (a human). I used AI to generate the formulae and the charts.