[R] How to eliminate this for loop ?
Greg Snow
Greg.Snow at imail.org
Tue Nov 9 18:34:14 CET 2010
OK, double oops. I first tested my code with length 100, then upped the number but forgot to up the preallocation part, I should have used a variable there instead so that only one place needed to be changed.
My version did have problems when I tried to do a vector of length 10,000, some values were NaN probably due to b^largenumber being essentially 0, then being in the denominator.
Though for really long vectors the round off error in any version could accumulate to the point of affecting the results. This could start a debate about whether the missing value could be seen as better than a potentially incorrect non-missing value. It would mainly depend on the purpose and I don't think there would be a general preference either way.
--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.snow at imail.org
801.408.8111
> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-
> project.org] On Behalf Of William Dunlap
> Sent: Tuesday, November 09, 2010 10:07 AM
> To: PLucas; r-help at r-project.org
> Subject: Re: [R] How to eliminate this for loop ?
>
> Note that for long vectors the OP's code would
> go much faster if he preallocated the output vector
> a to its eventual length. I.e., start with
> a <- numeric(N)
> instead of
> a <- c()
> I defined 2 functions that differed only in how
> a was initialized
> f0 <- function(b, c) {
> N <- length(c)
> a <- c()
> a[1] <- 1; # initial value
> for(i in 2:N) {
> a[i]<-a[i-1]*b - c[i-1] # b is a value, c is another vector
> }
> a
> }
>
> f1 <- function(b, c) {
> N <- length(c)
> a<- numeric(N)
> a[1] <- 1; # initial value
> for(i in 2:N) {
> a[i]<-a[i-1]*b - c[i-1] # b is a value, c is another vector
> }
> a
> }
> and timed them for a 100,000 long vector:
> > c <- rnorm(1e5)
> > system.time(a0 <- f0(b=0.5, c=c))
> user system elapsed
> 17.270 1.410 18.704
> > system.time(a1 <- f1(b=0.5, c=c))
> user system elapsed
> 0.400 0.000 0.401
> > identical(a0, a1)
> [1] TRUE
> If that is not fast enough then you have to
> start thinking harder. E.g., look at functions
> like filter() and/or do some algebra.
>
> (Greg's code also had an error of that sort,
> preallocating 100 entries where the eventual
> length was 1000).
>
> Bill Dunlap
> Spotfire, TIBCO Software
> wdunlap tibco.com
>
> > -----Original Message-----
> > From: r-help-bounces at r-project.org
> > [mailto:r-help-bounces at r-project.org] On Behalf Of Greg Snow
> > Sent: Tuesday, November 09, 2010 8:31 AM
> > To: Greg Snow; PLucas; r-help at r-project.org
> > Subject: Re: [R] How to eliminate this for loop ?
> >
> > Oops, my version added cc instead of subtracted, it still
> > works if you multiply cc by -1 (except the initial 1).
> >
> > --
> > Gregory (Greg) L. Snow Ph.D.
> > Statistical Data Center
> > Intermountain Healthcare
> > greg.snow at imail.org
> > 801.408.8111
> >
> >
> > > -----Original Message-----
> > > From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-
> > > project.org] On Behalf Of Greg Snow
> > > Sent: Monday, November 08, 2010 1:15 PM
> > > To: PLucas; r-help at r-project.org
> > > Subject: Re: [R] How to eliminate this for loop ?
> > >
> > > If you are willing to shift the c vector by 1 and have 1
> > (the initial
> > > value) as the start of c, then you can just do:
> > >
> > > cumsum( cc * b^( (n-1):0 ) ) / b^( (n-1):0 )
> > >
> > > to compare:
> > >
> > > cc <- c(1, rnorm(999) )
> > > b <- 0.5
> > > n <- length(cc)
> > >
> > > a1 <- numeric(100)
> > > a1[1] <- 1
> > >
> > > system.time(for(i in 2:n ) {
> > > a1[i] <- b*a1[i-1] + cc[i]
> > > })
> > >
> > > system.time(a2 <- cumsum( cc * b^( (n-1):0 ) ) / b^( (n-1):0 ))
> > >
> > > all.equal(a1,a2)
> > >
> > > Though you could have problems with the b^ part if the
> > length gets too
> > > long.
> > >
> > > --
> > > Gregory (Greg) L. Snow Ph.D.
> > > Statistical Data Center
> > > Intermountain Healthcare
> > > greg.snow at imail.org
> > > 801.408.8111
> > >
> > >
> > > > -----Original Message-----
> > > > From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-
> > > > project.org] On Behalf Of PLucas
> > > > Sent: Monday, November 08, 2010 2:26 AM
> > > > To: r-help at r-project.org
> > > > Subject: [R] How to eliminate this for loop ?
> > > >
> > > >
> > > > Hi, I would like to create a list recursively and eliminate my
> for
> > > loop
> > > > :
> > > >
> > > > a<-c()
> > > > a[1] <- 1; # initial value
> > > > for(i in 2:N) {
> > > > a[i]<-a[i-1]*b - c[i-1] # b is a value, c is another vector
> > > > }
> > > >
> > > >
> > > > Is it possible ?
> > > >
> > > > Thanks
> > > > --
> > > > View this message in context:
> > http://r.789695.n4.nabble.com/How-to-
> > > > eliminate-this-for-loop-tp3031667p3031667.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.
> > >
> > > ______________________________________________
> > > 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.
> >
> > ______________________________________________
> > 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.
> >
>
> ______________________________________________
> 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