[R] how to make simulation faster
stats12
skarmv at gmail.com
Fri Oct 26 15:08:17 CEST 2012
Hi,
Thank you for your reply. I updated my post with the code. Also, about
posting from Nabble, since I am a new user I didn't know about that problem.
If I post to the mailing list ( r-help at r-project.org), would it get rid of
that problem?
output1<-vector("numeric",length(1:r))
output2<-vector("numeric",length(1:r))
output3<-vector("numeric",length(1:r))
output4<-vector("numeric",length(1:r))
for (m in 1:r){
ll<-function(p){
cumhaz<-(time*exp(Zi*p[3]))*p[1]
cumhaz<-aggregate(cumhaz,by=list(i),FUN=sum)
lnhaz<-delta*log(exp(Zi*p[3])*p[1])
lnhaz<-aggregate(lnhaz,by=list(i),FUN=sum)
lnhaz<-lnhaz[2:(r+1)]
cumhaz<-cumhaz[2:(r+1)]
lik<-r[m]*log(p[2])-
sum((di[,m]+1/p[2])*log(1+cumhaz[,m]*p[2]))+sum(lnhaz[,m])-
n*log(gamma(1/p[2]))+sum(log(gamma(di[,m]+1/p[2])))
-lik
}
initial<-c(1,0.5,0.5)
estt<-nlm(ll,initial,hessian=T)
lambda<-estt$estimate[1]
theta<-estt$estimate[2]
beta<-estt$estimate[3]
hessian<-t$hessian
cov<-solve(hessian)
vtheta<-cov[2,2]
vbeta<-cov[3,3]
output1[m]<-theta
output2[m]<-beta
output3[m]<-vtheta
output4[m]<-vbeta
}
theta<-as.vector(output1)
beta<-as.vector(output2)
vtheta<-as.vector(output3)
vbeta<-as.vector(output4)
Code would be nice: I struggle to think of an basic MLE fitting that
would take ~36s per iteration and scale so badly if you are writing
idiomatic R:
http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example
Finally, I note you're posting from Nabble. Please include context in
your reply -- I don't believe Nabble does this automatically, so
you'll need to manually include it. Most of the regular respondents on
this list don't use Nabble -- it is a _mailing list_ after all -- so
we don't get the forum view you do, only emails of the individual
posts. Combine that with the high volume of posts, and it's quite
difficult to trace a discussion if we all don't make sure to include
context.
--
View this message in context: http://r.789695.n4.nabble.com/how-to-make-simulation-faster-tp4647492p4647541.html
Sent from the R help mailing list archive at Nabble.com.
More information about the R-help
mailing list