[R] replace loops with matrix

[Nikhil Kaza]

what  is nt? is that a typo for ns?
I  don't see why you need to calculate lia within the loop.

Also
library(fBasics)
ccl <-rowprod(lia)



Nikhil Kaza
Asst. Professor,
City and Regional Planning
University of North Carolina

nikhil.list at gmail.com

On Aug 17, 2010, at 6:22 PM, Hey Sky wrote:

> Hey, R users
>
> I am using numerical method for my research paper and the  
> computation burden is
> very heavy. first I tried to do it with loops, example code as  
> following, and it
> take hours to converge for only 200 obs. and my real data has 4000  
> obs. and the
> optimization command that I use is:
>
>
> optim(guess,myfunc1,data=mydata, method="BFGS",hessian=T))
>
> then I tried matrix form computation, it takes only 1/10 of the time  
> the loop
> method takes. it may still have room to improve it. at least, the  
> following
> part looks ugly.
> ccl[,m]<-lia[,1]*lia[,2]*lia[,3]*lia[,4]*lia[,5]
>
> any suggestion are appreciated.
>
>
> The Loop code:
> for(m in 1:ns){
>  for(i in 1:nt){
>     vbar2[,i]=a[1]+     eta[m]+acedu[,i]*a[2]+acwrk[,i]*a[3]
>     vbar3[,i]=b[1]+b[2]*eta[m]+acedu[,i]*b[3]+acwrk[,i]*b[4]
>
>     v8[,i]=1+exp(vbar2[,i])+exp(vbar3[,i])
>
>       for(j in 1:n){
>           if (edu[j,i]==1) lia[j,i]=1/v8[j,i]
>           if (wrk[j,i]==1) lia[j,i]=exp(vbar2[j,i])/v8[j,i]
>           if (home[j,i]==1) lia[j,i]=exp(vbar3[j,i])/v8[j,i]
>      }
>    ccl[,m]<-lia[,i]*ccl[,m]
>  }
> }
>
> The Matrix code:
> for(m in 1:ns){
>     vbar2[,1:nt]=a[1]+     eta[m]+acedu[,1:nt]*a[2]+acwrk[,1:nt]*a[3]
>     vbar3[,1:nt]=b[1]+b[2]*eta[m]+acedu[,1:nt]*b[3]+acwrk[,1:nt]*b[4]
>
>     v8[,1:nt]=1+exp(vbar2[,1:nt])+exp(vbar3[,1:nt])
>
>     lia[1:n,]<-ifelse(edu[1:n,]==1,1/v8[1:n,],
>                            ifelse(wrk[1:n,]==1,exp(vbar2[1:n,])/ 
> v8[1:n,],
>
> exp(vbar3[1:n,])/v8[1:n,]))
>
>    ccl[,m]<-lia[,1]*lia[,2]*lia[,3]*lia[,4]*lia[,5]
> }
>
> Nan
> from Montreal
>
>
>
> ______________________________________________
> 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.