From ?integrate:
When integrating over infinite intervals do so explicitly, rather
than just using a large number as the endpoint. This increases
the chance of a correct answer - any function whose integral over
an infinite interval is finite must be near zero for most of that
interval.
so the explanation was right there on the help page.
Your fuction is effectively 0 from x=2e3 and always very small.
integrate(f, 0, 2e3)
0.0001653797 with absolute error < 3.3e-06
If you multiply the result of f by 1e6 you get
integrate(f, 0, 2e3)
165.3797 with absolute error < 1.4e-08
integrate(f, 0, Inf)
165.3797 with absolute error < 0.0076
so this was a scaling issue.
Real names and proper signatures are preferred here.
On Thu, 21 Aug 2008, soneone ashamed of her real name wrote:
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.
--
Brian D. Ripley, [EMAIL PROTECTED]
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
______________________________________________
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.
