On 02-10-2012, at 17:23, Naser Jamil <jamilnase...@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@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.
