Re: [Rd] minor flaw in integrate()

From: Duncan Murdoch <murdoch_at_stats.uwo.ca>
Date: Tue, 03 Jul 2007 13:21:19 -0400

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 Received on Tue 03 Jul 2007 - 17:24:57 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 - 18:35:55 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.