[R] Vectorizing integrate()

[arun]

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]]
> > >
> > > ______________________________________________
> > > R-help at 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.
> >
> > David Winsemius, MD
> > Alameda, CA, USA
> 
> ______________________________________________
> R-help at 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.