BUG: stats.yulesimon: incorrect kurtosis values

[tk-yoshimura]

Describe your issue.

The results do not match wolfram's results for kurtosis of the yule simon distribution.

I think the following is correct:
_discrete_distns.py L1660

        g2 = np.where(alpha > 4,
              alpha + 3 + ((11 * alpha**3 - 49 * alpha - 22) /
                           (alpha * (alpha - 4) * (alpha - 3))),

input for wolfram alpha :
WaringYuleDistribution[5]

Reproducing Code Example

import scipy.stats as stats

dist = stats.yulesimon(5)

print(dist.stats("mvsk"))

Error message

result: 
-6.199999999999999

expected:
118.8

SciPy/NumPy/Python version and system information

1.13.0 1.23.0 sys.version_info(major=3, minor=9, micro=5, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: /c/opt/64/include
    lib directory: /c/opt/64/lib
    name: openblas
    openblas configuration: USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= ZEN MAX_THREADS=24
    pc file directory: c:/opt/64/lib/pkgconfig
    version: 0.3.26.dev
  lapack:
    detection method: pkgconfig
    found: true
    include directory: /c/opt/64/include
    lib directory: /c/opt/64/lib
    name: openblas
    openblas configuration: USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= ZEN MAX_THREADS=24
    pc file directory: c:/opt/64/lib/pkgconfig
    version: 0.3.26.dev
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 2.12.0
Compilers:
  c:
    commands: cc
    linker: ld.bfd
    name: gcc
    version: 10.3.0
  c++:
    commands: c++
    linker: ld.bfd
    name: gcc
    version: 10.3.0
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.10
  fortran:
    commands: gfortran
    linker: ld.bfd
    name: gcc
    version: 10.3.0
  pythran:
    include directory: C:\Users\runneradmin\AppData\Local\Temp\pip-build-env-z8a5b3_f\overlay\Lib\site-packages/pythran
    version: 0.15.0
Machine Information:
  build:
    cpu: x86_64
    endian: little
    family: x86_64
    system: windows
  cross-compiled: false
  host:
    cpu: x86_64
    endian: little
    family: x86_64
    system: windows
Python Information:
  path: C:\Users\runneradmin\AppData\Local\Temp\cibw-run-0mr_3znw\cp39-win_amd64\build\venv\Scripts\python.exe
  version: '3.9'

[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.

[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.