@@ -1366,23 +1366,19 @@ def nankurt(
13661366 m2 = adjusted2 .sum (axis , dtype = np .float64 )
13671367 m4 = adjusted4 .sum (axis , dtype = np .float64 )
13681368
1369+ # #57972: tolerance to consider the central moment equals to zero.
1370+ # We adapted the tolerance from scipy:
1371+ # https://github.com/scipy/scipy/blob/04d6d9c460b1fed83f2919ecec3d743cfa2e8317/scipy/stats/_stats_py.py#L1429
1372+ constant_tolerance2 = (np .finfo (m2 .dtype ).eps * total ) ** 2 # match order of m2
1373+ constant_tolerance4 = constant_tolerance2 ** 2 # match order of m4
1374+ m2 = _zero_out_fperr (m2 , constant_tolerance2 )
1375+ m4 = _zero_out_fperr (m4 , constant_tolerance4 )
1376+
13691377 with np .errstate (invalid = "ignore" , divide = "ignore" ):
13701378 adj = 3 * (count - 1 ) ** 2 / ((count - 2 ) * (count - 3 ))
13711379 numerator = count * (count + 1 ) * (count - 1 ) * m4
13721380 denominator = (count - 2 ) * (count - 3 ) * m2 ** 2
13731381
1374- # floating point error
1375- #
1376- # #18044 in _libs/windows.pyx calc_kurt follow this behavior
1377- # to fix the fperr to treat denom <1e-14 as zero
1378- # #57972 arbitrary <1e-14 tolerance leads to problematic behaviour on low variance.
1379- # We adapted the tolerance to use one similar to scipy:
1380- # https://github.com/scipy/scipy/blob/04d6d9c460b1fed83f2919ecec3d743cfa2e8317/scipy/stats/_stats_py.py#L1429
1381- constant_tolerance2 = (np .finfo (m2 .dtype ).eps * total ) ** 2 # match order of m2
1382- constant_tolerance4 = constant_tolerance2 ** 2 # match order of m4
1383- numerator = _zero_out_fperr (numerator , constant_tolerance2 )
1384- denominator = _zero_out_fperr (denominator , constant_tolerance4 )
1385-
13861382 if not isinstance (denominator , np .ndarray ):
13871383 # if ``denom`` is a scalar, check these corner cases first before
13881384 # doing division
0 commit comments