[R] errors in integrate function?

[Ole Christensen]

Dear Richard

I think the cause of your problem with integrate() can be seen from
plotting the function you want integrate :

gamma.by.normal <- function(y,x,shape,scale,stdev) {
   return(
dgamma(y,shape=shape,scale=scale)*dnorm(x=x-y,mean=0,sd=stdev) ) 
}

curve(gamma.by.normal(x,x=10,shape=1,scale=4,stdev=0.2),from=8, to=12)

shows that most of the probability mass is between 9.3 and 10.7. You are
integrating this function from zero to infinity by a procedure that
choses the discretisation points in an automatic way (which may not give
any points in [9.3;10.7] ).

In my experience, one should be suspicuous when using the integrate()
function. Always plot the function to see how it looks like.

You are probably aware that for integer values of the shape parameter,
you should not use numeric integration to calculate the integral. Closed
form expression exist. 

Cheers Ole

> I have been trying the integrate function in R, a function which would be 
> very useful for a current project of mine. But I am encountering errors 
> integrating the one function I have tried. The function to be integrated is 
> a product of a gamma demsity and a normal density: 
> 
> gamma.by.normal_function(y,x,shape,scale,stdev) 
>              return( dgamma(y,shape=shape,scale=scale)* 
>                          dnorm(x=x-y,mean=0,sd=stdev) ) 
> 
> > integrate(gamma.by.normal,lower=0,upper=Inf, 
>     rel.tol=1e-6,x=10,shape=1,scale=4,stdev=0.2) 
> > 4.589933e-11 with absolute error < 8.4e-11 
> 
> > integrate(gamma.by.normal,lower=0,upper=Inf, 
>     rel.tol=1e-12,x=10,shape=1,scale=4,stdev=0.2) 
> > 0.02054692 with absolute error < 4.5e-13 
> 
> The correct answer is the second. But with rel.tol=1e-6, the answer 
> returned is off by 1e-2 even though the error reported is 8e-11. 


-- 
Ole F. Christensen
Department of Mathematics and Statistics
Fylde College, Lancaster University 
Lancaster, LA1 4YF, England
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._