BUG: more floating point error robustness in rolling mean. close #2527

Author: wesm

RELEASE.rst

@@ -412,7 +412,8 @@ pandas 0.9.1
   - Fix variety of cut/qcut string-bin formatting bugs (GH1978_, GH1979_)
   - Raise Exception when xs view not possible of MultiIndex'd DataFrame (GH2117_)
   - Fix groupby(...).first() issue with datetime64 (GH2133_)
-  - Better floating point error robustness in some rolling_* functions (GH2114_)
+  - Better floating point error robustness in some rolling_* functions
+    (GH2114_, GH2527_)
   - Fix ewma NA handling in the middle of Series (GH2128_)
   - Fix numerical precision issues in diff with integer data (GH2087_)
   - Fix bug in MultiIndex.__getitem__ with NA values (GH2008_)
@@ -483,6 +484,7 @@ pandas 0.9.1
 .. _GH2117: https://github.com/pydata/pandas/issues/2117
 .. _GH2133: https://github.com/pydata/pandas/issues/2133
 .. _GH2114: https://github.com/pydata/pandas/issues/2114
+.. _GH2527: https://github.com/pydata/pandas/issues/2114
 .. _GH2128: https://github.com/pydata/pandas/issues/2128
 .. _GH2008: https://github.com/pydata/pandas/issues/2008
 .. _GH2179: https://github.com/pydata/pandas/issues/2179

pandas/algos.pyx

@@ -43,6 +43,7 @@ cdef inline int int_min(int a, int b): return a if a <= b else b
 cdef extern from "math.h":
     double sqrt(double x)
     double fabs(double)
+    int signbit(double)
 
 from . import lib
 
@@ -948,12 +949,12 @@ def roll_sum(ndarray[double_t] input, int win, int minp):
 
 def roll_mean(ndarray[double_t] input,
                int win, int minp):
-    cdef double val, prev, sum_x = 0
-    cdef Py_ssize_t nobs = 0, i
-    cdef Py_ssize_t N = len(input)
+    cdef:
+        double val, prev, result, sum_x = 0
+        Py_ssize_t nobs = 0, i, neg_ct = 0
+        Py_ssize_t N = len(input)
 
     cdef ndarray[double_t] output = np.empty(N, dtype=float)
-
     minp = _check_minp(win, minp, N)
 
     for i from 0 <= i < minp - 1:
@@ -963,6 +964,8 @@ def roll_mean(ndarray[double_t] input,
         if val == val:
             nobs += 1
             sum_x += val
+            if signbit(val):
+                neg_ct += 1
 
         output[i] = NaN
 
@@ -972,15 +975,27 @@ def roll_mean(ndarray[double_t] input,
         if val == val:
             nobs += 1
             sum_x += val
+            if signbit(val):
+                neg_ct += 1
 
         if i > win - 1:
             prev = input[i - win]
             if prev == prev:
                 sum_x -= prev
                 nobs -= 1
+                if signbit(prev):
+                    neg_ct -= 1
 
         if nobs >= minp:
-            output[i] = sum_x / nobs
+            result = sum_x / nobs
+            if neg_ct == 0 and result < 0:
+                # all positive
+                output[i] = 0
+            elif neg_ct == nobs and result > 0:
+                # all negative
+                output[i] = 0
+            else:
+                output[i] = result
         else:
             output[i] = NaN
 

pandas/stats/moments.py

@@ -278,7 +278,10 @@ def _rolling_moment(arg, window, func, minp, axis=0, freq=None,
     calc = lambda x: func(x, window, minp=minp, **kwargs)
     return_hook, values = _process_data_structure(arg)
     # actually calculate the moment. Faster way to do this?
-    result = np.apply_along_axis(calc, axis, values)
+    if values.ndim > 1:
+        result = np.apply_along_axis(calc, axis, values)
+    else:
+        result = calc(values)
 
     rs = return_hook(result)
     if center:

pandas/stats/tests/test_moments.py

@@ -306,6 +306,14 @@ def test_fperr_robustness(self):
         result = mom.rolling_var(arr, 2)
         self.assertTrue((result[1:] >= 0).all())
 
+        # #2527, ugh
+        arr = np.array([0.00012456, 0.0003, 0])
+        result = mom.rolling_mean(arr, 1)
+        self.assertTrue(result[-1] >= 0)
+
+        result = mom.rolling_mean(-arr, 1)
+        self.assertTrue(result[-1] <= 0)
+
     def _check_moment_func(self, func, static_comp, window=50,
                            has_min_periods=True,
                            has_center=True,

vb_suite/stat_ops.py

@@ -73,3 +73,12 @@
 
 stats_rank2d_axis0_average = Benchmark('df.rank()', setup,
                                        start_date=datetime(2011, 12, 12))
+
+# rolling functions
+
+setup = common_setup + """
+arr = np.random.randn(100000)
+"""
+
+stats_rolling_mean = Benchmark('rolling_mean(arr, 100)', setup,
+                               start_date=datetime(2011, 6, 1))