Finding the density of \(g(X)\)

If \(Y = g(X)\) where \(X\) is a random variable with a given density and \(g\) is invertible, then we can find the density of \(Y\) in two steps:

  1. Find the CDF of Y,
  2. Find the density of Y by taking the derivative of the CDF.

Example

If \(R \sim \mathrm{Unif}(0,2)\), and \(Y = \pi R^2\), then the CDF of \(Y\) is

\[\begin{align*} F_Y(x) &= P(Y \leq x) \\ &= P(\pi R^2 \leq x) \\ &= P(R \leq \sqrt{x/\pi}) \\ &= F_R(\sqrt{x/\pi}), \end{align*}\]

where \(F_R\) denotes the CDF of \(R\). We’ve computed the CDF of a \(\mathrm{Unif}(a,b)\) random variable to be \(F(x)=(x-a)/(b-a)\). This means \(F_R(x) = x/2\). Therefore \[F_Y(x) = F_R\left(\sqrt{x/\pi}\right) = \frac 12 \sqrt{x/\pi}.\]

The density of \(Y\) is then \[f_Y(x) = F_Y'(x) = \frac{d}{dx}\left(\frac 12 \sqrt{x/\pi}\right) = \frac{1}{2\sqrt{\pi}}\frac{d}{dx}\left(\sqrt{x}\right) = \frac{1}{4\sqrt{\pi x}}.\]

Simulation verification

We can check numerically whether our density is correct by generating lots of samples of \(Y\) and then plotting a histogram of our data:

n <- 1e6
r <- runif(n, 0, 2)
y <- pi*r^2
hist(y, probability = T, breaks = 50)
curve(1/(4*sqrt(pi*x)), from = 0, to = 4*pi, add = T, col = "red") # density for Y we computed

The above command overlays the predicted density curve on the histogram to show that our calculation was correct.

The inverse transform method

There is a dual question to the previous example: given a distribution, how do we generate data that fits that distribution? More concretely, how do we simulate a random variable \(Y\) whose CDF is \(F\), if we’re only allowed to use the built in runif command to generate uniformly random numbers? That is, can we find \(g\) so that \(Y = g(U)\) has CDF \(F\), where \(U\sim\mathrm{Unif}(0,1)\)?

Remark. This is a truly salient problem. R has built in tools to generate random numbers for the classical distributions, but in practice it’s common to build a model for which the desired distribution has no such built in command.

The following steps, known as the inverse transform method, give a method for this, assuming \(F\) is invertible:

  1. Generate \(U \sim \mathrm{Unif}(0,1)\),
  2. Let \(Y = F^{-1}(U)\).

Proof

Here is why the CDF of \(Y = F^{-1}(U)\) is \(F\): \[F_Y(x) = P(Y \leq x) = P(F^{-1}(U) \leq x) = P(U \leq F(x)) = F_U(F(x)) = F(x),\] where the last equality holds because \(F_U(x) = x\).

Note. The inverse transform method depends on \(F\) being invertible. There is an active field of probability theory that deals with developing algorithms for general CDFs (ie. without assuming invertibility) whose recent history goes back to the 1940s and the advent of modern numerical analysis and computation. Names you might see in this area include rejection sampling, Markov Chain Monte Carlo, and Gibbs sampling.

Example

Suppose we want to generate samples from the \(\mathrm{Exp}(\lambda)\) distribution. The CDF of \(\mathrm{Exp}(\lambda)\) is \[F(x) = 1-e^{-\lambda x}.\] Using standard algebra, we find that \[F^{-1}(x) = -\frac 1\lambda\ln(1-x).\] Thus, we could let \[Y = -\frac 1\lambda\ln(1-U)\] where \(U \sim \mathrm{Unif}(0,1)\), and the theory tells us that \(Y \sim \mathrm{Exp}(\lambda)\).

Simulation verification

We can check numerically whether the above method works. We’ll generate samples like we did in our last simulation, and then see if our data fits the desired density.

lambda <- 1
n <- 1e6
u <- runif(n, 0, 1)
y <- -1/lambda * log(1-u)
hist(y, probability = T, breaks = 50)
curve(lambda*exp(-lambda*x), from = 0, to = 12, add = T, col = "red") # density of exponential distribution