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

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.

• Could you precise what you mean by "ad hoc methods"? Jun 3, 2020 at 15:00
• On 'ad hoc', relatedly for the Gamma distribution random deviate generation is, for example, the Ahrens-Dieter acceptance–rejection method as cited here en.wikipedia.org/wiki/Gamma_distribution . Jun 4, 2020 at 14:22
• Thanks: Acceptance-rejection is an exact method so I do not understand in which sense it is ad hoc, especially when optimised. Devroye (1986, section IX.3) discusses the optimal strategies for different ranges of the Gamma scale, with acceptance probabilities being lower bounded. Jun 4, 2020 at 15:45
• True an exact method, but I would not consider it as intuitive as, for example, the path to solution via the CDF inversion method (proceed to solve directly or by numerical methods). Jun 4, 2020 at 17:00
• The other side of the coin is that the inverse CDF may cost much more to compute than running an accept-reject with optimised parameters. Jun 4, 2020 at 17:21

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.

• I have compared, for the special case of the Erlang-distribution (so a generalization of results to other distributions may not be valid), the inversion root method with NRM, to the procedure based on the sum of k natural logarithms of uniform variables. For the case k=3, the process based on the sum of k natural logs of k uniforms, base on routine sample mean and variance comparisons, appears to be superior, however, I did not look at higher moments, or extreme values occurrence. Jun 4, 2020 at 16:53