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.
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
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.Wouldn't it be better to return NA or NaN, for the same reason Inf/Inf doesn't return 1? Duncan Murdoch> > 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
If integrate is changed it would be nice at the same time to make it into an S3 generic. deriv already is an S3 generic but strangely integrate is not. Ryacas provides a deriv method but for integrate Ryacas inconsistently provides Integrate since integrate is not generic. On 6/28/07, Peter Ruckdeschel <Peter.Ruckdeschel at uni-bayreuth.de> 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. > > 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 >