[R] Fwd: ts.intersect bug?
Barry Rowlingson
B.Rowlingson at lancaster.ac.uk
Wed Sep 21 17:54:12 CEST 2005
Antonio, Fabio Di Narzo wrote:
> Ooops, bad example.
> Try this instead:
> a <- ts(runif(6500), start=0, freq=10)
> b <- lag(a, 1)
> c <- ts.intersect(a, b)
>
> Gives an error from .cbind.ts
or even:
> a=ts(numeric(2600),freq=10,start=0);b=lag(a,1);c=ts.intersect(a,b)
Error in "[<-"(`*tmp*`, , (1 + cs[i]):cs[i + 1], value = c(0, 0, 0, 0, :
number of items to replace is not a multiple of replacement length
seems to be sensitive to the lengths of the time series:
> a=ts(numeric(2601),freq=10,start=0);b=lag(a,1);c=ts.intersect(a,b)
> (works fine)
Digging into the code, the window() function is returning a different
length time series for each one in these failing cases. I reckon its a
floating-point precision situation, where the last time series point
should be included but the arithmetic precision of a 2600-long series at
separation of 1/10 is leaving it out.
> st;en
[1] 0
[1] 259.8
> a=ts(numeric(2600),freq=10,start=0);b=lag(a,1)
> length(window(a,st,en))
[1] 2599
> length(window(b,st,en))
[1] 2598
- ts.intersect is trying to put these two time series together, and so
fails. But:
> a=ts(numeric(2601),freq=10,start=0);b=lag(a,1)
> length(window(a,st,en))
[1] 2599
> length(window(b,st,en))
[1] 2599
- works.
Note that en is not precisely 259.8:
> en == 259.8
[1] FALSE
> en-259.8
[1] -5.684342e-14
I've computed 'en' as the .cbind.ts function does, and its not exactly
259.8. If it were, then it would work... Perhaps .cbind.ts should round
to the nearest true time point or something...
Note that it fails in plenty of smaller cases too:
> a=ts(numeric(13),freq=10,start=0);b=lag(a,1);c=ts.intersect(a,b)
Error in "[<-"(`*tmp*`, , (1 + cs[i]):cs[i + 1], value = c(0, 0, 0, 0, :
number of items to replace is not a multiple of replacement length
Seems to not like 13s and 10s and integer products thereof (2600,
6500). Are you superstitious?
Baz
PS without counting, how many letters are there in 'superstitious'?
More information about the R-help
mailing list