BUG: stats.yulesimon: incorrect kurtosis values

[tk-yoshimura]

BUG: stats.yulesimon: incorrect kurtosis values

[tk-yoshimura]

I am not familiar with the PR rules of this repository, so may I ask you to do this on my behalf?
It is already midnight in Japan and I am going to bed.

Thank you so much.

[mdhaber]

Yes, according to Wikipdia,

our implementation is missing that factor of 11.

It's probably missed by our generic tests because the generic calculation is inaccurate.

To generate your reference in Mathematica, are you using Kurtosis or ExcessKurtosis? I'm a surprised to see an explicit +3 in the formula if it's really the excess kurtosis.

Would you like to submit a PR @tk-yoshimura?

[mdhaber]

Just out of curiousity, I checked whether scipy.integrate._tanhsinh._nsum would have provided a better generic implementation here. There are probably acceleration methods that would achieve great accuracy with fewer function evaluations, but the integral approximation is nice since it will work for pretty much any monotonically decreasing sequence.

import scipy.stats as stats
stats.yulesimon._munp(4, 5)  # 51.50194329138698 and warning that sum doesn't converge
stats.yulesimon.moment(4, 5)  # 52.083333333333336 (after this correction)

import numpy as np
from scipy.integrate._tanhsinh import _nsum
from scipy._lib._util import _lazywhere

def f(x):
    return _lazywhere(np.isfinite(x), (x,), 
                      lambda x: stats.yulesimon._pmf(x, 5)*x**4, 
                      fillvalue=0)
    
_nsum(f, 1, np.inf)
#  success: True
#   status: 0
#      sum: 52.083333333330906
#    error: 2.792518871722942e-07
#     nfev: 33181

(actual relative error 5e-14)

[tk-yoshimura]

I am glad my code correction was correct.
Many thanks.

If you are interested in the stable distribution of half-integer shape parameter numbers not yet implemented in scipy, I have a C# implementation in my repository below.
https://github.com/tk-yoshimura/DoubleDoubleStatistic

[mdhaber]

Thanks @tk-yoshimura. I'll pass the link to the stable distribution to the right people. In the meantime, would you like to submit a PR to fix the kurtosis implementation, or would you like me to take care of it? If you do it, I can review and get it merged. If I do it, I'll need to wait for another maintainer to review and merge.