From hb at stat.berkeley.edu  Tue Dec  4 19:47:16 2007
From: hb at stat.berkeley.edu (Henrik Bengtsson)
Date: Tue, 4 Dec 2007 10:47:16 -0800
Subject: [Rd] seq() does not create proper numbers (PR#10489)
In-Reply-To: <47557050.6000603@stats.uwo.ca>
References: <20071204121012.9F5632834155@mail.pubhealth.ku.dk>
	<47557050.6000603@stats.uwo.ca>
Message-ID: <59d7961d0712041047s2dcc38ebu60c81e570005c02c@mail.gmail.com>

On 04/12/2007, Duncan Murdoch <murdoch at stats.uwo.ca> wrote:
> On 12/4/2007 7:10 AM, rihle at gwdg.de wrote:
> > Full_Name: Rico Ihle
> > Version: 2.6.1
> > OS: Windows XP Professional Version 2002 Pentium(R) 4 CPU 3.2 GHz
> > Submission from: (NULL) (134.76.183.24)
>
> This is not a bug.  See FAQ 7.31, "Why doesn't R think these numbers are
> equal?"
>
> >
> >
> > # Bug in seq() function:
> > x <- seq(-2,2,by=0.01)
>
> You've specified that the step size should be some number that's close
> to 0.01, but not exactly 0.01, since R doesn't know how to represent that.

Here is another example adding to the clarification (confusion?):

n <- 1e30;
print(n*(1/n) < 1);
## [1] TRUE

And another one:

c <- 100;
from <- -2;
to <- 2;
x <- seq(from=from, to=to, by=1/c);
y <- seq(from=from*c, to=to*c, by=1)/c;
print(identical(x,y));
## [1] FALSE
print(all.equal(x,y));
## [1] TRUE

Yet another example of problems with representation of certain values:

by <- 1/c;
z <- seq(from=from/by, to=to/by, by=1)*by;
print(identical(x,z));
## [1] FALSE
print(identical(y,z));  # <<< NOTE
## [1] FALSE
print(all.equal(x,z));
## [1] TRUE
print(all.equal(y,z));
## [1] TRUE

I guess the question is whether you prefer the 'x' or the 'y'
sequence.  Actually, if the representation of 'by=1/c' is, say,
strictly smaller than 1/c, for some simple tests it looks like the
deviation of d[k]=y[k]-x[k] increases as k grows.  Here is a plot
illustrating this:

c <- 10;
to <- 1e6;
x <- seq(from=0, to=to, by=1/c);
y <- seq(from=0, to=to*c, by=1)/c;
d <- y-x;
n <- length(d);

# Plot (t,d), where t is uniform on the logaritmic scale,
# because otherwise there are too many data points.
t <- seq(from=0, to=log10(n), length.out=1e3);
t <- unique(as.integer(10^t));
plot(t, d[t]);

/Henrik

>
> > which(x==0.05)# How that??
> > # although:
> > x# 0.05 seems to be at position 206 of x!!:
> > x[206]
> > # Why is this not equal to 0.05?
> >
> > # Reason:
> > x2 <- as.character(x);x2
> > x2[206]# Ooooh... It's really not equal to 0.05!! How that? (compare lines 5 and
> > 6!)
> >
> > # Remedy:
> > x3 <- round(as.numeric(x),2)
> > which(x3==0.05)
> >
> > # The (necessary) rounding is apparently and unfortunately NOT included in the
> > seq() function!!!
> > # But it should be!!!
>
> Why should it?  You didn't specify a round number as the step size.
>
> Duncan Murdoch
>
> > # Because if one doesn't know about the demonstrated "nice" feature of the seq()
> > function
> > # (and it is not visible in lines 5 or 6!!!)
> > # one gets mad that x[206] is not equal to 0.05 although x[206] is printed as
> > 0.05!!!
> >
> > # Similarly:
> > y <-  seq(-0.5,.5,by=0.01)
> > which(y == 0.05)# None? How that? Result should be 56!!
> > y[56]
> >
> > # but:
> > y2 <-as.character(y)
> > y2[56]
> > which(y2 == 0.05)
> >
> > # or rounding alternatively:
> > which(round(y,2) == 0.05)
> >
> > ______________________________________________
> > R-devel at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-devel
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>



[Rd] seq() does not create proper numbers (PR#10489)
Henrik Bengtsson