Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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.

share|improve this question
add comment

1 Answer 1

up vote 1 down vote accepted

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).

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.