[R] Integrate inside function

[R. Michael Weylandt]

I don't believe you did anything wrong here, just that integrate()
itself isn't (can't?) be vectorized.

E.g.,

x <- integrate(function(x) x^2, 0, 2)
y <- integrate(function(x) x^2, 0, c(2,3))

identical(x, y) # FALSE -- so there's a difference.

# Take a look
str(x)
str(y)

# It's just in the call so let's drop that

x$call <- y$call <- NULL

identical(x,y) # TRUE

So at some point you'll need to use Vectorize() if you want a
vectorized function, at at least one that looks vectorized.
(Vectorize() isn't actually magic: it's just a loop so it will be a
little slow on large vectors)

Michael

On Sat, Mar 17, 2012 at 8:46 AM, Mauro Rossi <mauro.rossi at irpi.cnr.it> wrote:
> Dear David and Micheal,
>    thanks for your suggestion. Vectorize does what I need.
> David you suggest me that I didn't built the function in a manner that would
> vectorize. Could you please explain me what's wrong?
> Thanks,
> Mauro
>
>
> On Mar 15, 2012, at 6:08 AM, Mauro Rossi wrote:
>
> Dear R users,
>    first I take this opportunity to greet all the R community for your
> continuous efforts.
> I wrote a function to calculate the pdf and cdf of a custom distribution
> (mixed gamma model).
> The function is the following:
>
> pmixedgamma3 <- function(y, shape1, rate1, shape2, rate2, prev)
>     {
>     density.function<-function(x)
>         {
>         shape1=shape1
>         rate1=rate1
>         shape2=shape2
>         rate2=rate2
>         prev=prev
>
> fun<-prev*((rate1^shape1)/gamma(shape1)*(x^(shape1-1))*exp(-rate1*x)) +
> (1-prev)*((rate2^shape2)/gamma(shape2)*(x^(shape2-1))*exp(-rate2*x))
>         return(fun)
>         }
>     den<-density.function(y)
>     p.mixedgamma<-integrate(f=density.function,lower=0,upper=y)
>     return(list(density=den,cumdensity=p.mixedgamma$value))
>     }
>
> Why doesn't this calculate cumdensity (p.mixedgamma) in case x is a vector?
>
>
> You did not build the function in a manner that would vectorize, but perhaps
> the convenience vuntion would help here:
>
> Vpmixedgamma3 <- Vectorize(pmixedgamma3)
>
>
>
> For instance try:
> pmixedgamma3(c(10,20),shape1=10,rate1=1,shape2=10,rate2=1,prev=0.5)
>
>
>
>> Vpmixedgamma3(c(10,20),shape1=10,rate1=1,shape2=10,rate2=1,prev=0.5)
>            [,1]      [,2]
> density    0.12511   0.002908153
> cumdensity 0.5420703 0.9950046
>
>
> As you can see cumdensity is calculated just for the first x value.
>
>
>
>
> --
>
> Mauro Rossi
>
> Istituto di Ricerca per la Protezione Idrogeologica
>
> Consiglio Nazionale delle Ricerche
>
> Via della Madonna Alta, 126
>
> 06128 Perugia
>
> Italia
>
> Tel. +39 075 5014421
>
> Fax +39 075 5014420
>
> Skype ID: mauro.rossi76