Continuous Distributions

Assume we want to generate a random variable \(X\) with cumulative distribution function (CDF) \(F_X\). The inverse transform sampling algorithm is simple:
1. Generate \(U \sim \text{Unif}(0,1)\)
2. Let \(X = F_X^{-1}(U)\).

Then, \(X\) will follow the distribution governed by the CDF \(F_X\), which was our desired result.

Note that this algorithm works in general but is not always practical. For example, inverting \(F_X\) is easy if \(X\) is an exponential random variable, but its harder if \(X\) is Normal random variable.

Discrete Distributions

Now we will consider the discrete version of the inverse transform method. Assume that \(X\) is a discrete random variable such that \(P(X = x_i) = p_i\). The algorithm proceeds as follows:
1. Generate \(U \sim \text{Unif}(0,1)\)
2. Determine the index \(k\) such that \(\sum_{j=1}^{k-1} p_j \leq U < \sum_{j=1}^k p_j\), and return \(X = x_k\).

Notice that the second step requires a search.

Why it works ?

Computationally, this method involves computing the quantile function of the distribution — in other words, computing the cumulative distribution function (CDF) of the distribution (which maps a number in the domain to a probability between 0 and 1) and then inverting that function many times. This is the source of the term “inverse” or “inversion” in most of the names for this method. Note that for a discrete distribution, computing the CDF is not in general too difficult: we simply add up the individual probabilities for the various points of the distribution. For a continuous distribution, however, we need to integrate the probability density function (PDF) of the distribution, which is impossible to do analytically for most distributions (including the normal distribution). As a result, this method may be computationally inefficient for many distributions that does not have intractable CDF and other methods are preferred; however, it is a useful method for building more generally applicable samplers such as those based on rejection sampling.

For the normal distribution, the lack of an analytical expression for the corresponding quantile function means that other methods (e.g. the Box–Muller transform) may be computationally favorable. It is often the case that, even for simple distributions, the inverse transform sampling method can be improved on.

Simple Implementation (Discrete Random Variable)

Assume our random variable \(X\) can take on any one of \(K\) values with probabilities \(\{p_1, \ldots, p_K\}\). We implement the algorithm below, assuming these probabilities are stored in a vector called prob

discrete_inverse_sampling <- function( prob ) {
  U  <- runif(1)
  if(U <= prob[1]){
    return(1)
  }
  for(state in 2:length(prob)) {
    if(sum(prob[1:(state-1)]) < U && U <= sum(prob[1:state]) ) {
      return(state)
    }
  }
}