[R] Vectorizing integrate()

[Berend Hasselman]

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

>      Has anyone suggested using the byte code compiler "compiler" package?  An analysis by John Nash suggested to me that it may be roughly equivalent to vectorization;  see "http://rwiki.sciviews.org/doku.php?id=tips:rqcasestudy&s=compiler".
> 
> 

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

> identical(eta1,eta2)
[1] TRUE
> identical(eta1,eta3)
[1] TRUE
> identical(eta1,eta4)
[1] TRUE
> identical(eta1,eta5)
[1] TRUE
> 
> 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"))
                       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