[R] Help Regarding 'integrate'

[Prof Brian Ripley]

As my post showed, it is a scaling issue.  The function has so small a 
peak that it is effectively 0 -- when scaled sensibly integrate then works 
out of the box.

On Thu, 21 Aug 2008, Thomas Lumley wrote:

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