Skip to content
Prev 44182 / 63421 Next

Fastest non-overlapping binning mean function out there?

On 10/02/2012 05:19 PM, Henrik Bengtsson wrote:
since x and bx are ordered (sorting x would be expensive), the C code seems 
pretty straight-forward and memory-efficient -- create a result vector as long 
as bx, then iterate through x accumulating n and it's sum until x[i] >= bx[i]. 
(I think R's implementation of mean() actually pays more attention to numerical 
error, but with this much data...)

library(inline)
binmean <- cfunction(signature(x="numeric", bx="numeric"),
"   int nx = Rf_length(x), nb = Rf_length(bx), i, j, n;
     SEXP ans = PROTECT(NEW_NUMERIC(nb));
     double sum, *xp = REAL(x), *bxp = REAL(bx), *ansp = REAL(ans);
     sum = j = n = 0;
     for (i = 0; i < nx; ++i) {
         while (xp[i] >= bxp[j]) {
              ansp[j++] = n > 0 ? sum / n : 0;
              sum = n = 0;
         }
         n += 1;
         sum += xp[i];
     }
     ansp[j] = n > 0 ? sum / n : 0;
     UNPROTECT(1);
     return ans;
")

with a test case

nx <- 4e7
nb <- 1e3
x <- sort(runif(nx))
bx <- do.call(seq, c(as.list(range(x)), length.out=nb))

leading to

 > bx1 <- c(bx[-1], bx[nb] + 1)
 > system.time(res <- binmean(x, bx1))
    user  system elapsed
   0.052   0.000   0.050

Martin