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))
---------------------------------------------------------------------------------------------------------------------------------
Thanks for your attention.
Kind regards,
Jamil.
[[alternative HTML version deleted]]
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
Hi
I am facing a problem of restricting an intercept of systems of equations.
Y1=f(X1,X2,X3)
Y2=(X1,X2,X4)
I want to restrict an intercept of equation 2 equal to coefficient of X2
of equation 1.
Please help
Dereje
________________________________
From: Naser Jamil <jamilnaser79@gmail.com>
To: r-help@r-project.org
Sent: Tuesday, October 2, 2012 10:23 AM
Subject: [R] Integration in R
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))
---------------------------------------------------------------------------------------------------------------------------------
Thanks for your attention.
Kind regards,
Jamil.
[[alternative HTML version deleted]]
______________________________________________
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.
[[alternative HTML version deleted]]
On 02-10-2012, at 20:50, Dereje Bacha <d_bacha at yahoo.com> wrote:> Hi > > I am facing a problem of restricting an intercept of systems of equations. > Y1=f(X1,X2,X3) > Y2=(X1,X2,X4) > > I want to restrict an intercept of equation 2 equal to coefficient of X2 of equation 1. >Please do not hijack a thread/conversation. Do not reply to a thread with a completely different subject. Start a new thread for a new subject. It is completely unclear what you want. Berend