I'm developing a library of routines for generating random numbers for simulations (it's on GitHub). I've included fast routines for normally distributed and exponentially distributed randoms, using Marsaglia's nice fast Ziggurat algorithm (with some improvements). I surmise from things I've read that Pareto-distributed randoms would also be useful, and I don't see anywhere on the net that does this with the Ziggurat, so I'd like to do that. I'm more of a codemonkey than a mathematician; I managed to calculate the constants, build the tables and make the code for the exponentials and normals, but I have questions about the "alpha" parameter of the Pareto distribution that are a little beyond my mathematical chops.
The Ziggurat code tables are based on a particular curve; the cases I have done areexp(-x) and exp(-0.5 * x * x), with no other parameters. If you want numbers from a curve with different parameters, this can be done afterward with reasonably fast multiply or add operations, still keeping the expensive exp() outside the loop. But I'm not sure about Pareto's alpha. If I build a Ziggurat based on, say, alpha = 1 (which would make f(x) = 1 / (x + 1) I believe; is that right?), could the user then obtain values distributed with a different alpha with fast operations, or would he have to do calculations that would make the Ziggurat a waste of effort? Or, alternatively, could Pareto-distributed numbers be calculated cheaply from the other distributions, making the Ziggurat unnecessary? Thirdly, perhaps there is a value of alpha that is so commonly used that a Ziggurat for that value alone would be worth the effort, even if one had to calculate the others expensively?
Best Answer
Suppose $U$ is an exponentially distributed random variable and $\Pr(U>u)=e^{-au}$ for $u>0$ (and $\Pr(U>u)=1$ for $u\le 0$). Then $$ \Pr(e^U > x) = \Pr(U > \log_e x) = e^{-a\log_e x} = x^{-a} \text{ for }x>1, $$ so $e^U$ is a Pareto-distributed random variable. If you want the cutoff to be a positive number other than $1$, then just multiply $e^U$ by whatever that number is.