> Hello,
>
> I wrote the function below to integrate polysplines and thought that
> it may be useful to others. Please consider this code released under
> the GPL2 or later.
>
> Thanks,
>
> Bill
> <<integrate.polySpline.txt>>
Notice: This e-mail message, together with any attachments, contains
information of Merck & Co., Inc. (One Merck Drive, Whitehouse Station,
New Jersey, USA 08889), and/or its affiliates (which may be known
outside the United States as Merck Frosst, Merck Sharp & Dohme or
MSD and in Japan, as Banyu - direct contact information for affiliates is
available at http://www.merck.com/contact/contacts.html) that may be
confidential, proprietary copyrighted and/or legally privileged. It is
intended solely for the use of the individual or entity named on this
message. If you are not the intended recipient, and have received this
message in error, please notify us immediately by reply e-mail and
then delete it from your system.
integrate.polySpline <- function(spl=NA, lower=NA, upper=NA,
oob.action=stop) {
## Verify that all input arguments are given
if (any(is.na(spl))) {
stop("The spline to integrate must be given")
} else if (any(is.na(lower))) {
stop("The lower bound(s) to integrate must be given")
} else if (any(is.na(upper))) {
stop("The upper bound(s) to integrate must be given")
}
## Verify that the inputs are within the interpolation limits
if (any(min(lower) < min(spl$knots) |
max(lower) > max(spl$knots) |
min(upper) < min(spl$knots) |
max(upper) > max(spl$knots))) {
oob.action("The bounds are not within the bounds of the spline")
}
## We won't require that the lower bound is smaller than the upper
## bound because there are times that it may be desired to have them
## inverted, but we will internally swap the arguments to simplify
## the integration.
negate <- lower > upper
tmp <- lower
tmp[negate] <- upper[negate]
upper[negate] <- lower[negate]
lower[negate] <- tmp[negate]
negate <- negate*-1+(!negate)*1
coef <- spl$coefficients
## divide the coefficients by their exponents to prevent the need at
## every iteration.
for (i in 2:length(coef[1,])) {
coef[,i] <- coef[,i]/i
}
r <- vector(mode="numeric", length=length(lower))
for (i in 2:length(spl$knots)) {
## Only work on parts of the integral that are requied
keep <- lower <= spl$knots[i] & upper >= spl$knots[i-1]
## Only do the integration if necessary
if (any(keep)) {
## Lower and upper bounds for this piecewise integration
lb <- pmax(lower[keep], spl$knots[i-1]) - spl$knots[i-1]
ub <- pmin(upper[keep], spl$knots[i]) - spl$knots[i-1]
## This is pulled out of the loop for efficiency (removing the ^1)
r[keep] <- r[keep] + coef[i-1,1]*(ub-lb)
for (j in 2:length(coef[i,])) {
r[keep] <- r[keep] + coef[i-1,j]*(ub^j-lb^j)
}
}
}
return(negate*r)
}
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
