[R] Alternatives for explicit for() loops

[Duncan Murdoch]

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