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