Raw File
polyfit.R
###
### p o l y f i t . R  Polynom
###

polyfit <- function(x, y, n = 1) {
    if (!is.numeric(x) || !is.numeric(y))
        stop("Arguments x and y must be numeric.")
    if (length(x) != length(y))
        stop("Vectors/matrices x and y must be of same length.")
    if (is.null(n) || n < 0 || ceiling(n) != floor(n))
        stop("Degree n must be a non-negative integer.")

    x <- x[1:length(x)]; y <- y[1:length(y)]
    A <- outer(x, seq(n, 0), "^")
    p <- qr.solve(A, y)
    return(p)
}

polyfit2 <- function(x, y, n = 1, p0 = NULL) {
    if (!is.numeric(x) || !is.numeric(y))
        stop("Arguments 'x' and 'y' must be numeric.")
    if (length(x) != length(y))
        stop("Vectors/matrices 'x' and 'y' must be of same length.")
    if (is.null(n) || n <= 0 || ceiling(n) != floor(n))
        stop("Argument 'n', order of fit, must be a positive integer.")
    if (is.null(p0))
        return(polyfit(x, y, n = n))
    else if (!is.numeric(p0) || length(p0) != 2)
        stop("Argument 'p0' must be a numeric vector of length 2.")

    x0 <- p0[1];  y0 <- p0[2]
    xx <- x - x0; yy <- y - y0

    M <- matrix(0, length(x), n)
    M[, n] <- xx
    for (i in (n-1):1) {
        M[, i] <- xx * M[, i+1]
    }
    pt <- qr.solve(M, yy)
    pt <- c(pt, y0)
    p <- numeric(n+1)
    for (i in (n+1):1) {
        p[i] <- polyval(pt, -x0)
        pt   <- polyder(pt)/(n-i+2)
    }
    return(p)
}
back to top