[R] Vectorizing integrate()

[Spencer Graves]

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


       If that has been suggested, please excuse me for entering this 
thread late without reading all the previous posts.


       Spencer


On 12/7/2012 8:54 AM, David L Carlson wrote:
> 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.


-- 
Spencer Graves, PE, PhD
President and Chief Technology Officer
Structure Inspection and Monitoring, Inc.
751 Emerson Ct.
San José, CA 95126
ph:  408-655-4567
web:  www.structuremonitoring.com