[R] errors in integrate function?

[Richard Condit]

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.

The error can be further demonstrated by checking the integral for a 
sequence of values, x=seq(-5,25,by=0.1).

dgamma.norm_function(x,shape,scale,stdev)
{
  prob=numeric()

  for(i in 1:length(x))
    prob[i]=integrate(gamma.by.normal,lower=0,upper=Inf,rel.tol=1e-6,
                      x=x[i],shape=shape,scale=scale,stdev=stdev)$val

  return(prob)
}

 > y=dgamma.norm(x,shape=1,scale=4,stdev=.2)
 > plot(x,y)

The errors will be readily apparent in the plot, which should be a smooth 
probability density.

For some parameter combinations, integrate works properly for all x, eg 
shape=1, scale=2, stdev=2.

I need a fast integration routine, but right now I can't use this one. I 
would appreciate any suggestions.




Richard Condit
Smithsonian Tropical Research Institute
Unit 0948
APO AA 34002-0948

fax 507-212-8148 (Panama)
ph  507-212-8045 (Panama)
condit at ctfs.si.edu
conditctfs at yahoo.com

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