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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
On Mon, 25 Feb 2002, Ole Christensen wrote:> Dear Richard > > I think the cause of your problem with integrate() can be seen from > plotting the function you want integrate :Yes, but integrate() still seems to be unduly optimistic about its accuracy -- eight orders of magnitude is a lot to be wrong by. -thomas Thomas Lumley Asst. Professor, Biostatistics tlumley at u.washington.edu University of Washington, Seattle -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._