Extending Random Deviate Generation Using CDF Inversion Method with Iterative Newton-Raphson, Improved Results?

Question

As a test example for my assertion, I applied a Monte Carlo inversion method to generate random deviates for the case of the Erlang-distribution, whose cumulative distribution can be expressed as:

$F(x) = 1- \sum\limits^{k-1}_{n = 0} \frac{1}{n!} exp(-λx)(λx)^n$

Or, for example, in the case k=3, one can write:

$(1- F(x))exp(λx) = 1 + λx + 1/2(λx)^2$

Following the standard process for the Monte Carlo inversion method to simulate an Erlang random deviate, one can set 1-F(x) = U, a generated Uniform random deviate. Further, taking the natural log of both side (log transform is apparently required for convergence):

$ Ln(U) + λx = ln(1 + λx + 1/2(λx)^2)$

Implementing a numerical analysis Newton-Raphson method (NRM) to derive a solution for 'x' yielding a Erlang distribution random deviate, I define a function g(x), which should approach zero with an exact solution value for ‘x’, as:

$g(x) = ln(U) + λx - ln(1 + λx + 1/2(λx)^2) -> 0$

One also needs the derivative of g(x) with respect to x for the implementing the NRM:

$g’(x) = λ - ( λ + λ^2 x)/(1 + λx + 1/2(λx)^2)$

With the NRM, the next iterative value is given by:

$ x_{i+1} = x_{i} - g(x_i)/g’(x_i) $

where I employed the random uniform deviate (U) as the initial starting guess value for 'x' and repeated six cycles where it was evident that $g(x_{i+1})$ did, in fact, approach zero.

Per a worksheet implementation, I witnessed successful convergence of g(x) to 0.0000xxxx or better for all 100 generated Erlang random deviates in test runs. In particular, when k=3 and λ=1.5, observed, for example, a sample mean of 2.45 (expected 2.00) and a variance of 1.14 (expected 1.333).

When I performed a related analysis for k=2, I generally observed more precision in estimates of mean and variance, and even a greater accuracy for the expected ratio of the sample mean to variance (which should be λ).

So, for other seemingly non-invertible CDF, this numerical analysis based methodology may be of assistance. Further, I suspect this method may actually perform better than other ad hoc methods in both efficiency and accuracy in generating random deviates.

Support for my supposition, per the cited Wikipedia source on the Erlang distribution, on generating Erlang-distributed random variates, the approximate suggested methodology is the sum of minus the natural log of k uniform random deviates, scaled by the inverse of λ, which first employs multiple uniform random deviates and clearly is but an approximation.

Answer 1

The Erlang(k) distribution is exactly the distribution of the sum of $k$ Exponential(1) random variables, which is exactly the distribution of the sum of $k$ natural logarithms of uniform variables. So for random number generation it's not an approximation and any difference between the empirical CDF and the Erlang CDF is sampling variation that should be there. (You can still improve on the sum-of-$k$-log-uniforms approach for speed when $k$ is large, and the case of non-integer $k$ doesn't have a simple exact algorithm).

Your approach of taking a starting value and doing Newton-Raphson to improve on it is useful (and sometimes used) for evaluating the quantile function as part of the inversion method in situations where there isn't a straightforward exact method and where there's no closed form for the quantile function. For example, R uses inversion (by default) to generate Normal random variates, and evaluates the quantile function using a starting approximation (a rational polynomial) followed by Newton-Raphson. See in qnorm.c

*  DESCRIPTION
 *
 *  Compute the quantile function for the normal distribution.
 *
 *  For small to moderate probabilities, algorithm referenced
 *  below is used to obtain an initial approximation which is
 *  polished with a final Newton step.
 *
 *  For very large arguments, an algorithm of Wichura is used.

For Erlang-type distributions non-integer $k$ or large $k$ (and for general Gamma distributions) there are other exact (up to rounding) approaches for random-number generation in this context. For example, R uses rejection sampling (an algorithm by Ahrens & Dieter), which appears to be faster than inversion.