[R] Help Regarding 'integrate'

[Thomas Lumley]

On Thu, 21 Aug 2008, Moshe Olshansky wrote:

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

Numerical integration over a long (or infinite) interval is fine. The 
problem here is that integrate appears not to find the peak of the 
integrand.  The integrand is very flat near zero, which makes things 
harder.

It seems to work well to split the integral somewhere not too far to 
the left of its peak. Splits at 20, 50,100,or 200 give good agreement.

> integrate(f,0,100)
0.0001600355 with absolute error < 6.5e-10
> integrate(f,100,Inf)
5.355248e-06 with absolute error < 3.4e-06

>  integrate(f,0,200)
0.0001653797 with absolute error < 3.1e-05
> integrate(f,200,Inf)
3.141137e-13 with absolute error < 1.2e-13

> integrate(f,0,50)
2.702318e-05 with absolute error < 1.3e-16
> integrate(f,50, Inf)
0.0001383181 with absolute error < 8e-05

Transforming the integral is only going to be helpful if you have a good 
idea of what it looks like and what the integrate() function is expecting. 
An automated transformation to a bounded interval won't help since that is 
the first thing the underlying Fortran code does.

 	-thomas

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

Thomas Lumley			Assoc. Professor, Biostatistics
tlumley at u.washington.edu	University of Washington, Seattle