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 <[EMAIL PROTECTED]> wrote:
> From: A <[EMAIL PROTECTED]>
> Subject: [R] Help Regarding 'integrate'
> To: r-help@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@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.
______________________________________________
R-help@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.
