[Rd] minor flaw in integrate()
Martin Maechler
maechler at stat.math.ethz.ch
Tue Jul 3 10:45:33 CEST 2007
>>>>> "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.
Martin
>>
>> Quick fix:
>>
>> ### old code ###
>> ### [snip]
>> else {
>> if (is.na(lower) || is.na(upper))
>> stop("a limit is missing")
>> if (is.finite(lower)) {
>> inf <- 1
>> bound <- lower
>> }
>> else if (is.finite(upper)) {
>> inf <- -1
>> bound <- upper
>> }
>> else {
>> inf <- 2
>> bound <- 0
>> }
>> wk <- .External("call_dqagi", ff, rho = environment(),
>> as.double(bound), as.integer(inf), as.double(abs.tol),
>> as.double(rel.tol), limit = limit, PACKAGE = "base")
>> }
>> ### [snip]
>>
>> ### new code to replace the old one ###
>>
>> ### [snip]
>> else {
>> if (is.na(lower) || is.na(upper))
>> stop("a limit is missing")
>>
>> if (lower == upper){
>>
>> wk <- list("value" = 0, "abs.error" = 0,
>> "subdivisions" = subdivisions,
>> "ierr" = 0 )
>>
>> } else {
>> if (is.finite(lower)) {
>> inf <- 1
>> bound <- lower
>> }
>> else if (is.finite(upper)) {
>> inf <- -1
>> bound <- upper
>> }
>> else {
>> inf <- 2
>> bound <- 0
>> }
>> wk <- .External("call_dqagi", ff, rho = environment(),
>> as.double(bound), as.integer(inf),
>> as.double(abs.tol), as.double(rel.tol),
>> limit = limit, PACKAGE = "base")
>>
>> }
>> }
>> ### [snip]
>>
>> Best, Peter
>>
>> ______________________________________________
>> R-devel at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
DM> ______________________________________________
DM> R-devel at r-project.org mailing list
DM> https://stat.ethz.ch/mailman/listinfo/r-devel
More information about the R-devel
mailing list