[Rd] minor flaw in integrate()

[Prof Brian Ripley]

I think throwing an error is a better solution: this is rather unlikely to 
be deliberate and returning NaN might postpone the detection too long.


On Tue, 3 Jul 2007, Duncan Murdoch wrote:

> On 7/3/2007 11:55 AM, Martin Maechler wrote:
>>>>>>> "PetRd" == Peter Ruckdeschel <Peter.Ruckdeschel at uni-bayreuth.de>
>>>>>>>     on Tue, 03 Jul 2007 17:26:43 +0200 writes:
>>
>>     PetRd> Thanks Martin and Duncan for your
>>     PetRd> comments,
>>
>>     PetRd> Martin Maechler wrote:
>>    >>>>>>> "DM" == Duncan Murdoch <murdoch at stats.uwo.ca>
>>    >>>>>>> on Mon, 02 Jul 2007 21:56:23 -0400 writes:
>>    >>
>>     DM> On 28/06/2007 5:05 PM, Peter Ruckdeschel wrote:
>>    >>>> Hi,
>>    >>>>
>>    >>>> I noticed a minor flaw in integrate() from package stats:
>>    >>>>
>>    >>>> Taking up arguments lower and upper from integrate(),
>>    >>>>
>>    >>>> if (lower ==  Inf) && (upper ==  Inf)
>>    >>>>
>>    >>>> or
>>    >>>>
>>    >>>> if (lower == -Inf) && (upper == -Inf)
>>    >>>>
>>    >>>> integrate() calculates the value for (lower==-Inf) && (upper==Inf).
>>    >>>>
>>    >>>> Rather, it should return 0.
>>    >>
>>     DM> Wouldn't it be better to return NA or NaN, for the same reason Inf/Inf
>>     DM> doesn't return 1?
>>    >>
>>     DM> Duncan Murdoch
>>    >>
>>    >> Yes indeed, I think it should return NaN.
>>
>>     PetRd> not quite convinced --- or more precisely:
>>
>>     PetRd> [ Let's assume lower = upper = Inf here,
>>     PetRd> case lower = upper = -Inf is analogue ]
>>
>>     PetRd> I'd say it depends on whether the (Lebesgue-) integral
>>
>>     PetRd> integral(f, lower = <some finite value>, upper = Inf)
>>
>>     PetRd> is well defined. Then, by dominated convergence, the integral should
>>     PetRd> default to 0.
>>
>>     PetRd> But I admit that then a test
>>
>>     PetRd> is.finite(integrate(f, lower = <some finite value>, upper = Inf)$value)
>>
>>     PetRd> would be adequate, too, which makes evaluation a little more expensive :-(
>>
>> No, that's not the Duncan's point I agreed on.
>> The argument is different:
>>
>> consider       Int(f, x, x^2)
>> 	       Int(f, x, 2*x)
>> 	       Int(f, x, exp(x))
>> etc,
>> These could conceivably give very different values,
>> with different limits for  x --> Inf
>>
>> Hence, 	       Int(f, Inf, Inf)
>>
>> is mathematically undefined, hence NaN
>
> In the case Peter was talking about, those limits would all be zero.
> But I don't think we could hope for the integrate() function in R to
> recognize integrability.  For example,
>
> > integrate(function(x) 1/x, 1e8, Inf)
> 1.396208e-05 with absolute error < 2.6e-05
>
> where the correct answer is Inf, since the integral is divergent.
>
> So I'd be fairly strongly opposed to returning 0.  Whether we return NaN
> or NA is harder:  I suspect the reason Inf-Inf or Inf/Inf is NaN is
> because this is handled on most platforms by the floating point hardware
> or the C run-time, rather than because we've made a deliberate decision
> for that.  Do we have other cases where NaN is used to mean "unable to
> determine the answer"?
>
> Duncan Murdoch
>
>> Martin
>>
>>     PetRd> If, otoh
>>
>>     PetRd> integrate(f, lower = <some finite value>, upper = Inf)
>>
>>     PetRd> throws an error, I agree, there should be a NaN ...
>>     PetRd> Best, Peter
>>
>> ______________________________________________
>> R-devel at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

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