partial cumsum

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com  

-----Original Message-----
From: r-help-bounces at r-project.org 
[mailto:r-help-bounces at r-project.org] On Behalf Of smu
Sent: Wednesday, November 11, 2009 7:58 AM
To: r-help at r-project.org
Subject: [R] partial cumsum

Hello,

I am searching for a function to calculate "partial" cumsums.

For example it should calculate the cumulative sums until a 
NA appears,
and restart the cumsum calculation after the NA.

this:

x <- c(1, 2, 3, NA, 5, 6, 7, 8, 9, 10)

should become this:

1 3 6 NA 5  11  18  26  35  45

Perhaps

   > ave(x, rev(cumsum(rev(is.na(x)))), FUN=cumsum)

    [1]  1  3  6 NA  5 11 18 26 35 45

Nice simple function!  Here's a different approach I use that's faster for long vectors with many NA values.  Note however, that this approach can suffer from catastrophic round-off error because it does a cumsum over the whole vector (after replacing NA's with zeros) and then subtracting out the cumsum at the most recent NA values.  Most of the body of this function is devoted to allowing (an unreasonable degree of) flexibility in specification of where to reset.

cumsum.reset <- function(x, reset.at=which(is.na(x)), na.rm=F) {
    # compute the cumsum of x, resetting the cumsum to 0 at each element indexed by reset.at
    if (is.logical(reset.at)) {
        if (length(reset.at)>length(x)) {
            if ((length(reset.at) %% length(x))!=0)
                stop("length of reset.at must be a multiple of length of x")
            x <- rep(x, len=length(reset.at))
        } else if (length(reset.at)<length(x)) {
            if ((length(x) %% length(reset.at))!=0)
                stop("length of x must be a multiple of length of reset.at")
            reset.at <- rep(reset.at, len=length(x))
        }
        reset.at <- which(reset.at)
    } else if (!is.numeric(reset.at)) {
        stop("reset.at must be logical or numeric")
    }
    if (length(i <- which(reset.at<=1)))
        reset.at <- reset.at[-i]
    if (any(i <- is.na(x[reset.at])))
        x[reset.at[i]] <- 0
    if (na.rm && any(i <- which(is.na(x))))
        x[i] <- 0
    y <- cumsum(x)
    d <- diff(c(0, y[reset.at-1]))
    r <- rep(0, length(x))
    r[reset.at] <- d
    return(y - cumsum(r))
}

x <- c(1, 2, 3, NA, 5, 6, 7, 8, 9, 10)
cumsum.reset(x)

[1]  1  3  6  0  5 11 18 26 35 45

ave(x, rev(cumsum(rev(is.na(x)))), FUN=cumsum)

[1]  1  3  6 NA  5 11 18 26 35 45

The speedup from not breaking the input vector into smaller vectors is (to me) surprisingly small -- only a factor of 3:

x <- replace(rnorm(1e6), sample(1e6, 10000), NA)
all.equal(replace(ave(x, rev(cumsum(rev(is.na(x)))), FUN=cumsum), is.na(x), 0), cumsum.reset(x))

[1] TRUE

system.time(cumsum.reset(x))

user  system elapsed 
   0.31    0.03    0.35

system.time(ave(x, rev(cumsum(rev(is.na(x)))), FUN=cumsum))

user  system elapsed 
   0.99    0.05    1.15

So, I'd go with the ave() approach unless this type of cumsum is the core of a long computationally intensive job.  And if that's the case, it would make sense to code it in C.

-- Tony Plate

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com 
 

any ideas?

thank you and best regards,

    stefan

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide 
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.