Hi All,
I'm writing a long code that takes long time to execute. So I used the
Rprof() function and found out that the function that takes about 80% of
the time is the incomb () fucntion (below), and this is most probably
because of the many explicit for() loops I'm using.
n=18;m=4;p=0.3;alpha=0.2;beta=2
x=c(3,0,0)
LD<-list()
for (i in 1:(m-1)) {
LD[[i]]<-seq(0,x[i],1)
}
LD[[m]]<-seq(0,(n-m-sum(x)),1)
LED<-expand.grid (LD)
LED<-as.matrix(LED)
store1<-numeric(nrow(LED))
h<- numeric(m-1)
lm<- numeric(m-1)
for (j in 1:length(store1) )
{
incomb<-function(x,alpha,beta) {
g<-((-1)^(sum(LED[j,])))*(gamma((1/beta)+1))*((alpha)^(-(1/beta)))
for (i in 1:(m-1)) {
h[i]<- choose(x[i],LED[j,i])
}
ik<-prod(h)*choose((n-m-sum(x)),LED[j,m])
for (i in 1:(m-1)) {
lm[i]<-(sum(LED[j,1:i])) + i
}
plm<-prod(lm)
gil<-g*ik/(plm)
hlm<-numeric(sum(LED[j,])+(m-1))
dsa<-length(hlm)
for (i in 1:dsa)
{
ppp<- sum(LED[j,])+(m-1)
hlm[i]<-
(choose(ppp,i))*((-1)^(i))*((i+1)^((-1)*((1/beta)+1)))
}
shl<-gil*(sum(hlm)+1)
return (shl)
}
store1[j]<-incomb(x,alpha=0.2,beta=2)
}
I'm trying to use alternatives (for ex. to vectorize things) to the
explicit for() loops, but things don't work out.
Any suggestions that can help me to speed up the execution of the incomb()
function are much appreciated.
Thanks a lot in advance.
Maram Salem
[[alternative HTML version deleted]]
On 01/11/2015 7:31 AM, Maram SAlem wrote:> Hi All, > > I'm writing a long code that takes long time to execute. So I used the > Rprof() function and found out that the function that takes about 80% of > the time is the incomb () fucntion (below), and this is most probably > because of the many explicit for() loops I'm using. > > n=18;m=4;p=0.3;alpha=0.2;beta=2 > x=c(3,0,0) > LD<-list() > for (i in 1:(m-1)) { > LD[[i]]<-seq(0,x[i],1) > } > LD[[m]]<-seq(0,(n-m-sum(x)),1) > LED<-expand.grid (LD) > LED<-as.matrix(LED) > store1<-numeric(nrow(LED)) > h<- numeric(m-1) > lm<- numeric(m-1) > for (j in 1:length(store1) ) > { > incomb<-function(x,alpha,beta) { > > g<-((-1)^(sum(LED[j,])))*(gamma((1/beta)+1))*((alpha)^(-(1/beta))) > for (i in 1:(m-1)) { > h[i]<- choose(x[i],LED[j,i]) > } > ik<-prod(h)*choose((n-m-sum(x)),LED[j,m]) > for (i in 1:(m-1)) { > lm[i]<-(sum(LED[j,1:i])) + i > } > plm<-prod(lm) > gil<-g*ik/(plm) > hlm<-numeric(sum(LED[j,])+(m-1)) > dsa<-length(hlm) > for (i in 1:dsa) > { > ppp<- sum(LED[j,])+(m-1) > hlm[i]<- > (choose(ppp,i))*((-1)^(i))*((i+1)^((-1)*((1/beta)+1))) > } > shl<-gil*(sum(hlm)+1) > return (shl) > } > store1[j]<-incomb(x,alpha=0.2,beta=2) > } > > > I'm trying to use alternatives (for ex. to vectorize things) to the > explicit for() loops, but things don't work out. > > Any suggestions that can help me to speed up the execution of the incomb() > function are much appreciated. >Generally I'd suggest starting with your innermost loops, and vectorize those. Continue outwards from there if you can. For example, you have> for (i in 1:(m-1)) { > h[i]<- choose(x[i],LED[j,i]) > }which could be vectorized as i <- 1:(m-1) h[i] <- choose(x[i], LED[j, i]) and maybe as h <- choose(x, LED[j, ]) (but I haven't read your code closely enough to know if those are equivalent). If m is large, this will make a big difference to that part of the code. Duncan Murdoch
Why are you recreating the incomb function within the loop instead of defining it outside the loop? Also you are referencing several variables that are global (e.g., m & j); you should be passing these in as parameters to the function. Jim Holtman Data Munger Guru What is the problem that you are trying to solve? Tell me what you want to do, not how you want to do it. On Sun, Nov 1, 2015 at 7:31 AM, Maram SAlem <marammagdysalem at gmail.com> wrote:> Hi All, > > I'm writing a long code that takes long time to execute. So I used the > Rprof() function and found out that the function that takes about 80% of > the time is the incomb () fucntion (below), and this is most probably > because of the many explicit for() loops I'm using. > > n=18;m=4;p=0.3;alpha=0.2;beta=2 > x=c(3,0,0) > LD<-list() > for (i in 1:(m-1)) { > LD[[i]]<-seq(0,x[i],1) > } > LD[[m]]<-seq(0,(n-m-sum(x)),1) > LED<-expand.grid (LD) > LED<-as.matrix(LED) > store1<-numeric(nrow(LED)) > h<- numeric(m-1) > lm<- numeric(m-1) > for (j in 1:length(store1) ) > { > incomb<-function(x,alpha,beta) { > > g<-((-1)^(sum(LED[j,])))*(gamma((1/beta)+1))*((alpha)^(-(1/beta))) > for (i in 1:(m-1)) { > h[i]<- choose(x[i],LED[j,i]) > } > ik<-prod(h)*choose((n-m-sum(x)),LED[j,m]) > for (i in 1:(m-1)) { > lm[i]<-(sum(LED[j,1:i])) + i > } > plm<-prod(lm) > gil<-g*ik/(plm) > hlm<-numeric(sum(LED[j,])+(m-1)) > dsa<-length(hlm) > for (i in 1:dsa) > { > ppp<- sum(LED[j,])+(m-1) > hlm[i]<- > (choose(ppp,i))*((-1)^(i))*((i+1)^((-1)*((1/beta)+1))) > } > shl<-gil*(sum(hlm)+1) > return (shl) > } > store1[j]<-incomb(x,alpha=0.2,beta=2) > } > > > I'm trying to use alternatives (for ex. to vectorize things) to the > explicit for() loops, but things don't work out. > > Any suggestions that can help me to speed up the execution of the incomb() > function are much appreciated. > > Thanks a lot in advance. > > Maram Salem > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see > 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. >[[alternative HTML version deleted]]
Hi Jim,
Thanks a lot for replying.
In fact I'm trying to run a simulation study that enables me to calculate
the Bayes risk of a sampling plan selected from progressively type-II
censored Weibull model. One of the steps involves evaluating the expected
test time, which is a rather complicated formula that involves nested
multiple summations where the counters of the summation signs are
dependent, that's why I thought of I should create the incomb() function
inside the loop, or may be I didn't figure out how to relate its arguments
to the ones inside the loop had I created it outside it. I'm trying to
create a matrix of all the possible combinations involved in the summations
and then use the apply() function on each row of that matrix. The problem
is that the code I wrote works perfectly well for rather small values of
the sample size,n, and the censoring number, m (for example, n=8,m=4),but
when n and m are increased (say, n=25,m=15) the code keeps on running for
days with no output. That's why I thought I should try to avoid explicit
loops as much as possible, so I did my best in this regard but still the
code takes too long to execute,(more than three days), thus, i believe
there must be something wrong.
Here's the full code:
library(pbapply)
f1 <- function(n, m) {
stopifnot(n > m)
r0 <- t(diff(combn(n-1, m-1)) - 1L)
r1 <- rep(seq(from=0, len=n-m+1), choose( seq(to=m-2, by=-1, len=n-m+1),
m-2))
cbind(r0[, ncol(r0):1, drop=FALSE], r1, deparse.level=0)
}
simpfun<- function (x,n,m,p,alpha,beta)
{
a<-factorial(n-m)/(prod((factorial(x)))*(factorial((n-m)- sum(x))))
b <- ((m-1):1)
c<- a*((p)^(sum(x)))*((1-p)^(((m-1)*(n-m))- sum(x%*%(as.matrix(b)))))
d <- n - cumsum(x) - (1:(m-1))
e<- n*(prod(d))*c
LD<-list()
for (i in 1:(m-1)) {
LD[[i]]<-seq(0,x[i],1)
}
LD[[m]]<-seq(0,(n-m-sum(x)),1)
LED<-expand.grid (LD)
LED<-as.matrix(LED)
store1<-numeric(nrow(LED))
for (j in 1:length(store1) )
{
incomb<-function(x,alpha,beta) {
g<-((-1)^(sum(LED[j,])))*(gamma((1/beta)+1))*((alpha)^(-(1/beta)))
h <- choose(x, LED[j,-m])
ik<-prod(h)*choose((n-m-sum(x)),LED[j,m])
lm<-cumsum(LED[j,-m]) + (1:(m-1))
plm<-prod(lm)
gil<-g*ik/(plm)
hlm<-numeric(sum(LED[j,])+(m-1))
dsa<-length(hlm)
for (i in 1:dsa)
{
ppp<- sum(LED[j,])+(m-1)
hlm[i]<-
(choose(ppp,i))*((-1)^(i))*((i+1)^((-1)*((1/beta)+1)))
}
shl<-gil*(sum(hlm)+1)
return (shl)
}
store1[j]<-incomb(x,alpha=0.2,beta=2)
}
val1<- sum(store1)*e
return(val1)
}
va<-pbapply(s,1,simpfun,n=6,m=4,p=0.3,alpha=0.2,beta=2)
EXP<-sum(va)
Any help would be greatly appreciated.
Thanks a lot for your time.
Best Regards,
Maram Salem
On 2 November 2015 at 00:27, jim holtman <jholtman at gmail.com> wrote:
> Why are you recreating the incomb function within the loop instead of
> defining it outside the loop? Also you are referencing several variables
> that are global (e.g., m & j); you should be passing these in as
parameters
> to the function.
>
>
> Jim Holtman
> Data Munger Guru
>
> What is the problem that you are trying to solve?
> Tell me what you want to do, not how you want to do it.
>
> On Sun, Nov 1, 2015 at 7:31 AM, Maram SAlem <marammagdysalem at
gmail.com>
> wrote:
>
>> Hi All,
>>
>> I'm writing a long code that takes long time to execute. So I used
the
>> Rprof() function and found out that the function that takes about 80%
of
>> the time is the incomb () fucntion (below), and this is most probably
>> because of the many explicit for() loops I'm using.
>>
>> n=18;m=4;p=0.3;alpha=0.2;beta=2
>> x=c(3,0,0)
>> LD<-list()
>> for (i in 1:(m-1)) {
>> LD[[i]]<-seq(0,x[i],1)
>> }
>> LD[[m]]<-seq(0,(n-m-sum(x)),1)
>> LED<-expand.grid (LD)
>> LED<-as.matrix(LED)
>> store1<-numeric(nrow(LED))
>> h<- numeric(m-1)
>> lm<- numeric(m-1)
>> for (j in 1:length(store1) )
>> {
>> incomb<-function(x,alpha,beta) {
>>
>> g<-((-1)^(sum(LED[j,])))*(gamma((1/beta)+1))*((alpha)^(-(1/beta)))
>> for (i in 1:(m-1)) {
>> h[i]<- choose(x[i],LED[j,i])
>> }
>> ik<-prod(h)*choose((n-m-sum(x)),LED[j,m])
>> for (i in 1:(m-1)) {
>> lm[i]<-(sum(LED[j,1:i])) + i
>> }
>> plm<-prod(lm)
>> gil<-g*ik/(plm)
>> hlm<-numeric(sum(LED[j,])+(m-1))
>> dsa<-length(hlm)
>> for (i in 1:dsa)
>> {
>> ppp<- sum(LED[j,])+(m-1)
>> hlm[i]<-
>> (choose(ppp,i))*((-1)^(i))*((i+1)^((-1)*((1/beta)+1)))
>> }
>> shl<-gil*(sum(hlm)+1)
>> return (shl)
>> }
>> store1[j]<-incomb(x,alpha=0.2,beta=2)
>> }
>>
>>
>> I'm trying to use alternatives (for ex. to vectorize things) to the
>> explicit for() loops, but things don't work out.
>>
>> Any suggestions that can help me to speed up the execution of the
incomb()
>> function are much appreciated.
>>
>> Thanks a lot in advance.
>>
>> Maram Salem
>>
>> [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> 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.
>>
>
>
[[alternative HTML version deleted]]