cummean is more stable against FP errors. closes tidyverse#1387

NEWS.md

@@ -13,7 +13,9 @@
 * `lead` and `lag` correctly handle default values for string columns in
   hybrid (#1403).
 
-* `bind_rows` handles `POSIXct` stored as integer (#1402).   
+* `bind_rows` handles `POSIXct` stored as integer (#1402).
+
+* `cummean` is more stable against floating point errors (#1387).
 
 # dplyr 0.4.3
 

src/window.cpp

@@ -19,10 +19,10 @@ LogicalVector cumall(LogicalVector x) {
   if( current == NA_LOGICAL) return out ;
   if( current == FALSE){
     std::fill( out.begin(), out.end(), FALSE ) ;
-    return out ;    
+    return out ;
   }
   for (int i = 1; i < n; i++) {
-    current = x[i] ; 
+    current = x[i] ;
     if( current == NA_LOGICAL ) break ;
     if( current == FALSE ){
       std::fill( out.begin() + i, out.end(), FALSE ) ;
@@ -44,18 +44,18 @@ LogicalVector cumany(LogicalVector x) {
   if( current == NA_LOGICAL ) return out ;
   if( current == TRUE ){
     std::fill( out.begin(), out.end(), TRUE ) ;
-    return out ;      
+    return out ;
   }
   for (int i = 1; i < n; i++) {
     current = x[i] ;
     if( current == NA_LOGICAL ) break ;
     if( current == TRUE ){
       std::fill( out.begin() + i, out.end(), TRUE ) ;
-      break ;        
+      break ;
     }
     out[i] = current || out[i - 1];
   }
-  
+
   return out;
 }
 
@@ -64,11 +64,12 @@ LogicalVector cumany(LogicalVector x) {
 // [[Rcpp::export]]
 NumericVector cummean(NumericVector x) {
   int n = x.length();
-  NumericVector out(n);
+  NumericVector out = no_init(n);
 
-  out[0] = x[0];
-  for (int i = 1; i < n; i++) {
-    out[i] = out[i - 1] * ((i * 1.0) / (i + 1)) + x[i] / (i + 1);
+  double sum = out[0] = x[0];
+  for (int i = 1; i < n; i++ ) {
+    sum += x[i];
+    out[i] = sum / (i + 1.0);
   }
 
   return out;

tests/testthat/test-window.R

@@ -7,23 +7,23 @@ test_that("If n = 0, lead and lag return x", {
 
 test_that("If n = length(x), returns all missing", {
   miss <- rep(NA_integer_, 2)
-  
+
   expect_equal(lead(1:2, 2), miss)
   expect_equal(lag(1:2, 2), miss)
-  
+
 })
 
 test_that("cumany handles NA (#408)", {
   batman <- c(NA,NA,NA,NA,NA)
   expect_true(all(is.na(cumany(batman))))
   expect_true(all(is.na(cumall(batman))))
-  
+
   x <- c(FALSE,NA)
   expect_true( all( !cumall(x) ) )
-  
+
   x <- c(TRUE,NA)
   expect_true( all( cumany(x) ) )
-  
+
 })
 
 test_that("percent_rank ignores NAs (#1132)", {
@@ -34,3 +34,7 @@ test_that("cume_dist ignores NAs (#1132)", {
   expect_equal( cume_dist(c(1:3, NA)), c(1/3, 2/3, 1, NA) )
 })
 
+test_that( "cummean is not confused by FP error (#1387)", {
+  a <- rep(99, 9)
+  expect_true( all( cummean(a) == a) )
+})