From 5066d54db5c6339f1eec95e404f29712e476d22e Mon Sep 17 00:00:00 2001 From: phofl Date: Sun, 11 Oct 2020 15:13:59 +0200 Subject: [PATCH 1/9] Implement bedfords algorithm --- pandas/_libs/window/aggregations.pyx | 31 +++++++++++++++++----------- 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/pandas/_libs/window/aggregations.pyx b/pandas/_libs/window/aggregations.pyx index 937f7d8df7728..56c2346f5c6c5 100644 --- a/pandas/_libs/window/aggregations.pyx +++ b/pandas/_libs/window/aggregations.pyx @@ -302,7 +302,7 @@ cdef inline float64_t calc_var(int64_t minp, int ddof, float64_t nobs, result = ssqdm_x / (nobs - ddof) # Fix for numerical imprecision. # Can be result < 0 once Kahan Summation is implemented - if result < 1e-15: + if result < 0: #1e-15: result = 0 else: result = NaN @@ -311,10 +311,10 @@ cdef inline float64_t calc_var(int64_t minp, int ddof, float64_t nobs, cdef inline void add_var(float64_t val, float64_t *nobs, float64_t *mean_x, - float64_t *ssqdm_x) nogil: + float64_t *ssqdm_x, float64_t *m_x) nogil: """ add a value from the var calc """ cdef: - float64_t delta + float64_t delta, prev_mean # `isnan` instead of equality as fix for GH-21813, msvc 2017 bug if isnan(val): @@ -324,15 +324,18 @@ cdef inline void add_var(float64_t val, float64_t *nobs, float64_t *mean_x, # a part of Welford's method for the online variance-calculation # https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance delta = val - mean_x[0] + prev_mean = mean_x[0] mean_x[0] = mean_x[0] + delta / nobs[0] - ssqdm_x[0] = ssqdm_x[0] + ((nobs[0] - 1) * delta ** 2) / nobs[0] + m_x[0] = m_x[0] + (val - prev_mean) * (val - mean_x[0]) + ssqdm_x[0] = m_x[0] + # ssqdm_x[0] = ssqdm_x[0] + ((nobs[0] - 1) * delta ** 2) / nobs[0] cdef inline void remove_var(float64_t val, float64_t *nobs, float64_t *mean_x, - float64_t *ssqdm_x) nogil: + float64_t *ssqdm_x, float64_t *m_x) nogil: """ remove a value from the var calc """ cdef: - float64_t delta + float64_t delta, prev_mean if notnan(val): nobs[0] = nobs[0] - 1 @@ -340,11 +343,15 @@ cdef inline void remove_var(float64_t val, float64_t *nobs, float64_t *mean_x, # a part of Welford's method for the online variance-calculation # https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance delta = val - mean_x[0] + prev_mean = mean_x[0] mean_x[0] = mean_x[0] - delta / nobs[0] - ssqdm_x[0] = ssqdm_x[0] - ((nobs[0] + 1) * delta ** 2) / nobs[0] + m_x[0] = m_x[0] - (val - prev_mean) * (val - mean_x[0]) + # ssqdm_x[0] = ssqdm_x[0] - ((nobs[0] + 1) * delta ** 2) / nobs[0] + ssqdm_x[0] = m_x[0] else: mean_x[0] = 0 ssqdm_x[0] = 0 + m_x[0] = 0 def roll_var(ndarray[float64_t] values, ndarray[int64_t] start, @@ -353,7 +360,7 @@ def roll_var(ndarray[float64_t] values, ndarray[int64_t] start, Numerically stable implementation using Welford's method. """ cdef: - float64_t mean_x = 0, ssqdm_x = 0, nobs = 0, + float64_t mean_x = 0, ssqdm_x = 0, nobs = 0, m_x = 0, float64_t val, prev, delta, mean_x_old int64_t s, e Py_ssize_t i, j, N = len(values) @@ -375,7 +382,7 @@ def roll_var(ndarray[float64_t] values, ndarray[int64_t] start, if i == 0 or not is_monotonic_bounds: for j in range(s, e): - add_var(values[j], &nobs, &mean_x, &ssqdm_x) + add_var(values[j], &nobs, &mean_x, &ssqdm_x, &m_x) else: @@ -384,17 +391,17 @@ def roll_var(ndarray[float64_t] values, ndarray[int64_t] start, # calculate adds for j in range(end[i - 1], e): - add_var(values[j], &nobs, &mean_x, &ssqdm_x) + add_var(values[j], &nobs, &mean_x, &ssqdm_x, &m_x) # calculate deletes for j in range(start[i - 1], s): - remove_var(values[j], &nobs, &mean_x, &ssqdm_x) + remove_var(values[j], &nobs, &mean_x, &ssqdm_x, &m_x) output[i] = calc_var(minp, ddof, nobs, ssqdm_x) if not is_monotonic_bounds: for j in range(s, e): - remove_var(values[j], &nobs, &mean_x, &ssqdm_x) + remove_var(values[j], &nobs, &mean_x, &ssqdm_x, &m_x) return output From 9f83a4c0f5cdf53f944988cf0b5f4910fd092533 Mon Sep 17 00:00:00 2001 From: phofl Date: Sun, 11 Oct 2020 17:39:13 +0200 Subject: [PATCH 2/9] Use kahan summation for var --- pandas/_libs/window/aggregations.pyx | 43 ++++++++++++++-------------- 1 file changed, 22 insertions(+), 21 deletions(-) diff --git a/pandas/_libs/window/aggregations.pyx b/pandas/_libs/window/aggregations.pyx index 56c2346f5c6c5..1aa6c79aedb73 100644 --- a/pandas/_libs/window/aggregations.pyx +++ b/pandas/_libs/window/aggregations.pyx @@ -311,10 +311,10 @@ cdef inline float64_t calc_var(int64_t minp, int ddof, float64_t nobs, cdef inline void add_var(float64_t val, float64_t *nobs, float64_t *mean_x, - float64_t *ssqdm_x, float64_t *m_x) nogil: + float64_t *ssqdm_x, float64_t *compensation) nogil: """ add a value from the var calc """ cdef: - float64_t delta, prev_mean + float64_t delta, prev_mean, y, t # `isnan` instead of equality as fix for GH-21813, msvc 2017 bug if isnan(val): @@ -323,35 +323,36 @@ cdef inline void add_var(float64_t val, float64_t *nobs, float64_t *mean_x, nobs[0] = nobs[0] + 1 # a part of Welford's method for the online variance-calculation # https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance - delta = val - mean_x[0] + y = val - compensation[0] + t = y - mean_x[0] + compensation[0] = t + mean_x[0] - y + delta = t prev_mean = mean_x[0] mean_x[0] = mean_x[0] + delta / nobs[0] - m_x[0] = m_x[0] + (val - prev_mean) * (val - mean_x[0]) - ssqdm_x[0] = m_x[0] - # ssqdm_x[0] = ssqdm_x[0] + ((nobs[0] - 1) * delta ** 2) / nobs[0] + ssqdm_x[0] = ssqdm_x[0] + (val - prev_mean) * (val - mean_x[0]) cdef inline void remove_var(float64_t val, float64_t *nobs, float64_t *mean_x, - float64_t *ssqdm_x, float64_t *m_x) nogil: + float64_t *ssqdm_x, float64_t *compensation) nogil: """ remove a value from the var calc """ cdef: - float64_t delta, prev_mean + float64_t delta, prev_mean, y, t if notnan(val): nobs[0] = nobs[0] - 1 if nobs[0]: # a part of Welford's method for the online variance-calculation # https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance - delta = val - mean_x[0] + y = val - compensation[0] + t = y - mean_x[0] + compensation[0] = t + mean_x[0] - y + delta = t prev_mean = mean_x[0] mean_x[0] = mean_x[0] - delta / nobs[0] - m_x[0] = m_x[0] - (val - prev_mean) * (val - mean_x[0]) - # ssqdm_x[0] = ssqdm_x[0] - ((nobs[0] + 1) * delta ** 2) / nobs[0] - ssqdm_x[0] = m_x[0] + ssqdm_x[0] = ssqdm_x[0] - (val - prev_mean) * (val - mean_x[0]) else: mean_x[0] = 0 ssqdm_x[0] = 0 - m_x[0] = 0 def roll_var(ndarray[float64_t] values, ndarray[int64_t] start, @@ -360,7 +361,7 @@ def roll_var(ndarray[float64_t] values, ndarray[int64_t] start, Numerically stable implementation using Welford's method. """ cdef: - float64_t mean_x = 0, ssqdm_x = 0, nobs = 0, m_x = 0, + float64_t mean_x = 0, ssqdm_x = 0, nobs = 0, compensation_add = 0, compensation_remove = 0, float64_t val, prev, delta, mean_x_old int64_t s, e Py_ssize_t i, j, N = len(values) @@ -382,26 +383,26 @@ def roll_var(ndarray[float64_t] values, ndarray[int64_t] start, if i == 0 or not is_monotonic_bounds: for j in range(s, e): - add_var(values[j], &nobs, &mean_x, &ssqdm_x, &m_x) + add_var(values[j], &nobs, &mean_x, &ssqdm_x, &compensation_add) else: # After the first window, observations can both be added # and removed - # calculate adds - for j in range(end[i - 1], e): - add_var(values[j], &nobs, &mean_x, &ssqdm_x, &m_x) - # calculate deletes for j in range(start[i - 1], s): - remove_var(values[j], &nobs, &mean_x, &ssqdm_x, &m_x) + remove_var(values[j], &nobs, &mean_x, &ssqdm_x, &compensation_remove) + + # calculate adds + for j in range(end[i - 1], e): + add_var(values[j], &nobs, &mean_x, &ssqdm_x, &compensation_add) output[i] = calc_var(minp, ddof, nobs, ssqdm_x) if not is_monotonic_bounds: for j in range(s, e): - remove_var(values[j], &nobs, &mean_x, &ssqdm_x, &m_x) + remove_var(values[j], &nobs, &mean_x, &ssqdm_x, &compensation_remove) return output From 214d5637f27efc20306c820e793a29a719fe955a Mon Sep 17 00:00:00 2001 From: phofl Date: Sun, 11 Oct 2020 17:45:15 +0200 Subject: [PATCH 3/9] Add whatsnew and adjust comments to reflect new behavior --- doc/source/whatsnew/v1.2.0.rst | 1 + pandas/_libs/window/aggregations.pyx | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/doc/source/whatsnew/v1.2.0.rst b/doc/source/whatsnew/v1.2.0.rst index 2b4b10c39602a..8fe75a3a4cdc0 100644 --- a/doc/source/whatsnew/v1.2.0.rst +++ b/doc/source/whatsnew/v1.2.0.rst @@ -191,6 +191,7 @@ Other enhancements - :meth:`DatetimeIndex.searchsorted`, :meth:`TimedeltaIndex.searchsorted`, :meth:`PeriodIndex.searchsorted`, and :meth:`Series.searchsorted` with datetimelike dtypes will now try to cast string arguments (listlike and scalar) to the matching datetimelike type (:issue:`36346`) - Added methods :meth:`IntegerArray.prod`, :meth:`IntegerArray.min`, and :meth:`IntegerArray.max` (:issue:`33790`) - Where possible :meth:`RangeIndex.difference` and :meth:`RangeIndex.symmetric_difference` will return :class:`RangeIndex` instead of :class:`Int64Index` (:issue:`36564`) +- :meth:`Rolling.var()` and :meth:`Rolling.std()` use Kahan summation and Welfords Method to avoid numerical issues (:issue:`37051`) .. _whatsnew_120.api_breaking.python: diff --git a/pandas/_libs/window/aggregations.pyx b/pandas/_libs/window/aggregations.pyx index 1aa6c79aedb73..1a4f6839b5f78 100644 --- a/pandas/_libs/window/aggregations.pyx +++ b/pandas/_libs/window/aggregations.pyx @@ -321,7 +321,7 @@ cdef inline void add_var(float64_t val, float64_t *nobs, float64_t *mean_x, return nobs[0] = nobs[0] + 1 - # a part of Welford's method for the online variance-calculation + # Welford's method for the online variance-calculation using Kahan summation # https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance y = val - compensation[0] t = y - mean_x[0] @@ -341,7 +341,7 @@ cdef inline void remove_var(float64_t val, float64_t *nobs, float64_t *mean_x, if notnan(val): nobs[0] = nobs[0] - 1 if nobs[0]: - # a part of Welford's method for the online variance-calculation + # Welford's method for the online variance-calculation using Kahan summation # https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance y = val - compensation[0] t = y - mean_x[0] From 602443c9f21e7f9c701bd7c177028e4d7e8687a5 Mon Sep 17 00:00:00 2001 From: phofl Date: Sun, 11 Oct 2020 17:51:34 +0200 Subject: [PATCH 4/9] Add test --- pandas/tests/window/test_rolling.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/pandas/tests/window/test_rolling.py b/pandas/tests/window/test_rolling.py index 73831d518032d..16390ce397db3 100644 --- a/pandas/tests/window/test_rolling.py +++ b/pandas/tests/window/test_rolling.py @@ -868,3 +868,18 @@ def test_rolling_period_index(index, window, func, values): result = getattr(ds.rolling(window, closed="left"), func)() expected = pd.Series(values, index=index) tm.assert_series_equal(result, expected) + + +@pytest.mark.parametrize( + ("func", "values"), + [ + ("var", [5e33, 0, 0.5, 0.5, 2, 0]), + ("std", [7.071068e16, 0, 0.7071068, 0.7071068, 1.414214, 0]), + ], +) +def test_rolling_var_numerical_issues(func, values): + # GH: 37051 + ds = pd.Series([99999999999999999, 1, 1, 2, 3, 1, 1]) + result = getattr(ds.rolling(2), func)() + expected = pd.Series([np.nan] + values) + tm.assert_series_equal(result, expected) From a953f4669f4d9c3caffc286d24a854a5fc321ed8 Mon Sep 17 00:00:00 2001 From: phofl Date: Sun, 11 Oct 2020 17:54:50 +0200 Subject: [PATCH 5/9] Delete old comment --- pandas/_libs/window/aggregations.pyx | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/pandas/_libs/window/aggregations.pyx b/pandas/_libs/window/aggregations.pyx index 1a4f6839b5f78..f13801aec7f89 100644 --- a/pandas/_libs/window/aggregations.pyx +++ b/pandas/_libs/window/aggregations.pyx @@ -300,9 +300,7 @@ cdef inline float64_t calc_var(int64_t minp, int ddof, float64_t nobs, result = 0 else: result = ssqdm_x / (nobs - ddof) - # Fix for numerical imprecision. - # Can be result < 0 once Kahan Summation is implemented - if result < 0: #1e-15: + if result < 0: result = 0 else: result = NaN From 354693550099962e0b10bdf528772f80c260f829 Mon Sep 17 00:00:00 2001 From: phofl Date: Sun, 11 Oct 2020 18:05:16 +0200 Subject: [PATCH 6/9] Fix flake problems --- pandas/_libs/window/aggregations.pyx | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/pandas/_libs/window/aggregations.pyx b/pandas/_libs/window/aggregations.pyx index f13801aec7f89..a175520c930f1 100644 --- a/pandas/_libs/window/aggregations.pyx +++ b/pandas/_libs/window/aggregations.pyx @@ -319,7 +319,8 @@ cdef inline void add_var(float64_t val, float64_t *nobs, float64_t *mean_x, return nobs[0] = nobs[0] + 1 - # Welford's method for the online variance-calculation using Kahan summation + # Welford's method for the online variance-calculation + # using Kahan summation # https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance y = val - compensation[0] t = y - mean_x[0] @@ -339,7 +340,8 @@ cdef inline void remove_var(float64_t val, float64_t *nobs, float64_t *mean_x, if notnan(val): nobs[0] = nobs[0] - 1 if nobs[0]: - # Welford's method for the online variance-calculation using Kahan summation + # Welford's method for the online variance-calculation + # using Kahan summation # https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance y = val - compensation[0] t = y - mean_x[0] @@ -359,7 +361,8 @@ def roll_var(ndarray[float64_t] values, ndarray[int64_t] start, Numerically stable implementation using Welford's method. """ cdef: - float64_t mean_x = 0, ssqdm_x = 0, nobs = 0, compensation_add = 0, compensation_remove = 0, + float64_t mean_x = 0, ssqdm_x = 0, nobs = 0, compensation_add = 0, + float64_t compensation_remove = 0, float64_t val, prev, delta, mean_x_old int64_t s, e Py_ssize_t i, j, N = len(values) From d29e7c5ef424dc7dea978c0a81671de93abe0c8d Mon Sep 17 00:00:00 2001 From: phofl Date: Sun, 11 Oct 2020 19:03:37 +0200 Subject: [PATCH 7/9] Adjust kahan summation and add new tests --- pandas/_libs/window/aggregations.pyx | 11 ++++++----- pandas/tests/window/test_rolling.py | 12 +++++++----- 2 files changed, 13 insertions(+), 10 deletions(-) diff --git a/pandas/_libs/window/aggregations.pyx b/pandas/_libs/window/aggregations.pyx index a175520c930f1..72a6f654b1544 100644 --- a/pandas/_libs/window/aggregations.pyx +++ b/pandas/_libs/window/aggregations.pyx @@ -322,11 +322,11 @@ cdef inline void add_var(float64_t val, float64_t *nobs, float64_t *mean_x, # Welford's method for the online variance-calculation # using Kahan summation # https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance + prev_mean = mean_x[0] - compensation[0] y = val - compensation[0] t = y - mean_x[0] compensation[0] = t + mean_x[0] - y delta = t - prev_mean = mean_x[0] mean_x[0] = mean_x[0] + delta / nobs[0] ssqdm_x[0] = ssqdm_x[0] + (val - prev_mean) * (val - mean_x[0]) @@ -336,18 +336,17 @@ cdef inline void remove_var(float64_t val, float64_t *nobs, float64_t *mean_x, """ remove a value from the var calc """ cdef: float64_t delta, prev_mean, y, t - if notnan(val): nobs[0] = nobs[0] - 1 if nobs[0]: # Welford's method for the online variance-calculation # using Kahan summation # https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance + prev_mean = mean_x[0] - compensation[0] y = val - compensation[0] t = y - mean_x[0] compensation[0] = t + mean_x[0] - y delta = t - prev_mean = mean_x[0] mean_x[0] = mean_x[0] - delta / nobs[0] ssqdm_x[0] = ssqdm_x[0] - (val - prev_mean) * (val - mean_x[0]) else: @@ -393,7 +392,8 @@ def roll_var(ndarray[float64_t] values, ndarray[int64_t] start, # calculate deletes for j in range(start[i - 1], s): - remove_var(values[j], &nobs, &mean_x, &ssqdm_x, &compensation_remove) + remove_var(values[j], &nobs, &mean_x, &ssqdm_x, + &compensation_remove) # calculate adds for j in range(end[i - 1], e): @@ -403,7 +403,8 @@ def roll_var(ndarray[float64_t] values, ndarray[int64_t] start, if not is_monotonic_bounds: for j in range(s, e): - remove_var(values[j], &nobs, &mean_x, &ssqdm_x, &compensation_remove) + remove_var(values[j], &nobs, &mean_x, &ssqdm_x, + &compensation_remove) return output diff --git a/pandas/tests/window/test_rolling.py b/pandas/tests/window/test_rolling.py index 16390ce397db3..f792a362164a8 100644 --- a/pandas/tests/window/test_rolling.py +++ b/pandas/tests/window/test_rolling.py @@ -871,15 +871,17 @@ def test_rolling_period_index(index, window, func, values): @pytest.mark.parametrize( - ("func", "values"), + ("func", "third_value", "values"), [ - ("var", [5e33, 0, 0.5, 0.5, 2, 0]), - ("std", [7.071068e16, 0, 0.7071068, 0.7071068, 1.414214, 0]), + ("var", 1, [5e33, 0, 0.5, 0.5, 2, 0]), + ("std", 1, [7.071068e16, 0, 0.7071068, 0.7071068, 1.414214, 0]), + ("var", 2, [5e33, 0.5, 0, 0.5, 2, 0]), + ("std", 2, [7.071068e16, 0.7071068, 0, 0.7071068, 1.414214, 0]), ], ) -def test_rolling_var_numerical_issues(func, values): +def test_rolling_var_numerical_issues(func, third_value, values): # GH: 37051 - ds = pd.Series([99999999999999999, 1, 1, 2, 3, 1, 1]) + ds = pd.Series([99999999999999999, 1, third_value, 2, 3, 1, 1]) result = getattr(ds.rolling(2), func)() expected = pd.Series([np.nan] + values) tm.assert_series_equal(result, expected) From 411ba688264fb0bed808ec33395405504fb77b4c Mon Sep 17 00:00:00 2001 From: phofl Date: Sun, 11 Oct 2020 19:05:40 +0200 Subject: [PATCH 8/9] Run black --- pandas/tests/window/test_rolling.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pandas/tests/window/test_rolling.py b/pandas/tests/window/test_rolling.py index f792a362164a8..ac0adeaf308bd 100644 --- a/pandas/tests/window/test_rolling.py +++ b/pandas/tests/window/test_rolling.py @@ -876,7 +876,7 @@ def test_rolling_period_index(index, window, func, values): ("var", 1, [5e33, 0, 0.5, 0.5, 2, 0]), ("std", 1, [7.071068e16, 0, 0.7071068, 0.7071068, 1.414214, 0]), ("var", 2, [5e33, 0.5, 0, 0.5, 2, 0]), - ("std", 2, [7.071068e16, 0.7071068, 0, 0.7071068, 1.414214, 0]), + ("std", 2, [7.071068e16, 0.7071068, 0, 0.7071068, 1.414214, 0]), ], ) def test_rolling_var_numerical_issues(func, third_value, values): From 95d652c852b668f88203ea9fe86ffa8e072a9143 Mon Sep 17 00:00:00 2001 From: phofl Date: Sun, 11 Oct 2020 20:05:01 +0200 Subject: [PATCH 9/9] Add old comment --- pandas/_libs/window/aggregations.pyx | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pandas/_libs/window/aggregations.pyx b/pandas/_libs/window/aggregations.pyx index 72a6f654b1544..b05f887a269df 100644 --- a/pandas/_libs/window/aggregations.pyx +++ b/pandas/_libs/window/aggregations.pyx @@ -300,7 +300,9 @@ cdef inline float64_t calc_var(int64_t minp, int ddof, float64_t nobs, result = 0 else: result = ssqdm_x / (nobs - ddof) - if result < 0: + # Fix for numerical imprecision. + # Can be result < 0 once Kahan Summation is implemented + if result < 1e-15: result = 0 else: result = NaN