BUG: fix polluted window in skewness computation

[Alvaro-Kothe]

• closes BUG: Rolling skew error #47461 (Replace xxxx with the GitHub issue number)
• xref BUG: df.rolling.{std, skew, kurt} gives unexpected value #61416
• xref QST: Consistently apply Welford Method and Kahan Summation in roll_xxx functions #59715
• Tests added and passed if fixing a bug or adding a new feature
• All code checks passed.
• Added type annotations to new arguments/methods/functions.
• Added an entry in the latest doc/source/whatsnew/vX.X.X.rst file if fixing a bug or adding a new feature.

---

This PR changes the algorithm used to compute the skewness. Previously, the algorithm used computed the summation of $X^{3}$, $X^{2}$, and $X$, and used kahan summation for stability. The problem is that $EX^3$ - $(EX)^3$ can suffer heavilly from catastrophic cancellation.

Now, it uses a Welford-like algorithm to compute up to the 3rd central moment, and continuously check for stability when updating the 3rd central moment, by checking for severe loss of significant digits when updating its value. If it detects a stability problem, it recomputes the window.

Benchmark

asv continuous -f 1.1 -E virtualenv:3.13 upstream/main HEAD -b rolling.Methods

Change	Before [ea75dd7]	After [488e99e] <fix/polluted-window-skew>	Ratio	Benchmark (Parameter)
-	1.41±0.01ms	1.21±0.01ms	0.86	rolling.Methods.time_method('DataFrame', ('expanding', {}), 'float', 'skew')
-	1.45±0.01ms	1.24±0.01ms	0.86	rolling.Methods.time_method('DataFrame', ('expanding', {}), 'int', 'skew')
-	1.41±0ms	1.22±0.03ms	0.86	rolling.Methods.time_method('Series', ('expanding', {}), 'int', 'skew')
-	1.39±0.01ms	1.18±0.08ms	0.85	rolling.Methods.time_method('Series', ('expanding', {}), 'float', 'skew')

Note

The performance increase is mainly due to the fact that it didn't have to recompute the window. This change can have a worse case scenario of $O(n * \text{window\_size})$

[Alvaro-Kothe]

Reproduction of #47461

import pandas as pd data = [ -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, 2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, 7729964027, 7729965641, 7729965103, 7729966179, 7729968330, 7729969406, 7729971020, 7729969406, 7729969944, 7729969406, 7729969406, 7729969944, 7729971558, 7729972633, 7729973171, 7729973171, 7729973171, 7729969406, ] df = pd.DataFrame(data, columns=["data"]) print(df.data.rolling(10).skew()[-20:]) # 34 NaN # 35 NaN # 36 3.162278e+00 # 37 1.778781e+00 # 38 1.035098e+00 # 39 4.841229e-01 # 40 2.426355e-13 # 41 -4.841229e-01 # 42 -1.035098e+00 # 43 -1.778781e+00 # 44 -3.162278e+00 # 45 -3.928221e-01 # 46 -6.794734e-01 # 47 -1.272001e+00 # 48 -1.034738e+00 # 49 8.839369e-01 # 50 9.294772e-01 # 51 4.561840e-01 # 52 1.688976e-01 # 53 1.688976e-01 # Name: data, dtype: float64 print(df.data.tail(10).skew()) # 0.16889762763904356

[eicchen]

Could I have you review my PR for #61416? The work is done so should save some repeated effort.

PR: #61481

[Alvaro-Kothe]

The window from the example #62863 (comment) was recomputed 3 times.

1. at index 23, removing 2. The 3rd moment dropped from 46.08 to -7e-15
2. at index 40, adding 7729968330. The 3rd moment dropped from 2e29 to 1e17.
3. at index 45, removing -2. The 3rd moment increased from -3e29 to 4e14.

[mroeschke]

Thanks you for using descriptive names for these variables!

[mroeschke]

How did you come up with the 1e10 threshold?

[Alvaro-Kothe]

It loses signifcant digits when $|A| + |B| >> |A - B|$, so I thought that a change in 10 significant digits in base 10 should be considered a catastrophic cancellation.

[mroeschke]

Might be interesting to use fused add multiply fma where applicable

[Alvaro-Kothe]

I couldn't figure out an effective way to reformulate the computation of term1 to use fma. But I managed to change the update functions for the 3rd moment. Now I am using fma(term1, n + 2.0, -3.0 * m2[0])

[Alvaro-Kothe]

I feel like this PR is changing a lot of the current behavior in the skew computation. I will minimize the changes to focus on numerical stability.