Re: [Rd] minor flaw in integrate()

From: Prof Brian Ripley <ripley_at_stats.ox.ac.uk>
Date: Tue, 03 Jul 2007 18:27:39 +0100 (BST)

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

______________________________________________
R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Tue 03 Jul 2007 - 17:32:20 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Tue 03 Jul 2007 - 19:35:54 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.