On Mon, May 10, 2010 at 8:37 PM,

I went googling and found a new interpretation

numpy.random.pareto is actually the Lomax distribution also known as Pareto 2, Pareto (II) or Pareto Second Kind distribution

Great!

So, from this it looks like numpy.random does not have a Pareto distribution, only Lomax, and the confusion came maybe because somewhere in the history the (II) (second kind) got dropped in the explanations.

and actually it is in scipy.stats.distributions, but without rvs

# LOMAX (Pareto of the second kind.) # Special case of Pareto of the first kind (location=-1.0)

I understand the point with this last comment, but I think it can be confusing in that the Pareto (of the first kind) has no "location" parameter and people might think you are referring to the Generalized Pareto distribution. I think its much clearer to say: # Special case of the Pareto of the first kind, but shifted to the left by 1. x --> x + 1

2) Modify numpy/random/mtrand/distributions.c in the following way:

double rk_pareto(rk_state *state, double a) { //return exp(rk_standard_exponential(state)/a) - 1; return 1.0 / rk_double(state)**(1.0 / a); }

I'm not an expert on random number generator, but using the uniform distribution as in http://en.wikipedia.org/wiki/Pareto_distribution#Generating_a_random_sample_... and your Devroy reference seems better, than based on the relationship to the exponential distribution http://en.wikipedia.org/wiki/Pareto_distribution#Relation_to_the_exponential...

Correct. The exp relationship was for the existing implementation (which corresponds to the Lomax). I commented that line out and just used 1/U^(1/a).

I think without changing the source we can rewrite the docstring that this is Lomax (or Pareto of the Second Kind), so that at least the documentation is less misleading.

But I find calling it Pareto very confusing, and I'm not the only one anymore, (and I don't know if anyone has used it assuming it is classical Pareto), so my preferred solution would be

* rename numpy.random.pareto to numpy.random.lomax * and create a real (classical, first kind) pareto distribution (even though it's just adding or subtracting 1, ones we know it)

I personally have used numpy.random.pareto thinking it was the Pareto distribution of the first kind---which led to this post in the first place. So, I'm in strong agreement. While doing this, perhaps we should increase functionality and allow users the ability to specify the scale of the distribution (instead of just the shape)? I can make a ticket for this and give a stab at creating the necessary patch.

What's the backwards compatibility policy with very confusing names in numpy?

It seems reasonable that we might have to follow the deprecation route, but I'd be happier with a "faster" fix. 1.5 - Provide numpy.random.lomax. Make numpy.random.pareto raise a DeprecationWarning and then call lomax. 2.0 (if there is no 1.6) - Make numpy.random.pareto behave as Pareto distribution of 1st kind. Immediately though, we can modify the docstring that is currently in there to make the situation clear, instructing users how they can generate samples from the "standard" Pareto distribution. This is the first patch I'll submit. Perhaps it is better to only change the docstring and then save all changes in functionality for 2.0. Deferring to others on this one...