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.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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?

share|cite|improve this question
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.

share|cite|improve this answer
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
The Ziggurat algorithm generates normals and exponentials with a rejection-sampling algorithm that over 99% of the time generates one uniform random, compares that to a value in a table, and then multiplies it by a value in another table--and that's it. The other 1% of the time it does the math. I'd like to apply this to Paretos as well as exps and norms. If the single exp() is a lot faster than generating a Pareto from a uniform from scratch, then that might be a useful way to do it. But I'm really interested in whether trying to do it with the Ziggurat is better. – Lee Daniel Crocker May 20 '13 at 4:13
To give you some real numbers: doing Paretos the standard way (1/(u**(1/lambda))) takes about 22 seconds for 100,000,000 variates (slow laptop, no optimization). Exponentials from the Ziggurat take about 9 seconds (10 with a lambda). Generating Paretos with exp(expvariate()) is about 20 seconds. So maybe the best bang for the buck is to optimize my exponentials and generate Paretos from them. I'd sure like to apply the Ziggurat to something else, though, just for the exercise. Thanks a lot for your help. – Lee Daniel Crocker May 20 '13 at 8:11

Your Answer


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.