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.s...@imail.org 801.408.8111 > -----Original Message----- > From: r-help-boun...@r-project.org [mailto:r-help-boun...@r- > project.org] On Behalf Of William Dunlap > Sent: Tuesday, November 09, 2010 10:07 AM > To: PLucas; r-help@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-boun...@r-project.org > > [mailto:r-help-boun...@r-project.org] On Behalf Of Greg Snow > > Sent: Tuesday, November 09, 2010 8:31 AM > > To: Greg Snow; PLucas; r-help@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.s...@imail.org > > 801.408.8111 > > > > > > > -----Original Message----- > > > From: r-help-boun...@r-project.org [mailto:r-help-boun...@r- > > > project.org] On Behalf Of Greg Snow > > > Sent: Monday, November 08, 2010 1:15 PM > > > To: PLucas; r-help@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.s...@imail.org > > > 801.408.8111 > > > > > > > > > > -----Original Message----- > > > > From: r-help-boun...@r-project.org [mailto:r-help-boun...@r- > > > > project.org] On Behalf Of PLucas > > > > Sent: Monday, November 08, 2010 2:26 AM > > > > To: r-help@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@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. > > > > ______________________________________________ > > 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. ______________________________________________ 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.
