Bremsstrahlung

I complained earlier about trying to generate random values with a Planck distribution. This morning I had the good fortune to locate a most informative book which has helped me solve my problem.

The book in question is Non-Uniform Random Variate Generation. It turns out that the author, Luc Devroye, is apparently most dissatisfied with his publisher, and so is distributing the book for free. I would encourage anyone with an interesting the topic to take a look; it's quite comprehensive, and does a nice job of covering both the mathematical underpinnings and the practical implementation aspects.

For my problem, it turns out that the Planck distribution can be handily constructed from a gamma distribution and a Zipf distribution. For a Planck distribution given by:

A random variable fitting this distribution is given by:

Where G is a gamma distribution and Z a Zipf distribution. All that now remains is to implement these two more basic distributions.

In general, a gamma distribution requires a good deal of work, but turns out that one can get a gamma distribution with integer shape parameter (and the scale parameter implicitly taken to be one) simply by adding up exponential distributions with scale parameter one:

where E is the aforementioned exponential random variable. Naturally this formula is best used for small values of a, but since I only need a Planck distribution with a=3, I only need a gamma distribution with a=4, so no trouble there.

The Zipf distribution is more esoteric, and I admittedly don't understand it all that thoroughly, so I'll just quote Devroye's algorithm. For a Zipf distribution:

use the rejection algorithm:

b = 2^(a+1)
do
    x = U^(-1/(a-1))
    t = (1+(1/x))^(a-1)
while V*x*(t-1)/(b-1) <= (t/b)
x is the result

where ^ denotes exponentiation, and U and V are uniform random variates in [0,1).

Finally, here's an example of the result from the Planck distribution¹ with a=3, b=0.5: