[R] Integration in R

[Rui Barradas]

Hello,

Yes, it's possible to remove the loop. Since the loop is used to compute 
a running product and all we want is the final result, use the 
vectorized behavior of R and a final ?prod().
Seedup: another 2x. And 4x2 == 8 == 1 [decimal] order of magnitude.


lf2 <-function (x) {
    v<-1
    x1 <- x[1]
    x2 <- x[2]
    x3 <- x[3]
    x4 <- x[4]
    z1 <- exp(x1+x2*dose)
    z2 <- exp(x3+x4*dose)
    psi0<-1/((1+z1)*(1+z2))
    psi1<-z1*psi0
    v <- (psi0^y0)*(psi1^y1)*((1-psi0-psi1)^y2)
    return( prod(v) )
}

lf2.c <- cmpfun(lf2)

Hope this helps,

Rui Barradas
Em 02-10-2012 18:21, Berend Hasselman escreveu:
> 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
>
> ______________________________________________
> 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.