Skip to content
Prev 313205 / 398506 Next

Vectorizing integrate()

On 07-12-2012, at 18:12, Spencer Graves wrote:

            
Not yet.
But here are some results for alternative ways of doing what  the OP wanted.


# Initial parameters

N <- 1000
B <- c(0,1)
sem1 <- runif(N, 1, 2)
x <- rnorm(N)
X <- cbind(1, x)

# load compiler package

library(compiler)

# functions

# Original loop solution with function fun defined outside loop

f1 <- function(X, B, x, sem1) {
    eta <- numeric(nrow(X))
    fun <- function(u, m, s) 1/ (1 + exp(- (B[1] + B[2] * (m + u)))) * dnorm(u, 0, s)
    for(j in 1:nrow(X)){
    	eta[j] <- integrate(fun, -Inf, Inf, m=x[j], s=sem1[j])$value
    }
    eta
}
f2 <- cmpfun(f1)

# sapply solution with fun defined outside function

fun <- function(u, m, s) 1/ (1 + exp(- (B[1] + B[2] * (m + u)))) * dnorm(u, 0, s)

f3 <- function(X, B, x, sem1) sapply(1:nrow(X), function(i) integrate(fun, -Inf, Inf, m=x[i], s=sem1[i])$value)
f4 <- cmpfun(f3)

# sapply solution with fun defined within function

f5 <- function(X, B, x, sem1) {
    fun <- function(u, m, s) 1/ (1 + exp(- (B[1] + B[2] * (m + u)))) * dnorm(u, 0, s)
    sapply(1:nrow(X), function(i) integrate(fun, -Inf, Inf, m=x[i], s=sem1[i])$value)
}
f6 <- cmpfun(f5)

# Testing

eta1 <- f1(X, B, x, sem1)
eta2 <- f2(X, B, x, sem1)
eta3 <- f3(X, B, x, sem1)
eta4 <- f4(X, B, x, sem1)
eta5 <- f5(X, B, x, sem1)
eta6 <- f6(X, B, x, sem1)

identical(eta1,eta2)
identical(eta1,eta3)
identical(eta1,eta4)
identical(eta1,eta5)

library(rbenchmark)

benchmark(eta1 <- f1(X, B, x, sem1), eta2 <- f2(X, B, x, sem1), eta3 <- f3(X, B, x, sem1),
          eta4 <- f4(X, B, x, sem1), eta5 <- f5(X, B, x, sem1), eta6 <- f6(X, B, x, sem1),
          replications=10, columns=c("test","elapsed","relative"))


# Results
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
+           eta4 <- f4(X, B, x, sem1), eta5 <- f5(X, B, x, sem1), eta6 <- f6(X, B, x, sem1),
+           replications=10, columns=c("test","elapsed","relative"))
                       test elapsed relative
1 eta1 <- f1(X, B, x, sem1)   1.873    1.207
2 eta2 <- f2(X, B, x, sem1)   1.552    1.000
3 eta3 <- f3(X, B, x, sem1)   1.807    1.164
4 eta4 <- f4(X, B, x, sem1)   1.841    1.186
5 eta5 <- f5(X, B, x, sem1)   1.852    1.193
6 eta6 <- f6(X, B, x, sem1)   1.601    1.032

As you can see using the compiler package is beneficial speedwise.
f2 and f6, both the the result of using the compiler package, are the quickest.
It's quite likely that more can be eked out of this.

Berend