On 24/06/2010 4:08 PM, Lasse Kliemann wrote:
> a <- 0 ; for(i in (1:200000000)) a <- a + 1/i
> b <- 0 ; for(i in (200000000:1)) b <- b + 1/i
> c <- sum(1/(1:200000000))
> d <- sum(1/(200000000:1))
> order(c(a,b,c,d))
[1] 1 2 4 3
> b<c
[1] TRUE
> c==d
[1] FALSE
I'd expected b being the largest, since we sum up the smallest
numbers first. Instead, c is the largest, which is sum() applied
to the vector ordered with largest numbers first.
Can anyone shed some light on this?
I don't know why you'd expect b to be larger than the others. When
you're using sum(), you should expect the most accurate answer, but it
won't necessarily be the biggest: some rounding errors will make the
answer bigger than it should be. For example, if we only had whole
number arithmetic, 0.6 + 0.6 + 0.6 might come out to 3, if each value
was rounded to 1 before adding.
As to whether c or d comes out better: you'd need to program it
yourself if you care about getting the last bit right.
What is the best way in R to compute a sum while avoiding
cancellation effects?
Use sum(). If it's not good enough, then do it in C, accumulating in
extended precision (which is what sum() does), from smallest to largest
if all terms are positive. If there are mixed signs it's a lot harder
to optimize, but I believe the optimal order would be something like the
one that keeps each partial sum as close as possible to zero. It would
rarely be practical to implement that, so summing from smallest absolute
value to largest would be my recommendation.
By the way, sum() in the above example is much faster than the
loops, so it would be nice if we could utilize it.
Why do you think you can't?
Duncan Murdoch
------------------------------------------------------------------------
______________________________________________
R-help@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@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.
