[R] Help Regarding 'integrate'

[Prof Brian Ripley]

>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 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.
>

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
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