[R] Vectorizing integrate()

[David Winsemius]

On Dec 7, 2012, at 7:14 AM, Doran, Harold wrote:

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

I had been thinking (before you presented the data case and allowed  
the problem to be seen more fully)  that X[ n, ] was being passed to  
that expression. It's not. You are using the index of one object to  
pass positions of a vector, so there is really only one value. Using  
mapply would reduce to using sapply as Carlson illustrated.

-- 
David
>
>
>> -----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
>

David Winsemius, MD
Alameda, CA, USA