DOC: error in the description of `numpy.random.pareto`

[ViktorBezb]

Issue with current documentation:

The 'Notes' part of the documentation on Pareto distribution is wrong. The density in both pages reads as $p(x) = \frac{am^a}{x^{a+1}}$, but this is obviously not even a density, because the function $p(x)$ is not integrable for $a>0$.

https://numpy.org/doc/stable/reference/random/generated/numpy.random.Generator.pareto.html#numpy.random.Generator.pareto

https://numpy.org/doc/stable/reference/random/generated/numpy.random.pareto.html

Idea or request for content:

No response

[Mukulikaa]

Hey @Sinusxdx, I looked it up on Wikipedia and the docs seem correct - Pareto distribution.

[ziyaj]

Hey @Sinusxdx,

The probability density $p(x$) function is correct. It is integrable. When $x < m$, $p(x)=0$, so we only need to integral from $m$ to $+\infty$.

$\int p(x)dx = \int_{m}^{+\infty} \frac{am^a}{x^{a+1}} dx = a m^a \int_{m}^{+\infty} \frac{1}{x^{a+1}} dx = a m^a \left. (-\frac{1}{a} x^{-a})\right|_{m}^{+\infty} = m^a (m^{-a} - +\infty^{-a}) = m^0= 1$.

Note that $a$ and $m$ are constants, so they can be pulled out of the integral in the beginning. So we are pretty much left with $\int \frac{1}{x^{a+1}} dx$, which is one of the most basic integral. Note that $+\infty^{-a} = \frac{1}{+\infty^a} = 0$.

[ViktorBezb]

Hello @Mukulikaa and @ziyaj,

Indeed $p(x)$ is a density if considered only for $x > m$; however this is not reflected in the documentation. It should clearly specify that $p(x) = 0$ for $x < m$. Furthermore, the value of $m$ is not given anywhere: neither in the description, nor it is available as an argument. What is the value of $m$ used in the generation of a random sample?

Running this command

import numpy as np
min(np.random.pareto(7, size=10000))

repeatedly I get small numbers like 2.766544932031323e-06. Now, Pareto random variables cannot take values below $m$. Since I am getting very close to $0$, it seems that the distribution is supported on $[0, +\infty)$. However $p(x) = \frac{a m ^a}{x^{a+1}}$ is not integrable on $[0, +\infty)$ whatever $m > 0$ is. Something is wrong here.

[Mukulikaa]

Perhaps @melissawm can help out with this!

[melissawm]

It feels like something like this could be enough:

However, if we do have some code path on forbidden intervals maybe it's worth including a warning?

[ViktorBezb]

It feels like something like this could be enough:

However, if we do have some code path on forbidden intervals maybe it's worth including a warning?

I don't think it is sufficient. Two issues remain:

1. It is not indicated anywhere what numerical value of $m$ is used.
2. The distribution seems to be supported on $[0, +\infty)$, which should imply $m = 0$. But $m = 0$ is impossible.

I tried to find the code for numpy.random.pareto, and it seems that the source written in another language. I found this file. Is this how Pareto random variables generated? If so, this is the code:

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

Hence a Pareto random variable with parameter $a$ is generated as $e^{\xi/a}-1$, where $\xi$ is a standard unit exponential. From here we get the density $f(x) = \frac{a}{(x+1)^{a+1}}$, $x \geq 0$. Note an important difference with the current notes: $x+1$ instead of $x$ in the denominator, hence everything is integrable on $[0, +\infty)$.

Suggestion to modify the notes. Given the parameter $a > 0$, the probability density for the generated (pseudo-)random number is

$$ f(x) = \frac{a}{(x+1)^{a+1}}, \ \ \ x \geq 0,$$

This corresponds to the Pareto II distribution and can be seen as Pareto I distribution shape $a$ and scale $m=1$ shifted by $-1$.

Disclaimer: I am a mathematician and I do not have much experience with writing libraries.

[melissawm]

Maybe @rkern can take a look - as far as I can see the definition you mention now lives in this file

[rkern]

I implemented the distribution based on the algorithm in Luc Devroye's Non-Uniform Random Variate Generation on pg 262, which does indeed match @ViktorBezb's suggested text. Most of the Notes in the docstrings were written later by other people at a time when we had fewer tools for reviews. Given the plethora of minor variants of "Pareto distributions" that are out there, those volunteers just found the wrong variant to document by mistake.

[Matt711]

I'll work on this issue