[R] cumsum vs. sum
Duncan Murdoch
murdoch at stats.uwo.ca
Wed Feb 18 19:28:26 CET 2009
On 18/02/2009 12:41 PM, Stavros Macrakis wrote:
> Hmm. Why not use the same method to guarantee the same result? Or at
> least document the possibility that cumsum(x)[length(x)] != sum(x)...
> that seems like an easy trap to fall into.
Assuming equality of floating point numbers computed by two different
paths is always a trap.
R doesn't try to obtain results that are equal to the last bit in other
circumstances; why should it do so here? For example, one somewhat
controversial choice in R is to use 64 bit precision in intermediate
computations when available, rather than rounding everything to 52 bits
as it does when stored to memory in doubles. This means that the value
you get is likely to be closer to the truth than if you did the rounding
earlier, but it is also subject to change according to optimization
level, compiler version, etc.
Duncan Murdoch
>
> -s
>
> On Wed, Feb 18, 2009 at 11:39 AM, Martin Maechler
> <maechler at stat.math.ethz.ch> wrote:
>>>>>>> "SM" == Stavros Macrakis <macrakis at alum.mit.edu>
>>>>>>> on Wed, 18 Feb 2009 10:00:40 -0500 writes:
>> SM> Nice! Glad to hear it. It sounds as though it is still possible for
>> SM> cumsum(x)[length(x)] to not be exactly equal to sum, though?
>>
>> Well, possible, probably yes, platform-dependently;
>> However I vaguely remember that I didn't see one such case in the few
>> experiments I did.
>>
>> Martin
>>
>> SM> On Wed, Feb 18, 2009 at 8:03 AM, Martin Maechler
>> SM> <maechler at stat.math.ethz.ch> wrote:
>> SM> ...
>> >> o cumsum(x) and cumprod(x) for double precision x now use a long
>> >> double accumulator where available and so more closely match
>> >> sum() and prod() in potentially being more accurate.
>>
>
> ______________________________________________
> 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.
More information about the R-help
mailing list