Vectorizing integrate()

I have written a program to solve a particular logistic regression problem using IRLS. In one step, I need to integrate something out of the linear predictor. The way I'm doing it now is within a loop and it is as you would expect slow to process, especially inside an iterative algorithm.

I'm hoping there is a way this can be vectorized, but I have not found it so far. The portion of code I'd like to vectorize is this

for(j in 1:nrow(X)){
 fun <- function(u) 1/ (1 + exp(- (B[1] + B[2] * (x[j] + u)))) * dnorm(u, 0, sd[j])
               eta[j] <- integrate(fun, -Inf, Inf)$value
}

Here X is an n x p model matrix for the fixed effects, B is a vector with the estimates of the fixed effects at iteration t, x is a predictor variable in the jth row of X, and sd is a variable corresponding to x[j].

Is there a way this can be done without looping over the rows of X?

Thanks,
Harold

	[[alternative HTML version deleted]]

v <- 1:5
w <- rep(1, 5)
fun <- function(x, m, n) n*exp(-m*x)
eta <- sapply(1:5, function(i) integrate(fun, 0, Inf, m=v[i],

eta

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-
project.org] On Behalf Of David Winsemius
Sent: Thursday, December 06, 2012 12:59 PM
To: Doran, Harold
Cc: r-help at r-project.org
Subject: Re: [R] Vectorizing integrate()


On Dec 6, 2012, at 10:10 AM, Doran, Harold wrote:

I have written a program to solve a particular logistic regression

problem using IRLS. In one step, I need to integrate something out of
the linear predictor. The way I'm doing it now is within a loop and it
is as you would expect slow to process, especially inside an iterative
algorithm.

I'm hoping there is a way this can be vectorized, but I have not

found it so far. The portion of code I'd like to vectorize is this

for(j in 1:nrow(X)){
 fun <- function(u) 1/ (1 + exp(- (B[1] + B[2] * (x[j] + u)))) *

dnorm(u, 0, sd[j])

               eta[j] <- integrate(fun, -Inf, Inf)$value
}

The Vectorize function is just a wrapper to mapply. If yoou are able to
get that code to execute properly for your un-posted test cases, then
why not use mapply?

Here X is an n x p model matrix for the fixed effects, B is a vector

with the estimates of the fixed effects at iteration t, x is a
predictor variable in the jth row of X, and sd is a variable
corresponding to x[j].

Is there a way this can be done without looping over the rows of X?

Thanks,
Harold

	[[alternative HTML version deleted]]

guide.html

David Winsemius, MD
Alameda, CA, USA

-----Original Message-----
From: David Winsemius [mailto:dwinsemius at comcast.net]
Sent: Thursday, December 06, 2012 1:59 PM
To: Doran, Harold
Cc: r-help at r-project.org
Subject: Re: [R] Vectorizing integrate()


On Dec 6, 2012, at 10:10 AM, Doran, Harold wrote:

I have written a program to solve a particular logistic regression problem

using IRLS. In one step, I need to integrate something out of the linear
predictor. The way I'm doing it now is within a loop and it is as you would
expect slow to process, especially inside an iterative algorithm.

I'm hoping there is a way this can be vectorized, but I have not found
it so far. The portion of code I'd like to vectorize is this

for(j in 1:nrow(X)){
 fun <- function(u) 1/ (1 + exp(- (B[1] + B[2] * (x[j] + u)))) * dnorm(u, 0,

sd[j])

               eta[j] <- integrate(fun, -Inf, Inf)$value }

The Vectorize function is just a wrapper to mapply. If yoou are able to get
that code to execute properly for your un-posted test cases, then why not
use mapply?

Here X is an n x p model matrix for the fixed effects, B is a vector with the

estimates of the fixed effects at iteration t, x is a predictor variable in the jth
row of X, and sd is a variable corresponding to x[j].

Is there a way this can be done without looping over the rows of X?

Thanks,
Harold

	[[alternative HTML version deleted]]

David Winsemius, MD
Alameda, CA, USA

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

eta <- sapply(1:nrow(X), function(i) integrate(fun, 

eta

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-
project.org] On Behalf Of Doran, Harold
Sent: Friday, December 07, 2012 9:14 AM
To: David Winsemius
Cc: r-help at r-project.org
Subject: Re: [R] Vectorizing integrate()

David et al

Thanks, I should have made the post more complete. I routinely use
apply functions, but often avoid mapply() as I find it so non-
intuitive. In this instance, I think the situation demands I change
that position, though.

Reproducible code for the current implementation of the function is

B <- c(0,1)
sem1 = runif(10, 1, 2)
x <- rnorm(10)
X <- cbind(1, x)
eta <- numeric(10)

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

I can't get my head around how mapply() would work here. It accepts as
its first argument a function. But, in my case I have two functions:
the user defined integrand, fun(), an then of course calling the R
function integrate().

I was thinking maybe along these lines, but this is obviously wrong.

mapply(integrate(function(u) 1/ (1 + exp(- (B[1] + B[2] * (x + u)))) *
dnorm(u, 0, sem1), -Inf, Inf)$value, MoreArgs = list(B, x, sem1))

-----Original Message-----
From: David Winsemius [mailto:dwinsemius at comcast.net]
Sent: Thursday, December 06, 2012 1:59 PM
To: Doran, Harold
Cc: r-help at r-project.org
Subject: Re: [R] Vectorizing integrate()


On Dec 6, 2012, at 10:10 AM, Doran, Harold wrote:

I have written a program to solve a particular logistic regression

problem

using IRLS. In one step, I need to integrate something out of the

linear

predictor. The way I'm doing it now is within a loop and it is as you

would

expect slow to process, especially inside an iterative algorithm.

I'm hoping there is a way this can be vectorized, but I have not

found

it so far. The portion of code I'd like to vectorize is this

for(j in 1:nrow(X)){
 fun <- function(u) 1/ (1 + exp(- (B[1] + B[2] * (x[j] + u)))) *

dnorm(u, 0,

sd[j])

               eta[j] <- integrate(fun, -Inf, Inf)$value }

The Vectorize function is just a wrapper to mapply. If yoou are able

to get

that code to execute properly for your un-posted test cases, then why

not

use mapply?

Here X is an n x p model matrix for the fixed effects, B is a

vector with the

estimates of the fixed effects at iteration t, x is a predictor

variable in the jth

row of X, and sd is a variable corresponding to x[j].

Is there a way this can be done without looping over the rows of X?

Thanks,
Harold

	[[alternative HTML version deleted]]

David Winsemius, MD
Alameda, CA, USA

-----Original Message-----
From: David Winsemius [mailto:dwinsemius at comcast.net]
Sent: Thursday, December 06, 2012 1:59 PM
To: Doran, Harold
Cc: r-help at r-project.org
Subject: Re: [R] Vectorizing integrate()


On Dec 6, 2012, at 10:10 AM, Doran, Harold wrote:

I have written a program to solve a particular logistic regression problem

using IRLS. In one step, I need to integrate something out of the linear
predictor. The way I'm doing it now is within a loop and it is as you would
expect slow to process, especially inside an iterative algorithm.

I'm hoping there is a way this can be vectorized, but I have not found
it so far. The portion of code I'd like to vectorize is this

for(j in 1:nrow(X)){
? fun <- function(u) 1/ (1 + exp(- (B[1] + B[2] * (x[j] + u)))) * dnorm(u, 0,

sd[j])

? ? ? ? ? ? ? ? eta[j] <- integrate(fun, -Inf, Inf)$value }

The Vectorize function is just a wrapper to mapply. If yoou are able to get
that code to execute properly for your un-posted test cases, then why not
use mapply?

Here X is an n x p model matrix for the fixed effects, B is a vector with the

estimates of the fixed effects at iteration t, x is a predictor variable in the jth
row of X, and sd is a variable corresponding to x[j].

Is there a way this can be done without looping over the rows of X?

Thanks,
Harold

??? [[alternative HTML version deleted]]

David Winsemius, MD
Alameda, CA, USA

B <- c(0,1)
sem1 = runif(10, 1, 2)
x <- rnorm(10)
X <- cbind(1, x)
eta <- numeric(10)
for(j in 1:nrow(X)){

eta

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

eta2 <- sapply(1:nrow(X), function(i) integrate(fun, 

eta2

identical(eta, eta2)

res<-mapply(function(i)

res

identical(eta, res)

-----Original Message-----
From: arun [mailto:smartpink111 at yahoo.com]
Sent: Friday, December 07, 2012 10:36 AM
To: Doran, Harold
Cc: R help; David L Carlson; David Winsemius
Subject: Re: [R] Vectorizing integrate()



Hi,
Using David's function:
fun <- function(u, m, s) 1/ (1 + exp(- (B[1] + B[2] *
????? (m + u)))) * dnorm(u, 0, s)

?res<-mapply(function(i) integrate(fun,-
Inf,Inf,m=x[i],s=sem1[i])$value,1:nrow(X))
res
# [1] 0.5212356 0.6214989 0.5306124 0.5789414 0.3429795 0.6972879
0.5952949
?#[8] 0.7531899 0.4740685 0.7576587
identical(res,eta)
#[1] TRUE
A.K.

----- Original Message -----
From: "Doran, Harold" <HDoran at air.org>
To: David Winsemius <dwinsemius at comcast.net>
Cc: "r-help at r-project.org" <r-help at r-project.org>
Sent: Friday, December 7, 2012 10:14 AM
Subject: Re: [R] Vectorizing integrate()

David et al

Thanks, I should have made the post more complete. I routinely use
apply functions, but often avoid mapply() as I find it so non-
intuitive. In this instance, I think the situation demands I change
that position, though.

Reproducible code for the current implementation of the function is

B <- c(0,1)
sem1 = runif(10, 1, 2)
x <- rnorm(10)
X <- cbind(1, x)
eta <- numeric(10)

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

I can't get my head around how mapply() would work here. It accepts as
its first argument a function. But, in my case I have two functions:
the user defined integrand, fun(), an then of course calling the R
function integrate().

I was thinking maybe along these lines, but this is obviously wrong.

mapply(integrate(function(u) 1/ (1 + exp(- (B[1] + B[2] * (x + u)))) *
dnorm(u, 0, sem1), -Inf, Inf)$value, MoreArgs = list(B, x, sem1))

-----Original Message-----
From: David Winsemius [mailto:dwinsemius at comcast.net]
Sent: Thursday, December 06, 2012 1:59 PM
To: Doran, Harold
Cc: r-help at r-project.org
Subject: Re: [R] Vectorizing integrate()


On Dec 6, 2012, at 10:10 AM, Doran, Harold wrote:

I have written a program to solve a particular logistic regression

problem

using IRLS. In one step, I need to integrate something out of the

linear

predictor. The way I'm doing it now is within a loop and it is as you

would

expect slow to process, especially inside an iterative algorithm.

I'm hoping there is a way this can be vectorized, but I have not

found

it so far. The portion of code I'd like to vectorize is this

for(j in 1:nrow(X)){
? fun <- function(u) 1/ (1 + exp(- (B[1] + B[2] * (x[j] + u)))) *

dnorm(u, 0,

sd[j])

? ? ? ? ? ? ? ? eta[j] <- integrate(fun, -Inf, Inf)$value }

The Vectorize function is just a wrapper to mapply. If yoou are able

to get

that code to execute properly for your un-posted test cases, then why

not

use mapply?

Here X is an n x p model matrix for the fixed effects, B is a

vector with the

estimates of the fixed effects at iteration t, x is a predictor

variable in the jth

row of X, and sd is a variable corresponding to x[j].

Is there a way this can be done without looping over the rows of X?

Thanks,
Harold

??? [[alternative HTML version deleted]]

David Winsemius, MD
Alameda, CA, USA

B <- c(0,1)
sem1 = runif(10, 1, 2)
x <- rnorm(10)
X <- cbind(1, x)
eta <- numeric(10)
for(j in 1:nrow(X)){

eta

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

eta2 <- sapply(1:nrow(X), function(i) integrate(fun, 

eta2

identical(eta, eta2)

res<-mapply(function(i)

res

identical(eta, res)

-----Original Message-----
From: arun [mailto:smartpink111 at yahoo.com]
Sent: Friday, December 07, 2012 10:36 AM
To: Doran, Harold
Cc: R help; David L Carlson; David Winsemius
Subject: Re: [R] Vectorizing integrate()



Hi,
Using David's function:
fun <- function(u, m, s) 1/ (1 + exp(- (B[1] + B[2] *
????? (m + u)))) * dnorm(u, 0, s)

?res<-mapply(function(i) integrate(fun,-
Inf,Inf,m=x[i],s=sem1[i])$value,1:nrow(X))
res
# [1] 0.5212356 0.6214989 0.5306124 0.5789414 0.3429795 0.6972879
0.5952949
?#[8] 0.7531899 0.4740685 0.7576587
identical(res,eta)
#[1] TRUE
A.K.

----- Original Message -----
From: "Doran, Harold" <HDoran at air.org>
To: David Winsemius <dwinsemius at comcast.net>
Cc: "r-help at r-project.org" <r-help at r-project.org>
Sent: Friday, December 7, 2012 10:14 AM
Subject: Re: [R] Vectorizing integrate()

David et al

Thanks, I should have made the post more complete. I routinely use
apply functions, but often avoid mapply() as I find it so non-
intuitive. In this instance, I think the situation demands I change
that position, though.

Reproducible code for the current implementation of the function is

B <- c(0,1)
sem1 = runif(10, 1, 2)
x <- rnorm(10)
X <- cbind(1, x)
eta <- numeric(10)

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

I can't get my head around how mapply() would work here. It accepts as
its first argument a function. But, in my case I have two functions:
the user defined integrand, fun(), an then of course calling the R
function integrate().

I was thinking maybe along these lines, but this is obviously wrong.

mapply(integrate(function(u) 1/ (1 + exp(- (B[1] + B[2] * (x + u)))) *
dnorm(u, 0, sem1), -Inf, Inf)$value, MoreArgs = list(B, x, sem1))

-----Original Message-----
From: David Winsemius [mailto:dwinsemius at comcast.net]
Sent: Thursday, December 06, 2012 1:59 PM
To: Doran, Harold
Cc: r-help at r-project.org
Subject: Re: [R] Vectorizing integrate()


On Dec 6, 2012, at 10:10 AM, Doran, Harold wrote:

I have written a program to solve a particular logistic regression

problem

using IRLS. In one step, I need to integrate something out of the

linear

predictor. The way I'm doing it now is within a loop and it is as you

would

expect slow to process, especially inside an iterative algorithm.

I'm hoping there is a way this can be vectorized, but I have not

found

it so far. The portion of code I'd like to vectorize is this

for(j in 1:nrow(X)){
? fun <- function(u) 1/ (1 + exp(- (B[1] + B[2] * (x[j] + u)))) *

dnorm(u, 0,

sd[j])

? ? ? ? ? ? ? ? eta[j] <- integrate(fun, -Inf, Inf)$value }

The Vectorize function is just a wrapper to mapply. If yoou are able

to get

that code to execute properly for your un-posted test cases, then why

not

use mapply?

Here X is an n x p model matrix for the fixed effects, B is a

vector with the

estimates of the fixed effects at iteration t, x is a predictor

variable in the jth

row of X, and sd is a variable corresponding to x[j].

Is there a way this can be done without looping over the rows of X?

Thanks,
Harold

??? [[alternative HTML version deleted]]

David Winsemius, MD
Alameda, CA, USA

David et al

Thanks, I should have made the post more complete. I routinely use  
apply functions, but often avoid mapply() as I find it so non- 
intuitive. In this instance, I think the situation demands I change  
that position, though.

Reproducible code for the current implementation of the function is

B <- c(0,1)
sem1 = runif(10, 1, 2)
x <- rnorm(10)
X <- cbind(1, x)
eta <- numeric(10)

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

I can't get my head around how mapply() would work here. It accepts  
as its first argument a function. But, in my case I have two  
functions: the user defined integrand, fun(), an then of course  
calling the R function integrate().

I was thinking maybe along these lines, but this is obviously wrong.

mapply(integrate(function(u) 1/ (1 + exp(- (B[1] + B[2] * (x + u))))  
* dnorm(u, 0, sem1), -Inf, Inf)$value, MoreArgs = list(B, x, sem1))

Must be a cut and paste issue. All three agree on the results but they are
different from those in arun's message:

B <- c(0,1)
sem1 = runif(10, 1, 2)
x <- rnorm(10)
X <- cbind(1, x)
eta <- numeric(10)
for(j in 1:nrow(X)){

+ fun <- function(u) 1/ (1 + exp(- (B[1] + B[2] * (x[j] + u)))) * dnorm(u,
0, sem1[j])
+ eta[j] <- integrate(fun, -Inf, Inf)$value
+ }

eta

  [1] 0.6586459 0.5565611 0.5226462 0.6824012 0.6908079 0.7578524 0.4715976
  [8] 0.2165210 0.2378657 0.3492133

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

+       (m + u)))) * dnorm(u, 0, s)

eta2 <- sapply(1:nrow(X), function(i) integrate(fun,

+       -Inf, Inf, m=x[i], s=sem1[i])$value)

eta2

  [1] 0.6586459 0.5565611 0.5226462 0.6824012 0.6908079 0.7578524 0.4715976
  [8] 0.2165210 0.2378657 0.3492133

identical(eta, eta2)

[1] TRUE

res<-mapply(function(i)

integrate(fun,-Inf,Inf,m=x[i],s=sem1[i])$value,1:nrow(X))

res

  [1] 0.6586459 0.5565611 0.5226462 0.6824012 0.6908079 0.7578524 0.4715976
  [8] 0.2165210 0.2378657 0.3492133

identical(eta, res)

[1] TRUE

-------
David C

-----Original Message-----
From: arun [mailto:smartpink111 at yahoo.com]
Sent: Friday, December 07, 2012 10:36 AM
To: Doran, Harold
Cc: R help; David L Carlson; David Winsemius
Subject: Re: [R] Vectorizing integrate()



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

  res<-mapply(function(i) integrate(fun,-
Inf,Inf,m=x[i],s=sem1[i])$value,1:nrow(X))
res
# [1] 0.5212356 0.6214989 0.5306124 0.5789414 0.3429795 0.6972879
0.5952949
  #[8] 0.7531899 0.4740685 0.7576587
identical(res,eta)
#[1] TRUE
A.K.

----- Original Message -----
From: "Doran, Harold" <HDoran at air.org>
To: David Winsemius <dwinsemius at comcast.net>
Cc: "r-help at r-project.org" <r-help at r-project.org>
Sent: Friday, December 7, 2012 10:14 AM
Subject: Re: [R] Vectorizing integrate()

David et al

Thanks, I should have made the post more complete. I routinely use
apply functions, but often avoid mapply() as I find it so non-
intuitive. In this instance, I think the situation demands I change
that position, though.

Reproducible code for the current implementation of the function is

B <- c(0,1)
sem1 = runif(10, 1, 2)
x <- rnorm(10)
X <- cbind(1, x)
eta <- numeric(10)

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

I can't get my head around how mapply() would work here. It accepts as
its first argument a function. But, in my case I have two functions:
the user defined integrand, fun(), an then of course calling the R
function integrate().

I was thinking maybe along these lines, but this is obviously wrong.

mapply(integrate(function(u) 1/ (1 + exp(- (B[1] + B[2] * (x + u)))) *
dnorm(u, 0, sem1), -Inf, Inf)$value, MoreArgs = list(B, x, sem1))

-----Original Message-----
From: David Winsemius [mailto:dwinsemius at comcast.net]
Sent: Thursday, December 06, 2012 1:59 PM
To: Doran, Harold
Cc: r-help at r-project.org
Subject: Re: [R] Vectorizing integrate()


On Dec 6, 2012, at 10:10 AM, Doran, Harold wrote:

I have written a program to solve a particular logistic regression

problem

using IRLS. In one step, I need to integrate something out of the

linear

predictor. The way I'm doing it now is within a loop and it is as you

would

expect slow to process, especially inside an iterative algorithm.

I'm hoping there is a way this can be vectorized, but I have not

found

it so far. The portion of code I'd like to vectorize is this

for(j in 1:nrow(X)){
   fun <- function(u) 1/ (1 + exp(- (B[1] + B[2] * (x[j] + u)))) *

dnorm(u, 0,

sd[j])

                 eta[j] <- integrate(fun, -Inf, Inf)$value }

The Vectorize function is just a wrapper to mapply. If yoou are able

to get

that code to execute properly for your un-posted test cases, then why

not

use mapply?

Here X is an n x p model matrix for the fixed effects, B is a

vector with the

estimates of the fixed effects at iteration t, x is a predictor

variable in the jth

row of X, and sd is a variable corresponding to x[j].

Is there a way this can be done without looping over the rows of X?

Thanks,
Harold

     [[alternative HTML version deleted]]

David Winsemius, MD
Alameda, CA, USA

     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".

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),

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

On 12/7/2012 9:40 AM, Berend Hasselman wrote:

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.

     So the compiler (f2, f4, f6) provided a slight improvement over f1 and f3 but not f2, and in any event, the improvement was not great.

On 07-12-2012, at 19:37, Spencer Graves wrote:

On 12/7/2012 9:40 AM, Berend Hasselman wrote:

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.

      So the compiler (f2, f4, f6) provided a slight improvement over f1 and f3 but not f2, and in any event, the improvement was not great.

I don't understand the "but not f2".
And I don't understand the conclusion for (f2,f4,f6). f4 is a compiled version of f3 and is slower than its non compiled version.
f2 and f6 are the quickest compiled versions.
Indeed the improvement is not earth shattering but it does demonstrate what you can achieve by using the compiler package.

Berend

HI,

I can see ?runif(), ?rnorm() without a set seed.

A.K.

----- Original Message -----
From: David L Carlson <dcarlson at tamu.edu>
To: 'arun' <smartpink111 at yahoo.com>; "'Doran, Harold'" <HDoran at air.org>
Cc: 'R help' <r-help at r-project.org>; 'David Winsemius' <dwinsemius at comcast.net>
Sent: Friday, December 7, 2012 11:54 AM
Subject: RE: [R] Vectorizing integrate()

Must be a cut and paste issue. All three agree on the results but they are
different from those in arun's message:

B <- c(0,1)
sem1 = runif(10, 1, 2)
x <- rnorm(10)
X <- cbind(1, x)
eta <- numeric(10)
for(j in 1:nrow(X)){

+ fun <- function(u) 1/ (1 + exp(- (B[1] + B[2] * (x[j] + u)))) * dnorm(u,
0, sem1[j])
+ eta[j] <- integrate(fun, -Inf, Inf)$value
+ }

eta

[1] 0.6586459 0.5565611 0.5226462 0.6824012 0.6908079 0.7578524 0.4715976
[8] 0.2165210 0.2378657 0.3492133

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

+       (m + u)))) * dnorm(u, 0, s)

eta2 <- sapply(1:nrow(X), function(i) integrate(fun, 

+       -Inf, Inf, m=x[i], s=sem1[i])$value)

eta2

[1] 0.6586459 0.5565611 0.5226462 0.6824012 0.6908079 0.7578524 0.4715976
[8] 0.2165210 0.2378657 0.3492133

identical(eta, eta2)

[1] TRUE

res<-mapply(function(i)

integrate(fun,-Inf,Inf,m=x[i],s=sem1[i])$value,1:nrow(X))

res

[1] 0.6586459 0.5565611 0.5226462 0.6824012 0.6908079 0.7578524 0.4715976
[8] 0.2165210 0.2378657 0.3492133

identical(eta, res)

[1] TRUE

-------
David C

-----Original Message-----
From: arun [mailto:smartpink111 at yahoo.com]
Sent: Friday, December 07, 2012 10:36 AM
To: Doran, Harold
Cc: R help; David L Carlson; David Winsemius
Subject: Re: [R] Vectorizing integrate()



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

 res<-mapply(function(i) integrate(fun,-
Inf,Inf,m=x[i],s=sem1[i])$value,1:nrow(X))
res
# [1] 0.5212356 0.6214989 0.5306124 0.5789414 0.3429795 0.6972879
0.5952949
 #[8] 0.7531899 0.4740685 0.7576587
identical(res,eta)
#[1] TRUE
A.K.

----- Original Message -----
From: "Doran, Harold" <HDoran at air.org>
To: David Winsemius <dwinsemius at comcast.net>
Cc: "r-help at r-project.org" <r-help at r-project.org>
Sent: Friday, December 7, 2012 10:14 AM
Subject: Re: [R] Vectorizing integrate()

David et al

Thanks, I should have made the post more complete. I routinely use
apply functions, but often avoid mapply() as I find it so non-
intuitive. In this instance, I think the situation demands I change
that position, though.

Reproducible code for the current implementation of the function is

B <- c(0,1)
sem1 = runif(10, 1, 2)
x <- rnorm(10)
X <- cbind(1, x)
eta <- numeric(10)

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

I can't get my head around how mapply() would work here. It accepts as
its first argument a function. But, in my case I have two functions:
the user defined integrand, fun(), an then of course calling the R
function integrate().

I was thinking maybe along these lines, but this is obviously wrong.

mapply(integrate(function(u) 1/ (1 + exp(- (B[1] + B[2] * (x + u)))) *
dnorm(u, 0, sem1), -Inf, Inf)$value, MoreArgs = list(B, x, sem1))

-----Original Message-----
From: David Winsemius [mailto:dwinsemius at comcast.net]
Sent: Thursday, December 06, 2012 1:59 PM
To: Doran, Harold
Cc: r-help at r-project.org
Subject: Re: [R] Vectorizing integrate()


On Dec 6, 2012, at 10:10 AM, Doran, Harold wrote:

I have written a program to solve a particular logistic regression

problem

using IRLS. In one step, I need to integrate something out of the

linear

predictor. The way I'm doing it now is within a loop and it is as you

would

expect slow to process, especially inside an iterative algorithm.

I'm hoping there is a way this can be vectorized, but I have not

found

it so far. The portion of code I'd like to vectorize is this

for(j in 1:nrow(X)){
  fun <- function(u) 1/ (1 + exp(- (B[1] + B[2] * (x[j] + u)))) *

dnorm(u, 0,

sd[j])

                eta[j] <- integrate(fun, -Inf, Inf)$value }

The Vectorize function is just a wrapper to mapply. If yoou are able

to get

that code to execute properly for your un-posted test cases, then why

not

use mapply?

Here X is an n x p model matrix for the fixed effects, B is a

vector with the

estimates of the fixed effects at iteration t, x is a predictor

variable in the jth

row of X, and sd is a variable corresponding to x[j].

Is there a way this can be done without looping over the rows of X?

Thanks,
Harold

    [[alternative HTML version deleted]]

David Winsemius, MD
Alameda, CA, USA