[R] Help Regarding 'integrate'

[Moshe Olshansky]

The phenomenon is most likely caused by numerical errors. I do not know how 'integrate' works but numerical integration over a very long interval does not look a good idea to me.
I would do the following:

f1<-function(x){
   return(dchisq(x,9,77)*((13.5/x)^5)*exp(-13.5/x))
}

f2<-function(y){
   return(dchisq(1/y,9,77)*((13.5*y)^5)*exp(-13.5*y)/y^2)
}
# substitute y = 1/x

I1 <- integrate(f1,0,1,abs.tol=abs.tol=1.0e-20*.Machine$double.eps^0.25)

I2 <- integrate(f2,0,1,abs.tol=abs.tol=1.0e-20*.Machine$double.eps^0.25)

--- On Thu, 21/8/08, A <born.to.b.wyld at gmail.com> wrote:

> From: A <born.to.b.wyld at gmail.com>
> Subject: [R] Help Regarding 'integrate'
> To: r-help at r-project.org
> Received: Thursday, 21 August, 2008, 4:44 PM
> I have an R function defined as:
> 
> f<-function(x){
>   return(dchisq(x,9,77)*((13.5/x)^5)*exp(-13.5/x))
> }
> 
> Numerically integrating this function, I observed a couple
> of things:
> 
> A) Integrating the function f over the entire positive real
> line gives an
> error:
> > integrate(f,0,Inf)
> Error in integrate(f, 0, Inf) : the integral is probably
> divergent
> 
> B)  Increasing the interval of integration actually
> decreases the value of
> the integral:
> > integrate(f,0,10^5)
> 9.813968e-06 with absolute error < 1.9e-05
> > integrate(f,0,10^6)
> 4.62233e-319 with absolute error < 4.6e-319
> 
> 
> Since the function f is uniformly positive, B) can not
> occur. Also, the
> theory tells me that the infinite integral actually exists
> and is finite, so
> A) can not occur. That means there are certain problems
> with the usage of
> function 'integrate' which I do not understand. The
> help document tells me
> that 'integrate' uses quadrature approximation to
> evaluate integrals
> numerically. Since I do not come from the numerical methods
> community, I
> would not know the pros and cons of various methods of
> quadrature
> approximation. One naive way that I thought of evaluating
> the above integral
> was by first locating the maximum of the function (may be
> by using
> 'optimize' in R) and then applying the
> Simpson's rule to an interval around
> the maximum. However, I am sure that the people behind the
> R project and
> other users have much better ideas, and I am sure the
> 'correct' method has
> already been implemented in R. Therefore, I would
> appreciate if someone can
> help me find it.
> 
> 
> Thanks
> 
> 	[[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.