Hi, Today I was flabbergasted to see something that looks like a rounding error in the very basic seq function in R.> a = seq(0.1,0.9,by=0.1) > a[1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9> a[1] == 0.1[1] TRUE> a[2] == 0.2[1] TRUE> a[3] == 0.3[1] FALSE It turns out that the alternative> a = (1:9)/10works just fine. Are there any good guides out there on how to deal with issues like this? I am normally aware of rounding errors, but it really surprised me to see that an elementary function like seq would behave in this way. Thanks, Michael Knudsen -- Michael Knudsen micknudsen at gmail.com http://sites.google.com/site/micknudsen/
On 9/30/2009 2:40 PM, Michael Knudsen wrote:> Hi, > > Today I was flabbergasted to see something that looks like a rounding > error in the very basic seq function in R. > >> a = seq(0.1,0.9,by=0.1) >> a > [1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 >> a[1] == 0.1 > [1] TRUE >> a[2] == 0.2 > [1] TRUE >> a[3] == 0.3 > [1] FALSE > > It turns out that the alternative > >> a = (1:9)/10 > > works just fine. Are there any good guides out there on how to deal > with issues like this? I am normally aware of rounding errors, but it > really surprised me to see that an elementary function like seq would > behave in this way.Why? You asked for an increment of 1 in the second case (which is exactly represented in R), then divided by 10, so you'll get the same as 0.3 gives you. In the seq() case you asked for an increment of a number close to but not equal to 1/10 (because 1/10 is not exactly representable in R), so you got something different. Duncan Murdoch
On Wed, Sep 30, 2009 at 8:40 PM, Michael Knudsen <micknudsen at gmail.com> wrote:>> a = seq(0.1,0.9,by=0.1) >> a > [1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 >> a[1] == 0.1 > [1] TRUE >> a[2] == 0.2 > [1] TRUE >> a[3] == 0.3 > [1] FALSEA friend of mine just pointed out a possible solution:> a=seq(0.1,0.9,by=0.1) > a = seq(0.1,0.9,by=0.1) > a[3]==0.3[1] FALSE> all.equal(a[3],0.3)[1] TRUE The all.equal function checks if two objects are "nearly" equal. -- Michael Knudsen micknudsen at gmail.com http://sites.google.com/site/micknudsen/
See: http://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-doesn_0027t-R-think-these-numbers-are-equal_003f On Wed, Sep 30, 2009 at 2:40 PM, Michael Knudsen <micknudsen at gmail.com> wrote:> Hi, > > Today I was flabbergasted to see something that looks like a rounding > error in the very basic seq function in R. > >> a = seq(0.1,0.9,by=0.1) >> a > [1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 >> a[1] == 0.1 > [1] TRUE >> a[2] == 0.2 > [1] TRUE >> a[3] == 0.3 > [1] FALSE > > It turns out that the alternative > >> a = (1:9)/10 > > works just fine. Are there any good guides out there on how to deal > with issues like this? I am normally aware of rounding errors, but it > really surprised me to see that an elementary function like seq would > behave in this way. > > Thanks, > Michael Knudsen > > -- > Michael Knudsen > micknudsen at gmail.com > http://sites.google.com/site/micknudsen/ > > ______________________________________________ > 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. >
For my own edification more than anything (I never took computer science): is> a = seq(0.1,0.9,by=0.1) > a <- as.character(a) > a[3] == "0.3"[1] TRUE safe? -Ista On Wed, Sep 30, 2009 at 3:46 PM, cls59 <chuck at sharpsteen.net> wrote:> > > Martin Batholdy wrote: >> >> hum, >> >> can you explain that a little more detailed? >> Perhaps I miss the background knowledge - but it seems just absurd to >> me. >> >> 0.1+0.1+0.1 is 0.3 - there is no rounding involved, is there? >> >> > > Unfortunately this comes as an utter shock to many people who never take a > Computer Science course. I watch it nail engineering students all the time. > > Basically, if you have a fraction and the denominator is not equal to 2^n > for some integer n, that fraction will NEVER be stored as an exact "floating > point" number-- instead it will contain some error due to concessions that > must be made in order to use an efficient binary number scheme. > > These errors are generally small, but they do propagate-- especially if you > are carrying the same numbers through a large computation. A good example is > large-scale numerical solutions to nonlinear problems where iterative > algorithms are employed repetitively at each solution step. As the > calculation progresses the roundoff error can rot away the computational > soundness of the algorithm. > > If this concerns you, I would suggest reading up on common internal > representations of floating point numbers as well as the propagation of > roundoff error. > > At the very least I hope this revelation will instill an appropriate sense > of paranoia concerning the numbers calculated by those magic boxes sitting > on our desks. > > -Charlie > > ----- > Charlie Sharpsteen > Undergraduate > Environmental Resources Engineering > Humboldt State University > -- > View this message in context: http://www.nabble.com/Rounding-error-in-seq%28...%29-tp25686630p25687626.html > Sent from the R help mailing list archive at Nabble.com. > > ______________________________________________ > 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. >-- Ista Zahn Graduate student University of Rochester http://yourpsyche.org
>> can you explain that a little more detailed?>> Perhaps I miss the background knowledge - but it seems just absurd to me. >> >> 0.1+0.1+0.1 is 0.3 - there is no rounding involved, is there? >> >> why is >> x <- 0.1 + 0.1 +0.1 >> not equal to >> y <- 0.3 > > Remember that this is in BINARY arithmetic. It's really not any stranger > than the fact that 1/3 + 1/3 != 2/3 in finite accuracy decimal arithmetic > (0.33333 + 0.33333 = 0.66666 != 0.66667). A nice description of this is in section 2.2.4 of Braun and Murdoch's book. They give a striking example. > x <- 1:11 > mean(x) [1] 6 > var(x) [1] 11 > sum( (x - mean(x))^2) / 10 # two-pass formula for mean [1] 11 > (sum(x^2) - 11 * mean(x)^2) / 10 # one-pass formula for mean [1] 11 > A <- 1.e10 # add large value to each x_i > x <- 1:11 + A > var(x) # R function survives [1] 11 > sum( (x - mean(x))^2) / 10 # two-pass formula survives [1] 11 > (sum(x^2) - 11 * mean(x)^2) / 10 # one-pass formula suffers catastrophic loss of precision! [1] 0