Random Variate Generation by Numerical Inversion when only the Density Is Known

We present a numerical inversion method for generating random variates
from continuous distributions when only the density function is
given. The algorithm is based on polynomial interpolation of the
inverse CDF and Gauss-Lobatto integration. The user can select the
required precision, which may be close to machine precision for
smooth, bounded densities; the necessary tables have moderate
size. Our computational experiments with the classical standard
distributions (normal, beta, gamma, t-distributions) and with the
noncentral chi-square, hyperbolic, generalized hyperbolic, and stable
distributions showed that our algorithm always reaches the required
precision. The setup time is moderate and the marginal execution time
is very fast and nearly the same for all distributions. Thus for the
case that large samples with fixed parameters are required the
proposed algorithm is the fastest inversion method known. Speed-up
factors up to 1000 are obtained when compared to inversion algorithms
developed for the specific distributions. This makes our algorithm
especially attractive for the simulation of copulas and for
quasi--Monte Carlo applications.