BUG: stats.zipf: incorrect pmf values about scipy HOT 4 OPEN

Thanks for the report. Could you elaborate with a minimal working example for which parameter values the error occurs? I have a hard time interpreting the Wolfram Alpha computarion. The same holds for the other issue you opened.

from scipy.

Thanks for adding the MWE, that helps a lot.

from scipy.

For example, if the following code is executed:

import numpy as np
import scipy.stats as stats

x = np.arange(0, 16).astype(np.int32)
dist = stats.zipf(9.0)
pmf = dist.pmf(x)

print(pmf)

The values are as expected.

[0.00000000e+00 9.97995633e-01 1.94921022e-03 5.07034310e-05
3.80705121e-06 5.10973764e-07 9.90301387e-08 2.47312621e-08
7.43564690e-09 2.57600117e-09 9.97995633e-10 4.23247571e-10
1.93418240e-10 9.41105843e-11 4.83032464e-11 2.59601567e-11]

However, if the shape parameter is an integer, an incorrect value is obtained for x>=11, as shown below.

import numpy as np
import scipy.stats as stats

x = np.arange(0, 16).astype(np.int32)
dist = stats.zipf(9)
pmf = dist.pmf(x)

print(pmf)

[0.00000000e+00 9.97995633e-01 1.94921022e-03 5.07034310e-05
3.80705121e-06 5.10973764e-07 9.90301387e-08 2.47312621e-08
7.43564690e-09 2.57600117e-09 9.97995633e-10 0.00000000e+00
1.15400158e-09 4.95390192e-10 0.00000000e+00 0.00000000e+00]

When executing the following lines
_discrete_distns.py L1309

def _pmf(self, k, a):
    # zipf.pmf(k, a) = 1/(zeta(a) * k**a)
    Pk = 1.0 / special.zeta(a, 1) / k**a
    return Pk

The input value k and the shape parameter a are both integer types, which may be causing the overflow.

from scipy.

That is indeed the issue. The maximum value for int32 is $2147483647$ but $11^9=2357947691$. We could cast k to float, then do *k**-a instead. Would you be interested in drafting a PR?

from scipy.