# Fast generation of Pareto-distributed randoms.

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?

-
Maybe you should specify what you mean by $f(x)=1/(x+1)$. Often in probability theory one uses capital $F$ for a cumulative probability distribution function, so that if the random variable is (capital) $X$, then $\Pr(X\le x)=F(x)$, and a lower-case $f$ for the probability density function. For the Pareto distribution, one has $\Pr(X>x)=(x/c)^{-\alpha}$ for $x\ge c$, so $F(x)=1-(x/c)^{-\alpha}$ for $x\ge c$, and then $f(x) = \alpha(x/c)^{-\alpha-1}$ for $x\ge c$ and $f(x)=0$ for $x<c$. – Michael Hardy May 20 '13 at 4:03
Yeah, I could be clearer. Let me try this: I have code that can generate an exponential variate X much faster than the usual X = -log(U) (where U is a uniform variate). If we want exps of rate lambda, we just multiply these by 1/lambda, and we're still faster. The standard way to generate Paretos is X = 1 / (U ^ (1 / alpha)), which is slow. You remind me below that I could instead do X = exp(Y / alpha), where Y is one of my fast exponentials. That might be a little faster. But I know I could get rid of the exp() for any one value of alpha; my question is how useful that would be. – Lee Daniel Crocker May 20 '13 at 5:36

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.
Not sure how that applies to my question. If I understand correctly (again, codemonkey, not mathematician), you're saying I can generate Pareto-distributed numbers from the exponentially-distributed ones with another exp() operation, but that rather defeats the point: I could get Pareto-distributed numbers simply but expensively anyway. – Lee Daniel Crocker May 20 '13 at 3:37
So you're saying "$\exp$" is considered expensive and you're looking for a way to avoid it? – Michael Hardy May 20 '13 at 4:04
@LeeDanielCrocker : Why use rejection sampling to generate exponentially distributed random variables? Let $U$ be uniformly distributed in $[0,1]$; then $-\log_e U$ is exponentially distributed. $\qquad$ – Michael Hardy Jun 23 at 16:17