How accurate is it to generate random number using the inverse-cdf-method?

Question

I am reading the following approach of generating a random number from a given probability distribution.

https://blogs.sas.com/content/iml/2013/07/22/the-inverse-cdf-method.html

I understand that cdf is 1-1 mapping, so there is a guarantee that the generated number will be within the domain of the pdf. However, I am wondering what's the accuracy of this approach? Is there a mathematical prove?

I applied this method on the exponential distribution, and found the the average of the generated numbers are smaller than expected ... I am wondering is this because I just didn't draw enough samples or this method is not very accurate? Thanks!

Answer 1

Proof

In what follows, $F(x)$ is a CDF and $f(x)$ is a PDF.

Here is a proof for your case, if $$U \sim U([0,1])$$ Then let's prove that (assuming that the CDF of $X$ is bijective) $$X = F^{-1}(U)$$ is distributed according to $f_X(x)$. This is what you're looking for So \begin{equation} Pr(F^{-1}(u) \leq x) = Pr(u \leq F(x)) \end{equation} where we have inverted the CDF to the other side, now we get \begin{equation} Pr(F^{-1}(u) \leq x)= \int\limits_{-\infty}^{\infty}f_U(u) \ du = \int\limits_{0}^{F(x)} \ du \end{equation} This is so because $F(x) \geq 0$. Now, we get \begin{equation} Pr(F^{-1}(u) \leq x)= F(x) \end{equation} We have just proved that \begin{equation} Pr(F^{-1}(u) \leq x) = F(x) \end{equation} which is the CDF of $X$.

Empirical results

The approach is accurate and this is what I got on MATLAB

This is the MATLAB code

M = 10000 % number of samples to generate using MATLAB's exprnd

l = 3; % lambda parameter

mu = 0; % mean

f =  @(x) (-1/l)*log(1-x) % using your approach

u =rand(M,1); % generate M uniform random samples between 0 and 1

z = f(u); % z are exponentially distributed with your approach
% just plotting
histogram(exprnd(1/l,M,1),'Normalization','probability')
hold on
histogram(z,'Normalization','probability')
legend('MATLAB Exponential Random Number Generator','Simulated using your approach')
title('Exp. PDF with \lambda = 3')
grid on
grid minor

Answer 2

About a decade ago I was involved in providing a reference implementation for the random distributions/generators package in teh C++11 standard. The issue of how to create random distributions that are "faithful" to the last machine epsilon (or a small number times that) is subtle. Your method using CDF mapping is very general but suffers from the following headache:

The inverse CDF needs to be represented as a function, and except in certain special cases, that function will not have a closed form in terms of existing other functions (like sines and exponentials) available in the standard. So you are forced to approximate it, normally using a table and interpolation method (rescaled cubic splines is pretty close to optimal for this purpose). Even though the approximation gives 4-th order accuracy, in general for table sizes of a few K bytes (which are acceptable memory costs) you get faithful representations at the float level ($2^{-32}$) but not at the double level ($2^{-53}$).

To address this issue, you need to use a rejection method (these are explained well in Numerical Recipes). Rejection methods utilize the supplied cdf function, which you are entitled to presume is accurate. Combining cubic splines with rejection costs at least one extra uniform random number generation all the time, but if done with care the average number of extra randoms needed is very close to just one.