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.
