# Generate random numbers from a family of PDFs

For a part of a simulation task, I need to generate (lots of) random numbers from the distribution

$$P(E_k | N) = \frac{1}{E_k}\left(\frac{E_k}{k_BT}\right)^{N-1}\frac{1}{(N-2)!} e^{-E_k/k_BT}$$

For $N>2$ and constant $k_BT$. (Depending on the case, N is usually Poisson distributed about N=5 or N=10, but it shouldn't break on larger N).

My current method is to use a trivial rejection algorithm for $N<10$ and a Normal approximation for $N>=10$. The problem with this technique is that it

1. truncates the distribution
2. for the $N<10$ cases, uses an average of 29 uniform random numbers per successful $E_k$ (I could trade this off for issue (1)).
3. for the $N>10$ cases, Normal is still a somewhat poor approximation

Is there some way of making a generalized deterministic transformation from a uniformly distributed RV to this distribution? I would rather not have to make a set of tables to use with Ziggurat or something, but the current state of this generator needs work.

If you would like a current state estimate, in a simple 4-minute test run, just this generation alone (65M cases) ate 16% of my execution time, of which 7.2% was generating 870M uniform randoms, 1.8% was 35M normal randoms, and the rest was internal computation: primarily the testing of the rejection part. Even a modest improvement in execution speed and accuracy would be quite helpful.

-

I can see the probability distribution for $E_k$ given $N$ is of a Gamma random variate with scale parameter $k_BT$ and shape parameter $N-1$. You can use any of the standard Statistical packages to generate $E_k$ for a given $N>2$ and a constant $k_BT$. For example, function rgamma generates Gamma random variates for given Gamma parameters in R that uses modified rejection algorithm developed by Ahrens and Dieter (1982).