[R] area under the curve

[Comcast]

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.