[R] Vectorizing integrate()

Berend Hasselman bhh at xs4all.nl
Fri Dec 7 20:31:55 CET 2012 

On 07-12-2012, at 18:01, arun wrote:

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

Good eyesight! Correct.
Oops.

Inserting set.seed(413) at the start I now get this

# > 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)   2.248    1.392
# 2 eta2 <- f2(X, B, x, sem1)   1.738    1.076
# 3 eta3 <- f3(X, B, x, sem1)   2.024    1.253
# 4 eta4 <- f4(X, B, x, sem1)   2.072    1.283
# 5 eta5 <- f5(X, B, x, sem1)   2.038    1.262
# 6 eta6 <- f6(X, B, x, sem1)   1.615    1.000

Relative conclusions stay the same. The compiled versions f2, f6 are the quickest.

Berend
 

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

More information about the R-help mailing list