[R] Integration in R

[Berend Hasselman]

On 02-10-2012, at 17:23, Naser Jamil <jamilnaser79 at gmail.com> wrote:

> Dear R-users,
> I am facing problem with integrating in R a likelihood function which is a
> function of four parameters. It's giving me the result at the end but
> taking more than half an hour to run. I'm wondering is there any other
> efficient way deal with. The following is my code. I am ready to provide
> any other description of my function if you need to move forward.
> 
> ------------------------------------------------------------------------------------------------------------------------------------
> 
> library(cubature)
> dose<-c(2,3,5)
> y0<-c(2,1,0)
> y1<-c(1,1,1)
> y2<-c(0,1,2)
> 
> lf<-function (x) {
>    v<-1
>    for (i in 1:length(dose)) {
>        psi0<-1/((1+exp(x[1]+x[2]*dose[i]))*(1+exp(x[3]+x[4]*dose[i])))
> 
> psi1<-exp(x[1]+x[2]*dose[i])/((1+exp(x[1]+x[2]*dose[i]))*(1+exp(x[3]+x[4]*dose[i])))
>       v<-v*(psi0^y0[i])*(psi1^y1[i])*((1-psi0-psi1)^y2[i])
>        }
>       return(v)
>        }
> 
> adaptIntegrate(lf, lowerLimit = c(-20, 0,-20,0), upperLimit = c(0,10,0,10))
> 

There are several things you could do.

1. Use the compiler package to make a compiled version of your function.
2. rewrite your function by putting x[1] etc. in separate variables x1, x2,... so avoiding the [..] indexing. You can do the same for dose[i].
And also compiling this version of the function.
3. do not recompute expressions such as exp(x1+x2*dose.i) several times. Store the result in a temporary variable and use that variable.

With these functions

library(compiler)

lf.c <- cmpfun(lf)

lf1 <-function (x) {
   v<-1 
   x1 <- x[1]
   x2 <- x[2]
   x3 <- x[3]
   x4 <- x[4]
   for (i in 1:length(dose)) {                     
       dose.i <- dose[i]
       z1 <- exp(x1+x2*dose.i)
       z2 <- exp(x3+x4*dose.i)
       psi0<-1/((1+z1)*(1+z2))
       psi1<-z1*psi0
       v<-v*(psi0^y0[i])*(psi1^y1[i])*((1-psi0-psi1)^y2[i])
   }
   return(v)
}

lf1.c <- cmpfun(lf1)

I tested relative speeds with this code (small tolerance and max. function evaluations)

library(rbenchmark)

f1 <- function() adaptIntegrate(lf   , lowerLimit = c(-20, 0,-20,0), upperLimit = c(0,10,0,10), tol=1e-3,maxEval=50000)
f2 <- function() adaptIntegrate(lf.c , lowerLimit = c(-20, 0,-20,0), upperLimit = c(0,10,0,10), tol=1e-3,maxEval=50000)
f3 <- function() adaptIntegrate(lf1  , lowerLimit = c(-20, 0,-20,0), upperLimit = c(0,10,0,10), tol=1e-3,maxEval=50000)
f4 <- function() adaptIntegrate(lf1.c, lowerLimit = c(-20, 0,-20,0), upperLimit = c(0,10,0,10), tol=1e-3,maxEval=50000)
benchmark(z1 <- f1(),z2 <- f2(), z3 <- f3(),z4 <- f4(),replications=1)

Result:

> benchmark(z1 <- f1(),z2 <- f2(), z3 <- f3(),z4 <- f4(),replications=1)
        test replications elapsed relative user.self sys.self user.child
1 z1 <- f1()            1   3.197    4.535     3.177    0.008          0
2 z2 <- f2()            1   1.240    1.759     1.235    0.003          0
3 z3 <- f3()            1   2.171    3.079     2.167    0.002          0
4 z4 <- f4()            1   0.705    1.000     0.700    0.004          0

As you can see you should be able to get at least a fourfold speedup by using the compiled version of the optimized function.
I would certainly  set tol and maxEval to something reasonable initially.
Checking that z1, z2, z3, and z4 are equal is left to you.

Finally, it may even be possible to eliminate the for loop in your function but I'll leave that for someone else.

Berend