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.l...@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@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.
