HI there, I have been trying to use a code posted on R help to be able to
calculate area under the curve for complicated data points and there seems to be
an issue with the code: no "b" object found. I am not a good R user
and can''t find were the problem is. Any help? Thanks!!
This is the code ( as a test run I gave it this info because I know the answer:
x<-seq(1:50)
y<-seq(1:50)
simp <- function (y, a = NULL, b = NULL, x = NULL, n = 200)
{
if (is.null(a) | is.null(b)) {
if (is.null(x))
stop("No x values provided to integrate over.\n")
}
else {
x <- c(a, b)
}
fff <- 1
if (length(x) == 2) {
if (x[1] == x[2])
return(0)
if (x[2] < x[1]) {
fff <- -1
x <- rev(x)
}
x <- seq(x[1], x[2], length = n)
if (is.function(y))
y <- y(x)
else {
cat("y must be a function when x is\n")
cat("of length equal to 2.\n")
stop("Bailing out.\n")
}
equisp <- TRUE
}
else {
if (is.function(y))
y <- y(x)
else if (length(y) != length(x))
stop("Mismatch in lengths of x and y.\n")
s <- order(x)
x <- x[s]
ddd <- diff(x)
if (any(ddd == 0))
stop("Gridpoints must be distinct.\n")
equisp <- isTRUE(all.equal(diff(ddd), rep(0, length(ddd) - 1)))
y <- y[s]
}
n <- length(x) - 1
if (equisp) {
old.op <- options(warn = -1)
on.exit(options(old.op))
M <- matrix(y, nrow = n + 2, ncol = 4)[1:(n - 2), ]
h <- x[2] - x[1]
fc <- h * c(-1, 13, 13, -1)/24
aa <- apply(t(M) * fc, 2, sum)
a1 <- h * sum(y[1:3] * c(5, 8, -1))/12
an <- h * sum(y[(n - 1):(n + 1)] * c(-1, 8, 5))/12
return(fff * sum(c(a1, aa, an)))
}
m <- n%/%2
i <- 1:(m + 1)
a <- x[2 * i] - x[2 * i - 1]
i <- 1:m
b <- x[2 * i + 1] - x[2 * i]
o <- (a[i] * b + 2 * a[i] * a[i] - b * b)/(6 * a[i])
p <- (a[i] + b)^3/(6 * a[i] * b)
q <- (a[i] * b + 2 * b * b - a[i] * a[i])/(6 * b)
k <- numeric(n + 1)
k[1] <- o[1]
i <- 1:(m - 1)
k[2 * i] <- p[i]
k[2 * i + 1] <- q[i] + o[-1]
if (n > 2 * m) {
aa <- a[m + 1]
bb <- b[m]
den <- 6 * bb * (bb + aa)
k[2 * m] <- p[m] - (aa^3)/den
k[2 * m + 1] <- q[m] + (aa^3 + 4 * bb * aa^2 + 3 * aa *
bb^2)/den
k[2 * m + 2] <- (2 * bb * aa^2 + 3 * aa * bb^2)/den
}
else {
k[2 * m] <- p[m]
k[2 * m + 1] <- q[m]
}
fff * sum(k * y)
}
--
Mariana Varela, DVM, PhD
MRC-University of Glasgow Centre for Virus Research
Institute of Infection, Immunity and Inflammation
College of Medical, Veterinary and Life Sciences
Garscube Estate, 464 Bearsden Road
Henry Wellcome Building, room 436
Glasgow, G61 1QH
Scotland (UK)
Phone: +44 (0) 141 330 2196
[[alternative HTML version deleted]]
Why re-invent? Use lrm in rms. -- David Sent from my iPhone On Aug 15, 2011, at 10:07 AM, Mariana Varela <Mariana.Varela at glasgow.ac.uk> wrote:> HI there, I have been trying to use a code posted on R help to be able to calculate area under the curve for complicated data points and there seems to be an issue with the code: no "b" object found. I am not a good R user and can't find were the problem is. Any help? Thanks!! > This is the code ( as a test run I gave it this info because I know the answer: > x<-seq(1:50) > y<-seq(1:50) > > > > simp <- function (y, a = NULL, b = NULL, x = NULL, n = 200) > { > if (is.null(a) | is.null(b)) { > if (is.null(x)) > stop("No x values provided to integrate over.\n") > } > else { > x <- c(a, b) > } > fff <- 1 > if (length(x) == 2) { > if (x[1] == x[2]) > return(0) > if (x[2] < x[1]) { > fff <- -1 > x <- rev(x) > } > x <- seq(x[1], x[2], length = n) > if (is.function(y)) > y <- y(x) > else { > cat("y must be a function when x is\n") > cat("of length equal to 2.\n") > stop("Bailing out.\n") > } > equisp <- TRUE > } > else { > if (is.function(y)) > y <- y(x) > else if (length(y) != length(x)) > stop("Mismatch in lengths of x and y.\n") > s <- order(x) > x <- x[s] > ddd <- diff(x) > if (any(ddd == 0)) > stop("Gridpoints must be distinct.\n") > equisp <- isTRUE(all.equal(diff(ddd), rep(0, length(ddd) - 1))) y <- y[s] > } > n <- length(x) - 1 > if (equisp) { > old.op <- options(warn = -1) > on.exit(options(old.op)) > M <- matrix(y, nrow = n + 2, ncol = 4)[1:(n - 2), ] > h <- x[2] - x[1] > fc <- h * c(-1, 13, 13, -1)/24 > aa <- apply(t(M) * fc, 2, sum) > a1 <- h * sum(y[1:3] * c(5, 8, -1))/12 > an <- h * sum(y[(n - 1):(n + 1)] * c(-1, 8, 5))/12 > return(fff * sum(c(a1, aa, an))) > } > m <- n%/%2 > i <- 1:(m + 1) > a <- x[2 * i] - x[2 * i - 1] > i <- 1:m > b <- x[2 * i + 1] - x[2 * i] > o <- (a[i] * b + 2 * a[i] * a[i] - b * b)/(6 * a[i]) > p <- (a[i] + b)^3/(6 * a[i] * b) > q <- (a[i] * b + 2 * b * b - a[i] * a[i])/(6 * b) > k <- numeric(n + 1) > k[1] <- o[1] > i <- 1:(m - 1) > k[2 * i] <- p[i] > k[2 * i + 1] <- q[i] + o[-1] > if (n > 2 * m) { > aa <- a[m + 1] > bb <- b[m] > den <- 6 * bb * (bb + aa) > k[2 * m] <- p[m] - (aa^3)/den > k[2 * m + 1] <- q[m] + (aa^3 + 4 * bb * aa^2 + 3 * aa * > bb^2)/den > k[2 * m + 2] <- (2 * bb * aa^2 + 3 * aa * bb^2)/den > } > else { > k[2 * m] <- p[m] > k[2 * m + 1] <- q[m] > } > fff * sum(k * y) > } > > > -- > Mariana Varela, DVM, PhD > MRC-University of Glasgow Centre for Virus Research > Institute of Infection, Immunity and Inflammation > College of Medical, Veterinary and Life Sciences > Garscube Estate, 464 Bearsden Road > Henry Wellcome Building, room 436 > Glasgow, G61 1QH > Scotland (UK) > Phone: +44 (0) 141 330 2196 > > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at 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.
On 16/08/11 05:07, Mariana Varela wrote:> HI there, I have been trying to use a code posted on R help to be able to calculate area under the curve for complicated data points and there seems to be an issue with the code: no "b" object found. I am not a good R user and can't find were the problem is. Any help? Thanks!! > This is the code ( as a test run I gave it this info because I know the answer: > x<-seq(1:50) > y<-seq(1:50)There is nothing wrong with the code, except for a few typographical errors in the code *as you posted it*. (A missing "^" and a missing line feed.) How did you actually call simp()? Your posting does not show this. With your x and y (why wrap the 1:50 up in seq()? just x <- 1:50 is what is required) I get > simp(y,x=x) [1] 1249.5 (which is the right answer! :-) ) cheers, Rolf Turner