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

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!

• Assuming a continuous CDF, the approach itself is exact assuming you are able to simulate draws with infinite precision. But you're working with a computer that's limited by discretization. Your question should be more about the limits of accuracy to implementing things on computers, not about the inverse CDF approach. Sep 14, 2018 at 21:51
• @stats_model Thanks. Is there a mathematical prove existing? Sep 14, 2018 at 21:56
• You're going to need to state more precisely what you want to prove. The proof that the CDF is exactly what it should be is just mechanical: $F(t) = P(X \leq t) = P(U \leq F(t)) = P(F^{-1}(U) \leq t)$. But I see you want to know something about the pdf. As a general rule, the CDF is more fundamental to the distribution. If you are able to implement the CDF exactly, then you have implemented the whole distribution because $P(a \leq x \leq b) = P(x \leq b) - P(x \leq a) = F(b) - F(a)$. The PDF is just something that exists when your CDF happens to be sufficiently regular. Sep 14, 2018 at 22:02

## 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 $$Pr(F^{-1}(u) \leq x) = Pr(u \leq F(x))$$ where we have inverted the CDF to the other side, now we get $$Pr(F^{-1}(u) \leq x)= \int\limits_{-\infty}^{\infty}f_U(u) \ du = \int\limits_{0}^{F(x)} \ du$$ This is so because $F(x) \geq 0$. Now, we get $$Pr(F^{-1}(u) \leq x)= F(x)$$ We have just proved that $$Pr(F^{-1}(u) \leq x) = F(x)$$ 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

• Thanks. Other than the empirical experiments, is there a mathematical prove existing? Sep 14, 2018 at 21:57
• welcome. do you want me to prove it for the exponential case ? Sep 14, 2018 at 21:57
• Yes, that would be great! Thank you! Sep 14, 2018 at 22:01
• I think I should provide a proof for all PDF's .. I'm editing Sep 14, 2018 at 22:01
• Thanks for sharing the proof. I am studying it, but don't understand the 4th equation: Pr(F−1(u)≤x)=∫−∞∞fU(u) du Could you please explain why? Thank you! Sep 17, 2018 at 17:02

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.