# Simulating Random Variables with Inverse Transform Sampling¶

In studying the transformation of random variables in All of Statistics and working on a few related exercises I've been in search of bigger picture motivation and a very cool result finally clicked: that we can simulate samples from any distribution by applying its inverse CDF to samples taken from a uniform random variable. (Note: technically this only works when the CDF has a closed form inverse function). Or in plainer terms, how can we transform results from math.random() to simulate draws from another distribution?

## Transformation of Random Variables¶

Let's consider how to take the transformation of a random variable $X$ with cumulative distribution function $F_X(x)$. Let $Y = r(X)$, that is, $Y$ is the transformation of $X$ via function $r(X)$.

To get the CDF of $Y$ we use the definition of CDFs:

$F_Y(y) = P(Y \leq y) = P(r(X) \leq y)$

We have $F_X(x)$ and want to know how to compute $F_Y(y)$ in terms of $F_X(x)$. To get there we can take the inverse of $r(x)$ on both sides of the inequality:

$F_Y(y) = P(Y \leq y) = P(r(X) \leq y) = P(X \leq r^{-1}(y))$

hey, that kind of looks like the CDF of $X$:

$P(X \leq r^{-1}(y)) = F_X(r^{-1}(y))$

and that's how we get $F_Y(y)$ in terms of $F_X(x)$. We can compute the density function $f_y(y)$ by differentiating $F_Y(y)$, applying the chain rule:

$f_y(y) = f_y(r^{-1}(y)) \times \frac{d}{dy} r^{-1}(y)\mathop{dy}$

Note that it is only this simple if $r(x)$ is one-to-one and strictly monotone increasing; it gets more complicated to reason about the regions where $Y$ is defined otherwise, but this hopefully provides the background for why we use with the inverse of the CDF.

## Transforming a Uniform Random Variable¶

Getting back to the "cool result" I mentioned earlier, let's put this to use to show how we can take samples from other distributions if all we have available is a random number generator. Random number generation is a whole other topic, but you are likely familiar with a random function available in most languages that will give you a number between 0 and 1. This can be viewed as a function that will give you samples from a uniform random variable in the range 0 to 1.

So we have a uniform random variable $U$ and another distribution $X$ with CDF $F$ that we'd like to draw samples from. If we'd like to be able to get $X$ as a function of $F$ and $U$, we need to work backwards and use the inverse of $F$ as the transformation function of $U$, so that when we apply the inverse of the inverse we get back to $X$:

$F_X(x) = P(F^{-1}(U) \leq x) = P(U \leq F(x)) = F_U(F_X(x))$

$F_U(F(x)) = \begin{cases} 0 & F(x) < 0 \\\ F(x) & 0 \leq F(x) \leq 1 \\\ 1 & F(x) > 1 \end{cases}$

Note that $F(x)$, being a CDF, can only take ranges from 0 to 1, so this collapses to just:

$F_X(x) = F_U(F(x)) = F(x)$ wherever $F(x)$ is defined, so $X \sim F$.

Now that we've established that $X = F^{-1}(U)$, let's put it to use in code.

## Sampling from Uniform¶

First, let's plot a histogram of samples from the $U(0, 1)$:

In [1]:
%matplotlib inline
%config InlineBackend.figure_format = 'retina'

import random
import matplotlib.pyplot as plt

plt.hist([random.random() for i in range(10000)], 50, normed=True)
plt.title('Samples from Uniform(0,1)')
None


yep, that looks pretty uniform.

## Sampling from an Exponential Distribution¶

Now let's say we want to sample from an exponential distribution to simulate the interval between requests hitting a web server (or see plenty of other examples of where this function might come in handy).

The exponential distribution has CDF:

$F_X(x) = 1 - e^{-\lambda x}$

following the plan outlined above, we can create our own exponential sampling function based on the random number generator if we can compute the CDF's inverse:

$1 - e^{-\lambda x} = y$
$e^{-\lambda x} = 1 - y$
$-\lambda x = ln(1 - y)$
$x = \frac{-1}{\lambda}ln(1 - y)$

so

$F^{-1}_X(x) = \frac{-1}{\lambda}ln(1 - x)$

Applying this inverse function to a sample taken from our random number generator should do the trick. Let's write this function and plot its histogram.

In [2]:
import math

def my_exp(lmbda=1.0):
return (-1 / lmbda)*math.log(1 - random.random())

plt.hist([my_exp() for i in range(10000)], 50, normed=True)
plt.title('Samples from hand rolled exponential function')
None


and just to make sure this looks right, let's use numpy's exponential function and compare:

In [3]:
import numpy.random

plt.hist([numpy.random.exponential() for i in range(10000)], 50, normed=True)
plt.title('Samples from numpy.random.exponential')
None


Looking good.

## Why does this work?¶

To get some intuition as to why this works, picture a cumulative distribution function and imagine taking a random point on the Y axis and seeing where it intersects the CDF.

Here's a plot of where 0.7 intersects the CDF for the Exponential distribution:

In [4]:
import numpy as np
from scipy.stats import expon

x = np.arange(0, 5, 0.005)
plt.plot(x, expon.cdf(x))
plt.plot(x, [0.7 for el in x], 'r--')
plt.title('CDF of Exponential')
None


this is exactly what it means to take the inverse; where the dotted line intersects is $F^{-1}(0.7)$.

So with a random number generator in hand, we can randomly start on the Y axis and use a random variable's inverse CDF function to get back to a point sampled according to its distribution.