[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