Hi, Im new to R and having some trouble with my code - it works, its just
very slow! Ive tried lots of things, but nothing quite seems to work, so any
help and suggestions would be really appreciated!
I want to calculate the marginal likelihood for every element of a row of a
matrix and the corresponding element in every other row. Then sum these for
each row, so I get an upper triangular matrix which consists of the sum of
the marginal likelihoods for row i and row j.
X<-read.table()
X1<-as.matrix(X)
mlr<-function(X1,p,q,Vinv=1){
X2<-X1^2
P<-length(p)
Q<-length(q)
Res<-matrix(ncol=P,nrow=P)
for(i in 1:(P-1)){
lpi<-p[i]
Y1<-matrix(nrow=(P-i),ncol=Q)
Y2<-matrix(nrow=(P-i),ncol=Q)
tB<-matrix(nrow=(P-i),ncol=Q)
for(j in (i+1):P){
lpj<-p[j]
for(k in 1:Q){
Y1[j-i,k]<-sum(X1[c(i,j),k])
Y2[j-i,k]<-sum(X2[c(i,j),k])
}
tB[j-i,]<-(lpi+lpj)*q
}
V1<- (Vinv + tB)^(-1)
c<-tB/2
a1<-1 + c
b1<-1 + (Y2 - V1*(Y1^2))/2
z<-log(V1)/2 + lgamma(a1) - c*log.pi - a1*log(b1)
Res[i,]<-c(rep(NA,i),apply(z,1,sum))
}
return(Res)
}
Its fine for a matrix of 100x100, but the data Im working with is 3538x116,
and so this can take hours!
Any help would be really appreciated!
Thanks,
Anna
--
View this message in context:
http://old.nabble.com/Speed-up-R-code-tp26474898p26474898.html
Sent from the R help mailing list archive at Nabble.com.
Loops tend to dramatically increase computation time. You may re-write a vectorized version of your code if possible, i.e. use matrix algebra. Calculus is a lot faster if one can avoid loops (at least some of them) . Best, Stephane 2009/11/23 AnnaFowler <a.fowler09 at imperial.ac.uk>:> > Hi, Im new to R and having some trouble with my code - it works, its just > very slow! Ive tried lots of things, but nothing quite seems to work, so any > help and suggestions would be really appreciated! > > I want to calculate the marginal likelihood for every element of a row of a > matrix and the corresponding element in every other row. Then sum these for > each row, so I get an upper triangular matrix which consists of the sum of > the marginal likelihoods for row i and row j. > > > X<-read.table() > X1<-as.matrix(X) > > mlr<-function(X1,p,q,Vinv=1){ > X2<-X1^2 > ?P<-length(p) > ?Q<-length(q) > > ?Res<-matrix(ncol=P,nrow=P) > > ?for(i in 1:(P-1)){ > > ? ? ?lpi<-p[i] > ? ? ? ?Y1<-matrix(nrow=(P-i),ncol=Q) > ? ? ? ?Y2<-matrix(nrow=(P-i),ncol=Q) > ? ? ? ?tB<-matrix(nrow=(P-i),ncol=Q) > > ? ? ? ? for(j in (i+1):P){ > ? ? ? ? lpj<-p[j] > ? ? ? ? for(k in 1:Q){ > ? ? ? ? ? ?Y1[j-i,k]<-sum(X1[c(i,j),k]) > ? ? ? ? ? ?Y2[j-i,k]<-sum(X2[c(i,j),k]) > > ? ? ? ? ? ?} > ? ? ? ? ? tB[j-i,]<-(lpi+lpj)*q > ? ? ? ? ?} > > ? ? ? ?V1<- (Vinv + tB)^(-1) > ? ? ? ?c<-tB/2 > ? ? ? ?a1<-1 + c > ? ? ? ?b1<-1 + (Y2 - V1*(Y1^2))/2 > ? ? ? z<-log(V1)/2 + lgamma(a1) - c*log.pi - a1*log(b1) > ? ? ? Res[i,]<-c(rep(NA,i),apply(z,1,sum)) > ?} > ?return(Res) > } > > > Its fine for a matrix of 100x100, but the data Im working with is 3538x116, > and so this can take hours! > > Any help would be really appreciated! > Thanks, > Anna > -- > View this message in context: http://old.nabble.com/Speed-up-R-code-tp26474898p26474898.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. >
Thanks Stephane, Thats a great help! SL-16 wrote:> > Loops tend to dramatically increase computation time. You may re-write > a vectorized version of your code if possible, i.e. use matrix > algebra. Calculus is a lot faster if one can avoid loops (at least > some of them) . > > Best, > Stephane > > 2009/11/23 AnnaFowler <a.fowler09 at imperial.ac.uk>: >> >> Hi, Im new to R and having some trouble with my code - it works, its just >> very slow! Ive tried lots of things, but nothing quite seems to work, so >> any >> help and suggestions would be really appreciated! >> >> I want to calculate the marginal likelihood for every element of a row of >> a >> matrix and the corresponding element in every other row. Then sum these >> for >> each row, so I get an upper triangular matrix which consists of the sum >> of >> the marginal likelihoods for row i and row j. >> >> >> X<-read.table() >> X1<-as.matrix(X) >> >> mlr<-function(X1,p,q,Vinv=1){ >> X2<-X1^2 >> ?P<-length(p) >> ?Q<-length(q) >> >> ?Res<-matrix(ncol=P,nrow=P) >> >> ?for(i in 1:(P-1)){ >> >> ? ? ?lpi<-p[i] >> ? ? ? ?Y1<-matrix(nrow=(P-i),ncol=Q) >> ? ? ? ?Y2<-matrix(nrow=(P-i),ncol=Q) >> ? ? ? ?tB<-matrix(nrow=(P-i),ncol=Q) >> >> ? ? ? ? for(j in (i+1):P){ >> ? ? ? ? lpj<-p[j] >> ? ? ? ? for(k in 1:Q){ >> ? ? ? ? ? ?Y1[j-i,k]<-sum(X1[c(i,j),k]) >> ? ? ? ? ? ?Y2[j-i,k]<-sum(X2[c(i,j),k]) >> >> ? ? ? ? ? ?} >> ? ? ? ? ? tB[j-i,]<-(lpi+lpj)*q >> ? ? ? ? ?} >> >> ? ? ? ?V1<- (Vinv + tB)^(-1) >> ? ? ? ?c<-tB/2 >> ? ? ? ?a1<-1 + c >> ? ? ? ?b1<-1 + (Y2 - V1*(Y1^2))/2 >> ? ? ? z<-log(V1)/2 + lgamma(a1) - c*log.pi - a1*log(b1) >> ? ? ? Res[i,]<-c(rep(NA,i),apply(z,1,sum)) >> ?} >> ?return(Res) >> } >> >> >> Its fine for a matrix of 100x100, but the data Im working with is >> 3538x116, >> and so this can take hours! >> >> Any help would be really appreciated! >> Thanks, >> Anna >> -- >> View this message in context: >> http://old.nabble.com/Speed-up-R-code-tp26474898p26474898.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. > >-- View this message in context: http://old.nabble.com/Speed-up-R-code-tp26474898p26477645.html Sent from the R help mailing list archive at Nabble.com.