How can I speed up the following function in R?

Question

I wrote the following function in R. I want to iterate it, say 50000 times. I used "sapply" in my function but it runs slowly in R.My PC is still working about 20h now and I have no idea about the run time. Any ideas how to speed up this operation? Thanks.

data=matrix(c(0.01132162,1,0.04056053,1,0.11676735,0,0.12029087,1,
           0.16197702,1,0.17190980,1,0.20386841,1,0.21251687,0,
           0.36536492,0,0.40256414,1),ncol=2,byrow=T)

 GIBS=function(data,a,b,beta,R)
    {
       m=length(R)
       n1=sum(data[,2]==1)
       n2=m-n1
       n=m+sum(R)
       N=c(n1,n2)
       R1=c(0,R)
       nstar=c()
       for(i in 1:m) nstar[i]=n-(i-1)-sum(R1[1:i])
       Z=c(data[1,1],data[2:m,1]-data[1:(m-1),1])
       f=function(x)
        {
          A=0
          for(i in 1:m) A=A+x^i*nstar[i]*Z[i]
          FR=1
          for(j in 1:2) FR=FR*(A+b[j])^(N[j]+a[j])
          return(x^(m*(m+1)/2)*exp(-beta*(x-1))/FR)
        }
      INT=integrate(f,1,Inf)$value
      SG=function(it)
       {
        uu=runif(1)
        g0=function(t) integrate(f,1,t)$value/INT-uu
        aa=5
        if(g0(1)>0) {while(g0(aa)>0) aa=aa+1} else {while(g0(aa)<0) aa=aa+1}
        ra=uniroot(g0,c(1,aa))$root
        A1=sum(ra^(1:m)*nstar*Z)
        rl1=rgamma(1,n1+a[1],A1+b[1])
        rl2=rgamma(1,n2+a[2],A1+b[2])
        return(c(ra,rl1,rl2))
      }
     return(colMeans(t(sapply(1:10000,SG,simplify = "array"))))
  }
########
BGI=matrix(NA,ncol=3,nrow=50000)
for(iter in 1:50000)
  {
    BGI[iter,]=GIBS(data,c(2,1.6),c(2,2),5,c(10,rep(0,9)))
    cat(iter, "of 50000\r") 
   flush.console()
 }

The reason it's not completing is that there are bugs in the code. Object lengths are not conformable in a number of instances. Since the function is named `GIBS` may I ask if this is supposed to be a Gibbs sampler? If so, you may find a pre-existing package/function to suit your needs. — Hack-R, Aug 05 '16 at 12:36
if you put a `print(iter)` in your for loop then you'd know how far along it is and by extension could guess how much longer it would take. As to the problem itself, R isn't very fast with loops and between the `sapply`, `while`, and `for` loops there's a lot in your code. Your best bet, assuming this is important enough to justify the work, would be to use Rcpp which would unfortunately require you to rewrite much of your code in cpp. — Dean MacGregor, Aug 05 '16 at 12:47
@K.Ahmadi: Concernings speed [this](http://stackoverflow.com/questions/2908822/speed-up-the-loop-operation-in-r) perhaps helps — Christoph, Aug 05 '16 at 14:24

Hack-R · Answer 1 · 2016-08-05T12:59:52.930

Regarding the updated version of the function I would suggest to leverage multicore/HPC functions to run it if you need that many iterations.

Below are the comments on the original version of the function:

The problem is with this part of your function (and its inputs, f and the data objects):

integrate(f,1,Inf)

Beginning with the line of f that says

FR=prod((A+b)^(N+a))

It gives the following warnings, which lead to subsequent problems and effectively an infinite loop:

Warning messages:
1: In x^(1:m) : longer object length is not a multiple of shorter object length
2: In x^(1:m) * nstar :
  longer object length is not a multiple of shorter object length
3: In x^(1:m) * nstar * Z :
  longer object length is not a multiple of shorter object length
4: In A + b : longer object length is not a multiple of shorter object length
5: In N + a : longer object length is not a multiple of shorter object length
6: In x^(1:m) : longer object length is not a multiple of shorter object length
7: In x^(1:m) * nstar :
  longer object length is not a multiple of shorter object length
8: In x^(1:m) * nstar * Z :
  longer object length is not a multiple of shorter object length
9: In A + b : longer object length is not a multiple of shorter object length
10: In N + a : longer object length is not a multiple of shorter object length
11: In x^(1:m) : longer object length is not a multiple of shorter object length
12: In x^(1:m) * nstar :
  longer object length is not a multiple of shorter object length
13: In x^(1:m) * nstar * Z :
  longer object length is not a multiple of shorter object length
14: In A + b : longer object length is not a multiple of shorter object length
15: In N + a : longer object length is not a multiple of shorter object length

I don't understand what this function is trying to accomplish so that I may help you fix it, but if you make that part more clear I will try.

@K.Ahmadi OK, great. If this solved your problem please mark it as the solution. If not, please update your question with the annotated changes. — Hack-R, Aug 05 '16 at 12:56

How can I speed up the following function in R?

1 Answers1