Big Green CD 8/BGCD_8_Dev.iso/NEXTSTEP/UNIX/Educational/R-0.49-MI/R-0.49-I/library/base

Text File | 1997-09-14 | 175.6 KB | 6,323 lines

abline <-
function(a=NULL, b=NULL, h=NULL, v=NULL, reg=NULL, coef=NULL,
    col=par("col"), lty=par("lty"), ...)
{
    if(!is.null(reg)) a <- reg
    if(!is.null(a) && is.list(a)) {
        temp <- as.vector(coefficients(a))
        if(length(temp) == 1) {
            a <- 0
            b <- temp
        }
        else {
            a <- temp[1]
            b <- temp[2]
        }
    }
    if(!is.null(coef)) {
        a <- coef[1]
        b <- coef[2]
    }
    .Internal(abline(a, b, h, v, col, lty, ...))
    invisible()
}
aperm <- function(a, perm, resize=TRUE) {
    if (missing(perm))
        perm<-(length(dim(a)):1)
    else {
        if(length(perm) != length(dim(a)))
            stop("perm has incorrect length")
        if(!all(sort(perm)==1:length(perm)))
            stop("perm is not a permutation")
    }
    .Internal(aperm(a,perm,resize))
}
append <- function (x, values, after = length(x))
{
        lengx <- length(x)
        if (after <= 0)
                c(values, x)
        else if (after >= lengx)
                c(x, values)
        else c(x[1:after], values, x[(after + 1):lengx])
}
"apply"<-
function(X, MARGIN, FUN, ...)
{
    # ENSURE THAT FUN IS A FUNCTION
    if(is.character(FUN))
        FUN <- get(FUN, mode = "function")
    else if(mode(FUN) != "function") {
        f <- substitute(FUN)
        if(is.name(f))
            FUN <- get(as.character(f), mode = "function")
        else stop(paste("\"", f, "\" is not a function", sep = ""))
    }
    # ENSURE THAT X IS AN ARRAY OBJECT
    d <- dim(X)
    dl <- length(d)
    ds <- 1:length(d)
    if(dl == 0)
        stop("dim(X) must have a positive length")
    if(length(class(X)) > 0)
        X <- if(dl == 2) as.matrix(X) else as.array(X)
    dn <- dimnames(X)
    # EXTRACT THE MARGINS AND ASSOCIATED DIMNAMES
    s.call <- (1:length(d))[-MARGIN]
    s.ans <- (1:length(d))[MARGIN]
    d.call <- d[-MARGIN]
    d.ans <- d[MARGIN]
    dn.call <- dn[-MARGIN]
    dn.ans <- dn[MARGIN]
    # dimnames(X) <- NULL
    # DO THE CALLS
    newX <- aperm(X, c(s.call, s.ans))
    dim(newX) <- c(prod(d.call), prod(d.ans))
    d2 <- dim(newX)[2]
    ans <- vector("list", d2)
    for(i in 1:d2)
        ans[[i]] <- FUN(array(newX[,i], d.call, dn.call), ...)
    # ANSWER DIMS AND DIMNAMES
    ans.names <- names(ans[[1]])
    ans.list <- is.recursive(ans[[1]])
    ans.length <- length(ans[[1]])
    if(!ans.list)
        ans.list <- any(unlist(lapply(ans, length)) != ans.length)
    if(!ans.list)
        ans <- unlist(ans, recursive = F)
    if(length(MARGIN) == 1 && length(ans) == d2) {
        if(length(dn.ans[[1]]) > 0)
            names(ans) <- dn.ans[[1]]
        else names(ans) <- NULL
        return(ans)
    }
    else if(length(ans) == d2)
        return(array(ans, d.ans, dn.ans))
    else if(length(ans) > 0 && length(ans) %% d2 == 0) {
        if(is.null(dn.ans))
            return(array(ans, c(length(ans)/d2, d[MARGIN])))
        else return(array(ans, c(length(ans)/d2, d.ans),
                c(list(ans.names), dn.ans)))
    }
    else return(ans)
}
approx <- function(x, y, xout, method="lines", n=50, rule=1) {
    if( !is.numeric(x) || !is.numeric(y) )
        stop("approx: x and y must be numeric")
    if(length(x) != length(y)) stop("x and y must have equal lengths")
    ok <- !(is.na(x) | is.na(y))
    x <- x[ok]
    y <- y[ok]
    if(length(x) < 2) stop("approx requires at least two values to interpolate")
    o <- order(x)
    x <- x[o]
    y <- y[o]
    if(missing(xout)) {
        if(n <= 0) stop("approx requires n >= 1")
        xout <- seq(x[1], x[length(x)], length=n)
    }
    if(rule == 1) {
        low <- y[1]
        high <- y[length(x)]
    }
    else if(rule == 2){
        low <- NA
        high <- low
    }
    else stop("invalid extrapolation rule in approx")
    y<-.C("approx", as.double(x), as.double(y), length(x), xout=as.double(xout), length(xout), as.double(low), as.double(high))$xout
    return(list(x=xout,y=y))
}
approxfun <- function(x, y, method="lines", rule=1) {
    if(length(x) != length(y)) stop("x and y must have equal lengths")
    if( !is.numeric(x) || !is.numeric(y) )
        error("approxfun: x and y must be numeric")
    ok <- !(is.na(x) | is.na(y))
    x <- x[ok]
    y <- y[ok]
    if(length(x) < 2) stop("approx requires at least two values to interpolate")
    o <- order(x)
    x <- x[o]
    y <- y[o]
    if(rule == 1) {
        low <- y[1]
        high <- y[length(x)]
    }
    else if(rule == 2){
        low <- as.double(NA)
        high <- low
    }
    rm(method, rule)
    function(v)
        .C("approx", as.double(x), as.double(y), length(x), xout=as.double(v), length(v), low, high)$xout
}
array <- function(data = NA, dim = length(data), dimnames = NULL)
{
    data <- as.vector(data)
    vl <- prod(dim)
    if( length(data) != vl  ) {
        t1 <- ceiling(vl/length(data))
        data <- rep(data,t1)
        if( length(data) != vl )
            data <- data[1:vl]
    }
    dim(data) <- dim
    if(is.list(dimnames))
        dimnames(data) <- dimnames
    data
}
arrows <- function(x0, y0, x1, y1, length=0.25, angle=30, code=2,
        col=par("fg"), lty=NULL, xpd=FALSE) {
    .Internal(arrows(
        x0,
        y0,
        x1,
        y1,
        length=length,
        angle=angle,
        code=code,
        col=col,
        lty=lty,
        xpd=xpd))
}
as.logical <- function(x) .Internal(as.vector(x,"logical"))
as.integer <- function(x) .Internal(as.vector(x,"integer"))
as.real <- function(x) .Internal(as.vector(x,"real"))
as.complex <- function(x) .Internal(as.vector(x, "complex"))
as.double <- function(x) .Internal(as.vector(x,"real"))
as.single <- function(x)
{
    warning("type single is not supported in R")
    .Internal(as.vector(x,"real"))
}
as.character <- function(x) .Internal(as.vector(x,"character"))
as.list <- function(x) .Internal(as.vector(x,"list"))
as.vector <- function(x, mode="any") .Internal(as.vector(x,mode))
as.matrix <- function(x)
{
    UseMethod("as.matrix")
}
as.matrix.default <- function(x)
{
    if( is.matrix(x) )
        x
    else
        array(x, c(length(x),1), if(!is.null(names(x))) list(names(x), NULL) else NULL)
}
as.matrix.data.frame <- function(x)
{
    y <- .Internal(as.matrix.data.frame(x))
    dimnames(y) <- dimnames(x)
    y
}
as.null <- function(x) NULL
as.function <- function(x) stop("mode function cannot be assigned")
as.array <- function(x)
{
    if( is.array(x) )
        return(x)
    dim(x) <-length(x)
    return(x)
}
as.name <- function(x) .Internal(as.name(x))
# as.call <- function(x) stop("type call cannot be assigned")
as.numeric <- as.double
as.qr <- function(x) stop("you cannot be serious")
as.ts <- function(x) if(is.ts(x)) x else ts(x)
as.formula <- function(object)
    if(inherits(object, "formula")) object else formula(object)
assign <-
function(x, value, envir=sys.frame(sys.parent()), inherits=FALSE,
immediate=TRUE)
.Internal(assign(x,value,envir,inherits))
attach <- function(x, pos=2)
    .Internal(attach(x, pos, deparse(substitute(x))))
detach <- function(name, pos=2)
{
    if(!missing(name)) {
        name <- substitute(name)
        if(!is.character(name))
            name <- deparse(name)
        pos <- match(name, search())
        if(is.na(pos))
            stop("invalid name")
    }
    .Internal(detach(pos))
}
"objects" <-
function (name, pos = -1, envir = NULL, all.files = FALSE, pattern)
{
        if (!missing(name)) {
                name <- substitute(name)
                if (!is.character(name))
                        name <- deparse(name)
                pos <- match(name, search())
                if (is.na(pos))
                        stop("invalid name")
        }
        else if (!missing(pos)) {
                if (pos < 1 || pos > length(search()))
                        stop("invalid pos value")
        }
        else if (!missing(envir)) {
                pos <- 0
        }
        else {
                pos <- -1
        }
        all.files <- .Internal(ls(pos, envir, all.files))
    if(!missing(pattern))
        grep(pattern, all.files, value = TRUE)
    else all.files
}
ls <- objects
# Average a vector over the levels of a factor.
ave  <-  function (x, ...)
{
    l <- list(...)
    if (is.null(l)) {
        x[] <- mean(x)
    }
    else {
        g <- 1
        nlv <- 1
        for (i in 1:length(l)) {
            l[[i]] <- as.factor(l[[i]])
            g <- g + nlv * (as.numeric(l[[i]]) - 1)
            nlv <- nlv * length(levels(l[[i]]))
        }
        x[] <- lapply(split(x, g), mean)[g]
    }
    x
}
"axis" <-
function (which, at, labels = TRUE, ...)
{
    if (which%%2 == 1) {
        axp <- par("xaxp")
        usr <- par("usr")[1:2]
        log <- par("xlog")
    }
    else {
        axp <- par("yaxp")
        usr <- par("usr")[3:4]
        log <- par("ylog")
    }
    if (missing(at)) {
        if (log) {
            if (usr[2] < usr[1] + 1) {
                at <- seq(axp[1], axp[2], length = axp[3] + 1)
            }
            else {
                p1 <- ceiling(min(usr))
                p2 <- floor(max(usr))
                if (p2 - p1 < 2) {
                 at <- c(1, 2, 5) * 10^rep(p1:p2, rep(3, p2 - p1 + 1))
                 at <- at[10^usr[1] <= at & at < 10^usr[2]]
                }
                else at <- 10^seq(p1, p2, by = 1)
            }
        }
        else at <- seq(axp[1], axp[2], length = axp[3] + 1)
        ind <- rep(TRUE, length(at))
    }
    else {
        ind <- (usr[1] <= at & at <= usr[2])
    }
    if (is.logical(labels)) {
        if (labels) {
            if (!log)
                at[abs(at/(max(at) - min(at))) < 0.001] <- 0
            labels <- format(at, trim = T)
        }
        else labels <- rep("", length(at))
    }
    else labels <- format(labels, trim = T)
    .Internal(axis(which, as.double(at[ind]), labels[ind], ...))
}
backsolve <-
function(r, x, k=ncol(r))
{
    r <- as.matrix(r)
    x <- as.matrix(x)
    if(k <= 0 || nrow(x) != k) stop("invalid parameters in backsolve")
    z <- .Fortran("bkslv",
        as.double(r),
        nrow(r),
        as.integer(k),
        as.double(x),
        as.integer(k),
        y=matrix(0, k, ncol(x)),
        as.integer(1),
        info=integer(1),
        DUP=FALSE)
    if(z$info != 0) stop("singular matrix in backsolve")
    z$y
}
"barplot" <-
function (height, names.arg, col=NULL, border=par("fg"),
    beside=FALSE, space=0.2, legend.text,
    main=NULL, xlab=NULL, ylab=NULL, xlim, ylim, ...)
{
    opar <- par(yaxs="i", xpd=TRUE)
    on.exit(par(opar))
    if (is.matrix(height)) {
        if (beside) {
            delta <- 0.5 * (1 - space)
            if (missing(xlim))
                xlim <- c(0, ncol(height)) + 0.5
            if (missing(ylim))
                ylim <- range(-0.01, height)
            plot.new()
                        plot.window(xlim, ylim, log = "")
            for (i in 1:ncol(height)) {
                xx <- seq(i-delta, i+delta, length=nrow(height)+1)
                xl <- xx[1:nrow(height)]
                xr <- xx[1:nrow(height)+1]
                rect(xl, 0, xr, height[,i], col=col, xpd=TRUE)
            }
        }
        else {
            delta <- 0.5 * (1 - space)
            nheight <- rbind(0, apply(height, 2, cumsum))
            if (missing(xlim))
                xlim <- c(0, ncol(height)) + 0.5
            if (missing(ylim))
                ylim <- range(-0.01, nheight)
            plot.new()
            plot.window(xlim, ylim, log = "")
            for (i in 1:ncol(height))
                rect(i - delta, nheight[-1, i],
                     i + delta, nheight[1:nrow(height), i],
                     col = col, xpd=TRUE)
        }
        if(missing(names.arg))
            names.arg <- dimnames(height)[[2]]
        if(!is.null(names.arg)) {
            if(length(names.arg) != ncol(height))
                stop("incorrect number of names")
            for(i in 1:length(names.arg))
            axis(1, at=1:length(names.arg),
                labels=names.arg, lty=0)
        }
    }
    else {
        delta <- 0.5 * (1 - space)
        if (missing(xlim))
            xlim <- c(0, length(height)) + 0.5
        if (missing(ylim))
            ylim <- range(-0.01, height)
        plot.new()
        plot.window(xlim, ylim, log = "")
        rect(1:length(height) - delta, 0,
             1:length(height) + delta, height, col, xpd=TRUE)
        if(missing(names.arg))
            names.arg <- names(height)
        if(!is.null(names.arg))
            for(i in 1:length(names.arg))
            axis(1,1:length(height), labels=names.arg, lty=0)
    }
    if(!missing(legend.text) && !missing(col)) {
        xy <- par("usr")
        legend(xy[2]-xinch(0.1),xy[4]-yinch(0.1),
            legend=rev(legend.text), fill=rev(col),
            xjust=1, yjust=1)
    }
    title(main=main, xlab=xlab, ylab=ylab, ...)
    axis(2)
}
box <- function(lty="solid", ...)
    .Internal(box(lty=lty, ...))
boxplot <- function(x, ..., range=1.5, width=NULL, varwidth=FALSE,
    notch=FALSE, names, data=sys.frame(sys.parent()),
    plot=TRUE, border=par("fg"), col=NULL, log="", pars=NULL)
{
    if(is.language(x)) {
        if(length(x) == 3 && deparse(x[[1]]) == '~') {
            groups <- eval(x[[3]], data)
            x <- eval(x[[2]], data)
            groups <- split(x, groups)
        }
        else stop("invalid first argument")
        apars <- list(...)
        pars <- c(apars[named.elements(apars)], pars)
    }
    else {
        groups <- list(x, ...)
        n <- named.elements(groups)
        pars <- c(groups[n], pars)
        groups[n] <- NULL
        if(length(groups)==1 && is.list(x))
            groups <- x
    }
    n <- length(groups)
    if(!missing(names)) attr(groups, "names") <- names
    else if(is.null(attr(groups, "names"))) attr(groups, "names") <- 1:n
    for(i in 1:n)
        groups[i] <- list(boxplot.stats(groups[[i]], range))
    if(plot) {
        bxp(groups, width, varwidth=varwidth, notch=notch,
            border=border, col=col, log=log, pars=pars)
        invisible(groups)
    }
    else groups
}
boxplot.stats <- function(x, coef)
{
    nna <- !is.na(x)
    n <- length(nna)
    stats <- fivenum(x, na.rm=TRUE)
    iqr <- diff(stats[c(2, 4)])
    out <- x < (stats[2]-coef*iqr) | x > (stats[4]+coef*iqr)
    if(coef > 0) stats[c(1, 5)] <- range(x[!out], na.rm=TRUE)
    conf <- stats[3]+c(-1.58, 1.58)*diff(stats[c(2, 4)])/sqrt(n)
    list(stats=stats, n=n, conf=conf, out=x[out&nna])
}
bxp <- function(z, notch=FALSE, width=NULL, varwidth=FALSE,
    border=par("fg"), col=NULL, log="", pars=NULL, ...)
{
    bplt <- function(x, wid, stats, out, conf, notch, border, col)
    {
        if(!any(is.na(stats))) {
            wid <- wid/2
            if(notch) {
                xx <- x+wid*c(-1,1,1,0.5,1,1,-1,-1,-0.5,-1)
                yy <- c(stats[c(2,2)],conf[1],stats[3],conf[2],
                    stats[c(4,4)],conf[2],stats[3],conf[1])
                polygon(xx, yy, col=col, border=border)
                segments(x-wid/2, stats[3], x+wid/2, stats[3], col=border)
            }
            else {
                xx <- x+wid*c(-1,1,1,-1)
                yy <- stats[c(2,2,4,4)]
                polygon(xx, yy, col=col, border=border)
                segments(x-wid,stats[3],x+wid,stats[3],col=border)
            }
            segments(rep(x,2),stats[c(1,5)], rep(x,2), stats[c(2,4)], lty="dashed",col=border)
            segments(rep(x-wid/2,2),stats[c(1,5)],rep(x+wid/2,2), stats[c(1,5)],col=border)
            points(rep(x,length(out)), out, col=border)
        }
    }
    n <- length(z)
    limits <- numeric(0)
    nmax <- 0
    for(i in 1:n) {
        nmax <- max(nmax,z[[i]]$n)
        limits <- range(limits, z[[i]]$stats, z[[i]]$out)
    }
    if(!is.null(width)) {
        if(length(width) != n | any(is.na(width)) | any(width <= 0))
            stop("invalid boxplot widths")
        width <- 0.8*width/max(width)
    }
    else if(varwidth) {
        width <- 0.8*sqrt(z[[i]]$n/nmax)
    }
    if(n == 1) width <- 0.4
    else width <- rep(0.8,n)
    plot.new()
    plot.window(xlim=c(0.5,n+0.5), ylim=limits, log=log)
    for(i in 1:n) {
        if(missing(border) || length(border)==0)
            border <- par("fg")
        bplt(i,width[i],z[[i]]$stats,z[[i]]$out,
            z[[i]]$conf,notch=notch,
            border=border[(i-1)%%length(border)+1],
            col=if(is.null(col)) col
            else col[(i-1)%%length(col)+1])
    }
    if(n > 1) axis(1, at=1:n, labels=names(z))
    axis(2)
    do.call("title", c(pars, list(...)))
    box()
    invisible(1:n)
}
builtins <- function(internal=FALSE)
.Internal(builtins(internal))
cbind <- function(...) {
    if(any.data.frame(...))
        data.frame(...)
    else
        .Internal(cbind(...))
}
cat <- function(...,file="",sep=" ", fill=FALSE, labels=NULL,append=FALSE)
    .Internal(cat(list(...),file,sep,fill,labels,append))
#nchar <- function(x) {
#    x<-as.character(x)
#    .Internal(nchar(x))
#}
substr <- function(x,start,stop) {
    x<-as.character(x)
    .Internal(substr(x,as.integer(start),as.integer(stop)))
}
strsplit <- function(x,split) {
    x<-as.character(x)
    split<-as.character(split)
    .Internal(strsplit(x,split))
}
substring <- function(text,first,last=1000000)
{
        storage.mode(text) <- "character"
        n <- max(length(text), length(first), length(last))
        text <- rep(text, length = n)
        first <- rep(first, length = n)
        last <- rep(last, length = n)
        substr(text, first, last)
}
abbreviate<-function(names.arg, minlength = 4, use.classes = T, dot = F)
{
        #we just ignore use.classes
        if(minlength<=0)
                return(rep("",length(names.arg)))
        names.arg<-as.character(names.arg)
        dups<-duplicated(names.arg)
        old<-names.arg
        if(any(dups))
                names.arg<-names.arg[!dups]
        dup2<-rep(T,length(names.arg))
        x<-these<-names.arg
        repeat {
                ans<-.Internal(abbreviate(these,minlength,use.classes))
                x[dup2]<-ans
                dup2<-duplicated(x)
                if(!any(dup2))
                        break
                minlength<-minlength+1
                dup2 <- dup2 | match(x, x[duplicated(x)], 0)
                these<-names.arg[dup2]
        }
        if(any(dups))
                x<-x[match(old,names.arg)]
        if(dot)
                x<-paste(x,".",sep="")
        names(x)<-old
        x
}
chisquare.test <- function (x, y=NULL, correct=TRUE,
    p = rep(1/length(x), length(x)))
{
    if (is.matrix(x))
        if ((nrow(x) == 1) || (ncol(x) == 1))
            x <- as.vector(x)
    if (!is.null(y) && !is.matrix(x))  {
        x <- factor(x)
        y <- factor(y)
        x <- table(x,y)
    }
    if (is.matrix(x)) {
        row.totals <- apply(x, 1, sum)
        col.totals <- apply(x, 2, sum)
        E <- outer(row.totals, col.totals, "*")/sum(x)
        df <- (nrow(x) - 1) * (ncol(x) - 1)
        if (correct && nrow(x) == 2 && ncol(x) == 2)
            yates <- .5
        else yates <- 0
        chi <- (abs(x - E) - yates)^2/E
        dimnames(E) <- dimnames(x)
        dimnames(chi) <- dimnames(x)
    }
    else {
        if (length(x) != length(p))
            stop("Arguments must be same length")
        E <- sum(x) * p
        df <- length(x) - 1
        chi <- (x - E)^2/E
        names(chi) <- names(x)
        names(E) <- names(x)
    }
    cat("X =", round(sum(chi), 4), " df =", df, "  p-value =",
        round(1 - pchisq(sum(chi), df), 4), "\n")
    invisible(list(E = E, chi = chi))
}
chol <- function(x)
{
    if(!is.numeric(x))
        stop("non-numeric argument to chol")
    if(is.matrix(x)) {
        if(nrow(x) != ncol(x))
            stop("non-square matrix in chol")
        n <- nrow(x)
    }
    else {
        if(length(x) != 1)
            stop("non-matrix argument to chol")
        n <- as.integer(1)
    }
    if(!is.double(x)) storage.mode(x) <- "double"
    z <- .Fortran("chol",
        x=x,
        n,
        n,
        v=matrix(0, nr=n, nc=n),
        info=integer(1),
        DUP=FALSE)
    if(z$info)
        stop("singular matrix in chol")
    z$v
}
chol2inv <- function(x, size=ncol(x))
{
    if(!is.numeric(x))
        stop("non-numeric argument to chol2inv")
    if(is.matrix(x)) {
        nr <- nrow(x)
        nc <- ncol(x)
    }
    else {
        nr <- length(x)
        nc <- as.integer(1)
    }
    size <- as.integer(size)
    if(size <= 0 || size > nr || size > nc)
        stop("invalid size argument in chol2inv")
    if(!is.double(x)) storage.mode(x) <- "double"
    z <- .Fortran("ch2inv",
        x=x,
        nr,
        size,
        v=matrix(0, nr=size, nc=size),
        info=integer(1),
        DUP=FALSE)
    if(z$info)
        stop("singular matrix in chol2inv")
    z$v
}
colnames <- function(x) {
    dn <- dimnames(x)
    if(is.null(dn)) dn else dn[[2]]
}
"colnames<-" <- function(x, value) {
    dn <- dimnames(x)
    if(is.null(dn)) dimnames(x) <- list(dn, value)
    else dimnames(x) <- list(dn[[1]], value)
    x
}
rgb <- function(r, g, b, names=NULL)
.Internal(rgb(r, g, b, names))
hsv <- function(h=1,s=1,v=1,gamma=1)
.Internal(hsv(h,s,v,gamma))
# This is a quick little ``rainbow'' function.
rainbow <-
function (n, s=1, v=1, start=0, end=(n-1)/n, gamma=1)
{
    hsv(seq(start, end, length=n), s, v, gamma)
}
"topo.colors" <-
function (n)
{
    j <- round(n/3)
    k <- round(n/3)
    i <- n - j - k
    rval <- rainbow(i, start = 43/60, end = 31/60)
    rval <- c(rval, rainbow(j, start = 23/60, end = 10/60))
    rval <- c(rval, hsv(seq(from = 10/60, to = 6/60, length.out = k),
        s = seq(from = 1, to = 0.3, length.out = k), 1))
    rval
}
"terrain.colors" <-
function (n)
{
    j <- round(n/3)
    k <- round(n/3)
    i <- n - j - k
    rval <- hsv(23/60, 1, v = seq(0.6, 0.85, len = i))
    rval <- c(rval, hsv(seq(23/60, 10/60, length = j), s = 1,
        v = seq(0.85 , 1, length = j)))
    rval <- c(rval, hsv(seq(from = 10/60, to = 6/60, length.out = k),
        s = seq(from = 1 , to = 0.3, length.out = k), 1))
    rval
}
"heat.colors" <-
function (n)
{
    j <- round(n/4)
    i <- n - j
    rval <- rainbow(i, start = 0, end = 1/6)
    if (i>0) rval <- c(rval, hsv(1/6, seq(from = 1, to = 1/(2*j),
            length.out = j), 1))
}
complete.cases <- function(...) .Internal(complete.cases(...))
pi <- 4*atan(1)
letters <- c(
"a","b","c","d","e","f","g","h","i","j","k","l", "m",
"n","o","p","q","r","s","t","u","v","w","x","y","z")
LETTERS <- c(
"A","B","C","D","E","F","G","H","I","J","K","L", "M",
"N","O","P","Q","R","S","T","U","V","W","X","Y","Z")
month.name <- c(
"January", "February", "March", "April", "May", "June",
"July", "August", "September", "October", "November", "December")
month.abb <- c(
"Jan", "Feb", "Mar", "Apr", "May", "Jun",
"Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
contour <- function(x=seq(0,1,len=nrow(z)), y=seq(0,1,len=ncol(z)), z,
    nlevels=10, levels=pretty(range(z,na.rm=TRUE), nlevels), labcex=0,
    col=par("fg"), lty=par("lty"), add=FALSE)
{
    if(is.list(x)) {
        x <- x$x
        y <- x$y
    }
    if(any(diff(x) <= 0) || any(diff(y) <= 0))
        stop("increasing x and y values expected")
    if(!add) {
        plot.new()
        plot.window(range(x), range(y), "")
    }
    .Internal(contour(x, y, z, levels, col=col, lty=lty))
    if(!add) {
        axis(1)
        axis(2)
        box()
    }
    invisible()
}
contrasts <-
function(x, contrasts=TRUE)
{
    if (!is.factor(x))
        stop("contrasts apply only to factors")
    ctr <- attr(x,"contrasts")
    if(is.null(ctr)) {
        if(is.ordered(x))
            ctr <- get(options("contrasts")[[1]][[2]])(levels(x), contrasts=contrasts)
        else
            ctr <- get(options("contrasts")[[1]][[1]])(levels(x), contrasts=contrasts)
        dimnames(ctr) <- list(levels(x), dimnames(ctr)[[2]])
    }
    else if(is.character(ctr))
        ctr <- get(ctr)(levels(x), contrasts=contrasts)
    ctr
}
"contrasts<-" <-
function(x, ctr)
{
    if(!is.factor(x))
        stop("contrasts apply only to factors")
    if(is.numeric(ctr)) {
        ctr <- as.matrix(ctr)
        nlevs <- nlevels(x)
        if(nrow(ctr) != nlevs || ncol(ctr) >= nlevs)
            stop("invalid contrast matrix extents")
        cm <- qr(cbind(1,ctr))
        if(cm$rank != ncol(ctr)+1) stop("singular contrast matrix")
        cm <- qr.qy(cm, diag(nlevs))[,2:nlevs]
        cm[,1:ncol(ctr)] <- ctr
        dimnames(cm) <- list(levels(x),NULL)
    }
    else if(is.character(ctr))
        cm <- ctr
    else if(is.null(ctr))
        cm <- NULL
    else stop("numeric contrasts or contrast name expected")
    attr(x, "contrasts") <- cm
    x
}
contr.poly <-
function(n, contrasts=TRUE)
{
    normalize <- function(x) x/sqrt(sum(x^2))
    if(is.numeric(n) && length(n) == 1)
        levs <- 1:n
    else {
        levs <- n
        n <- length(n)
    }
    if(n < 2)
        stop(paste("Contrasts not defined for", n - 1,
            "degrees of freedom"))
    contr <- matrix(0, n, n)
    x <- 1:n
    d <- x - mean(x)
    contr[,1] <- rep(1/sqrt(n),n)
    contr[,2] <- normalize(d)
    if(n > 2)
        for(i in 3:n) {
            a1 <- sum(d*contr[,i-1]*contr[,i-1])
            a2 <- sum(d*contr[,i-1]*contr[,i-2])
            contr[,i] <- normalize((d-a1)*contr[,i-1]-a2*contr[,i-2])
        }
    dimnames(contr) <- list(levs, paste("^", 0:(n-1), sep=""))
    if(contrasts) {
        contr[, -1, drop=FALSE]
    }
    else {
        contr[, 1] <- 1
        contr
    }
}
contr.helmert <-
function (n, contrasts=TRUE)
{
    if (length(n) <= 1) {
        if(is.numeric(n) && length(n) == 1 && n > 1) levels <- 1:n
        else stop("contrasts are not defined for 0 degrees of freedom")
    }
    else levels <- n
    lenglev <- length(levels)
    if (contrasts) {
        cont <- array(-1, c(lenglev, lenglev-1), list(levels, NULL))
        cont[col(cont) <= row(cont) - 2] <- 0
        cont[col(cont) == row(cont) - 1] <- 1:(lenglev-1)
    }
    else {
        cont <- array(0, c(lenglev, lenglev), list(levels, levels))
        cont[col(cont) == row(cont)] <- 1
    }
    cont
}
contr.treatment <-
function(n, contrasts = TRUE)
{
    if(is.numeric(n) && length(n) == 1)
        levs <- 1:n
    else {
        levs <- n
        n <- length(n)
    }
    contr <- array(0, c(n, n), list(levs, levs))
    contr[seq(1, n^2, n + 1)] <- 1
    if(contrasts) {
        if(n < 2)
            stop(paste("Contrasts not defined for", n - 1,
                "degrees of freedom"))
        contr <- contr[, -1, drop = FALSE]
    }
    contr
}
contr.sum <-
function (n, contrasts=TRUE)
{
    if (length(n) <= 1) {
        if (is.numeric(n) && length(n) == 1 && n > 1)
            levels <- 1:n
        else stop("Not enough degrees of freedom to define contrasts")
    }
    else levels <- n
    lenglev <- length(levels)
    if (contrasts) {
        cont <- array(0, c(lenglev, lenglev - 1), list(levels, NULL))
        cont[col(cont) == row(cont)] <- 1
        cont[lenglev, ] <- -1
    }
    else {
        cont <- array(0, c(lenglev, lenglev), list(levels, levels))
        cont[col(cont) == row(cont)] <- 1
    }
    cont
}
"co.intervals" <-
function (x, number = 6, overlap = 0.5)
{
        x <- sort(x[!is.na(x)])
        n <- length(x)
        # "from the record"
        r <- n/(number * (1 - overlap) + overlap)
        l <- round(1 + 0:(number - 1) * (1 - overlap) * r)
        u <- round(r + 0:(number - 1) * (1 - overlap) * r)
        cbind(x[l], x[u])
}
panel.smooth <-
function(x, y, col, pch, f=2/3, iter=3, ...)
{
    points(x, y, pch=pch, col=col)
    lines(lowess(x, y, f=f, iter=iter), ...)
}
"coplot" <-
function (formula, data, given.values, panel=points, rows, columns, show.given = TRUE,
    col = par("fg"), pch=par("pch"), ...)
{
    deparen <- function(expr) {
        while (is.language(expr) && !is.name(expr) && deparse(expr[[1]]) == "(") expr <- expr[[2]]
        expr
    }
    bad.formula <- function() stop("invalid conditioning formula")
    bad.lengths <- function() stop("incompatible variable lengths")
    #  parse and check the formula
    formula <- deparen(formula)
    if (deparse(formula[[1]]) != "~")
        bad.formula()
    y <- deparen(formula[[2]])
    rhs <- deparen(formula[[3]])
    if (deparse(rhs[[1]]) != "|")
        bad.formula()
    x <- deparen(rhs[[2]])
    rhs <- deparen(rhs[[3]])
    if (is.language(rhs) && !is.name(rhs)
        && (deparse(rhs[[1]]) == "*" || deparse(rhs[[1]]) == "+")) {
        have.b <- TRUE
        a <- deparen(rhs[[2]])
        b <- deparen(rhs[[3]])
    }
    else {
        have.b <- FALSE
        a <- rhs
    }
    #  evaluate the formulae components to get the data values
    if (missing(data))
        data <- sys.frame(sys.parent())
    x.name <- deparse(x)
    x <- eval(x, data)
    nobs <- length(x)
    y.name <- deparse(y)
    y <- eval(y, data)
    if(length(y) != nobs) bad.lengths()
    a.name <- deparse(a)
    a <- eval(a, data)
    if(length(a) != nobs) bad.lengths()
    if (have.b) {
        b.name <- deparse(b)
        b <- eval(b, data)
        if(length(b) != nobs) bad.lengths()
    }
    else b <- NULL
    #  generate the given value intervals
    bad.givens <- function() stop("invalid given.values")
    if(missing(given.values)) {
        if(is.factor(a)) {
            a.intervals <- cbind(1:nlevels(a), 1:nlevels(a))
            a <- codes(a)
        }
        else a.intervals <- co.intervals(a)
        b.intervals <- NULL
        if (have.b)  {
            if(is.factor(b)) {
                b.intervals <- cbind(1:nlevels(b), 1:nlevels(b))
                b <- codes(b)
            }
            else b.intervals <- co.intervals(b)
        }
    }
    else {
        if(!is.list(given.values))
            given.values <- list(given.values)
        if(length(given.values) != (if(have.b) 2 else 1))
            bad.givens()
        a.intervals <- given.values[[1]]
        if(is.factor(a)) {
            if(is.character(a.intervals))
                a.levels <- match(a.levels, levels(a))
            else a.levels <- cbind(a.levels, a.levels)
            a <- codes(a)
        }
        else if(is.numeric(a)) {
            if(!is.numeric(a)) bad.givens()
            if(!is.matrix(a.intervals) || ncol(a.intervals) != 2)
                a.intervals <- cbind(a.intervals, a.intervals)
        }
        if(have.b) {
            b.intervals <- given.values[[2]]
            if(is.factor(b)) {
                if(is.character(b.intervals))
                    b.levels <- match(b.levels, levels(b))
                else b.levels <- cbind(b.levels, b.levels)
                b <- codes(b)
            }
            else if(is.numeric(b)) {
                if(!is.numeric(b)) bad.givens()
                if(!is.matrix(b.intervals) || ncol(b.intervals) != 2)
                    b.intervals <- cbind(b.intervals, b.intervals)
            }
        }
    }
    if(any(is.na(a.intervals))) bad.givens()
    if(have.b)
        if(any(is.na(b.intervals))) bad.givens()
    #  compute the page layout
    if (have.b) {
        rows <- nrow(b.intervals)
        columns <- nrow(b.intervals)
        nplots <- rows * columns
        total.rows <- rows + if (show.given)
            1
        else 0
        total.columns <- columns + if (show.given)
            1
        else 0
    }
    else {
        nplots <- nrow(a.intervals)
        if (missing(rows)) {
            if (missing(columns)) {
                rows <- ceiling(round(sqrt(nplots)))
                columns <- ceiling(nplots/rows)
            }
            else rows <- ceiling(nplots/columns)
        }
        else if (missing(columns))
            columns <- ceiling(nplots/rows)
        if (rows * columns < nplots)
            stop("rows * columns too small")
        total.rows <- rows + if (show.given)
            1
        else 0
        total.columns <- columns
    }
    #  plot that sucker!
    if(have.b) oma <- rep(5, 4) else oma <- c(5, 6, 5, 4)
    opar <- par(mfrow = c(total.rows, total.columns),
            oma = oma,
            mar = if (have.b) rep(0, 4) else c(0.5, 0, 0.5, 0),
            new = FALSE)
    on.exit(par(opar))
    plot.new()
    xlim <- range(x, na.rm = TRUE)
    ylim <- range(y, na.rm = TRUE)
    pch <- rep(pch, length=nobs)
    col <- rep(col, length=nobs)
    do.panel <- function(index) {
        istart <- (total.rows - rows) + 1
        i <- total.rows - ((index - 1)%/%columns)
        j <- (index - 1)%%columns + 1
        par(mfg = c(i, j, total.rows, total.columns), new = TRUE)
        plot.new()
        plot.window(xlim, ylim, log = "")
        if(any(id)) {
            grid(lty="solid")
            panel(x[id], y[id], col = col[id], pch=pch[id], ...)
        }
        if ((i == total.rows) && (j%%2 == 0))
            axis(1)
        if ((i == istart || index + columns > nplots) && (j%%2 == 1))
            axis(3)
        if ((j == 1) && ((total.rows - i)%%2 == 0))
            axis(2)
        if ((j == columns || index == nplots) && ((total.rows - i)%%2 == 1))
            axis(4)
#        if (i == total.rows)
#            axis(1, labels = (j%%2 == 0))
#        if (i == istart || index + columns > nplots)
#            axis(3, labels = (j%%2 == 1))
#        if (j == 1)
#            axis(2, labels = ((total.rows - i)%%2 == 0))
#        if (j == columns || index == nplots)
#            axis(4, labels = ((total.rows - i)%%2 == 1))
        box()
    }
    if(have.b) {
        count <- 1
        for(i in 1:rows) {
            for(j in 1:columns) {
                id <- ((a.intervals[j,1] <= a)
                    & (a <= a.intervals[j,2])
                    & (b.intervals[i,1] <= b)
                    & (b <= b.intervals[i,2]))
                do.panel(count)
                count <- count + 1
            }
        }
    }
    else {
        for (i in 1:nplots) {
            id <- ((a.intervals[i,1] <= a)
                & (a <= a.intervals[i,2]))
            do.panel(i)
        }
    }
    mtext(x.name, side=1, at=0.5*(columns/total.columns),
        outer=TRUE, line=5, xpd=TRUE)
    mtext(y.name, side=2, at=0.5*(rows/total.rows),
        outer=TRUE, line=4, xpd=TRUE)
    if(show.given) {
        mar <- par("mar")
        nmar <- mar + c(4,0,0,0)
        par(fig = c(0, columns/total.columns, rows/total.rows, 1),
            mar=nmar, new=TRUE)
        plot.new()
        nint <- nrow(a.intervals)
        plot.window(range(a.intervals, na.rm=T),
            c(0.5, nint+0.5), log="")
        rect(a.intervals[,1], 1:nint-0.3,
            a.intervals[,2], 1:nint+0.3, col=gray(0.9))
        axis(3)
        axis(1, labels=FALSE)
        box()
        mtext(paste("Given :", a.name),
            side=3, at=mean(par("usr")[1:2]), line=3, xpd=T)
        if(have.b) {
                nmar <- mar + c(0, 4, 0, 0)
            par(fig = c(columns/total.columns, 1,
                0, rows/total.rows), mar=nmar, new=TRUE)
            plot.new()
            nint <- nrow(b.intervals)
            plot.window(c(0.5, nint+0.5),
                range(b.intervals, na.rm=T), log="")
            rect(1:nint-0.3, b.intervals[,1],
                1:nint+0.3, b.intervals[,2], col=gray(0.9))
            axis(4)
            axis(2, labels=FALSE)
            box()
            mtext(paste("Given :", b.name),
                side=4, at=mean(par("usr")[3:4]), line=3, xpd=T)
        }
    }
}
cor <- function (x, y=NULL, use="all.obs")
{
    na.method <- pmatch(use, c("all.obs", "complete.obs", "pairwise.complete.obs"))
    .Internal(cor(as.matrix(x), if(is.null(y)) y else as.matrix(y), na.method))
}
cov <- function (x, y=NULL, use="all.obs")
{
    na.method <- pmatch(use, c("all.obs", "complete.obs", "pairwise.complete.obs"))
    .Internal(cov(as.matrix(x), if(is.null(y)) y else as.matrix(y), na.method))
}
curve <- function(expr, from, to, n=100, add=FALSE, type="l", ...) {
    expr <- substitute(expr)
    lims <- par("usr")
    if(missing(from)) from <- lims[1]
    if(missing(to)) to <- lims[2]
    x <- seq(from,to,length=n)
    y <- eval(expr)
    if(add)
        lines(x, y, ...)
    else
        plot(x, y, type="l", ...)
}
cut <- function(x, breaks, labels) {
    if (!is.numeric(x))
        stop("cut: x must be numeric")
    if (length(breaks) == 1) {
        if (is.na(breaks) | breaks < 2)
            stop("invalid number of breaks")
        nbreaks <- breaks
        breaks <- range(x,na.rm=TRUE)
        #breaks <- pretty(breaks + c(0, diff(breaks)/1000), nbreaks)
        breaks <- seq(breaks[1]-diff(breaks)/1000,breaks[2]+diff(breaks)/1000, len=nbreaks+1)
    }
    breaks <- sort(breaks)
    if(any(duplicated(breaks)))
        stop("cut: breaks are not unique")
    if(missing(labels)) {
        labels <- paste("Range",1:(length(breaks)-1))
    }
    else {
        if(length(labels)!=length(breaks)-1)
            stop("labels/breaks length conflict")
    }
    codes <- .C("bincode",
        as.double(x),
        length(x),
        as.double(breaks),
        length(breaks),
        integer(length(x)),
        TRUE)[[5]]
    factor(codes,1:max(codes),labels)
}
data <- function(..., list=character(0))
{
    names <- c(as.character(substitute(list(...))[-1]), list)
    if(length(names) == 0) {
        system(paste("$RHOME/cmd/pager", system.file("data",
            "index.doc")))
    }
    else
        for(name in names) {
            file <- system.file("data", name)
            if(file == "") stop(paste("no data set called", name))
            else source(file)
    }
    invisible(names)
}
data.matrix <-
function(frame)
{
    if(!is.data.frame(frame))
        return(as.matrix(frame))
    log <- unlist(lapply(frame, is.logical))
    num <- unlist(lapply(frame, is.numeric))
    fac <- unlist(lapply(frame, is.factor))
    if(!all(log|fac|num))
        stop("non-numeric data type in frame")
    d <- dim(frame)
    x <- matrix(nr=d[1],nc=d[2],dimnames=dimnames(frame))
    for(i in 1:length(frame)) {
        xi <- frame[[i]]
        if(is.logical(xi)) x[,i] <- as.numeric(xi)
        else if(is.numeric(xi)) x[,i] <- xi
        else x[,i] <- codes(xi)
    }
    x
}
frame.cvt <-
function(x, as.is=FALSE)
{
    if(!as.is && (is.character(x) || is.logical(x)))
        x <- factor(x)
    x
}
data.frame <-
function (..., row.names=NULL, col.names=NULL, as.is=FALSE)
{
    frame <- list(...)
    n <- length(frame)
    as.is <- rep(as.is, length=n)
    for(i in 1:n) {
        as.is.i <- attr(frame[[i]], "AsIs")
        if(!is.null(as.is.i))
            as.is[i] <- as.is.i
    }
    if (is.null(col.names)) {
        v <- substitute(list(...))[-1]
        for (i in 1:length(v)) if (!is.symbol(v[[i]]))
            v[[i]] <- paste("X", i, sep="")
        arg.names <- as.character(v)
        col.names <- names(frame)
        if (is.null(col.names))
            col.names <- arg.names
        else {
            nameless <- (nchar(col.names) == 0)
            col.names[nameless] <- arg.names[nameless]
        }
    }
    names(frame) <- as.character(col.names)
    for (i in 1:length(frame)) {
        if (is.list(frame[[i]]))
            for (j in 1:length(frame[[i]])){
                if (!is.numeric(frame[[i]][[j]]) && !is.factor(frame[[i]][[j]]))
                    frame[[i]][[j]] <- frame.cvt(frame[[i]][[j]], as.is=as.is[i])
        }
        else {
            if (!is.numeric(frame[[i]]) && !is.factor(frame[[i]]))
                frame[[i]] <- frame.cvt(frame[[i]], as.is=as.is[i])
        }
    }
    .Internal(data.frame(frame, as.character(row.names), as.logical(as.is)))
}
row.names <- function(x) attr(x,"row.names")
"row.names<-" <- function(x,value) {
    if( !is.data.frame(x) )
        return(data.frame(x, row.names=value))
    else
        attr(x,"row.names") <- as.character(value)
    x
}
"is.na.data.frame" <-
function (x)
{
        y <- do.call("cbind", lapply(x, "is.na"))
        rownames(y) <- row.names(x)
        y
}
I <- function(x) {
    attr(x,"AsIs") <- TRUE
    x
}
plot.data.frame <-
function (x, ...)
{
    if (!is.data.frame(x))
        stop("plot.data.frame applied to non data frame")
    xm <- data.matrix(x)
    if (ncol(xm) == 1) {
        stripplot(x, ...)
    }
    else if (ncol(xm) == 2) {
        plot(xm)
    }
    else {
        pairs(xm, ...)
    }
}
t.data.frame<- function(x)
{
    x <- as.matrix(x)
    NextMethod("t")
}
de.ncols <- function(inlist)
{
    ncols <- matrix(0, nrow=length(inlist), ncol=2)
    i <- 1
    for( telt in inlist ) {
        if( is.matrix(telt) ) {
                ncols[i, 1] <- ncol(telt)
                ncols[i, 2] <- 2
        }
        else if( is.list(telt) ) {
            for( telt2 in telt )
                if( !is.vector(telt2) ) stop("wrong argument to dataentry")
            ncols[i, 1] <- length(telt)
            ncols[i, 2] <- 3
        }
        else if( is.vector(telt) ) {
            ncols[i, 1] <- 1
            ncols[i, 2] <- 1
        }
        else stop("wrong argument to dataentry")
        i <- i+1
    }
    return(ncols)
}
de.setup <- function(ilist, list.names, incols)
{
    ilen <- sum(incols)
    ivec <- vector("list", ilen)
    inames <- vector("list", ilen)
    i <- 1
    k <- 0
    for( telt in ilist ) {
        k <- k+1
        if( is.list(telt) ) {
            y <- names(telt)
            for( j in 1:length(telt) ) {
                ivec[[i]] <- telt[[j]]
                if( is.null(y) || y[j]=="" )
                    inames[[i]] <- paste("var", i, sep="")
                else inames[[i]] <- y[j]
                i <- i+1
            }
        }
        else if( is.vector(telt) ) {
            ivec[[i]] <- telt
            inames[[i]] <- list.names[[k]]
            i <- i+1
        }
        else if( is.matrix(telt) ) {
            y <- dimnames(telt)[[2]]
            for( j in 1:ncol(telt) ) {
                ivec[[i]] <- telt[, j]
                if( is.null(y) || y[j]=="" )
                    inames[[i]] <- paste("var", i, sep="")
                else inames[[i]] <- y[j]
                i <- i+1
            }
        }
        else stop("wrong argument to dataentry")
    }
    names(ivec) <- inames
    return(ivec)
}
# take the data in inlist and restore it to the format described by ncols and coltypes
de.restore <- function(inlist, ncols, coltypes, argnames, args)
{
    rlist <- vector("list", length=length(ncols))
    rnames <- vector("character", length=length(ncols))
    j <- 1
    lnames <- names(inlist)
    for( i in 1:length(ncols) ) {
        if(coltypes[i]==2) {
            tlen <- length(inlist[[j]])
            x <- matrix(0, nrow=tlen, ncol=ncols[i])
            cnames <- vector("character", ncol(x))
            for( ind1 in 1:ncols[i]) {
                if(tlen != length(inlist[[j]]) ) {
                    warning("could not restore type information")
                    return(inlist)
                }
                x[, ind1] <- inlist[[j]]
                cnames[ind1] <- lnames[j]
                j <- j+1
            }
            if( dim(x) == dim(args[[i]]) )
                rn <- dimnames(args[[i]])[[1]]
            else rn <- NULL
            if( any(cnames!="") )
                dimnames(x) <- list(rn, cnames)
            rlist[[i]] <- x
            rnames[i] <- argnames[i]
        }
        else if(coltypes[i]==3) {
            x <- vector("list", length=ncols[i])
            cnames <- vector("character", ncols[i])
            for( ind1 in 1:ncols[i]) {
                x[[ind1]] <- inlist[[j]]
                cnames[ind1] <- lnames[j]
                j <- j+1
            }
            if( any(cnames!="") )
                names(x) <- cnames
            rlist[[i]] <- x
            rnames[i] <- argnames[i]
        }
        else {
            rlist[[i]] <- inlist[[j]]
            j <- j+1
            rnames[i] <- argnames[i]
        }
    }
    names(rlist) <- rnames
    return(rlist)
}
de <- function(..., Modes=NULL, Names=NULL)
{
    sdata <- list(...)
    snames <- as.character(substitute(list(...))[-1])
    if( is.null(sdata) ) {
        if( is.null(Names) ) {
            if( !is.null(Modes) ) {
                odata <- vector("list", length=length(Modes))
            }
            else odata <- vector("list", length=1)
        }
        else {
            if( (length(Names) != length(Modes)) && !is.null(Modes) ) {
                warning("modes argument ignored")
                Modes <- NULL
            }
            odata <- vector("list", length=length(Names))
            names(odata) <- Names
        }
        ncols <- rep(1, length(odata))
        coltypes <- rep(1, length(odata))
    }
    else {
        ncols <- de.ncols(sdata)
        coltypes <- ncols[, 2]
        ncols <- ncols[, 1]
        odata <- de.setup(sdata, snames, ncols)
        if( !is.null(Names) )
            if( length(Names) != length(odata) )
                warning("names argument ignored")
            else names(odata) <- Names
        if( !is.null(Modes) )
            if( length(Modes) != length(odata) ) {
                warning("modes argument ignored")
                Modes <- NULL
            }
    }
    rdata <- dataentry(odata, Modes)
    t1 <- length(rdata)==sum(ncols)
    if( t1 && any(coltypes!=1) )
        rdata <- de.restore(rdata, ncols, coltypes, snames, sdata)
    else if( any(coltypes!=1) ) warning("could not restore data types properly")
    return(rdata)
}
data.entry <- function(..., Modes=NULL, Names=NULL)
{
    tmp1 <- de(..., Modes=Modes, Names=Names)
    j <- 1
    for(i in names(tmp1) ) {
        assign(i, tmp1[[j]], env=.GlobalEnv)
        j <- j+1
    }
    invisible(NULL)
}
"density" <-
function (x, bw, adjust=1, kernel="gaussian", n=512, width, from,
    to, cut = 3, plot.graph = FALSE)
{
    if (!is.numeric(x))
        stop("argument must be numeric")
    name <- deparse(substitute(x))
    x <- x[!is.na(x)]
        k.list <- c("gaussian", "rectangular", "triangular", "cosine")
        method <- pmatch(kernel, k.list)
        if(is.na(method))
          stop(paste("kernel must be a 'pmatch' of",
                     paste(k.list,collapse=', ')))
        if(n > 512) n <- 2^ceiling(log2(n)) #- to be fast with FFT
        if (missing(bw)) {
      if(missing(width))
        bw <- adjust * 1.06 * min(sd(x), IQR(x)/1.34) * length(x)^-0.2
      else
        bw <- 0.25 * width
        }
    if (missing(from))
        from <- min(x) - cut * bw
    if (missing(to))
        to <- max(x) + cut * bw
    y <- .C("massdist",
                x = as.double(x),
                nx= length(x),
                xlo = as.double(from),
                xhi = as.double(to),
                y = double(2 * n),
                ny= as.integer(n)) $ y
    xords <- seq(from, by = (to - from)/(n - 1), length = 2 * n)
    kords <- xords - from
    kords[(n + 2):(2 * n)] <- -kords[n:2]
    if (method == 1) {
        bw <- bw
        kords <- dnorm(kords, sd = bw)
    }
    else if (method == 2) {
        a <- bw/0.2886751
        kords <- ifelse(abs(kords) < 0.5 * a, 1/a, 0)
    }
    else if (method == 3) {
        a <- bw/0.4082483
        kords <- ifelse(abs(kords) < a, (1 - abs(kords)/a)/a, 0)
    }
    else if (method == 4) {
        a <- bw/1.135724
        kords <- ifelse(abs(kords) < a * pi, (1 + cos(kords/a))/(2*pi*a), 0)
    }
    else stop("unknown density estimation kernel")
    kords <- convolve(y, kords)[1:n]
    xords <- seq(from, by = (to - from)/(n - 1), length = n)
    structure(list(x = xords, y = kords, bw = bw,
        call=match.call(), name=name),
        class="density")
}
plot.density <-
function(s, main="", xlab, ylab="Density", type="l", ...)
{
    if(missing(xlab)) xlab <- paste("Bandwidth =", s$bw)
    plot.default(s, main=main, xlab=xlab, ylab=ylab, type=type, ...)
}
bw.ucv <-
function(x, samples=100)
{
        fucv <- function(h)
                .C("ucv", length(x), x,  as.double(h), u=1)$u
        n <- length(x)
        if(samples > 0 && n > samples) x <- sample(x, samples)
        hmax <- 1.144 * sqrt(var(x)) * length(x)^(-1/5) * 4
        storage.mode(x) <- "double"
        0.25 * optimize(fucv, c(0.1*hmax, hmax), tol=0.01*hmax)$minimum * (length(x)/n)^0.2
}
bw.bcv <- function(x, samples=100)
{
    fbcv <- function(h)
        .C("bcv", length(x), x,  as.double(h), u=1)$u
    n <- length(x)
    if(samples > 0 && n > samples) x <- sample(x, samples)
    hmax <- 1.144 * sqrt(var(x)) * length(x)^(-1/5) * 4
    storage.mode(x) <- "double"
    0.25 * optimize(fbcv, c(0.1*hmax, hmax), tol=0.01*hmax)$minimum * (length(x)/n)^0.2
}
bw.sj <- function(x, samples=100)
{
    SDh <- function(x, h) .C("phi4", length(x), x,  as.double(h), u=double(1))$u
    TDh <- function(x, h) .C("phi6", length(x), x, as.double(h), u=double(1))$u
    fSD <- function(h, x, alph2, c1) (c1/SDh(x, alph2 * h^(5/7)))^(1/5) - h
    lambda <- IQR(x)
    n1 <- length(x)
    if(samples > 0 && n1 > samples) x <- sample(x, samples)
    storage.mode(x) <- "double"
    n <- length(x)
    hmax <- 1.144 * sqrt(var(x)) * n^(-1/5)
    a <- 0.92 * lambda * n^(-1/7)
    b <- 0.912 * lambda * n^(-1/9)
    c1 <- 1/(2*sqrt(pi)*n)
    TD  <- -TDh(x, b)
    alph2 <- 1.357*(SDh(x,a)/TD)^(1/7)
    res <- uniroot(fSD, c(0.1*hmax, hmax), tol=0.01*hmax,
        x=x, alph2=alph2, c1=c1)$root
    res * (n/n1)^0.2
}
diag <-
function(x = 1, nrow, ncol = n)
{
    if(is.matrix(x) && nargs() == 1)
        return(as.matrix(x)[1 + 0:(min(dim(x)) - 1) * (dim(x)[1] + 1)])
    if(missing(x))
        n <- nrow
    else if(length(x) == 1 && missing(nrow) && missing(ncol)) {
        n <- as.integer(x)
        x <- 1
    }
    else n <- length(x)
    if(!missing(nrow))
        n <- nrow
    p <- ncol
    y <- array(0, c(n, p))
    y[1 + 0:(min(n, p) - 1) * (n + 1)] <- x
    y
}
"diag<-" <-
function(x, value)
{
    dx <- dim(x)
    if(length(dx) != 2 || prod(dx) != length(x))
        stop("only matrix diagonals can be replaced")
    i <- 1:min(dx)
    if(length(value) != 1 && length(value) != length(i))
        stop("replacement diagonal has wrong length")
    x[cbind(i, i)] <- value
    x
}
"diff" <-
function (x, lag = 1, differences = 1)
{
    ismat <- is.matrix(x)
    if (ismat)
        xlen <- dim(x)[1]
    else xlen <- length(x)
    if (lag < 1 | differences < 1)
        stop("Bad value for lag or differences")
    if (lag * differences >= xlen)
        return(x[0])
    r <- x
    s <- 1:lag
    if (is.matrix(r)) {
        for (i in 1:differences) {
            rlen <- dim(r)[1]
            r <- r[-s, , drop = FALSE] - r[-(rlen + 1 - s), , drop = FALSE]
        }
    }
    else for (i in 1:differences) {
        r <- r[-s] - r[-(length(r) + 1 - s)]
    }
    xtsp <- attr(x, "tsp")
    if (is.null(xtsp)) r
    else ts(r, end = xtsp[2], freq = xtsp[3])
}
dexp <- function(x, rate=1) .Internal(dexp(x, 1/rate))
pexp <- function(q, rate=1) .Internal(pexp(q, 1/rate))
qexp <- function(p, rate=1) .Internal(qexp(p, 1/rate))
rexp <- function(n, rate=1) .Internal(rexp(n, 1/rate))
dunif <- function(x, min=0, max=1) .Internal(dunif(x, min, max))
punif <- function(q, min=0, max=1) .Internal(punif(q, min, max))
qunif <- function(p, min=0, max=1) .Internal(qunif(p, min, max))
runif <- function(n, min=0, max=1) .Internal(runif(n, min, max))
dnorm <- function(x, mean=0, sd=1) .Internal(dnorm(x, mean, sd))
pnorm <- function(q, mean=0, sd=1) .Internal(pnorm(q, mean, sd))
qnorm <- function(p, mean=0, sd=1) .Internal(qnorm(p, mean, sd))
rnorm <- function(n, mean=0, sd=1) .Internal(rnorm(n, mean, sd))
dcauchy <- function(x, location=0, scale=1) .Internal(dcauchy(x, location, scale
))
pcauchy <- function(q, location=0, scale=1) .Internal(pcauchy(q, location, scale
))
qcauchy <- function(p, location=0, scale=1) .Internal(qcauchy(p, location, scale
))
rcauchy <- function(n, location=0, scale=1) .Internal(rcauchy(n, location, scale
))
dgamma <- function(x, shape, scale=1) .Internal(dgamma(x, shape, scale))
pgamma <- function(q, shape, scale=1) .Internal(pgamma(q, shape, scale))
qgamma <- function(p, shape, scale=1) .Internal(qgamma(p, shape, scale))
rgamma <- function(n, shape, scale=1) .Internal(rgamma(n, shape, scale))
dlnorm <- function(x, meanlog=0, sdlog=1) .Internal(dlnorm(x, meanlog, sdlog))
plnorm <- function(q, meanlog=0, sdlog=1) .Internal(plnorm(q, meanlog, sdlog))
qlnorm <- function(p, meanlog=0, sdlog=1) .Internal(qlnorm(p, meanlog, sdlog))
rlnorm <- function(n, meanlog=0, sdlog=1) .Internal(rlnorm(n, meanlog, sdlog))
dlogis <- function(x, location=0, scale=1) .Internal(dlogis(x, location, scale))
plogis <- function(q, location=0, scale=1) .Internal(plogis(q, location, scale))
qlogis <- function(p, location=0, scale=1) .Internal(qlogis(p, location, scale))
rlogis <- function(n, location=0, scale=1) .Internal(rlogis(n, location, scale))
dweibull <- function(x, shape, scale=1) .Internal(dweibull(x, shape, scale))
pweibull <- function(q, shape, scale=1) .Internal(pweibull(q, shape, scale))
qweibull <- function(p, shape, scale=1) .Internal(qweibull(p, shape, scale))
rweibull <- function(n, shape, scale=1) .Internal(rweibull(n, shape, scale))
"dotplot" <-
function (x, labels = NULL, groups = NULL, gdata = NULL, cex = par("cex"),
        pch = 21, gpch = 21, bg = par("bg"), color = par("fg"),
        gcolor = par("fg"), ...)
{
        opar <- par("mar", "cex", "yaxs")
        on.exit(par(opar))
        par(cex = cex, yaxs = "i")
        if (is.matrix(x)) {
                if (is.null(labels))
                        labels <- rownames(x)
                if (is.null(labels))
                        labels <- as.character(1:nrow(x))
                labels <- rep(labels, length = length(x))
                if (is.null(groups))
                        groups <- col(x, as.factor = TRUE)
                glabels <- levels(groups)
        }
        else {
                if (is.null(labels))
                        labels <- names(x)
                if (!is.null(groups))
                        glabels <- levels(groups)
                else glabels <- NULL
        }
        linch <- 0
        ginch <- 0
        goffset <- 0
        if (!is.null(labels))
                linch <- max(strwidth(labels, "inch"), na.rm = TRUE)
        if (!is.null(glabels)) {
                ginch <- max(strwidth(glabels, "inch"), na.rm = TRUE)
                goffset <- 0.4
        }
        if (!(is.null(labels) && is.null(glabels))) {
                nmai <- par("mai")
                nmai[2] <- nmai[4] + max(linch + goffset, ginch) +
                        0.1
                par(mai = nmai)
        }
        if (is.null(groups)) {
                o <- 1:length(x)
                y <- 1:length(x)
                ylim <- c(0, length(x) + 1)
        }
        else {
                o <- rev(order(codes(groups)))
                x <- x[o]
                groups <- groups[o]
                offset <- cumsum(c(0, diff(codes(groups)[o])))
                y <- 1:length(x) + 2 * offset
                ylim <- range(0, y + 2)
        }
        plot.new()
        plot.window(xlim = range(x, na.rm = T), ylim = ylim,
                log = "")
        box()
        xmin <- par("usr")[1]
        if (!is.null(labels)) {
                luser <- max(strwidth(labels, "user"), na.rm = TRUE)
                loffset <- luser + xinch(0.1)
                text(rep(xmin - loffset, length(x)), y, labels[o],
                        xpd = TRUE, adj = 0, col = color, ...)
        }
        abline(h = y, lty = "dotted")
        points(x, y, pch = pch, col = color, bg = bg)
        if (!is.null(groups)) {
                gpos <- rev(cumsum(tapply(groups, groups, length) +
                        2) - 1)
                guser <- max(strwidth(glabels, "user"), na.rm = TRUE)
                goffset <- max(luser + xinch(goffset), guser,
                        na.rm = TRUE) + xinch(0.1)
                text(rep(xmin - goffset, nlevels(groups)), gpos,
                        glabels, xpd = TRUE, adj = 0, col = gcolor,
                        ...)
                if (!is.null(gdata)) {
                        abline(h = gpos, lty = "dotted")
                        points(gdata, gpos, pch = gpch, col = gcolor,
                                bg = bg, ...)
                }
        }
        axis(1)
        invisible()
}
dput <- function(x,file="")
        .Internal(dput(x,file))
dump <- function(list, fileout="dumpdata")
.Internal(dump(list, fileout))
dyn.load <- function(x)
{
    x <- as.character(x)
    y <- substr(x, 1, 1)
    if (y == "/") {
        .Internal(dyn.load(x))
    }
    else {
        .Internal(dyn.load(
        paste(system("pwd", intern = T), x, sep = "/", collapse="")))
    }
}
edit <- function(name=NULL, file="", editor=options()$editor)
    .Internal(edit(name,file, editor))
vi <- function(name=NULL, file="") edit(name, file, editor="vi")
emacs <- function(name=NULL, file="") edit(name, file, editor="emacs")
xemacs <- function(name=NULL, file="") edit(name, file, editor="xemacs")
xedit <- function(name=NULL, file="") edit(name, file, editor="xedit")
eigen <- function(x) {
    if(!is.matrix(x) | nrow(x) != ncol(x))
        stop("non-square matrix in eigen")
    n <- nrow(x)
    z <- .C("eigen",
        n,
        as.double(x),
        vectors=matrix(0,n,n),
        values=double(n),
        double(n),
        ierr=integer(1))
    if(z$ierr)
        stop(paste("error code ",z$ierr," in eigen"))
    z[c("values", "vectors")]
}
environment <- function(fun=NULL) .Internal(environment(fun))
.GlobalEnv <- environment()
eval <-
function(expr, envir=sys.frame(sys.parent()))
{
    if(is.expression(expr))
        y<-for(i in 1:length(expr))
            .Internal(eval(expr[[i]],envir))
    else y<-.Internal(eval(expr, envir))
    y
}
quote <- function(x) substitute(x)
Recall <- function(...) .Internal(Recall(...))
exists <- function(x, where=NULL, envir=NULL, frame=NULL, mode="any", inherits=TRUE) {
    if(missing(envir) && !missing(frame))
        envir<-frame
    if(missing(envir) && !missing(where))
        envir<-where
    .Internal(exists(x,envir,mode,inherits))
}
factor <-
function(x, levels=sort(unique(x)), labels, exclude=NA, ordered=FALSE)
{
    if(length(x) == 0) return(character(0))
    if(length(exclude) > 0) {
        exclude <- as.vector(exclude, typeof(x))
        levels <- levels[is.na(match(levels, exclude))]
    }
    x <- .Internal(factor(match(x, levels), length(levels), ordered))
    if(missing(labels)) levels(x) <- levels
    else levels(x) <- labels
    x
}
as.factor <-
function(x, ordered=FALSE)
{
    test <- if(ordered) is.ordered else is.factor
        if(!test(x)) {
                levs <- sort(unique(x))
                x <- .Internal(factor(match(x, levs), length(levs), ordered))
                levels(x) <- levs
        }
        x
}
ordered <-
function(x, levels=sort(unique(x)), labels, exclude=NA, ordered=TRUE)
{
    if(length(exclude) > 0) {
        exclude <- as.vector(exclude, typeof(x))
        levels <- levels[is.na(match(levels, exclude))]
    }
    x <- .Internal(factor(match(x, levels), length(levels), ordered))
    if(missing(labels)) levels(x) <- levels
    else levels(x) <- labels
    x
}
"family" <-
function(x, ...)
UseMethod("family")
"print.family" <-
function(x, ...)
{
    cat("\nFamily:", x$family, "\n")
    cat("Link function:", x$link, "\n\n")
}
"power" <-
function(lambda = 1)
{
    if(lambda <= 0)
        return("log")
    return(lambda)
}
# this function is used with the glm function
# given a link it returns a link function, an inverse link
# function and the derivative dmu/deta
# Written by Simon Davies Dec 1995
## Modified by Thomas Lumley 26 Apr 97
## added valideta(eta) function returning TRUE if all of eta
## is in the domain of linkinv
"make.link" <-
function (link)
{
    recognise <- FALSE
    if (link == "logit") {
        linkfun <- function(mu) log(mu/(1 - mu))
        linkinv <- function(eta) exp(eta)/(1 + exp(eta))
        mu.eta <- function(eta) exp(eta)/(1 + exp(eta))^2
        valideta <- function(eta) TRUE
        recognise <- TRUE
    }
    if (link == "probit") {
        linkfun <- function(mu) qnorm(mu)
        linkinv <- pnorm
        mu.eta <- function(eta) 0.3989422 * exp(-0.5 * eta^2)
        valideta <- function(eta) TRUE
        recognise <- TRUE
    }
    if (link == "cloglog") {
        linkfun <- function(mu) log(-log(1 - mu))
        linkinv <- function(eta) 1 - exp(-exp(eta))
        mu.eta <- function(eta) exp(eta) * exp(-exp(eta))
        valideta <- function(eta) TRUE
        recognise <- TRUE
    }
    if (link == "identity") {
        linkfun <- function(mu) mu
        linkinv <- function(eta) eta
        mu.eta <- function(eta) rep(1, length(eta))
        valideta <- function(eta) TRUE
        recognise <- TRUE
    }
    if (link == "log") {
        linkfun <- function(mu) log(mu)
        linkinv <- function(eta) exp(eta)
        mu.eta <- function(eta) exp(eta)
        valideta <- function(eta) TRUE
        recognise <- TRUE
    }
    if (link == "sqrt") {
        linkfun <- function(mu) mu^0.5
        linkinv <- function(eta) eta^2
        mu.eta <- function(eta) 2 * eta
        valideta <- function(eta) all(eta>0)
        recognise <- TRUE
    }
    if (link == "1/mu^2") {
        linkfun <- function(mu) 1/mu^2
        linkinv <- function(eta) 1/eta^0.5
        mu.eta <- function(eta) -1/(2 * eta^1.5)
        valideta <- function(eta) all(eta>0)
        recognise <- TRUE
    }
    if (link == "inverse") {
        linkfun <- function(mu) 1/mu
        linkinv <- function(eta) 1/eta
        mu.eta <- function(eta) -1/(eta^2)
        valideta <- function(eta) all(eta!=0)
        recognise <- TRUE
    }
    if (!is.na(as.numeric(link))) {
        lambda <- as.numeric(link)
        linkfun <- function(mu) mu^lambda
        linkinv <- function(eta) eta^(1/lambda)
        mu.eta <- function(eta) (1/lambda) * eta^(1/lambda - 1)
        valideta <- function(eta) all(eta>0)
        recognise <- TRUE
    }
    if (!recognise)
        stop(paste(link, "link not recognised"))
    return(list(linkfun = linkfun,
        linkinv = linkinv,
        mu.eta = mu.eta,
        valideta=valideta))
}
"poisson" <-
function (link = "log")
{
    linktemp <- substitute(link)
    #this is a function used in the glm function
    #it holds everything personal to the family
    #converts link into character string
    if (!is.character(linktemp)) {
        linktemp <- deparse(linktemp)
        if (linktemp == "link")
            linktemp <- eval(link)
    }
    if (any(linktemp == c("log", "identity", "sqrt")))
        stats <- make.link(linktemp)
    else stop(paste(linktemp, "link not available for poisson",
            "family, available links are \"identity\", ",
            "\"log\" and \"sqrt\""))
    variance <- function(mu) mu
    validmu <- function(mu) all(mu>0)
    dev.resids <- function(y, mu, wt)
        2 * wt * (y * log(ifelse(y == 0, 1, y/mu)) - (y - mu))
    initialize <- expression({
        if (any(y < 0))
            stop(paste("Negative values not allowed for",
                "the Poisson family"))
        n <- rep(1, nobs)
        mustart <- y + 0.1
    })
    family <- list(family = "poisson",
            link = linktemp,
            linkfun = stats$linkfun,
            linkinv = stats$linkinv,
            variance = variance,
            dev.resids = dev.resids,
            mu.eta = stats$mu.eta,
            initialize = initialize,
            validmu = validmu,
            valideta = stats$valideta)
    class(family) <- "family"
    return(family)
}
"gaussian" <-
function ()
{
    stats <- make.link("identity")
    # this is a function used in the glm function
    # it holds everything personal to the family
    variance <- function(mu) rep(1, length(mu))
    dev.resids <- function(y, mu, wt) wt * ((y - mu)^2)
    initialize <- expression({
        n <- rep(1, nobs)
        mustart <- y
    })
    validmu <- function(mu) TRUE
    family <- list(family = "gaussian",
            link = "identity",
            linkfun = stats$linkfun,
            linkinv = stats$linkinv,
            variance = variance,
            dev.resids = dev.resids,
            mu.eta = stats$mu.eta,
            initialize = initialize,
            validmu = validmu,
            valideta = stats$valideta)
    class(family) <- "family"
    return(family)
}
"binomial" <-
function (link = "logit")
{
    linktemp <- substitute(link)
    # this is a function used in the glm function
    # it holds everything personal to the family
    # converts link into character string
    if (!is.character(linktemp)) {
        linktemp <- deparse(linktemp)
        if (linktemp == "link")
            linktemp <- eval(link)
    }
    if (any(linktemp == c("logit", "probit", "cloglog")))
        stats <- make.link(linktemp)
    else stop(paste(linktemp, "link not available for binomial",
        "family, available links are \"logit\", ",
        "\"probit\" and \"cloglog\""))
    variance <- function(mu) mu * (1 - mu)
    validmu <- function(mu) all(mu>0) && all(mu<1)
    dev.resids <- function(y, mu, wt)
        2 * wt * (y * log(ifelse(y == 0, 1, y/mu)) +
        (1 - y) * log(ifelse(y == 1, 1, (1 - y)/(1 - mu))))
    initialize <- expression({
        if (NCOL(y) == 1) {
            n <- rep(1, nobs)
            if (any(y < 0 | y > 1))
                stop("y values must be 0 <= y <= 1")
        }
        else if (NCOL(y) == 2) {
            n <- y[, 1] + y[, 2]
            y <- y[, 1]/n
            weights <- weights * n
        }
        else stop(paste("For the binomial family, y must be",
            "a vector of 0 and 1\'s or a 2 column",
            "matrix where col 1 is no. successes",
            "and col 2 is no. failures"))
        mustart <- (n * y + 0.5)/(n + 1)
    })
    family <- list(family = "binomial",
            link = linktemp,
            linkfun = stats$linkfun,
            linkinv = stats$linkinv,
            variance = variance,
            dev.resids = dev.resids,
            mu.eta = stats$mu.eta,
            initialize = initialize,
            validmu = validmu,
            valideta = stats$valideta)
    class(family) <- "family"
    return(family)
}
"Gamma" <-
function (link = "inverse")
{
    linktemp <- substitute(link)
    #this is a function used in the glm function
    #it holds everything personal to the family
    #converts link into character string
    if (!is.character(linktemp)) {
        linktemp <- deparse(linktemp)
        if (linktemp == "link")
            linktemp <- eval(link)
    }
    if (any(linktemp == c("inverse", "log", "identity")))
        stats <- make.link(linktemp)
    else stop(paste(linktemp, "link not available for gamma",
        "family, available links are \"inverse\", ",
        "\"log\" and \"identity\""))
    variance <- function(mu) mu^2
    validmu <- function(mu) all(mu>0)
    dev.resids <- function(y, mu, wt)
        -2 * wt * (log(ifelse(y == 0, 1, y/mu)) - (y - mu)/mu)
    initialize <- expression({
        if (any(y < 0))
            stop(paste("Negative values not",
                "allowed for the gamma family"))
        n <- rep(1, nobs)
        mustart <- y
    })
    family <- list(family = "Gamma",
            link = linktemp,
            linkfun = stats$linkfun,
            linkinv = stats$linkinv,
            variance = variance,
            dev.resids = dev.resids,
            mu.eta = stats$mu.eta,
            initialize = initialize,
            validmu = validmu,
            valideta = stats$valideta)
    class(family) <- "family"
    return(family)
}
"inverse.gaussian" <-
function()
{
    stats <- make.link("1/mu^2")
    variance <- function(mu) mu^3
    dev.resids <- function(y, mu, wt)  wt*((y - mu)^2)/(y*mu^2)
    initialize <- expression({
            if(any(y <= 0))
                stop(paste("Positive values only allowed for",
                    "the inverse.gaussian family"))
            n <- rep(1, nobs)
            mustart <- y
            })
    validmu <- function(mu) TRUE
    family <- list(family = "inverse.gaussian",
            link = "1/mu^2",
            linkfun = stats$linkfun,
            linkinv = stats$linkinv,
            variance = variance,
            dev.resids = dev.resids,
            mu.eta = stats$mu.eta,
            initialize = initialize,
            validmu = validmu,
            valideta = stats$valideta)
    class(family) <- "family"
    return(family)
}
"quasi" <-
function (link = "identity", variance = "constant")
{
    linktemp <- substitute(link)
    #this is a function used in the glm function
    #it holds everything personal to the family
    #converts link into character string
    if (is.expression(linktemp))
        linktemp <- eval(linktemp)
    if (!is.character(linktemp)) {
        linktemp <- deparse(linktemp)
        if (linktemp == "link")
            linktemp <- eval(link)
    }
    stats <- make.link(linktemp)
    #converts variance into character string
    variancetemp <- substitute(variance)
    if (!is.character(variancetemp)) {
        variancetemp <- deparse(variancetemp)
        if (linktemp == "variance")
            variancetemp <- eval(variance)
    }
    if (!any(variancetemp == c("mu(1-mu)",
        "mu", "mu^2", "mu^3", "constant")))
        stop(paste(variancetemp, "not recognised, possible variances",
            "are \"mu(1-mu)\", \"mu\", \"mu^2\", \"mu^3\" and",
            "\"constant\""))
    if (variancetemp == "constant") {
        variance <- function(mu) rep(1, length(mu))
        dev.resids <- function(y, mu, wt) wt * ((y - mu)^2)
        validmu <-function(mu) TRUE
    }
    if (variancetemp == "mu(1-mu)") {
        variance <- function(mu) mu * (1 - mu)
        validmu <-function(mu) all(mu>0) && all(mu<1)
        dev.resids <- function(y, mu, wt)
            2 * wt * (y * log(ifelse(y == 0, 1,
            y/mu)) + (1 - y) * log(ifelse(y == 1, 1, (1 - y)/(1 - mu))))
    }
    if (variancetemp == "mu") {
        variance <- function(mu) mu
        validmu<-function(mu) all(mu>0)
        dev.resids <- function(y, mu, wt)
            2 * wt * (y * log(ifelse(y == 0, 1, y/mu)) - (y - mu))
    }
    if (variancetemp == "mu^2") {
        variance <- function(mu) mu^2
        validmu<-function(mu) all(mu!=0)
        dev.resids <- function(y, mu, wt)
            -2 * wt * (log(y/mu) - (y - mu)/mu)
    }
    if (variancetemp == "mu^3") {
        variance <- function(mu) mu^3
        validmu <-function(mu) all(mu>0)
        dev.resids <- function(y, mu, wt) wt * ((y - mu)^2)/(y * mu^2)
    }
    initialize <- expression({
        n <- rep(1, nobs)
        mustart <- y
    })
    family <- list(family = "quasi",
            link = linktemp,
            linkfun = stats$linkfun,
            linkinv = stats$linkinv,
            variance = variance,
            dev.resids = dev.resids,
            mu.eta = stats$mu.eta,
            initialize = initialize,
            validmu = validmu,
            valideta = stats$valideta)
    class(family) <- "family"
    return(family)
}
fft <- function(z, inverse=FALSE)
.Internal(fft(z, inverse))
mvfft <- function(z, inverse=FALSE)
.Internal(mvfft(z, inverse))
nextn <- function(n, factors=c(2,3,5))
.Internal(nextn(n, factors))
convolve <-
function(x, y, conj=T) {
    if(length(x) != length(y))
        stop("length mismatch in convolution")
    if(conj)
        Re(fft(fft(x)*Conj(fft(y)),inv=T))/length(x)
    else
        Re(fft(fft(x)*fft(y),inv=T))/length(x)
}
fivenum <- function(x, na.rm=TRUE)
{
    xna <- is.na(x)
    if(na.rm) x <- x[!xna]
    else if(any(xna)) return(rep(NA,5))
    x <- sort(x)
    n <- length(x)
    if(n == 0) rep(NA,5)
    else {
        d <- c(1, 0.5*floor(0.5*(n+3)), 0.5*(n+1),
            n+1-0.5*floor(0.5*(n+3)), n)
        0.5*(x[floor(d)]+x[ceiling(d)])
    }
}
fix <- function(x) {
    subx <- substitute(x)
    if( is.name(subx) )
        subx<-deparse(subx)
    if (!is.character(subx) || length(subx) != 1)
        stop("fix requires a name")
    if(exists(subx, inherits=TRUE))
        x <- edit(get(subx))
    else
        stop(paste("no object named \"", subx, "\" to edit",sep=""))
    assign(subx, x, .GlobalEnv)
}
formals <- function(fun=sys.function(sys.parent())) {
    if(is.character(fun))
        fun <- get(fun, mode = "function")
    .Internal(formals(fun))
}
body <- function(fun=sys.function(sys.parent())) {
    if(is.character(fun))
        fun <- get(fun, mode = "function")
    .Internal(body(fun))
}
"formatC" <-
function (x, digits=NULL, width=max(0, digits) + 1, format=NULL,
    flag="", mode=NULL)
{
    # Copyright (C) Martin Maechler, 1994
    bl.string <- function(no) paste(rep(" ", no), collapse = "")
    if (is.null(x)) return("")
    n <- length(x)
    if (missing(mode))
        mode <- storage.mode(x)
    else if (any(mode == c("real", "integer")))
        storage.mode(x) <- mode
    else stop("\"mode\" must be \"real\" or \"integer\"")
    if (mode == "character" || (!is.null(format) && format == "s")) {
        if (mode != "character") {
            warning("should give \"character\" argument for format=\"s\" -- COERCE")
            x <- as.character(x)
        }
        nc <- nchar(x)
        if (width < 0) {
            flag <- "-"
            width <- -width
        }
        pad <- sapply(pmax(0, width - nc), bl.string)
        ## for R <= 0.16 (incompatibility to S):
        if(is.list(pad)) pad <- ""
        if (flag == "-")
            return(paste(x, pad, sep = ""))
        else    return(paste(pad, x, sep = ""))
    }
    some.special <- !all(Ok <- is.finite(x))
    if (some.special) {
        nQ <- nchar(rQ <- as.character(x[!Ok]))
        nX <- pmax(width - nQ, 0)
        #-- number of characters to add
        x[!Ok] <- 0
    }
    if (missing(format) || is.null(format))
        format <- if (mode == "integer")
            "d"
        else "g"
    else {
        if (any(format == c("f", "e", "E", "g", "G"))) {
            if (mode == "integer")
                mode <- storage.mode(x) <- "single"
        }
        else if (format == "d") {
            if (mode != "integer")
                mode <- storage.mode(x) <- "integer"
        }
        else stop("\"format\" must be in {\"f\",\"e\",\"E\",\"g\",\"G\", \"s\"}" )
    }
    if (missing(digits) || is.null(digits))
        digits <- if (mode == "integer")
            2
        else 4
    if (width == 0)
        stop("\"width\" must not be 0")
    r <- .C("str_signif",
        x = x,
        n = n,
        mode = as.character(mode),
        width = as.integer(width),
        digits = as.integer(digits),
        format = as.character(format),
        flag = as.character(flag),
        result = rep(bl.string(abs(width) + 10), n)
    )$result
    if (some.special)
        r[!Ok] <- rQ
    if (!is.null(x.atr <- attributes(x)))
        attributes(r) <- x.atr
    r
}
get <- function(x, envir=NULL, mode="any", inherits=TRUE)
.Internal(get(x,envir,mode,inherits))
# gl function of glim
gl <- function (n, k, length, labels=1:n, ordered=FALSE)
    factor(rep(rep(1:n,rep(k,n)), length=length),
        labels=labels, ordered=ordered)
# This function fits a generalized linear model via
# iteratively reweighted least squares for any family.
# Written by Simon Davies, Dec 1995
# glm.fit modified by Thomas Lumley, Apr 1997
glm <- function(formula, family=gaussian, data, weights=NULL,
    subset=NULL, na.action=na.fail, start=NULL, offset=NULL,
    control=glm.control(epsilon=0.0001, maxit=10, trace=FALSE),
    model=TRUE, method=glm.fit, x=FALSE, y=TRUE)
{
    call <- sys.call()
    # family
    if(is.character(family)) family <- get(family)
    if(is.function(family)) family <- family()
    if(is.null(family$family)) stop("family not recognised")
    # extract x, y, etc from the model formula and frame
    mt <- terms(formula)
    if(missing(data)) data <- sys.frame(sys.parent())
        mf <- match.call()
        mf$family <- mf$start <- mf$control <- mf$maxit <- NULL
    mf$model <- mf$method <- mf$x <- mf$y <- NULL
        mf$use.data <- TRUE
        mf[[1]] <- as.name("model.frame")
        mf <- eval(mf, sys.frame(sys.parent()))
    X <- model.matrix(mt, mf)
    Y <- model.response(mf, "numeric")
    weights <- model.weights(mf)
    offset <- model.offset(mf)
    # check weights and offset
    if( !is.null(weights) && any(weights<0) )
        stop("Negative wts not allowed")
    if(!is.null(offset) && length(offset) != NROW(Y))
        stop(paste("Number of offsets is", length(offset),
            ", should equal", NROW(Y), "(number of observations)"))
    # fit model via iterative reweighted least squares
    fit <- method(x=X, y=Y, weights=weights, start=start,
        offset=offset, family=family, control=control)
    # output
    if(model) fit$model <- mf
    if(x) fit$x <- X
    else fit$x <- NULL
    if(!y) fit$y <- NULL
    fit <- c(fit, list(call=call, formula=formula, x=x, terms=mt, data=data,
        offset=offset, control=control, method=method))
    class(fit) <- c("glm", "lm")
    return(fit)
}
glm.control <- function(epsilon = 0.0001, maxit = 10, trace = FALSE)
{
    if(epsilon <= 0)
        stop("value of epsilon must be > 0")
    if(maxit <= 0)
        stop("maximum number of iterations must be > 0")
    return(list(epsilon = epsilon, maxit = maxit, trace = trace))
}
## Modified by Thomas Lumley 26 Apr 97
## Added boundary checks and step halving
## Modified detection of fitted 0/1 in binomial
"glm.fit" <-
function (x, y, weights = rep(1, nobs), start = NULL, offset = rep(0, nobs),
family = gaussian(), control = glm.control(), intercept = TRUE)
{
    xnames <- dimnames(x)[[2]]
    ynames <- names(y)
    conv <- FALSE
    nobs <- NROW(y)
    nvars <- NCOL(x)
    # define weights and offset if needed
    if (is.null(weights))
        weights <- rep(1, nobs)
    if (is.null(offset))
        offset <- rep(0, nobs)
    # get family functions
    variance <- family$variance
    dev.resids <- family$dev.resids
    linkinv <- family$linkinv
    mu.eta <- family$mu.eta
    valideta<-family$valideta
    if (is.null(valideta)) valideta<-function(eta) TRUE
    validmu<-family$validmu
    if (is.null(validmu)) validmu<-function(mu) TRUE
    eval(family$initialize, sys.frame(sys.nframe()))
    if (NCOL(y) > 1)
        stop("y must be univariate unless binomial")
    # calculates the first estimate of eta and mu
    if (is.null(start)) {
      start<-c(0.5,rep(0,nvars-1))
      linkfun <- family$linkfun
      if (validmu(mustart)){
        etastart <- linkfun(mustart)
        if (valideta(etastart)){
          z <- etastart + (y - mustart)/mu.eta(etastart) - offset
          w <- ((weights * mu.eta(etastart)^2)/variance(mustart))^0.5
          fit <- qr(x * w)
          start <- qr.coef(fit, w * z)
          start[is.na(start)] <- 0
        }
      }
    }
      else if (length(start) != nvars)
        stop(paste("Length of start should equal", nvars, "and correspond to initial coefs for", deparse(xnames)))
    if (NCOL(x) == 1)
      eta <- as.vector(x * start)
      else eta <- as.vector(x %*% start)
    mu <- linkinv(eta + offset)
    if (!(validmu(mu) && valideta(eta)))
      stop("Can't find valid starting values: please specify with start=")
    # calculate initial deviance and coefficient
    devold <- sum(dev.resids(y, mu, weights))
    coefold <- start
    boundary<-FALSE
    # do iteration
    for (iter in 1:control$maxit) {
        mu.eta.val <- mu.eta(eta + offset)
        if (any(is.na(mu.eta.val))){
          mu.eta.val[is.na(mu.eta.val)]<-mu.eta(mu)[is.na(mu.eta.val)]
        }
        if (any(is.na(mu.eta.val)))
          stop("NAs in d(mu)/d(eta)")
        # calculate z and w using only values where mu.eta != 0
        good <- (mu.eta.val != 0)
        if (all(!good)) {
            conv <- FALSE
            warning("No observations informative at iteration",iter)
            break
        }
        z <- eta[good] + (y - mu)[good]/mu.eta.val[good]
        w <- ((weights * mu.eta.val^2)[good]/variance(mu)[good])^0.5
        x <- as.matrix(x)
        ngoodobs <- as.integer(nobs - sum(!good))
        ncols <- as.integer(1)
        # call linpack code
        fit <- .Fortran("dqrls",
            qr = x[good, ] * w,
            n = as.integer(ngoodobs),
            p = nvars,
            y = w * z,
            ny = ncols,
            tol = 1e-07,
            coefficients = mat.or.vec(nvars, 1),
            residuals = mat.or.vec(ngoodobs, 1),
            effects = mat.or.vec(ngoodobs, 1),
            rank = integer(1),
            pivot = as.integer(1:nvars),
            qraux = double(nvars),
            work = double(2 * nvars)
        )
        # stop if not enough parameters
        if (nobs < fit$rank)
            stop(paste("X matrix has rank", fit$rank, "but only", nobs, "observations"
            ))
        # if X matrix was not full rank then columns
        # were pivoted, hence we need to re-label the names
        if (fit$rank != nvars) {
            xnames <- xnames[fit$pivot]
            dimnames(fit$qr) <- list(NULL, xnames)
        }
        # calculate updated values of eta and mu with the new coef
        start <- coef <- fit$coefficients
        start[fit$pivot] <- coef
        if (nvars == 1)
            eta[good] <- x[good] * start
        else eta[good] <- as.vector(x[good, ] %*% start)
        mu <- linkinv(eta + offset)
        if (family$family == "binomial") {
            if (any(mu == 1) || any(mu == 0))
                warning("fitted probabilities of 0 or 1 occured")
            mu[mu == 1] <- 1 - 0.5 * control$epsilon/length(mu)
            mu[mu == 0] <- 0.5 * control$epsilon/length(mu)
        }
        if (family$family == "poisson") {
            if (any(mu == 0))
                warning("fitted rates of 0 occured")
            mu[mu == 0] <- 0.5 * control$epsilon/length(mu)^2
        }
        dev <- sum(dev.resids(y, mu, weights))
        if (control$trace)
            cat("Deviance =", dev, "Iterations -", iter, "\n")
        # check for divergence
        boundary<-FALSE
        if (any(is.na(dev)) || any(is.na(coef))) {
            warning("Step size truncated due to divergence")
            ii<-1
            while( (any(is.na(dev)) || any(is.na(start)))){
              if (ii>control$maxit) stop("Can't correct")
              ii<-ii+1
              start<-(start+coefold)/2
              if (nvars == 1)
                eta[good] <- x[good] * start
                else eta[good] <- as.vector(x[good, ] %*% start)
              mu <- linkinv(eta + offset)
              dev <- sum(dev.resids(y, mu, weights))
            }
            boundary<-TRUE
            coef<-start
            if (control$trace)
              cat("New Deviance =", dev, "\n")
              }
        ## check for fitted values outside domain.
        if (!(valideta(eta) && validmu(mu))) {
            warning("Step size truncated: out of bounds.")
            ii<-1
            while(!(valideta(eta) && validmu(mu))){
              if (ii>control$maxit) stop("Can't correct step size.")
              ii<-ii+1
              start<-(start+coefold)/2
              if (nvars == 1)
                eta[good] <- x[good] * start
                else eta[good] <- as.vector(x[good, ] %*% start)
              mu <- linkinv(eta + offset)
            }
            boundary<-TRUE
            coef<-start
            dev <- sum(dev.resids(y, mu, weights))
            if (control$trace)
              cat("New Deviance =", dev, "\n")
              }
        #check for convergence
        if (abs(dev - devold)/(0.1 + abs(dev)) < control$epsilon) {
          conv <- TRUE
          break
        }
          else {
            devold <- dev
            coefold <- coef
          }
          }
    if (!conv)
      warning("Algorithm did not converge")
    if (boundary)
      warning("Algorithm stopped at boundary value")
                    # calculate residuals
    residuals <- rep(NA, nobs)
    ##     residuals[good] <- z - eta
    residuals[good]<-z-eta[good]
                    # name output
    fit$qr <- as.matrix(fit$qr)
    Rmat <- fit$qr[1:nvars, 1:nvars]
    Rmat <- as.matrix(Rmat)
    Rmat[row(Rmat) > col(Rmat)] <- 0
    names(coef) <- xnames
    colnames(fit$qr) <- xnames
    dimnames(Rmat) <- list(xnames, xnames)
    names(residuals) <- ynames
    names(mu) <- ynames
    names(eta) <- ynames
    names(w) <- ynames
    names(weights) <- ynames
    names(y) <- ynames
                    # calculate null deviance
    if (intercept)
      wtdmu <- sum(weights * y)/sum(weights)
      else wtdmu <- linkinv(offset)
    nulldev <- sum(dev.resids(y, wtdmu, weights))
                    # calculate df
    nulldf <- nobs - as.numeric(intercept)
    resdf <- nobs - fit$rank - sum(weights == 0)
    qr <- list(qr = fit$qr, rank = fit$rank, qraux = fit$qraux)
    return(list(coefficients = coef, residuals = residuals, fitted.values = mu
            , effects = fit$effects, R = Rmat, rank = fit$rank, qr = qr, family = family
            , linear.predictors = eta, deviance = dev, null.deviance = nulldev, iter = iter
            , weights = w^2, prior.weights = weights, df.residual = resdf, df.null = nulldf
            , y = y,converged=conv,boundary=boundary))
}
print.glm <- function (x, digits = max(3, .Options$digits - 3),
    na.print="", ...)
{
    cat("\nCall: ", deparse(x$call), "\n\n")
    cat("Coefficients:\n")
    print.default(round(x$coefficients, digits), print.gap = 2)
    cat("\nDegrees of Freedom:", length(x$residuals), "Total;",
         x$df.residual, "Residual\n")
    cat("Null Deviance:", format(signif(x$null.deviance, digits)), "\n")
    cat("Residual Deviance:", format(signif(x$deviance, digits)), "\n")
    invisible(x)
}
anova.glm <- function(object, ..., test=NULL, na.action=na.omit)
{
    # check for multiple objects
    args <- function(...) nargs()
    if(args(...)) return(anova.glmlist(list(object, ...), test=test))
    # extract variables from model
    varlist <- attr(object$terms, "variables")
    if(!is.null(object$x) && !(is.logical(object$x) ||
        object$x==FALSE)) x <- object$x
    else {
        if(is.null(object$model)) {
            if(is.null(object$data))
                object$data <- sys.frame(sys.parent())
            object$model <- na.action(
                model.frame(eval(varlist, object$data),
                    as.character(varlist[-1]), NULL))
        }
        x <- model.matrix(object$terms, object$model)
    }
    varseq <- attr(x, "assign")
    nvars <- max(varseq)
    resdev <- resdf <- NULL
    # if there is more than one explanatory variable then
    # recall glm.fit to fit variables sequentially
    if(nvars > 1) {
        for(i in 1:(nvars-1)) {
            # explanatory variables up to i are kept in the model
            tempx <- x[, varseq <= i]
            # use method from glm to find residual deviance
            # and df for each sequential fit
            method <- object$method
            fit <- method(x=tempx, y=object$y,
                weights=object$prior.weights,
                start=object$start,
                offset=object$offset,
                family=object$family,
                control=object$control)
            resdev <- c(resdev, fit$deviance)
            resdf <- c(resdf, fit$df.residual)
        }
    }
    # add values from null and full model
    resdf <- c(object$df.null, resdf, object$df.residual)
    resdev <- c(object$null.deviance, resdev, object$deviance)
    # construct table and title
    table <- cbind(c(NA, -diff(resdf)), c(NA, -diff(resdev)), resdf, resdev)
    dimnames(table) <- list(c("NULL", attr(object$terms, "term.labels")),
                c("Df", "Deviance", "Resid. Df", "Resid. Dev"))
    title <- paste("Analysis of Deviance Table", "\n\nModel: ",
        object$family$family, ", link: ", object$family$link,
        "\n\nResponse: ", as.character(varlist[-1])[1],
        "\n\nTerms added sequentially (first to last)\n\n", sep="")
    # calculate test statistics if needed
    if(!is.null(test))
        table <- stat.anova(table=table, test=test, scale=sum(
            object$weights*object$residuals^2)/object$df.residual,
            df.scale=object$df.residual, n=NROW(x))
    # return output
    output <- list(title=title, table=table)
    class(output) <- "anova.glm"
    return(output)
}
anova.glmlist <- function(object, ..., test=NULL)
{
    # find responses for all models and remove
    # any models with a different response
    responses <- as.character(lapply(object, function(x) {
            as.character(x$formula[2])} ))
    sameresp <- responses==responses[1]
    if(!all(sameresp)) {
        object <- object[sameresp]
        warning(paste("Models with response", deparse(responses[
            !sameresp]), "removed because response differs from",
            "model 1"))
    }
    # calculate the number of models
    nmodels <- length(object)
    if(nmodels==1)    return(anova.glm(object[[1]], ..., test))
    # extract statistics
    resdf <- as.numeric(lapply(object, function(x) x$df.residual))
    resdev <- as.numeric(lapply(object, function(x) x$deviance))
    # construct table and title
    table <- cbind(resdf, resdev, c(NA, -diff(resdf)), c(NA, -diff(resdev)))
    variables <- as.character(lapply(object, function(x) {
            as.character(x$formula[3])} ))
    dimnames(table) <- list(variables, c("Resid. Df", "Resid. Dev", "Df",
                "Deviance"))
    title <- paste("Analysis of Deviance Table \n\nResponse: ", responses[1],
            "\n\n", sep="")
    # calculate test statistic if needed
    if(!is.null(test)) {
        bigmodel <- object[[(order(resdf)[1])]]
        table <- stat.anova(table=table, test=test, scale=sum(
            bigmodel$weights * bigmodel$residuals^2)/
            bigmodel$df.residual, df.scale=min(resdf),
            n=length(bigmodel$residuals))
    }
    # return output
    output <- list(table=table, title=title)
    class(output) <- "anova.glm"
    return(output)
}
stat.anova <- function(table, test, scale, df.scale, n)
{
    testnum <- match(test, c("Chisq", "F", "Cp"))
    cnames <- colnames(table)
    rnames <- rownames(table)
    if(is.na(testnum))
        stop(paste("Test \"", test, "\" not recognised", sep=""))
    if(testnum==1) {
        chisq <- 1-pchisq(abs(table[, "Deviance"]), abs(table[, "Df"]))
        table <- cbind(table, chisq)
        dimnames(table) <- list(rnames, c(cnames, "P(>|Chi|)"))
    }
    if(testnum==2) {
        Fvalue <- abs((table[, "Deviance"]/table[, "Df"])/scale)
        pvalue <- 1-pf(Fvalue, abs(table[, "Df"]), abs(df.scale))
        table <- cbind(table, Fvalue, pvalue)
        dimnames(table) <- list(rnames, c(cnames, "F", "Pr(>F)"))
    }
    if(testnum==3) {
        Cp <- table[, "Resid. Dev"] + 2*scale*(n - table[, "Resid. Df"])
        table <- cbind(table, Cp)
        dimnames(table) <- list(rnames, c(cnames, "Cp"))
    }
    return(table)
}
summary.glm <- function(object, dispersion = NULL,
    correlation = TRUE, na.action=na.omit)
{
    # calculate dispersion if needed
    if(is.null(dispersion)) {
        if(any(object$family$family == c("poisson", "binomial")))
            dispersion <- 1
        else {
            if(any(object$weights==0))
                warning(paste("observations with zero weight",
                "not used for calculating dispersion"))
            dispersion <- sum(object$weights*object$residuals^2)/
                    object$df.residual
        }
    }
    # extract x to get column names
    if(is.null(object$x)) {
        if(is.null(object$model)) {
            varlist <- attr(object$terms, "variables")
            if(is.null(object$data))
                object$data <- sys.frame(sys.parent())
            object$model <- na.action(model.frame(eval(varlist,
                object$data), as.character(varlist[-1]), NULL))
        }
        object$x <- model.matrix(object$terms, object$model)
    }
    # calculate scaled and unscaled covariance matrix
    p <- object$rank
    p1 <- 1:p
    covmat.unscaled <- chol2inv(object$qr$qr[p1,p1,drop=FALSE])
    dimnames(covmat.unscaled) <- list(names(object$coefficients)[p1], names(object$coefficients)[p1])
    covmat <- dispersion*covmat.unscaled
    dimnames(covmat) <- dimnames(covmat.unscaled)
    # calculate coef table
    nas <- is.na(object$coefficients)
    tvalue <- object$coefficients[p1]/diag(covmat)^0.5
    pvalue <- 2*(1-pt(abs(tvalue), object$df.residual))
    coef.table <- cbind(object$coefficients[p1], diag(covmat)^0.5, tvalue, pvalue)
    dimnames(coef.table) <- list(names(object$coefficients)[p1], c("Value", "Std.error",
                "t value", "P(>|t|)"))
    # return answer
    ans <- list(
        call=object$call,
        terms=object$terms,
        family=object$family,
        deviance.resid=residuals(object, type = "deviance"),
        coefficients=coef.table,
        dispersion=dispersion,
        df=c(object$rank, object$df.residual),
        deviance=object$deviance,
        df.residual=object$df.residual,
        null.deviance=object$null.deviance,
        df.null=object$df.null,
        iter=object$iter,
        cov.unscaled=covmat.unscaled,
        cov.scaled=covmat)
#        cov.scaled=covmat,
#        nas=nas)
    if(correlation) {
        ans$correlation <- as.matrix(covmat/(outer(diag(covmat),
            diag(covmat))^0.5))
    }
    class(ans) <- "summary.glm"
    return(ans)
}
print.summary.glm <- function (x, digits = max(3, .Options$digits - 3),
    na.print="", ...)
{
        cat("\nCall:\n")
        cat(paste(deparse(x$call), sep="\n", collapse="\n"), "\n\n", sep="")
    cat("Deviance Residuals: \n")
    if(x$df.residual > 5) {
        x$deviance.resid <- quantile(x$deviance.resid)
        names(x$deviance.resid) <- c("Min", "1Q", "Median", "3Q", "Max")
    }
    print.default(x$deviance.resid, digits=digits, na = "", print.gap = 2)
    cat("\nCoefficients:\n")
    print.default(x$coefficients, digits=digits, print.gap = 2)
    cat(paste("\n(Dispersion parameter for ", x$family$family,
        " family taken to be ", x$dispersion,
        ")\n\n    Null deviance: ", x$null.deviance,
        " on ", x$df.null, " degrees of freedom\n\n",
        "Residual deviance: ", x$deviance,
        " on ", x$df.residual, " degrees of freedom\n\n",
        "Number of Fisher Scoring iterations: ", x$iter,
        "\n\n", sep=""))
    correl <- x$correlation
    if(!is.null(correl)) {
        p <- dim(correl)[2]
        if(p > 1) {
            cat("Correlation of Coefficients:\n")
            correl[!lower.tri(correl)] <- NA
            print(correl[-1, -NCOL(correl), drop=FALSE], digits=digits,
                na="")
        }
        cat("\n")
    }
    invisible(x)
}
print.anova.glm <- function(x, digits = max(3, .Options$digits - 3),
    na.print = "", ...)
{
    cat("\n", x$title, sep="")
    print.default(x$table, digits=digits, na = "", print.gap = 2)
    cat("\n")
}
# Generic Functions
coefficients.glm <- function(object)
    object$coefficients
deviance.glm <- function(object)
    object$deviance
effects.glm <- function(object)
    object$effects
family.glm <- function(object) {
    family <- get(as.character(object$family$family), mode="function")
    family()
}
fitted.values.glm <- function(object)
    object$fitted.values
residuals.glm <- function(object, type="deviance")
{
    type <- match(type, c("deviance", "pearson", "working", "response"))
    y <- object$y
    mu <- object$fitted.values
    wts <- object$prior.weights
    switch(type,
        deviance={
            dev.resids <- object$family$dev.resids
            ifelse(y>mu, dev.resids(y, mu, wts)^0.5, -(dev.resids(
                y, mu, wts)^0.5))
        },
        pearson=object$residuals * object$weights^0.5,
        working=object$residuals,
        response=y - mu
        )
}
grep <-
function(pattern, x, ignore.case=FALSE, extended=TRUE, value=FALSE)
{
    .Internal(grep(pattern, x, ignore.case, extended, value))
}
sub <-
function(pattern, replacement, x, ignore.case=FALSE, extended=TRUE)
{
    .Internal(sub(pattern, replacement, x, ignore.case, extended))
}
gsub <-
function(pattern, replacement, x, ignore.case=FALSE, extended=TRUE)
{
    .Internal(gsub(pattern, replacement, x, ignore.case, extended))
}
"grid" <-
function (nx=3, ny=3, col="lightgray", lty="dotted")
{
    lims <- par("usr")
    if (nx > 1) {
        coord <- seq(lims[1], lims[2], len = nx + 2)[c(-1, -(nx + 2))]
        abline(v = coord, col = col, lty = lty)
    }
    if (ny > 1) {
        coord <- seq(lims[3], lims[4], len = ny + 2)[c(-1, -(ny + 2))]
        abline(h = coord, col = col, lty = lty)
    }
}
"help" <-
function (topic, data, library)
{
    if(!missing(topic)) {
        topic <- substitute(topic)
        if (is.character(topic) || is.name(topic)) {
            if (!is.character(topic))
                topic <- deparse(topic)
            if (!is.na(match(topic, c("+", "-", "*", "/", "^", "%%"))))
                topic <- "Arithmetic"
            else if (!is.na(match(topic, c("<", ">", "<=", ">=", "==", "!="))))
                topic <- "Comparison"
            else if (!is.na(match(topic, c("&", "&&", "|", "||", "!"))))
                topic <- "Logic"
            else if (!is.na(match(topic, c("[", "[[", "$"))))
                topic <- "Extract"
            system(paste("$RHOME/cmd/help", topic,
                paste(.Libraries, collapse = " "), "base"))
        }
        else {
            topic <- as.character(topic)
            if (topic[1] == "data") {
                file <- system.file("data", paste(topic[2], ".doc", sep = ""))
                if (file == "")
                    stop(paste("no documentation for dataset", topic[2]))
                else system(paste("$RHOME/cmd/pager", file))
            }
            else if (topic[1] == "library") {
                file <- system.file("help", paste(topic[2], "/INDEX", sep = ""))
                if (file == "")
                    stop(paste("no documentation for dataset", topic[2]))
                else system(paste("$RHOME/cmd/pager", file))
            }
            else stop("unimplemented help feature")
        }
    }
    else if(!missing(data)) {
        topic <- as.character(substitute(data))
        file <- system.file("data", paste(topic, ".doc", sep = ""))
        if (file == "")
            stop(paste("no documentation for dataset", topic))
        else system(paste("$RHOME/cmd/pager", file))
    }
    else if(!missing(library)) {
        topic <- as.character(substitute(library))
        file <- system.file("help", paste(topic, "/INDEX", sep = ""))
        if (file == "")
            stop(paste("no documentation for library", topic))
        else system(paste("$RHOME/cmd/pager", file))
    }
    else system("$RHOME/cmd/help help base")
}
# help.start <- function (gui = "irrelevant", browser="netscape -mono")
# {
#         system(paste(browser,"$RHOME/html/index.html&"))
# }
help.start <-
function(gui = "irrelevant", browser="netscape")
{
    file <- "$RHOME/html/index.html"
    system(paste(browser, " -remote \"openURL(", file,
        ")\" 2>/dev/null || ", browser, " ", file, " &", sep = ""))
}
"hist" <-
function (x, breaks, freq = TRUE, col = NULL, border = par("fg"),
    main = paste("Histogram of" , deparse(substitute(x))),
    xlim = range(breaks), ylim = range(counts, 0),
    xlab = deparse(substitute(x)), ylab, ...)
{
    if (!is.numeric(x))
        stop("hist: x must be numeric")
    if (missing(breaks)) {
        breaks <- range(x)
        breaks <- pretty(breaks + c(0, diff(breaks)/1000),
                    1 + log2(length(x)))
    }
    else if (length(breaks) == 1) {
        if (is.na(breaks) | breaks < 2)
            stop("invalid number of breaks")
        nbreaks <- breaks
        breaks <- range(x)
        breaks <- pretty(breaks + c(0, diff(breaks)/1000), nbreaks)
    }
    breaks <- sort(breaks)
    counts <- .C("bincount",
        as.double(x),
        length(x),
        as.double(breaks),
        length(breaks),
        counts = integer(length(breaks) - 1),
        TRUE
    )[[5]]
    if (any(counts < 0))
        browser()
    if (!freq) {
        counts <- counts/(sum(!is.na(x)) * diff(breaks))
        if (missing(ylab))
            ylab <- "Relative Frequency"
    }
    else if (missing(ylab))
        ylab <- "Frequency"
    plot.new()
    plot.window(xlim, ylim, "")
    title(main = main, xlab = xlab, ylab = ylab, ...)
    axis(1, ...)
    axis(2, ...)
    rect(breaks[-length(breaks)], 0, breaks[-1], counts,
        col = col, border = border)
    invisible(NULL)
}
print.htest<-function(x, digits = 4, quote = T, prefix = "")
{
        cat("\n\t", x$method, "\n\n")
        cat("data: ", x$data.name, "\n")
        if(!is.null(x$statistic))
                cat(names(x$statistic), " = ", format(round(x$statistic, 4)),
                        ", ", sep = "")
    if(!is.null (x$parameter))
        cat(paste(names(x$parameter), " = ", format(round(x$parameter,
            3)), ",", sep = ""), "")
        cat("p-value =", format(round(x$p.value, 4)), "\n")
        if(!is.null(x$alternative)) {
                if(!is.null(x$null.value)) {
                        if(length(x$null.value) == 1) {
                                if (x$alternative == "two.sided" )
                    alt.char <- "not equal to"
                else if( x$alternative == "less" )
                    alt.char <- "less than"
                else if( x$alternative == "greater" )
                                      alt.char <- "greater than"
                                cat("alternative hypothesis:", "true", names(x$
                                  null.value), "is", alt.char, x$null.value,
                                  "\n")
                        }
                        else {
                                cat("alternative hypothesis:", x$alternative,
                                  "\n")
                                cat("null values:\n")
                                print(x$null.value)
                        }
                }
                else cat("alternative hypothesis:", x$alternative, "\n")
        }
        if(!is.null(x$conf.int)) {
                cat(format(100 * attr(x$conf.int, "conf.level")),
                        "percent confidence interval:\n", format(c(x$conf.int[1
                        ], x$conf.int[2])), "\n")
        }
        if(!is.null(x$estimate)) {
                cat("sample estimates:\n")
                print(x$estimate)
        }
        cat("\n")
        invisible(x)
}
identify <- function(x, y=NULL, text=as.character(seq(x)), pos=FALSE, ...) {
    opar <- par(list(...))
    on.exit(par(opar))
    xy <- xy.coords(x, y)
    z <- .Internal(identify(xy$x,xy$y,as.character(text)))
    i <- seq(z[[1]])[z[[1]]]
    p <- z[[2]][z[[1]]]
    if(pos) list(ind=i,pos=p) else i
}
ifelse <-
function (test, yes, no)
{
        ans <- test
        test <- as.logical(test)
        nas <- is.na(test)
        ans[test] <- rep(yes, length = length(ans))[test]
        ans[!test] <- rep(no, length = length(ans))[!test]
        ans[nas] <- NA
        ans
}
"image" <-
function (x=seq(0,1,len=nrow(z)), y=seq(0,1,len=ncol(z)), z,
    zlim=range(z, na.rm=TRUE), col=heat.colors(12), ...)
{
    plot(0, 0, xlim=range(x,na.rm=TRUE), ylim=range(y,na.rm=TRUE),
        type="n", xaxs = "i", yaxs = "i", xlab="", ylab="", ...)
    .Internal(image(
            as.double(x),
            as.double(y),
            as.double(z),
            as.double(zlim),
            col))
}
"IQR" <-
function (x)
as.vector(diff(quantile(as.numeric(x), c(0.25, 0.75))))
is.vector <- function(x, mode="any") .Internal(is.vector(x,mode))
is.finite <- function(x) !is.na(x)
is.symbol <- function(x) typeof(x)=="symbol"
lapply <- function(x, FUN, ...)
{
    if(is.character(FUN))
        FUN <- get(FUN,mode="function")
    if(mode(FUN) != "function")
        stop(paste("\"",FUN,"\" is not a function",sep=" "))
    if(!is.list(x))
        stop("lapply can only be used for lists")
    rval <- vector("list",length(x))
    for(i in seq(along=x))
        rval[i] <- list(FUN(x[[i]],...))
    names(rval) <- names(x) # keep 'names' !
    return(rval)
}
legend <-
function(x, y, legend, fill, col="black", lty, pch, bty="o", bg=par("bg"),
    xjust=0, yjust=1, ...)
{
    xchar <- xinch(par("cin")[1])
    ychar <- yinch(par("cin")[2]) * 1.2
    xbox <- xinch(par("cin")[2] * 0.8)
    ybox <- yinch(par("cin")[2] * 0.8)
    yline <- 2*xchar
    w <- 2 * xchar + max(strwidth(legend))
    h <- (length(legend)+1)*ychar
    if(missing(y)) {
        if(is.list(x)) {
            y <- x$y
            x <- x$x
        }
    }
    if(!is.numeric(x) || !is.numeric(y))
        stop("non-numeric coordinates")
    if(length(x) <= 0 || length(x) != length(y))
        stop("differing coordinate lengths")
    if(length(x) != 1) {
        x <- mean(x)
        y <- mean(y)
        xjust <- 0.5
        yjust <- 0.5
    }
    if(!missing(fill)) {
        w <- w + xchar + xbox
    }
    if(!missing(pch)) {
        if(is.character(pch) && nchar(pch) > 1) {
            np <- nchar(pch)
            pch <- substr(rep(pch[1], np), 1:np, 1:np)
        }
        w <- w + 1.5 * xchar
    }
    if(!missing(lty))
        w <- w + 3 * xchar
    x <- x - xjust * w
    y <- y + (1 - yjust) * h
    xt <- rep(x, length(legend)) + xchar
    yt <- y - (1:length(legend)) * ychar
    if(bty != "n")
        rect(x, y, x+w, y-h, col=bg)
    x <- x + xchar
    if(!missing(fill)) {
        rect(xt, yt - 0.5 * ybox,
            xt + xbox, yt + 0.5 * ybox, col=fill)
        xt <- xt + xbox + xchar
    }
    if(!missing(pch)) {
        points(xt + 0.25 * xchar, yt, pch, col=col)
        xt <- xt + 1.5 * xchar
    }
    if(!missing(lty)) {
        segments(xt, yt, xt + 2 * xchar, yt, lty=lty, col=col)
        xt <- xt + 3 * xchar
    }
    text(xt, yt, text=legend, adj=c(0, 0.35))
}
library <- function(name)
{
    if (!exists(".Libraries", inherits=TRUE))
        assign(".Libraries", character(0), NULL)
    if(missing(name)) {
        cat("NAME\tDESCRIPTION\n")
        system(paste("cat", system.file("help","LibIndex")))
    }
    else {
        name <- substitute(name)
        if (!is.character(name))
            name <- deparse(name)
        if(is.na(match(name, .Libraries))) {
            file <- system.file("library", name)
            if(file == "") stop(paste("there is no library called", name))
            sys.source(file)
            assign(".Libraries", c(name, .Libraries), NULL)
        }
        invisible(.Libraries)
    }
}
library.dynam <- function(name)
{
    if (!exists(".Dyn.libs"))
        assign(".Dyn.libs", character(0), NULL)
    if(is.na(match(name, .Dyn.libs))) {
        .Internal(dyn.load(system.file("lib", name)))
        assign(".Dyn.libs", c(.Dyn.libs, name), NULL)
    }
    invisible(.Dyn.libs)
}
license <- function() {
cat("\nThis software is distributed under the terms of the GNU GENERAL\n")
cat("PUBLIC LICENSE Version 2, June 1991.  The terms of this license\n")
cat("are in a file called COPYING which you should have received with\n")
cat("this software.\n")
cat("\n")
cat("If you have not received a copy of this file, you can obtain one\n")
cat("by writing to:\n")
cat("\n")
cat("    The Free Software Foundation, Inc.\n")
cat("    675 Mass Ave, Cambridge, MA 02139, USA\n")
cat("\n")
cat("``Share and Enjoy.''\n\n")
}
lines <- function(x, ...)
UseMethod("lines")
lines.default <- function(x, y=NULL, type="l", col=par("col"), ...) {
    plot.xy(xy.coords(x, y), type=type, col=col, ...)
}
lm <- function(formula, data=NULL, subset=NULL, weights=NULL,
    na.action=na.fail, singular.ok=TRUE)
{
    mt <- terms(formula)
    if(is.null(data)) data <- sys.frame(sys.parent())
    mf <- match.call()
    mf$singular.ok <- NULL
    mf$use.data <- TRUE
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf, sys.frame(sys.parent()))
    x <- model.matrix(mt, mf);
    y <- model.response(mf, "numeric")
    w <- model.weights(mf)
    if(is.null(w)) {
        z <- lm.fit(x,y)
    }
    else {
        z <- lm.w.fit(x,y,w)
    }
    z$call <- match.call()
    z$terms <- mt
    z$model.frame <- mf
    class(z) <- if(is.matrix(y)) c("mlm","lm") else "lm"
    z
}
lm.fit <- function(x, y)
{
    n <- nrow(x)
    p <- ncol(x)
    ny <- NCOL(y)
    if(NROW(y) != n) stop("incompatible dimensions")
    z <- .Fortran("dqrls",
        qr=x,
        n=n,
        p=p,
        y=y,
        ny=ny,
        tol=1e-7,
        coefficients=mat.or.vec(p,ny),
        residuals=y,
        effects=y,
        rank=integer(1),
        pivot=as.integer(1:p),
        qraux=double(p),
        work=double(2*p))
    coef <- z$coefficients
    pivot <- z$pivot
    r1 <- 1:z$rank
    if(ny > 1) {
        coef[-r1,] <- NA
        coef[pivot,] <- coef
        dimnames(coef) <- list(dimnames(x)[[2]],dimnames(y)[[2]])
        dimnames(z$effects)[1] <- list(NULL)
    }
    else {
        coef[-r1] <- NA
        coef[pivot] <- coef
        names(coef) <- dimnames(x)[[2]]
        names(z$effects) <- NULL
    }
    z$coefficients <- coef
    c(z[c("coefficients","residuals","effects","rank")],
        list(fitted.values=y-z$residuals,
            assign=attr(x,"assign"),
            qr=z[c("qr","qraux","pivot","tol","rank")],
            df.residual=n-z$rank))
}
lm.w.fit <- function(x, y, w)
{
    n <- nrow(x)
    p <- ncol(x)
    ny <- NCOL(y)
    if(NROW(y) != n | length(w) != n)
        stop("incompatible dimensions")
    if(any(w < 0 | is.na(w)))
        stop("missing or negative weights not allowed")
    zero.weights <- FALSE
    if(any(w == 0)) {
        zero.weights <- TRUE
        save.r <- y
        save.f <- y
        save.w <- w
        ok <- w != 0
        nok <- !ok
        w <- w[ok]
        x0 <- x[!ok,]
        x <- x[ok,]
        y0 <- if(ny>1) y[!ok,,drop=FALSE] else y[!ok]
        y <- if(ny>1) y[ok,,drop=FALSE] else y[ok]
    }
    n <- nrow(x)
    p <- ncol(x)
    wts <- w^0.5
    z <- .Fortran("dqrls",
        qr=x*wts,
        n=n,
        p=p,
        y=y*wts,
        ny=ny,
        tol=1e-7,
        coefficients=mat.or.vec(p,ny),
        residuals=y,
        effects=mat.or.vec(n,ny),
        rank=integer(1),
        pivot=as.integer(1:p),
        qraux=double(p),
        work=double(2*p))
    coef <- z$coefficients
    pivot <- z$pivot
    r1 <- 1:z$rank
    if(ny > 1) {
        coef[-r1,] <- NA
        coef[pivot,] <- coef
        dimnames(coef) <- list(dimnames(x)[[2]],dimnames(y)[[2]])
        dimnames(z$residuals) <- dimnames(y)
        dimnames(z$effects)[[2]] <- dimnames(y)[[2]]
    }
    else {
        coef[-r1] <- NA
        coef[pivot] <- coef
        names(coef) <- dimnames(x)[[2]]
        names(z$residuals) <- names(y)
    }
    z$coefficients <- coef
    z$residuals <- z$residuals/wts
    z$fitted.values <- (y - z$residuals)
    z$weights <- w
    if(zero.weights) {
        coef[is.na(coef)] <- 0
        f0 <- x0 %*% coef
        if(ny > 1) {
            save.r[ok,] <- z$residuals
            save.r[ok,] <- y0 - f0
            save.f[ok,] <- fitted.values
            save.f[nok,] <- f0
        }
        else {
            save.r[ok] <- z$residuals
            save.r[ok] <- y0 - f0
            save.f[ok] <- fitted.values
            save.f[nok] <- f0
        }
        z$residuals <- save.r
        z$fitted.values <- save.f
        z$weights <- save.w
    }
    else {
        if(ny > 1) {
            dimnames(z$residuals) <- dimnames(y)
            dimnames(z$fitted.values) <- dimnames(y)
        }
        else {
            names(z$residuals) <- names(y)
            names(z$fitted.values) <- names(y)
        }
    }
    c(z[c("coefficients","residuals","fitted.values",
        "effects","weights","rank")], list(
            assign=attr(x,"assign"),
            qr=z[c("qr","qraux","pivot","tol","rank")],
            df.residual=n-z$rank))
}
update.lm <-
function(lm.obj, formula, data, weights, subset, na.action)
{
    call <- lm.obj$call
    if(!missing(formula))
        call$formula <- update.formula(call$formula, formula)
    if(!missing(data)) call$data <- substitute(data)
    if(!missing(subset)) call$subset <- substitute(subset)
    if(!missing(na.action)) call$na.action <- substitute(na.action)
    eval(call, sys.frame(sys.parent()))
}
residuals.lm <- function(z) z$residuals
fitted.values.lm <- function(z) z$fitted.values
coefficients.lm <- function(z) z$coefficients
weights.lm <- function(z) z$weights
df.residual.lm <- function(z) z$df.residual
deviance.lm <- function(z) sum((z$residuals)^2)
summary.lm <- function(z, correlation=FALSE)
{
    n <- NROW(z$qr$qr)
    p <- z$rank
    p1 <- 1:p
    r <- resid(z)
    f <- fitted(z)
    w <- weights(z)
    if (is.null(z$terms)) {
        stop("invalid 'lm' object:  no terms component")
    } else {
        if (attr(z$terms,"intercept")) {
            if(is.null(w)) {
                rss <- sum(r^2)
                mss <- sum((f-mean(f))^2)
            } else {
                wok <- (w!=0)
                u <- (sqrt(w)/sqrt(sum(w)))[wok]
                r <- sqrt(w)*r[wok]
                f <- sqrt(w)*f[wok]
                rss <- sum(r^2)
                mss <- sum((f - sum(f*u)*u)^2)
            }
        } else { #- no intercept
            rss <- sum(r^2)
            mss <- sum(f^2)
        }
    }
    resvar <- rss/(n-p)
    R <- chol2inv(z$qr$qr[p1,p1,drop=FALSE])
    se <- sqrt(diag(R)*resvar)
    est <- z$coefficients[z$qr$pivot[p1]]
    tval <- est/se
    ans <- z[c("call","terms")]
    ans$residuals <- r
    ans$coefficients <- cbind(est, se, tval, 2*(1-pt(abs(tval),n-p)))
    dimnames(ans$coefficients) <-
        list(names(z$coefficients)[z$qr$pivot[p1]],
            c("Estimate", "Std.Error","t Value", "Pr(>|t|)"))
    ans$sigma <- sqrt(resvar)
    ans$df <- c(p, n-p, NCOL(z$qr$qr))
    if(p != attr(z$terms,"intercept")) {
          df.int <- if(attr(z$terms,"intercept")) 1 else 0
        ans$r.squared <- mss/(mss+rss)
        ans$adj.r.squared <- 1-(1-ans$r.squared)*
          ((n - df.int) / (n - p))    #0.14 :    (n/(n-p))
        ans$fstatistic <- c((mss/(p-df.int))/(rss/(n-p)),p-df.int,n-p)
         #0.14: ans$fstatistic <- c((mss/(p-1))/(rss/(n-p)),p-1,n-p)
        names(ans$fstatistic) <- c("value","numdf","dendf")
    }
    ans$cov.unscaled <- R
    dimnames(ans$cov.unscaled) <- dimnames(ans$coefficients)[c(1,1)]
    if(correlation) {
        ans$correlation <- (R*resvar)/outer(se,se)
        dimnames(ans$correlation) <- dimnames(ans$cov.unscaled)
    }
    class(ans) <- "summary.lm"
    ans
}
print.lm <- function(x, digits = max(3, .Options$digits - 3), ...)
{
    cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
    cat("Coefficients:\n")
    print(coef(x))
    cat("\n")
}
print.summary.lm <- function(x, digits = max(3, .Options$digits - 3), ...)
{
    cat("\nCall:\n")
        cat(paste(deparse(x$call), sep="\n", collapse="\n"), "\n\n", sep="")
    resid <- x$residuals
    df <- x$df
    rdf <- df[2]
    if(rdf > 5) {
        cat("Residuals:\n")
        if(length(dim(resid)) == 2) {
            rq <- apply(t(resid), 1, quantile)
            dimnames(rq) <- list(c("Min", "1Q", "Median", "3Q", "Max"),
                dimnames(resid)[[2]])
        }
        else {
            rq <- quantile(resid)
            names(rq) <- c("Min", "1Q", "Median", "3Q", "Max")
        }
        print(rq, digits = digits, ...)
    }
    else if(rdf > 0) {
        cat("Residuals:\n")
        print(resid, digits = digits, ...)
    }
    if(nsingular <- df[3] - df[1])
        cat("\nCoefficients: (", nsingular,
            " not defined because of singularities)\n", sep = "")
    else cat("\nCoefficients:\n")
    print(x$coefficients, digits=digits, quote = FALSE, ...)
    cat("\nResidual standard error:", format(signif(x$sigma, digits)), "on",
        rdf, "degrees of freedom\n")
    if(!is.null(x$fstatistic)) {
        cat("Multiple R-Squared:",
            format(signif(x$r.squared, digits)))
        cat(",  Adjusted R-squared:",
            format(signif(x$adj.r.squared, digits)),"\n")
        cat("F-statistic:", format(signif(x$fstatistic[1], digits)),
            "on", x$fstatistic[2], "and", x$fstatistic[3],
            "degrees of freedom")
        cat(",  p-value:", format(signif(1 - pf(x$fstatistic[1],
            x$fstatistic[2], x$fstatistic[3]), digits)), "\n")
    }
    correl <- x$correlation
    if(!is.null(correl)) {
        p <- dim(correl)[2]
        if(p > 1) {
            cat("\nCorrelation of Coefficients:\n")
            correl[!lower.tri(correl)] <- NA
            print(correl[-1,-NCOL(correl)], digits=digits, na="")
        }
    }
    cat("\n")
    invisible(x)
}
anova.lm <- function(object, ...)
{
    if(length(list(object, ...)) > 1)
        return(anovalist.lm(object, ...))
    w <- weights(object)
    if(is.null(w)) ssr <- sum(resid(object)^2)
    else ssr <- sum(w*resid(object)^2)
    comp <- object$effects[1:object$rank]
    asgn <- object$assign[object$qr$pivot][1:object$rank]
    dfr <- df.residual(object)
    ss <- c(as.numeric(lapply(split(comp^2,asgn),sum)),ssr)
    df <- c(as.numeric(lapply(split(asgn,asgn),length)), dfr)
    if(attr(object$terms,"intercept")) {
        ss <- ss[-1]
        df <- df[-1]
    }
    ms <- ss/df
    f <- ms/(ssr/dfr)
    p <- 1-pf(f,df,dfr)
    table <- cbind(df,ss,ms,f,p)
    table[length(p),4:5] <- NA
    dimnames(table) <- list(c(attr(object$terms,"term.labels"),
        "Residual"), c("Df","Sum Sq", "Mean Sq", "F", "Pr(>F)"))
    result <- list(table=table, title="Analysis of Variance Table")
    class(result) <- "tabular"
    result
}
"anovalist.lm" <-
function (object, ..., test = NULL)
{
    objects <- list(object, ...)
    responses <- as.character(lapply(objects, function(x) {
        as.character(x$terms[[2]])
    }))
    sameresp <- responses == responses[1]
    if (!all(sameresp)) {
        objects <- objects[sameresp]
        warning(paste("Models with response",
            deparse(responses[!sameresp]),
            "removed because response differs from", "model 1"))
    }
    # calculate the number of models
    nmodels <- length(objects)
    if (nmodels == 1)
        return(anova.lm(object))
    models <- as.character(lapply(objects, function(x) x$terms))
    # extract statistics
    df.r <- unlist(lapply(objects, df.residual))
    ss.r <- unlist(lapply(objects, deviance))
    df <- c(NA, -diff(df.r))
    ss <- c(NA, -diff(ss.r))
    ms <- ss/df
    f <- p <- rep(NA,nmodels)
    for(i in 2:nmodels) {
        if(df[i] > 0) {
            f[i] <- ms[i]/(ss.r[i]/df.r[i])
            p[i] <- 1 - pf(f[i], df[i], df.r[i])
        }
        else {
            f[i] <- ms[i]/(ss.r[i-1]/df.r[i-1])
            p[i] <- 1 - pf(f[i], -df[i], df.r[i-1])
        }
    }
    table <- cbind(df.r,ss.r,df,ss,f,p)
    dimnames(table) <- list(1:nmodels, c("Res.Df", "Res.Sum-Sq", "Df",
        "Sum-Sq", "F", "Pr(>F)"))
    # construct table and title
    title <- "Analysis of Variance Table"
    topnote <- paste("Model ", format(1:nmodels),": ",
                models, sep="", collapse="\n")
    # calculate test statistic if needed
    output <- list(table = table, title = title, topnote=topnote)
    class(output) <- "tabular"
    return(output)
}
print.anova.lm <- function(x, digits = max(3, .Options$digits - 3), ...)
{
    class(x) <- NULL
    cat("\nAnalysis of Variance:\n\n")
    print.default(round(x, digits), na="", print.gap=2)
    cat("\n")
}
effects.lm <- function(z, term) {
    term <- deparse(substitute(term))
    k <- match(term,attr(z$terms,"term.labels"))
    if(is.na(k)) stop("effect not found")
    pattern <- attr(z$terms,"factors")[,k]
    factors <- as.logical(lapply(z$model.frame,is.factor))
    y <- model.response(z$model.frame,"numeric")
    k <- range(seq(length(z$assign))[z$assign==k])
    yhat0 <- if(k[1] > 1) qr.fitted(z$qr,y,k[1]-1) else 0
    yhat1 <- qr.fitted(z$qr,y,k[2])
    effects <- yhat1-yhat0
    tapply(effects,z$model.frame[factors & pattern!=0],mean,na.rm=TRUE)
}
lm.influence <-
function(z)
{
    n <- as.integer(nrow(z$qr$qr))
    k <- as.integer(z$qr$rank)
    .Fortran("lminfl",
        z$qr$qr, n, n, k,
        z$qr$qraux,
        z$coefficients,
        z$residuals,
        hat=double(n),
        coef=matrix(0,nr=n,nc=k),
        sigma=double(n),
        DUP=FALSE)[c("hat", "coef", "sigma")]
}
rstudent <-
function(z)
{
    infl <- lm.influence(z)
    residuals(z)/(infl$sigma*sqrt(1-infl$hat))
}
dfbetas <-
function(z)
{
    infl <- lm.influence(z)
    xxi <- chol2inv(z$qr$qr,z$qr$rank)
    d <- infl$coef/(outer(infl$sigma, sqrt(diag(xxi))))
    dn <- dimnames(z$qr$qr)
    dn[[2]] <- dn[[2]][1:z$qr$rank]
    dimnames(d) <- dn
    d
}
dffits <-
function(z)
{
    infl <- lm.influence(z)
    sqrt(infl$hat)*residuals(z)/(infl$sigma*(1-infl$hat))
}
covratio <-
function(z)
{
    infl <- lm.influence(z)
    n <- nrow(z$qr$qr)
    p <- z$rank
    e.star <- residuals(z)/(infl$sigma*sqrt(1-infl$hat))
    1/((((n - p - 1)+e.star^2)/(n - p))^p*(1-infl$hat))
}
save <- function(..., list = character(0), file = "", ascii = FALSE) {
    names <- as.character( substitute( list(...)))[-1]
    list<- c(list, names)
    invisible(.Internal(save( list, file, ascii)))
}
load <- function(file)
    .Internal(load(file))
locator <- function(n=1) {
    z <- .Internal(locator(n))
    x <- z[[1]]
    y <- z[[2]]
    n <- z[[3]]
    if(n==0) NULL else list(x=x[1:n],y=y[1:n])
}
log10 <- function(x) log(x,10)
log2 <- function(x) log(x,2)
lower.tri <- function(x, diag = FALSE)
{
    x <- as.matrix(x)
        if(diag) row(x) >= col(x)
        else row(x) > col(x)
}
lowess <- function(x, y=NULL, f=2/3, iter=3, delta=.01*diff(range(xy$x[o]))) {
    xy <- xy.coords(x,y)
    if(length(xy$x) != length(xy$y)) stop("x and y lengths differ")
    n <- length(xy$x)
    o <- order(xy$x)
    .C("lowess",
        x=as.double(xy$x[o]),
        as.double(xy$y[o]),
        n,
        as.double(f),
        as.integer(iter),
        as.double(delta),
        y=double(n),
        double(n),
        double(n))[c("x","y")]
}
lsfit <- function(x, y, wt=NULL, intercept=TRUE, tolerance=1e-07, yname=NULL)
{
    # find names of x variables (design matrix)
    x <- as.matrix(x)
    y <- as.matrix(y)
    xnames <- colnames(x)
    if( is.null(xnames) ) {
        if(ncol(x)==1) xnames <- "X"
        else xnames <- paste("X", 1:ncol(x), sep="")
    }
    if( intercept ) {
        x <- cbind(1, x)
        xnames <- c("Intercept", xnames)
    }
    # find names of y variables (responses)
    if(is.null(yname) && ncol(y) > 1) yname <- paste("Y", 1:ncol(y), sep="")
    # remove missing values
    good <- complete.cases(x, y, wt)
    dimy <- dim(as.matrix(y))
    if( any(!good) ) {
        warning(paste(sum(!good), "missing values deleted"))
        x <- as.matrix(x)[good, ]
        y <- as.matrix(y)[good, ]
        wt <- wt[good]
    }
    # check for compatible lengths
    nrx <- NROW(x)
    ncx <- NCOL(x)
    nry <- NROW(y)
    ncy <- NCOL(y)
    nwts <- length(wt)
    if(nry != nrx) stop(paste("X matrix has", nrx, "responses, Y",
            "has", nry, "responses."))
    if(nry < ncx) stop(paste(nry, "responses, but only", ncx, "variables"))
    # check weights if necessary
    if( !is.null(wt) ) {
        if(any(wt < 0)) stop("negative weights not allowed")
        if(nwts != nry) stop(paste("Number of weights =", nwts,
                ", should equal", nry, "(number of responses)"))
        wtmult <- wt^0.5
        if( any(wt==0) ) {
            xzero <- as.matrix(x)[wt==0, ]
            yzero <- as.matrix(y)[wt==0, ]
        }
        x <- x*wtmult
        y <- y*wtmult
        invmult <- 1/ifelse(wt==0, 1, wtmult)
    }
    # call linpack
    storage.mode(x) <- "double"
    storage.mode(y) <- "double"
    z <- .Fortran("dqrls",
        qr=x,
        n=nrx,
        p=ncx,
        y=y,
        ny=ncy,
        tol=tolerance,
        coefficients=mat.or.vec(ncx, ncy),
        residuals=mat.or.vec(nrx, ncy),
        effects=mat.or.vec(nrx, ncy),
        rank=integer(1),
        pivot=as.integer(1:ncx),
        qraux=double(ncx),
        work=double(2*ncx))
    # dimension and name output from linpack
    resids <- array(NA, dim=dimy)
    dim(z$residuals) <- c(nry, ncy)
    if(!is.null(wt)) {
        if(any(wt==0)) {
            if(ncx==1) fitted.zeros <- xzero * z$coefficients
            else fitted.zeros <- xzero %*% z$coefficients
            z$residuals[wt==0, ] <- yzero - fitted.zeros
        }
        z$residuals <- z$residuals*invmult
    }
    resids[good, ] <- z$residuals
    if(dimy[2] == 1 && is.null(yname)) {
        resids <- as.vector(resids)
        names(z$coefficients) <- xnames
    }
    else {
        colnames(resids) <- yname
        colnames(z$effects) <- yname
        dim(z$coefficients) <- c(ncx, ncy)
        dimnames(z$coefficients) <- list(xnames, yname)
    }
    z$qr <- as.matrix(z$qr)
    colnames(z$qr) <- xnames
    output <- list(coef=z$coefficients, residuals=resids)
    # if X matrix was collinear, then the columns would have been
    # pivoted hence xnames need to be corrected
    if( z$rank != ncx ) {
        xnames <- xnames[z$pivot]
        dimnames(z$qr) <- list(NULL, xnames)
        warning("X matrix was collinear")
    }
    # return weights if necessary
    if (!is.null(wt) ) {
        weights <- rep(NA, dimy[1])
        weights[good] <- wt
        output <- c(output, list(wt=weights))
    }
    # return rest of output
    rqr <- list(qt=z$effects, qr=z$qr, qraux=z$qraux, rank=z$rank,
        pivot=z$pivot, tol=z$tol)
    output <- c(output, list(intercept=intercept, qr=rqr))
    return(output)
}
ls.diag <- function(ls.out)
{
    resids <- as.matrix(ls.out$residuals)
    xnames <- colnames(ls.out$qr$qr)
    yname <- colnames(resids)
    # remove any missing values
    good <- complete.cases(resids, ls.out$wt)
    if( any(!good) ) {
        warning("missing observations deleted")
        resids <- as.matrix(resids)[good, ]
    }
    # adjust residuals if needed
    if( !is.null(ls.out$wt) ) {
        if( any(ls.out$wt[good] == 0) )
            warning(paste("Observations with 0 weight not used in",
                "calculating standard deviation"))
        resids <- resids * ls.out$wt[good]^0.5
    }
    # initialize
    p <- ls.out$qr$rank
    n <- nrow(resids)
    hatdiag <- rep(NA, NROW(ls.out$residuals))
    stats <- matrix(NA, nrow=NROW(ls.out$residuals), ncol=NCOL(resids))
    colnames(stats) <- yname
    stdres <- studres <- dfits <- Cooks <- stats
    # calculate hat matrix diagonals
    q <- qr.qy(ls.out$qr, rbind(diag(p), matrix(0, nrow=n-p, ncol=p)))
    hatdiag[good] <- apply(as.matrix(q^2), 1, sum)
    # calculate diagnostics
    stddev <- (apply(as.matrix(resids^2), 2, sum)/(n - p))^0.5
    stddevmat <- matrix(stddev, nrow=sum(good), ncol=ncol(resids), byrow=TRUE)
    stdres[good, ] <- resids/((1-hatdiag[good])^0.5 * stddevmat)
    studres[good, ] <- (stdres[good, ]*stddevmat)/(((n-p)*stddevmat^2 -
        resids^2/(1-hatdiag[good]))/(n-p-1))^0.5
    dfits[good, ] <- (hatdiag[good]/(1-hatdiag[good]))^0.5 * studres[good, ]
    Cooks[good, ] <- ((stdres[good, ]^2 * hatdiag[good])/p)/(1-hatdiag[good])
    if(ncol(resids)==1 && is.null(yname)) {
        stdres <- as.vector(stdres)
        Cooks <- as.vector(Cooks)
        studres <- as.vector(studres)
        dfits <- as.vector(dfits)
    }
    # calculate unscaled covariance matrix
    qr <- as.matrix(ls.out$qr$qr[1:p, 1:p])
    qr[row(qr)>col(qr)] <- 0
    qrinv <- solve(qr)
    covmat.unscaled <- qrinv%*%t(qrinv)
    dimnames(covmat.unscaled) <- list(xnames, xnames)
    # calculate scaled covariance matrix
    covmat.scaled <- sum(stddev^2) * covmat.unscaled
    # calculate correlation matrix
    cormat <- covmat.scaled/
        (outer(diag(covmat.scaled), diag(covmat.scaled))^0.5)
    # calculate standard error
    stderr <- outer(diag(covmat.unscaled)^0.5, stddev)
    dimnames(stderr) <- list(xnames, yname)
    return(list(std.dev=stddev, hat=hatdiag, std.res=stdres,
        stud.res=studres, cooks=Cooks, dfits=dfits,
        correlation=cormat, std.err=stderr,
        cov.scaled=covmat.scaled, cov.unscaled=covmat.unscaled))
}
ls.print <- function(ls.out, digits=4, print.it=TRUE)
{
    # calculate residuals to be used
    resids <- as.matrix(ls.out$residuals)
    if( !is.null(ls.out$wt) ) {
        if(any(ls.out$wt == 0))
            warning("Observations with 0 weights not used")
        resids <- resids * ls.out$wt^0.5
    }
    n <- apply(resids, 2, length)-apply(is.na(resids), 2, sum)
    p <- ls.out$qr$rank
    # calculate total sum sq and df
    if(ls.out$intercept) {
        if(is.matrix(ls.out$qr$qt))
            totss <- apply(ls.out$qr$qt[-1, ]^2, 2, sum)
        else totss <- sum(ls.out$qr$qt[-1]^2)
        degfree <- p - 1
    }
    else {
        totss <- apply(as.matrix(ls.out$qr$qt^2), 2, sum)
        degfree <- p
    }
    # calculate residual sum sq and regression sum sq
    resss <- apply(resids^2, 2, sum, na.rm=TRUE)
    resse <- (resss/(n-p))^.5
    regss <- totss - resss
    rsquared <- regss/totss
    fstat <- (regss/degfree)/(resss/(n-p))
    pvalue <- 1 - pf(fstat, degfree, (n-p))
    # construct summary
    summary <- cbind(format(round(resse, digits)), format(round(rsquared,
        digits)), format(round(fstat, digits)), format(degfree), format(
        n-p), format(round(pvalue, digits)))
    dimnames(summary) <- list(colnames(ls.out$residuals), c("Mean Sum Sq",
        "R Squared", "F-value", "Df 1", "Df 2", "Pr(>F)"))
    mat <- as.matrix(ls.out$qr$qr[1:p, 1:p])
    mat[row(mat)>col(mat)] <- 0
    qrinv <- solve(mat)
    # construct coef table
    coef.table <- as.list(1:ncol(ls.out$residuals))
    if(ncol(ls.out$residuals)==1) coef <- matrix(ls.out$coef, nc=1)
    else coef <- ls.out$coef
    for(i in 1:ncol(ls.out$residuals)) {
        covmat <- (resss[i]/(n[i]-p)) * (qrinv%*%t(qrinv))
        coef.table[[i]] <- cbind(coef[, i], diag(covmat)^.5,
            coef[, i]/diag(covmat)^.5,
            2*(1 - pt(abs(coef[, i]/diag(covmat)^.5), n[i]-p)))
        dimnames(coef.table[[i]]) <- list(colnames(ls.out$qr$qr),
            c("Estimate", "Std.Err", "t-value", "Pr(>|t|)"))
        #print results
        if(print.it) {
            if(ncol(ls.out$residuals)>1)
                cat("Response:", colnames(ls.out$residuals)[i],
                "\n\n")
            cat(paste("Residual Standard Error=", format(round(
                resse[i], digits)), "\nR-Square=", format(round(
                rsquared[i], digits)), "\nF-statistic (df=",
                format(degfree), ", ", format(n[i]-p), ")=",
                format(round(fstat[i], digits)), "\np-value=",
                format(round(pvalue[i], digits)), "\n\n", sep=""))
            print(round(coef.table[[i]], digits))
            cat("\n\n")
        }
    }
    names(coef.table) <- colnames(ls.out$residuals)
    invisible(list(summary=summary, coef.table=coef.table))
}
macintosh <- function() .Internal(device("Macintosh","",c(0,0,0)))
mad <- function(y, center, constant = 1.4826, na.rm = FALSE) {
    if(na.rm)
        y <- y[!is.na(y)]
    if(missing(center))
        constant * (median(abs(y - median(y))))
    else constant * (median(abs(y - center)))
}
match <- function(x, table, nomatch=NA)
        .Internal(match(x, table, nomatch))
match.call <- function(definition=NULL,call=sys.call(sys.parent()),
expand.dots=T)
        .Internal(match.call(definition,call,expand.dots))
pmatch <- function(x, table, nomatch=NA) {
    y<-.Internal(pmatch(x,table))
    y[is.na(y)]<-nomatch
    y
}
match.arg <- function(arg, choices) {
        if (missing(choices)) {
                formal.args <- formals(sys.function(sys.parent()))
                choices <- eval(formal.args[[deparse(substitute(arg))]])
        }
        rval <- choices[pmatch(arg, choices)]
        if (is.na(rval)) {
                stop(paste("argument should be one of",
            paste(choices,collapse=", "), sep = " "))
        }
    if( length(rval) > 1 ) {
        if(arg==choices)
            rval<-choices[1]
        else
            stop("there is more than one match in match.arg")
    }
        return(rval)
}
#just for compatiblity we have charmatch and char.expand
charmatch <- pmatch
char.expand <- function(input, target, nomatch = stop("no match"))
{
    if(length(input) != 1)
        stop("char.expand: input must have length 1")
    if(!(is.character(input) && is.character(target)))
        stop("char.expand: input must be character")
    y<-.Internal(pmatch(input,target))
    if(any(is.na(y))) eval(nomatch)
    target[y]
}
matrix <- function(data=NA, nrow=1, ncol=1, byrow=FALSE, dimnames=NULL) {
                if(missing(nrow)) nrow <- ceiling(length(data)/ncol)
                else if(missing(ncol)) ncol <- ceiling(length(data)/nrow)
                x <- .Internal(matrix(data, nrow, ncol, byrow))
                levels(x) <- levels(data)
                dimnames(x)<-dimnames
                x
}
nrow <- function(x) dim(x)[1]
ncol <- function(x) dim(x)[2]
NROW <- function(x) if(is.matrix(x)) nrow(x) else length(x)
NCOL <- function(x) if(is.matrix(x)) ncol(x) else as.integer(1)
rownames <- function(x) {
    dn <- dimnames(x)
    if(is.null(dn)) dn else dn[[1]]
}
"rownames<-" <- function(x, value) {
    dn <- dimnames(x)
    if(is.null(dn)) dimnames(x) <- list(value, dn)
    else dimnames(x) <- list(value, dn[[2]])
    x
}
colnames <- function(x) {
    dn <- dimnames(x)
    if(is.null(dn)) dn else dn[[2]]
}
"colnames<-" <- function(x, value) {
    dn <- dimnames(x)
    if(is.null(dn)) dimnames(x) <- list(dn, value)
    else dimnames(x) <- list(dn[[1]], value)
    x
}
row <- function(x, as.factor=FALSE) {
    if(as.factor) factor(.Internal(row(x)), labels=rownames(x))
    else .Internal(row(x))
}
col <- function(x, as.factor=FALSE) {
    if(as.factor) factor(.Internal(col(x)), labels=colnames(x))
    else .Internal(col(x))
}
crossprod <-
function(x, y=x)
.Internal(crossprod(x,y))
t <- function(x)
 UseMethod("t")
t.data.frame<- function(x)
{
        x <- as.matrix(x)
        NextMethod("t")
}
mean <- function(x, trim = 0, na.rm = FALSE) {
        if (na.rm)
                x<-x[!is.na(x)]
        trim <- trim[1]
    if(trim > 0) {
        if(trim >= 0.5) return(median(x, na.rm=FALSE))
        lo <- floor(length(x)*trim)+1
        hi <- length(x)+1-lo
        x <- sort(x, partial=unique(c(lo, hi)))[lo:hi]
    }
        sum(x)/length(x)
}
weighted.mean <- function(x, w, na.rm = FALSE ){
    if(missing(w)) w <- rep(1,length(x))
    if (na.rm) {
        w<-w[!is.na(x)]
        x<-x[!is.na(x)]
    }
    sum(x*w)/sum(w)
}
median <- function(x, na.rm = FALSE) {
    if(na.rm)
        x <- x[!is.na(x)]
    else if(any(is.na(x)))
        return(NA)
    n <- length(x)
    half <- (n + 1)/2
    if(n %% 2 == 1) {
        sort(x, partial = half)[half]
    }
    else {
        sum(sort(x, partial = c(half, half + 1))[c(half, half + 1)])/2
    }
}
menu<-function(x)
{
    xlen<-length(x)
    cat("\n")
    for(i in 1:xlen)
        cat(i,":",x[i],"\n",sep="")
    done<-0
    repeat {
        cat("Selection: ")
        ind<-.Internal(menu(as.character(x)))
        if(ind<=xlen)
            return(ind)
        cat("Enter an item from the menu, or 0 to exit\n")
    }
}
mode <- function(x) {
    if (is.expression(x)) return("expression")
    tx <- typeof(x)
    if (tx == "real" ) return("numeric")
    if (tx == "integer") return("numeric")
    if (tx == "closure" || tx == "builtin" ) return("function")
    if (tx == "language") {
        if(is.call(x)) return("call")
        if(is.name(x)) return("name")
    }
    tx
}
"mode<-" <- function(x, value)
{
    mde <- paste("as.",value,sep="")
    atr <- attributes(x)
    x <- eval(call(mde,x), sys.frame(sys.parent()))
    attributes(x) <- atr
    x
}
storage.mode <- function(x) {
        x <- typeof(x)
    if (x == "closure" || x == "builtin" ) return("function")
    x
}
"storage.mode<-" <- get("mode<-")
formula <- function(x, ...) UseMethod("formula")
formula.default <- function(x) {
    if (!is.null(x$formula))
        eval(x$formula)
    switch(typeof(x),
        NULL = structure(NULL, class="formula"),
        character = formula(eval(parse(text=x)[[1]])),
        call = eval(x),
        stop("invalid formula"))
}
formula.formula <- function(x) x
formula.terms <- function(x) {
    attributes(x) <- list(class="formula")
    x
}
print.formula <- function(x) print.default(unclass(x))
terms <- function(x, ...) UseMethod("terms")
print.terms <- function(x) print.default(unclass(x))
terms.default <- function(x) x$terms
terms.terms <- function(x) x
delete.response <-
function (termobj)
    terms(reformulate(attr(termobj, "term.labels"), NULL),
        specials=names(attr(termobj, "specials")))
reformulate <-
function (termlabels, response=NULL)
{
    if (is.null(response)){
        termtext <- paste("~", paste(termlabels, collapse="+"),collapse="")
        termobj <- eval(parse(text=termtext)[[1]])
    }
    else {
        termtext <- paste("response", "~", paste(termlabels, collapse="+"),
            collapse="")
        termobj <- eval(parse(text=termtext)[[1]])
        termobj[[2]] <- response
    }
    termobj
}
drop.terms <-
function(termobj, dropx=NULL, keep.response=FALSE)
{
    if (is.null(dropx))
        termobj
    else {
        newformula <- reformulate(attr(termobj, "term.labels")[-dropx], if (keep.response) termobj[[2]] else NULL)
        terms(newformula, specials=names(attr(termobj, "specials")))
    }
}
terms.formula <-
function (x, specials = NULL, abb = NULL, data = NULL, keep.order = FALSE)
{
    new.specials <- unique(c(specials, "offset"))
    terms <- .Internal(terms.formula(x, new.specials, abb, data, keep.order))
    offsets <- attr(terms,"specials")$offset
    if(!is.null(offsets)) {
        names <- dimnames(attr(terms,"factors"))[[1]][offsets]
        offsets <- match(names, dimnames(attr(terms,"factors"))[[2]])
        offsets <- offsets[!is.na(offsets)]
        if(length(offsets) > 0) {
            attr(terms, "factors") <- attr(terms,"factors")[,-offsets, drop=FALSE]
            attr(terms, "term.labels") <- attr(terms, "term.labels")[-offsets]
            attr(terms, "order") <- attr(terms, "order")[-offsets]
            attr(terms, "offset") <- attr(terms,"specials")$offset
        }
    }
    attr(terms, "specials")$offset <- NULL
    terms
}
coefficients <- function(x, ...)
UseMethod("coefficients")
coef <- coefficients
residuals <- function(x, ...)
UseMethod("residuals")
resid <- residuals
deviance <- function(x, ...)
UseMethod("deviance")
fitted.values <- function(x, ...)
UseMethod("fitted.values")
fitted <- fitted.values
anova <- function(x, ...)
UseMethod("anova")
effects <- function(x, ...)
UseMethod("effects")
weights <- function(x, ...)
UseMethod("weights")
df.residual <- function(x, ...)
UseMethod("df.residual")
offset <- function(x) x
na.action <- function(x, ...)
UseMethod("na.action")
na.action.default <- function(x) attr(x, "na.action")
na.fail <- function(frame)
{
    ok <- complete.cases(frame)
    if(all(ok)) frame else stop("missing values in data frame");
}
na.omit <- function(frame)
{
    ok <- complete.cases(frame)
    if (all(ok))
        return(frame)
    else return(frame[ok, ])
}
model.data.frame <- function(...) {
    cn <- as.character(substitute(list(...))[-1])
    data.frame(..., col.names=cn, as.is=TRUE)
}
"model.frame" <-
function (formula, data = sys.frame(sys.parent()), subset = NULL,
    na.action = eval(as.name(options("na.action")$na.action)),
    use.data = TRUE, process.offsets = TRUE, ...)
{
        if (!is.null(formula$model) && missing(data))
      return(formula$model)
    if (!missing(data) || is.null(formula$model.frame)) {
        dotsdata <- if (use.data)
            data
        else sys.frame(sys.parent())
        newframe <- substitute(list(...))
        dots <- eval(newframe, dotsdata)
        if (!is.null(dots)) {
            real.dots <- !unlist(lapply(dots, is.null))
            newframe <- as.call(newframe[c(T, real.dots)])
            dots <- dots[real.dots]
        }
        Terms <- terms(formula)
        frame <- attr(Terms, "variables")
        name.process <- function(x) paste("(", x, ")", sep = "")
        if (missing(data) && !is.null(formula$call)) {
            if (is.null(formula$call$data))
                data<-environment(NULL)
            else
                data <- eval(formula$call$data)
        }
        if (!(missing(subset) || exists(as.character(match.call()$subset), inherits = FALSE)))
            subset <- eval(match.call()$subset, data)
        if (is.null(dots))
            rval <- na.action(eval(frame, data)[subset, , drop = FALSE])
        else {
            dotnames <- sapply(names(eval(dots, data)), name.process)
            val <- eval(frame, data)
            newframe[[1]] <- as.name("model.data.frame")
            for (i in 1:length(dots)) newframe[[i + 1]] <- dots[[i]]
            dotsval <- eval(newframe, dotsdata)
            names(dotsval) <- dotnames
            if (dim(val)[1] == dim(dotsval)[1])
                newval <- c(val, dotsval)
            else stop("Mismatched dimensions in model.frame")
            class(newval) <- "data.frame"
            rval <- na.action(newval[subset, , drop = FALSE])
        }
        attr(rval, "terms") <- Terms
        offset.pos <- attr(Terms, "offset")
        if (process.offsets && (length(offset.pos) > 0)) {
            offset.total <- as.vector(as.matrix(rval[, offset.pos]) %*% rep(1, length(offset.pos
            )))
            rval[[offset.pos[1]]] <- offset.total
            names(rval)[offset.pos[1]] <- "(offset)"
        }
        rval
    }
    else formula$model.frame
}
model.weights <- function(x) x$"(weights)"
model.offset <- function(x) x$"(offset)"
model.matrix <- function (formula, data)
{
    t <- terms(formula)
    if (missing(data))
        data <- eval(attr(t, "variables"), sys.frame(sys.parent()))
    .Internal(model.matrix(t, data))
}
model.response <- function(data, type="numeric")
{
    if(attr(attr(data,"terms"), "response")) {
        if(is.list(data) | is.data.frame(data)) {
            v <- data[[1]]
            if(type == "numeric" | type == "double") {
                storage.mode(v) <- "double"
            }
            else stop("invalid response type")
            if(is.matrix(v) && ncol(v) == 1)
                dim(v) <- NULL
            return(v)
        }
        else stop("invalid data argument")
    }
    else
        return (NULL)
}
model.extract <- function(frame, component)
{
    component<-as.character(substitute(component))
    rval<-switch(component,
        response= model.response(frame),
        offset = model.offset(frame),
        weights = frame$"(weights)",
        start = frame$"(start)"
        )
    if(length(rval) == nrow(frame))
        names(rval)<-attr(frame, "row.names")
    else if(is.matrix(rval) && nrow(rval)==nrow(frame)) {
        t1<-dimnames(rval)
        dimnames(rval)<-list(attr(frame, "row.names"),t1[[2]])
    }
    return(rval)
}
update <-
function(x, ...)
    UseMethod("update")
mtext <- function(text, side=3, line=0, outer=FALSE, at=NULL, ...)
  .Internal(mtext(text, side, line, outer, at, ...))
named.elements <-
function(x)
{
    n <- names(x)
    if(is.null(n)) NULL
    else !is.na(n)
}
names <-
function(x, ...)
UseMethod("names")
names.default <-
function(x)
.Internal(names(x))
"names<-" <-
function(x, ...)
UseMethod("names<-")
"names<-.default" <-
function(x, n)
.Internal("names<-"(x, n))
nlm <-
function(f, p, hessian=FALSE, typsiz=rep(1,length(p)),
    fscale=1, print.level=0, ndigit=12, gradtl=1.e-6,
    stepmx=max(1000 * sqrt(sum((p/typsiz)^2)), 1000),
    steptl=1.e-6, itnlim=100)
{
    if(print.level == 0) msg <- msg <- 9
    else if(print.level == 1) msg <- 1
    else if(print.level == 2) msg <- 17
    .Internal(nlm(f, p, hessian, typsiz, fscale, msg, ndigit, gradtl,
        stepmx, steptl, itnlim))
}
optimize <-
function(f, interval, lower=min(interval), upper=max(interval),
    maximum=FALSE, tol=.Machine$double.eps^0.25, ...)
{
    if(maximum) {
        val <- .Internal(fmin(function(arg) -f(arg, ...), lower, upper, tol))
        list(maximum=val, objective=-f(val, ...))
    }
    else {
        val <- .Internal(fmin(function(arg) f(arg, ...), lower, upper, tol))
        list(minimum=val, objective=f(val, ...))
    }
}
uniroot <-
function(f, interval, lower=min(interval), upper=max(interval), tol=.Machine$double.eps^0.25, ...)
{
    if(f(interval[1], ...)*f(interval[2], ...) >= 0)
        stop("signs at end points not of opposite sign")
    val <- .Internal(zeroin(function(arg) f(arg, ...), lower, upper, tol))
    list(root=val, f.root=f(val, ...))
}
deriv <- function(x, ...)
UseMethod("deriv")
deriv.formula <- function(expr, namevec, function.arg=NULL, tag=".expr") {
if(length(expr) == 2) .Internal(deriv.default(expr[[2]], namevec, function.arg, tag))
else stop("invalid formula in deriv")
}
deriv.default <- function(expr, namevec, function.arg=NULL, tag=".expr")
.Internal(deriv.default(expr, namevec, function.arg, tag))
inherits <-
function(x, name)
{
    if(is.object(x)) any(!is.na(match(name,class(x))))
    else FALSE
}
NextMethod <-
function(generic=NULL, object=NULL, ...)
    .Internal(NextMethod(generic, object,...))
methods <-
function (generic.function, class)
{
    allnames <- unique(c(ls(.SystemEnv), ls(.GlobalEnv)))
    if (!missing(generic.function)) {
        if (!is.character(generic.function))
            generic.function <- deparse(substitute(generic.function))
        name <- paste("^", generic.function, ".", sep = "")
    }
    else if (!missing(class)) {
        if (!is.character(class))
            class <- paste(deparse(substitute(class)))
        name <- paste(".", class, "$", sep = "")
    }
    else stop("must supply generic.function or class")
    grep(gsub("\\.", "\\\\.", name), allnames, value = TRUE)
}
options <-
function(...) .Internal(options(...))
outer <- function(x, y, FUN="*", ...) {
        if(is.character(FUN))
                FUN <- get(FUN, mode="function", inherits=TRUE)
        nr <- length(x)
        nc <- length(y)
        matrix(
                FUN(matrix(x, nr, nc), matrix(y, nr, nc, byrow=TRUE), ...),
                nr, nc)
}
"%o%"<-outer
pairs <- function (x, labels, panel=points, main = NULL, font.main=par("font.main"),
    cex.main=par("cex.main"), ...)
{
    if(!is.matrix(x)) x <- data.matrix(x)
    if(!is.numeric(x)) stop("non-numeric argument to pairs")
    nc <- ncol(x)
    if(nc < 2) stop("only one column in the argument to pairs")
    if (missing(labels)) {
        labels <- dimnames(x)[[2]]
        if (is.null(labels))
            labels <- paste("var", 1:nc)
    }
    oma <- c(4, 4, 4, 4)
    if (!is.null(main))
        oma[3] <- 6
    opar <- par(mfrow = c(nc, nc), mar = rep(0.5, 4), oma = oma)
    on.exit(par(opar))
    for (i in 1:nc) for (j in 1:nc) {
        if (i == j) {
            plot(x[, j], x[, i], xlab = "", ylab = "", axes = FALSE, type = "n",
                ...)
            box()
            text(mean(par("usr")[1:2]), mean(par("usr")[3:4]), labels[i])
        }
        else {
            plot(x[, j], x[, i], type="n", xlab = "", ylab = "", axes = FALSE, ...)
            box()
            panel(x[, j], x[, i], ...)
        }
        if (j == 1 & 2 * floor(i/2) == i)
            axis(2)
        if (i == 1 & 2 * floor(j/2) == j)
            axis(3)
        if (j == nc & 2 * floor(i/2) != i)
            axis(4)
        if (i == nc & 2 * floor(j/2) != j)
            axis(1)
    }
    if (!is.null(main)) mtext(main, 3, 3, T, 0.5,
        cex=cex.main/par("cex"), font=font.main)
    invisible(NULL)
}
.Pars <- c(
"1em", "adj", "ask", "bty", "cex", "cin", "col", "cra", "crt", "csi",
"cxy", "din", "err", "exp", "fig", "fin", "font", "frm", "fty", "lab",
"las", "lty", "lwd", "mai", "mar", "mex", "mfg", "mgp", "new", "oma",
"omd", "omi", "pch", "pin", "plt", "pty", "rsz", "smo", "srt", "tck",
"uin", "usr", "xaxp", "xaxs", "xaxt", "xpd", "yaxp", "yaxs", "yaxt")
##>> R-alpha  0.16.1
##>>  pp_ par(); names(pp[sapply(pp,is.null)])
##>>  [1] "1em" "cxy" "din" "exp" "frm" "fty" "rsz" "uin"
##>>       ---   ---   ---   ---   ---   ---   ---   ---
##-- These are the ones used in 0.16.1 -- $RHOME/src/main/par.c  Query(..) :
.Pars <- c(
"adj", "ann", "ask", "bg", "bty",
"cex", "cex.axis", "cex.lab", "cex.main", "cex.sub", "cin",
"col", "col.axis", "col.lab", "col.main", "col.sub", "cra", "crt", "csi",
"err", "fg", "fig", "fin",
"font", "font.axis", "font.lab", "font.main", "font.sub", "lab", "las",
"lty", "lwd", "mai", "mar", "mex", "mfcol", "mfg", "mfrow", "mgp", "mkh",
"new", "oma", "omd", "omi", "pch", "pin", "plt", "ps", "pty",
"smo", "srt", "tck", "tmag", "type", "usr",
"xaxp", "xaxs", "xaxt", "xlog", "xpd",
"yaxp", "yaxs", "yaxt", "ylog")
par <-
function (...)
{
    single <- FALSE
    if (nargs() == 0) {
        args <- as.list(.Pars)
    }
    else {
        args <- list(...)
        if (length(args) == 1) {
            if (is.list(args[[1]]) | is.null(args[[1]]))
                args <- args[[1]]
            else
                if(is.null(names(args)))
                    single <- TRUE
        }
    }
    value <- if (single) .Internal(par(args))[[1]]
    else .Internal(par(args))
        if(!is.null(names(args))) invisible(value) else value
}
par2 <-
function (...)
{
    args <- list(...)
    if (length(args) == 1 && is.list(args[[1]]))
        args <- args[[1]]
    .Internal(par2(args))
}
#we don't use white; it's for compatibility
parse <- function(file="", n=NULL, text=NULL, prompt=NULL, white=FALSE)
as.expression(.Internal(parse(file, n, text, prompt)))
"[.expression"<-function(x,subs)
    structure(unclass(x)[subs],class="expression")
"[[.expression"<-function(x,subs)
    unclass(x)[[subs]]
is.expression<-function(x)
    inherits(x,"expression")
as.expression<-function(x)
{
    if(is.null(x) || is.expression(x))
        x
    else
        structure(as.list(x),class="expression")
}
expression <- function(...){
    y<-substitute(list(...))[-1]
    structure(y,class="expression")
}
print.expression<-function(x)
{
    print(as.call(c(as.name("expression"),x)))
}
paste <- function (..., sep = " ", collapse=NULL)
{
        args <- list(...)
    if(is.null(args)) ""
    else {
        for (i in 1:length(args)) args[[i]] <- as.character(args[[i]])
        .Internal(paste(args, sep, collapse))
    }
}
pictex <-
function(file="Rplots.tex", width=5, height=4, debug = FALSE,
        bg="white", fg="black")
{
        .Internal(device("pictex", as.character(c(file)),
                c(width, height, debug, bg, fg)))
        par(mar=c(5,4,2,4)+0.1)
}
piechart <- function (x, labels=names(x), edges=200,
    radius=0.8, col=NULL, main=NULL, ...)
{
    if (any(is.na(x) | x <= 0))
        stop("invalid values in pie")
    if (is.null(labels))
        labels <- as.character(1:length(x))
    x <- c(0, cumsum(x)/sum(x))
    dx <- diff(x)
    twopi <- 8 * atan(1)
    pin <- par("pin")
    xlim <- c(-1, 1)
    ylim <- c(-1, 1)
    if (pin[1] > pin[2]) xlim <- (pin[1]/pin[2]) * xlim
    else ylim <- (pin[2]/pin[1]) * ylim
    plot.new()
    plot.window(xlim, ylim, "")
    for (i in 1:length(dx)) {
        n <- floor(edges * dx[i])
        t <- seq(x[i], x[i + 1], length = n)
        xc <- c(cos(twopi * t), 0) * radius
        yc <- c(sin(twopi * t), 0) * radius
        polygon(xc, yc, col=col[(i-1)%%length(col)+1])
        t <- mean(x[i + 0:1])
        xc <- cos(twopi * t) * radius
        yc <- sin(twopi * t) * radius
        lines(c(1,1.05)*xc, c(1,1.05)*yc)
        text(1.1*xc, 1.1*yc, labels[i], xpd = TRUE, adj = ifelse(xc < 0, 1, 0))
    }
    title(main = main, ...)
    invisible(NULL)
}
xy.coords <- function(x, y, xlab=NULL, ylab=NULL) {
    if(is.null(y)) {
        ylab<- xlab
        if(is.language(x)) {
            if(length(x) == 3 && deparse(x[[1]]) == '~') {
                ylab <- deparse(x[[2]])
                xlab <- deparse(x[[3]])
                y <- eval(x[[2]], sys.frame(sys.parent()))
                x <- eval(x[[3]], sys.frame(sys.parent()))
            }
            else stop("invalid first argument")
        }
        else if(is.ts(x)) {
            if(is.matrix(x)) y <- x[,1]
            else y <- x
            x <- time(x)
            xlab <- "Time"
        }
        else if(is.complex(x)) {
            y <- Im(x)
            x <- Re(x)
            xlab <- paste("Re(", ylab, ")", sep="")
            ylab <- paste("Im(", ylab, ")", sep="")
        }
        else if(is.matrix(x) || is.data.frame(x)) {
            x <- data.matrix(x)
            if(ncol(x) == 1) {
                xlab <- "Index"
                y <- x[,1]
                x <- 1:length(y)
            }
            else {
                colnames <- dimnames(x)[[2]]
                if(is.null(colnames)) {
                    xlab <- paste(ylab,"[,1]",sep="")
                    ylab <- paste(ylab,"[,2]",sep="")
                }
                else {
                    xlab <- colnames[1]
                    ylab <- colnames[2]
                }
                y <- x[,2]
                x <- x[,1]
            }
        }
        else if(is.list(x)) {
            xlab <- paste(ylab,"$x",sep="")
            ylab <- paste(ylab,"$y",sep="")
            y <- x[["y"]]
            x <- x[["x"]]
        }
        else {
            if(is.factor(x)) x <- as.numeric(x)
            xlab <- "Index"
            y <- x
            x <- 1:length(x)
        }
    }
    else if(length(x) != length(y)) stop("x and y lengths differ")
    return(list(x=as.real(x), y=as.real(y), xlab=xlab, ylab=ylab))
}
plot <- function(x, ...)
UseMethod("plot")
plot.default <-
function (x, y=NULL, type="p", col=par("fg"), bg=NA, pch=par("pch"), xlim=NULL,
    ylim=NULL, log="", axes=TRUE, frame.plot=TRUE, panel.first=NULL,
    panel.last=NULL, ann=par("ann"), main=NULL, xlab=NULL, ylab=NULL,
    cex=par("cex"), lty=par("lty"), ...)
{
    xlabel <- if (!missing(x))
        deparse(substitute(x))
    else NULL
    ylabel <- if (!missing(y))
        deparse(substitute(y))
    else NULL
    xy <- xy.coords(x, y, xlabel, ylabel)
    xlab <- if (missing(xlab))
        xy$xlab
    else xlab
    ylab <- if (missing(ylab))
        xy$ylab
    else ylab
    xlim <- if (missing(xlim))
        range(xy$x, na.rm=TRUE)
    else xlim
    ylim <- if (missing(ylim))
        range(xy$y, na.rm=TRUE)
    else ylim
    plot.new()
    plot.window(xlim, ylim, log, ...)
    panel.first
    plot.xy(xy, type, col=col, pch=pch, cex=cex, bg=bg, lty=lty, ...)
    panel.last
    if (axes) {
        axis(1, ...)
        axis(2, ...)
        if (frame.plot)
            box(...)
    }
    if (ann)
        title(main=main, xlab=xlab, ylab=ylab, ...)
    invisible()
}
plot.factor <-
function(x, y, ...)
{
    if(missing(y))
        barplot(table(x), ...)
    else NextMethod("plot")
}
plot.xy <- function(xy, type, pch=1, lty="solid", col=par("fg"), bg=NA, cex=1, ...)
    .Internal(plot.xy(xy, type, pch, lty, col, bg=bg, cex, ...))
plot.new <- function(ask=NA)
    .Internal(plot.new(ask))
frame <- plot.new
pmax <-
function (..., na.rm = FALSE)
{
        elts <- list(...)
        maxmm <- as.vector(elts[[1]])
        for (each in elts[-1]) {
            work <- cbind(maxmm, as.vector(each))
            nas <- is.na(work)
            work[,1][nas[,1]] <- work[,2][nas[,1]]
            work[,2][nas[,2]] <- work[,1][nas[,2]]
            change <- work[,1] < work[,2]
            work[,1][change] <- work[,2][change]
            if (!na.rm) work[,1][nas[,1]+nas[,2] > 0] <- NA
            maxmm <- work[,1]
        }
        maxmm
}
pmin <-
function (..., na.rm = FALSE)
{
        elts <- list(...)
        minmm <- as.vector(elts[[1]])
        for (each in elts[-1]) {
            work <- cbind(minmm, as.vector(each))
            nas <- is.na(work)
            work[,1][nas[,1]] <- work[,2][nas[,1]]
            work[,2][nas[,2]] <- work[,1][nas[,2]]
            change <- work[,1] > work[,2]
            work[,1][change] <- work[,2][change]
            if (!na.rm) work[,1][nas[,1]+nas[,2] > 0] <- NA
            minmm <- work[,1]
        }
        minmm
}
points <- function(x, ...)
UseMethod("points")
points.default <- function(x, y=NULL, pch=1, col="black", ...) {
    plot.xy(xy.coords(x,y), type="p", pch=pch, col=col, ...)
}
polygon <-
function(x, y=NULL, border=par("fg"), ...)
{
    xy <- xy.coords(x, y)
    .Internal(polygon(xy$x, xy$y, border=border, ...))
}
postscript <- function(
    file="Rplots.ps",
    paper=options("papersize")$papersize,
    landscape=TRUE,
    width=0,
    height=0,
    family="Helvetica",
    pointsize=12,
    bg="white",
    fg="black",
    onefile,
    print.it,
    append)
{
    .Internal(device(
        "postscript",
        as.character(c(file, paper, family, bg, fg)),
        c(width, height, landscape, pointsize)))
}
ppoints <- function(x) {
    n <- length(x)
    if(n == 1) n <- x
    (1:n-0.5)/n
}
predict <- function(fit,...) UseMethod("predict")
predict.default <- function (object, ...) {
                 namelist <- list(...)
                 names(namelist) <- substitute(list(...))[-1]
                 m <- length(namelist)
                 X <- as.matrix(namelist[[1]])
                 if (m > 1)
                   for (i in (2:m)) X <- cbind(X, namelist[[i]])
                 if (object$intercept)
                   X <- cbind(rep(1, NROW(X)), X)
                 k <- NCOL(X)
                 if (length(object$coef) != k)
                   stop("Wrong number of predictors")
                 predictor <- X %*% object$coef
                 ip <- real(NROW(X))
                 for (i in (1:NROW(X))) ip[i] <- sum(X[i, ] *
                       (object$covmat %*% X[i, ]))
                 stderr1 <- sqrt(ip)
                 stderr2 <- sqrt(object$rms^2 + ip)
                 tt <- qt(0.975, object$df)
                 conf.l <- predictor - tt * stderr1
                 conf.u <- predictor + tt * stderr1
                 pred.l <- predictor - tt * stderr2
                 pred.u <- predictor + tt * stderr2
                 z <- cbind(predictor, conf.l, conf.u, pred.l, pred.u)
                 rownames(z) <- paste("P", 1:NROW(X), sep = "")
                 colnames(z) <- c("Predicted", "Conf lower", "Conf upper",
                                  "Pred lower", "Pred upper")
                 z
}
pretty <- function(x, n=5) {
    if(!is.numeric(x))
        stop("x must be numeric")
    if(is.na(n[1]) || n[1] < 1)
        stop("invalid n value")
    z <- .C("pretty",l=min(x),u=max(x),n=as.integer(n))
    seq(z$l,z$u,length=z$n+1)
}
print <- function(x, ...)
UseMethod("print")
print.default <- function(x,digits=NULL,quote=TRUE,na.print=NULL,print.gap=NULL)
{
    .Internal(print.default(x,digits,quote,na.print,print.gap))
}
print.atomic <- function(x,quote=TRUE,...) print.default(x,quote=quote)
prmatrix <- function(x, rowlab=character(0), collab=character(0),
    quote=TRUE, right=FALSE)
    .Internal(prmatrix(x,rowlab,collab,quote,right))
print.tabular <- function(table, digits = max(3, .Options$digits - 3),
        na.print = "", ...)
{
    if(!is.null(table$title)) cat("\n", table$title, "\n\n", sep="")
    if(!is.null(table$topnote))
        cat(paste(table$topnote, collapse="\n"), "\n\n", sep="")
        print.default(table$table, digits=digits, na = "", print.gap = 2)
    if(!is.null(table$botnote)) cat("\n",
        paste(table$botnote, collapse="\n"), sep="")
        cat("\n")
}
### -*- R -*-
prompt <- function(object, ...) UseMethod("prompt")
## Later, we may want  a data.frame method (as S).
prompt.default <- function(object, filename = paste0(name, ".man"),
               force.function = FALSE)
{
  paste0 <- function(...) paste(..., sep = "")
  name <- substitute(object)
  if(is.language(name) && !is.name(name)) name <- eval(name)
  name <- as.character(name)
  fn <- get(name)
  ##-- 'file' will contain the "lines" to be put in the  manual file
  if(is.function(fn) || force.function) {
    argls <- formals(fn)
    n <- length(argls)
    if(n > 0) { s <- 1:n; arg.names <- names(argls)
    }      else s <- integer(0)
    file <- paste0("TITLE(", name, " @@ ~~function to do ... ~~)")
    ##-- Construct the 'call':
    call <- paste0(name, "(")
    for(i in s) {
      n.i <- arg.names[i]
      call <- paste0(call, n.i, if(n.i == "...") "" else "=",
             if(mode(argls[[i]]) == "missing") ""
               else deparse(argls[[i]]))
      if(i != n) call <- paste0(call, ", ")
    }
    file <- c(file, "USAGE(", paste0(call, ")"),
          ")", paste0("ALIAS(", name, ")"))
    if(length(s))
      file <- c(file, "ARGUMENTS(",
        paste0("ARG(", arg.names, " @@",
               " ~~Describe ", arg.names, " here~~  )"),
        ")")
    fn.def <- deparse(fn)
    if(any(br <- substr(fn.def,1,1) == ")"))
      fn.def[br] <- paste(" ", fn.def[br])
    file <- c(file,
          "DESCRIPTION(",
          "~~ a precise description of what the function does. ~~",
          ")",
          "VALUE(",
          "~Describe the value returned",
          ")",
          "~~~~~~~ For  multiple values (list), use 'VALUES' INSTEAD !",
          "             ---------------              ------",
          "VALUES(",
          "A description of the LIST of values returned.  Use",
          "@@",
          "ARG(comp1 @@ Description of `comp1')",
          "ARG(comp2 @@ Description of `comp2')",
          "...",
          ")",
          "REFERENCES(",
          "~put references to the literature / web site here,...",
          ")",
          "NOTE(",
          "~further notes~",
          ")",
          "~make other sections like WARNING with SECTION(...)~",
          "SEEALSO(", "~ functions to SEE ALSO as  LANG(LINK(~~fun~~))",
          ")",
          "EXAMPLES(",
          "##---- Should be DIRECTLY executable !! ----",
          "##-- ==>  Define data, use random,",
          "##--      or do  help(data=index)  for the standard data sets.",
          "BLANK", "## The function is currently defined as",
          fn.def,
          ")"
          ##-- not yet: , "KEYWORD( ~keyword )"
          )
  } else { #-- not function --
    file <- c("NON_FUNCTION()",
          paste("TITLE(", name, "@@ ~~data-name / kind ...  )"),
          "DESCRIPTION(",
          "~~ a precise description of what the function does. ~~",
          ")")
  }
  cat(file, file = filename, sep = "\n")
  RHOME <- system("printenv RHOME", intern = TRUE)
  if(substr(RHOME,1,8) == "/tmp_mnt") RHOME <- substr(RHOME,9,1000)
  cat("created file named ", filename, " in the current directory.\n",
      " Edit the file and move it to the appropriate directory,\n",
      paste(RHOME,"src/manual/man/",sep="/") , "dropping the ending '.man'\n"
      )
  on.exit()  # DEBUG off
  invisible(file)
}
qqnorm <- function(y, ylim, main="Normal Q-Q Plot",
    xlab="Theoretical Quantiles", ylab="Sample Quantiles", ...) {
    y <- y[!is.na(y)]
    if(missing(ylim)) ylim <- c(min(y),max(y))
    x <- (1:length(y)-0.5)/length(y)
    plot(qnorm(x), sort(y), main=main ,xlab=xlab, ylab=ylab, ylim=ylim, ...)
}
qqline <-
function(y, ...)
{
    y <- quantile(y[!is.na(y)],c(0.25, 0.75))
    x <- qnorm(c(0.25, 0.75))
    slope <- diff(y)/diff(x)
    int <- y[1]-slope*x[1]
    abline(int, slope, ...)
}
qqplot <- function(x, y, plot.it = TRUE, xlab = deparse(substitute(x)),
ylab = deparse(substitute(y)), ...)
{
    sx<-sort(x)
    sy<-sort(y)
    lenx<-length(sx)
    leny<-length(sy)
    if( leny < lenx )
        sx<-approx(1:lenx, sx, n=leny)$y
    if( leny > lenx )
        sy<-approx(1:leny, sy, n=lenx)$y
    if(plot.it)
        plot(sx, sy, xlab = xlab, ylab = ylab, ...)
    invisible(list(x = sx, y = sy))
}
is.qr <- function(x) !is.null(x$qr)
qr <- function(x, tol= 1e-07)
{
    x <- as.matrix(x)
    p <- as.integer(ncol(x))
    n <- as.integer(nrow(x))
    if(!is.double(x))
        storage.mode(x) <- "double"
    .Fortran("dqrdc2",
        qr=x,
        n,
        n,
        p,
        as.double(tol),
        rank=integer(1),
        qraux = double(p),
        pivot = as.integer(1:p),
        double(2*p))[c(1,6,7,8)]
}
qr.coef <- function(qr, y)
{
    if( !is.qr(qr) )
        stop("first argument must be a QR decomposition")
    n <- nrow(qr$qr)
    p <- ncol(qr$qr)
    k <- as.integer(qr$rank)
    y <- as.matrix(y)
    ny <- as.integer(ncol(y))
    storage.mode(y) <- "double"
    if( nrow(y) != n )
        stop("qr and y must have the same number of rows")
    z <- .Fortran("dqrcf",
        as.double(qr$qr),
        n, k,
        as.double(qr$qraux),
        y,
        ny,
        coef=matrix(0,nr=k,nc=ny),
        info=integer(1))
    if(z$info != 0) stop("exact singularity in qr.coef")
    if(k < p) {
        coef <- matrix(as.double(NA),nr=p,nc=ny)
        coef[qr$pivot[1:k],] <- z$coef
    }
    else coef <- z$coef
    if(ncol(y) == 1)
        dim(coef) <- NULL
    return(coef)
}
qr.qy <- function(qr, y)
{
    if(!is.qr(qr)) stop("argument is not a QR decomposition")
    n <- as.integer(nrow(qr$qr))
    p <- as.integer(ncol(qr$qr))
    k <- as.integer(qr$rank)
    y <- as.matrix(y)
    ny <- as.integer(ncol(y))
    storage.mode(y) <- "double"
    if( nrow(y) != n )
        stop("qr and y must have the same number of rows")
    .Fortran("dqrqy",
        as.double(qr$qr),
        n, k,
        as.double(qr$qraux),
        y,
        ny,
        qy=mat.or.vec(n,ny))$qy
}
qr.qty <- function(qr, y)
{
    if(!is.qr(qr)) stop("argument is not a QR decomposition")
    n <- as.integer(nrow(qr$qr))
    p <- as.integer(ncol(qr$qr))
    k <- as.integer(qr$rank)
    y <- as.matrix(y)
    ny <- as.integer(ncol(y))
    storage.mode(y) <- "double"
    if( nrow(y) != n )
        stop("qr and y must have the same number of rows")
    .Fortran("dqrqty",
        as.double(qr$qr),
        n, k,
        as.double(qr$qraux),
        y,
        ny,
        qty=mat.or.vec(n,ny))$qty
}
qr.resid <- function(qr, y)
{
    if(!is.qr(qr)) stop("argument is not a QR decomposition")
    n <- as.integer(nrow(qr$qr))
    p <- as.integer(ncol(qr$qr))
    k <- as.integer(qr$rank)
    y <- as.matrix(y)
    ny <- as.integer(ncol(y))
    storage.mode(y) <- "double"
    if( nrow(y) != n )
        stop("qr and y must have the same number of rows")
    .Fortran("dqrrsd",
        as.double(qr$qr),
        n, k,
        as.double(qr$qraux),
        y,
        ny,
        rsd=mat.or.vec(n,ny))$rsd
}
qr.fitted <- function(qr, y, k=qr$rank)
{
    if(!is.qr(qr)) stop("argument is not a QR decomposition")
    n <- as.integer(nrow(qr$qr))
    p <- as.integer(ncol(qr$qr))
    k <- as.integer(k)
    if(k > qr$rank) stop("k is too large")
    y <- as.matrix(y)
    ny <- as.integer(ncol(y))
    storage.mode(y) <- "double"
    if( nrow(y) != n )
        stop("qr and y must have the same number of rows")
    .Fortran("dqrxb",
        as.double(qr$qr),
        n, k,
        as.double(qr$qraux),
        y,
        ny,
        xb=mat.or.vec(n,ny))$xb
}
## qr.solve is defined in 'solve'
quantile <- function (x, probs = seq(0, 1, 0.25), na.rm = FALSE)
{
    if (na.rm)
        x <- x[!is.na(x)]
    else if (any(!is.finite(x)))
        stop("Missing values, NaN\'s and Inf\'s not allowed if na.rm=FALSE")
    n <- length(x)
    if(any(probs < 0 | probs > 1))
        stop("probs outside [0,1]")
    index <- 1 + (n - 1) * probs
    lo <- floor(index)
    hi <- ceiling(index)
    x <- sort(x, partial=unique(c(lo,hi)))
    qs <- x[lo] + (x[hi] - x[lo]) * (index - lo)
    names(qs) <- paste(round(100*probs), "%", sep="")
    qs
}
quit <- function(save = "ask")
.Internal(quit(save))
q<-quit
range <- function(..., na.rm=FALSE)
c(min(..., na.rm=na.rm),max(..., na.rm=na.rm))
"count.fields" <-
function(file, sep = "", skip = 0)
{
    .Internal(count.fields(file, sep, skip))
}
"read.table" <-
function (file, header=FALSE, sep="", row.names, col.names, as.is=FALSE,
    na.strings="NA", skip=0)
{
    type.convert <-
    function(x, na.strings="NA", as.is=FALSE)
        .Internal(type.convert(x, na.strings, as.is))
        #  basic column counting and header determination
        #  we record whether it looks like we column names
        #  rlabp == 1
    row.lens <- count.fields(file, sep, skip)
    nlines <- length(row.lens)
    rlabp <- 0
    if (nlines > 1 && (row.lens[2] - row.lens[1]) == 1) {
        if (missing(header))
            header <- TRUE
        rlabp <- 1
    }
        # read in the header
    if (header) {
        col.names <- scan(file, what="", sep=sep, nlines=1, quiet=TRUE, skip=skip)
        skip <- skip + 1
        row.lens <- row.lens[-1]
        nlines <- nlines - 1
    }
    else if (missing(col.names))
        col.names <- paste("V", 1:row.lens[1], sep="")
        #  check that all rows have equal lengths
    cols <- unique(row.lens)
    if (length(cols) != 1)
        stop("all rows must have the same length",
            paste(row.lens, sep=","))
        #  set up for the scan of the file
        #  we read all values as character strings
        #  and convert later.
    what <- rep(list(""), cols)
    if (rlabp == 1)
        col.names <- c("row.names", col.names)
    names(what) <- col.names
    data <- scan(file=file, what=what, sep=sep, n=nlines*cols, skip=skip,
            na.strings=na.strings, quiet=TRUE)
        #  ok, now we have the data
        #  we now we convert to numeric or factor variables
        #  (depending on the specifies value of "as.is")
        #  we do this here so that columns match up
    cols <- length(data)
    if(is.logical(as.is)) {
        as.is <- rep(as.is, length=cols)
    }
    else if(is.numeric(as.is)) {
        if(any(as.is < 1 | as.is > cols))
            stop("invalid numeric as.is expression")
        i <- rep(FALSE, cols)
        i[as.is] <- TRUE
        as.is <- i
    }
    if (length(as.is) != cols)
        stop("as.is has the wrong length")
    for (i in 1:cols) {
        if (!as.is[i])
            data[[i]] <- type.convert(data[[i]], na.strings = na.strings)
    }
        #  now we determine row names
    if (missing(row.names)) {
        if (rlabp == 1) {
            row.names <- data[[1]]
            data <- data[-1]
        }
        else row.names <- as.character(1:nlines)
    }
    else if (is.null(row.names))
        row.names <- as.character(1:nlines)
    else if (is.character(row.names)) {
        if (length(row.names) == 1) {
            rowvar <- (1:cols)[match(col.names, row.names, 0) == 1]
            row.names <- data[[rowvar]]
            data <- data[-rowvar]
        }
    }
    else if (is.numeric(row.names) && length(row.names) == 1) {
        rlabp <- row.names
        row.names <- data[[rlabp]]
        data <- data[-rlabp]
    }
    else stop("invalid row.names specification")
        #  this is extremely underhanded
        #  we should use the constructor function ...
        #  don't try this at home kids
    class(data) <- "data.frame"
    row.names(data) <- row.names
    return(data)
}
rect <- function(xleft, ybottom, xright, ytop, col=NULL, border=par("fg"), lty=NULL, xpd=FALSE) {
    .Internal(rect(
        as.double(xleft),
        as.double(ybottom),
        as.double(xright),
        as.double(ytop),
        col=col,
        border=border,
        lty=lty,
        xpd=xpd))
}
rm <- function(...,list=character(0), envir=NULL,inherits=FALSE) {
    names<- as.character(substitute(list(...)))[-1]
    list<-c(list,names)
    .Internal(remove(list,envir,inherits))
}
rep <- function(x, times, length.out)
{
    if (length(x) == 0)
        return(x)
    if (missing(times))
        times <- ceiling(length.out/length(x))
    r<-.Internal(rep(x,times))
    if (!missing(length.out))
        return(r[1:length.out])
    return(r)
}
replace <-
function (x, list, values)
{
        x[list] <- values
        x
}
require <-
function (name, quietly = FALSE)
{
    if (missing(name))
        return(TRUE)
    if (!exists(".Libraries", inherits = TRUE))
        assign(".Libraries", character(0), NULL)
    if (!exists(".Provided", inherits=TRUE))
        assign(".Provided", character(0), NULL)
    name <- substitute(name)
    if (!is.character(name))
        name <- deparse(name)
    if (is.na(match(name, .Libraries))&& is.na(match(name, .Provided))) {
        file <- system.file("library", name)
        if (file == "") {
            if (!quietly)
    warning(paste("Required library ", name, " not found.\n"))
            return(FALSE)
        }
        if (!quietly)
            cat("Autoloading required library:", name, "\n")
        sys.source(file)
        assign(".Libraries", c(name, .Libraries), NULL)
    }
    return(TRUE)
}
provide <-
function(name)
{
    if (!exists(".Libraries", inherits = TRUE))
        assign(".Libraries", character(0), NULL)
    if (!exists(".Provided",inherits=TRUE))
        assign(".Provided", character(0), NULL)
    if (missing(name))
        return(list(provide=.Provided, library=.Libraries))
    name<-substitute(name)
    if (!is.character(name))
        name<-deparse(name)
    if (is.na(match(name, .Libraries)) && is.na(match(name, .Provided))){
        assign(".Provided",c(name,.Provided),NULL)
        return(TRUE)
    }
    else
        return(FALSE)
}
rev <- function(x) x[length(x):1]
rm <- function(...,list=character(0), envir=NULL,inherits=FALSE) {
    names<- as.character(substitute(list(...)))[-1]
    list<-c(list,names)
    .Internal(remove(list,envir,inherits))
}
rownames <- function(x) {
    dn <- dimnames(x)
    if(is.null(dn)) dn else dn[[1]]
}
"rownames<-" <- function(x, value) {
    dn <- dimnames(x)
    if(is.null(dn)) dimnames(x) <- list(value, dn)
    else dimnames(x) <- list(value, dn[[2]])
    x
}
sample <- function(x, size, replace=FALSE)
{
    if(length(x) == 1 && x >= 1) {
        if(missing(size)) size <- x
        .Internal(sample(x, size, replace))
    }
    else {
        if(missing(size)) size <- length(x)
        x[.Internal(sample(length(x), size, replace))]
    }
}
sapply <- function(X, FUN, ..., simplify = TRUE)
{
  if(is.character(FUN))
    FUN <- get(FUN, mode = "function")
    else if(mode(FUN) != "function") {
      farg <- substitute(FUN)
      if(mode(farg) == "name")
    FUN <- get(farg, mode = "function")
    else stop(paste("\"", farg, "\" is not a function", sep = ""))
    }
  answer <- lapply(as.list(X), FUN, ...)
  if(simplify && length(answer) &&
     length(common.len <- unique(unlist(lapply(answer, length)))) == 1) {
    if(common.len == 1)
      unlist(answer, recursive = FALSE)
    else if(common.len > 1)
      array(unlist(answer, recursive = FALSE),
        c(common.len, length(X)), list(NULL, names(answer)))
    else answer
  } else answer
}
scan <- function(file="", what=0, nmax=-1, sep="", skip=0, nlines=0,
    na.strings="NA", flush=FALSE, strip.white=FALSE, quiet=FALSE) {
    if( !missing(sep) )
        na.strings<-c(na.strings,"")
    .Internal(scan(file, what, nmax, sep, skip, nlines, na.strings,flush,strip.white, quiet))
}
sd <- function(x, na.rm=FALSE) sqrt(var(x, na.rm=na.rm))
segments <- function(x0, y0, x1, y1, col=par("fg"), lty=par("lty"))
    .Internal(segments(x0, y0, x1, y1, col=col, lty=lty))
seq <- function(from = 1, to = 1, by = ((to - from)/(length.out - 1)), length.out = NULL, along.with = NULL) {
    if(!missing(along.with))
        length.out <- length(along.with)
    else if(!missing(length.out))
        length.out <- ceiling(length.out)
    if(nargs() == 1 && !missing(from)) {
        if(mode(from) == "numeric" && length(from) == 1)
            1:from
        else seq(along.with = from)
    }
    else if(is.null(length.out))
        if(missing(by))
            from:to
        else {
            n <- (to - from)/by
            if(n < 0)
                stop("Wrong sign in by= argument")
            from + (0:n) * by
        }
    else if(length.out < 0)
        stop("Length cannot be negative")
    else if(length.out == 0)
        numeric(0)
    else if(missing(by)) {
        if(from == to || length.out < 2)
            by <- 1
        if(missing(to))
            to <- from + length.out - 1
        if(missing(from))
            from <- to - length.out + 1
        if(length.out > 2)
            if(from == to)
                rep(from, length.out)
            else as.vector(c(from, from + (1:(length.out - 2)) *
                    by, to))
        else as.vector(c(from, to))[1:length.out]
    }
    else if(missing(to))
        from + (0:(length.out - 1)) * by
    else if(missing(from))
        to - ((length.out - 1):0) * by
    else stop("Too many arguments")
}
sequence <- function(nvec)
{
    sequence <- NULL
    for(i in nvec) sequence<-c(sequence,seq(1:i))
    return(sequence)
}
qr.solve <- function(a,b, tol = 1e-7)
{
    if( !is.qr(a) )
        a <- qr(a, tol = tol)
    nc <- ncol(a$qr)
    if( a$rank != nc )
        stop("singular matrix 'a' in solve")
    if( missing(b) ) {
        if( nc != nrow(a$qr) )
            stop("only square matrices can be inverted")
        b<-diag(1,nc)
    }
    b<-as.matrix(b)
    return(qr.coef(a,b))
}
solve.qr <- qr.solve
solve <- qr.solve
sort <- function(x, partial=NULL, na.last=NA) {
    nas <- x[is.na(x)]
    x <- x[!is.na(x)]
    if(!is.null(partial))
        y <- .Internal(psort(x, partial))
    else {
        nms <- names(x)
        if(!is.null(nms)) {
            o <- order(x)
            y <- x[o]
            names(y) <- nms[o]
        }
        else
            y <- .Internal(sort(x))
    }
    if(!is.na(na.last)) {
        if(!na.last) y <- c(nas, y)
        else if (na.last) y <- c(y, nas)
    }
    y
}
source <- function (file, local=FALSE) {
    if(local) envir <- sys.frame(sys.parent())
    else if(!local) envir <- .GlobalEnv
    exprs <- parse(n = -1, file = file)
    if (length(exprs) == 0) return(invisible())
    for (i in exprs) {
        yy <- eval(i, envir)
    }
    invisible(yy)
}
sys.source <- function (file) {
    exprs <- parse(n = -1, file = file)
    if (length(exprs) == 0) return(invisible())
    for (i in exprs) {
        yy <- eval(i, NULL)
    }
    invisible(yy)
}
spline <- function(x, y, n=3*length(x), method="fmm", xmin=min(x), xmax=max(x)) {
    method <- match(method, c("periodic", "natural", "fmm"))
    if(is.na(method))
        stop("invalid interpolation method")
    if(length(x) != length(y))
        stop("x and y lengths differ in spline")
    if( !is.numeric(x) || !is.numeric(y) )
        stop("spline: x and y must be numeric")
    if(any(diff(x) <= 0))
        stop("invalid x array in spline")
    if(method == 1 && y[1] != y[length(y)]) {
        warning("first and last y values differ in spline - using y[1] for both")
        y[length(y)] <- y[1]
    }
    z <- .C("spline_coef",
        method=as.integer(method),
        n=length(x),
        x=as.double(x),
        y=as.double(y),
        b=double(length(x)),
        c=double(length(x)),
        d=double(length(x)),
        e=if(method == 1) double(length(x)) else double(0))
    u <- seq(xmin, xmax, length.out=n)
    .C("spline_eval",
        z$method,
        length(u),
        x=u,
        y=double(n),
        z$n,
        z$x,
        z$y,
        z$b,
        z$c,
        z$d)[3:4]
}
splinefun <- function(x, y, method="fmm") {
    method <- match(method, c("periodic", "natural", "fmm"))
    if(is.na(method))
        stop("invalid interpolation method")
    if(length(x) != length(y))
        stop("x and y lengths differ in spline")
    if( !is.numeric(x) || !is.numeric(y) )
        stop("splinefun: both x and y must be numeric")
    if(any(diff(x) <= 0))
        stop("invalid x array in spline")
    if(method == 1 && y[1] != y[length(y)]) {
        warning("first and last y values differ in spline - using y[1] for both")
        y[length(y)] <- y[1]
    }
    z <- .C("spline_coef",
        method=as.integer(method),
        n=length(x),
        x=as.double(x),
        y=as.double(y),
        b=double(length(x)),
        c=double(length(x)),
        d=double(length(x)),
        e=if(method == 1) double(length(x)) else double(0))
    rm(x,y,method)
    function(x) {
        .C("spline_eval",
            z$method,
            length(x),
            x=as.double(x),
            y=double(length(x)),
            z$n,
            z$x,
            z$y,
            z$b,
            z$c,
            z$d)$y
    }
}
split <- function(x, f) .Internal(split(x, as.factor(f)))
stem <- function(x,scale=1, width=80, atom=0.00000001) {
    if( !is.numeric(x) )
        stop("stem: x must be numeric")
    x <- x[!is.na(x)]
    if(length(x)==0) stop("no non-missing values")
    .C("stemleaf", as.double(x), length(x), as.double(scale), as.integer(width), as.double(atom))
    invisible(NULL)
}
# Dotplots a la Box, Hunter and Hunter
stripplot <- function(x, method="overplot", jitter=0.1, offset=1/3,
        vertical=FALSE, group.names,
        xlim, ylim, main="", ylab="", xlab="",
        pch=0, col=par("fg"), cex=par("cex"))
{
    method <- pmatch(method, c("overplot", "jitter", "stack"))[1]
    if(is.na(method) || method==0)
        error("invalid plotting method")
    if(is.language(x)) {
        if(length(x) == 3 && deparse(x[[1]]) == '~') {
            groups <- eval(x[[3]], sys.frame(sys.parent()))
            x <- eval(x[[2]], sys.frame(sys.parent()))
            groups <- split(x, groups)
        }
        else stop("invalid first argument")
    }
    else if(is.list(x)) {
        groups <- x
    }
    else if(is.numeric(x)) {
        groups <- list(x)
    }
    n <- length(groups)
    if(!missing(group.names)) attr(groups, "names") <- group.names
    else if(is.null(attr(groups, "names"))) attr(groups, "names") <- 1:n
    dlim <- rep(NA, 2)
    for(i in groups)
        dlim <- range(dlim, i, na.rm=T)
    glim <- c(1, n)
    if(method == 2) {
        if(n == 1) glim <- glim + c(-5, 5) * jitter
        else glim <- glim + c(-2, 2) * jitter
    }
    if(method == 3) {
        if(n == 1) glim <- glim + c(-1,1)
        else glim <- glim + c(0, 0.5)
    }
    if(missing(xlim)) {
        if(vertical) xlim <- glim
        else xlim <- dlim
    }
    if(missing(ylim)) {
        if(vertical) ylim <- dlim
        else ylim <- glim
    }
    plot(xlim, ylim, type="n", ann=FALSE, axes=FALSE)
    box()
    if(vertical) {
        if(n > 1) axis(1, at=1:n, lab=names(groups))
        axis(2)
    }
    else {
        axis(1)
        if(n > 1) axis(2, at=1:n, lab=names(groups))
    }
    if(vertical) csize <- cex*xinch(par("cin")[1])
    else csize <- cex*yinch(par("cin")[2])
    f <- function(x) seq(length(x))
    for(i in 1:length(groups)) {
        x <- groups[[i]]
        y <- rep(i,length(x))
        if(method == 2) y <- y + runif(length(y), -jitter, jitter)
        if(method == 3) {
            xg <- split(x, factor(x))
            xo <- lapply(xg, f)
            x <- unlist(xg)
            y <- y + (unlist(xo) - 1) * offset * csize
        }
        if(vertical) points(y, x, col=col[(i - 1)%%length(col) + 1],
            pch=pch[(i - 1)%%length(pch) + 1], cex=cex)
        else points(x, y, col=col[(i - 1)%%length(col) + 1],
            pch=pch[(i - 1)%%length(pch) + 1], cex=cex)
    }
    title(main=main, xlab=xlab, ylab=ylab)
}
"structure" <-
function (.Data, ...)
{
    specials <- c(".Dim", ".Dimnames", ".Names", ".Tsp", ".Label")
    replace <- c("dim", "dimnames", "names", "tsp", "levels")
    attrib <- list(...)
    if(!is.null(attrib)) {
        m <- match(names(attrib), specials)
        ok <- (!is.na(m) & m > 0)
        names(attrib)[ok] <- replace[m[ok]]
        if(any(names(attrib) == "tsp"))
            attrib$class <- unique(c("ts", attrib$class))
        if(is.numeric(.Data) && any(names(attrib) == "levels"))
            .Data <- factor(.Data)
        attributes(.Data) <- c(attributes(.Data), attrib)
    }
    return(.Data)
}
strwidth <- function(s, units="user", cex=NULL)
    .Internal(strwidth(s, pmatch(units, c("user", "figure", "inches")), cex))
sum <- function(..., na.rm=FALSE)
.Internal(sum(...,na.rm=na.rm))
min <- function(..., na.rm=FALSE)
.Internal(min(...,na.rm=na.rm))
max <- function(..., na.rm=FALSE)
.Internal(max(...,na.rm=na.rm))
prod <- function(...,na.rm=FALSE)
.Internal(prod(...,na.rm=na.rm))
summary <-
function (x, ...)
    UseMethod("summary")
summary.default <-
function(object, ..., digits = max(options()$digits - 3, 3))
{
    if(is.factor(object))
        return(summary.factor(object, ...))
    else if(is.matrix(object))
        return(summary.matrix(object, ...))
    value <- if(is.numeric(object)) {
        nas <- is.na(object)
        object <- object[!nas]
        qq <- quantile(object)
        qq <- signif(c(qq[1:3], mean(object), qq[4:5]), digits)
        names(qq) <- c("Min.", "1st Qu.", "Median", "Mean", "3rd Qu.", "Max.")
        if(any(nas))
            c(qq, "NA's" = sum(nas))
        else qq
    } else if(is.recursive(object) && !is.language(object) &&
        (n <- length(object))) {
        sumry <- array("", c(n, 3), list(names(object),
            c("Length", "Class", "Mode")))
        ll <- numeric(n)
        for(i in 1:n) {
            ii <- object[[i]]
            ll[i] <- length(ii)
            sumry[i, 2] <- if(is.object(ii)) class(ii) else "-none-"
            sumry[i, 3] <- mode(ii)
        }
        sumry[, 1] <- format(as.integer(ll))
        class(sumry) <- "table"
        sumry
    }
    else c(Length = length(object), Class = class(object), Mode = mode(object))
    class(value) <- "table"
    value
}
summary.factor <-
function(x, maxsum = 100, ...)
{
    nas <- is.na(x)
    ll <- levels(x)
    if(any(nas)) maxsum <- maxsum - 1
    tbl <- table(x)
    tt <- c(tbl) # names dropped ...
    names(tt) <- dimnames(tbl)[[1]]
    if(length(ll) > maxsum) {
        drop <- maxsum:length(ll)
        o <- rev(order(tt))
        tt <- c(tt[o[ - drop]], "(Other)" = sum(tt[o[drop]]))
    }
    if(any(nas)) c(tt, "NA's" = sum(nas)) else tt
}
summary.matrix <-
function(x, ...)
    summary.data.frame(data.frame(x))
summary.data.frame <-
function(x, maxsum = 7, ...)
{
    z <- lapply(as.list(x), summary, maxsum = maxsum)
    nv <- length(x)
    nm <- names(x)
    lw <- numeric(nv)
    nr <- max(unlist(lapply(z, length)))
    for(i in 1:nv) {
        sms <- z[[i]]
        lbs <- format(names(sms))
        sms <- paste(lbs, ":", format(sms), "  ", sep = "")
        lw[i] <- nchar(lbs[1])
        length(sms) <- nr
        z[[i]] <- sms
    }
    z <- unlist(z, use.names=FALSE)
    dim(z) <- c(nr, nv)
     blanks <- paste(character(max(lw) + 2), collapse = " ")
     pad <- floor(lw-nchar(nm)/2)
     nm <- paste(substring(blanks, 1, pad), nm, sep = "")
     dimnames(z) <- list(rep("", nr), nm)
    attr(z, "class") <- c("table") #, "matrix")
    z
}
print.table <-
function(x, digits= .Options$digits, quote = FALSE, na.print='', ...)
{
 print.default(unclass(x), digits=digits, quote=quote, na.print=na.print, ...)
}
svd <- function(x, nu=min(n,p), nv=min(n,p)) {
    if(!is.numeric(x))
        stop("argument to svd must be numeric")
    x <- as.matrix(x)
    dx <- dim(x)
    n <- dx[1]
    p <- dx[2]
    if(nu == 0) {
        job <- 0
        u <- double(0)
    }
    else if(nu == n) {
        job <- 10
        u <- matrix(0, n, n)
    }
    else if(nu == p) {
        job <- 20
        u <- matrix(0, n, p)
    }
    else
        stop("nu must be 0, nrow(x) or ncol(x)")
    if(nv == 0)
        job <- job + 0
    else if(nv == p || nv == n)
        job <- job + 1
    else
        stop("nv must be 0 or ncol(x)")
    if(job == 0)
        v <- double(0)
    else
        v <- matrix(0, p, p)
    mn <- min(n,p)
    mm <- min(n+1,p)
    z <- .Fortran("dsvdc",
        x,
        n,
        n,
        p,
        d=double(mm),
        double(p),
        u=u,
        n,
        v=v,
        p,
        double(n),
        as.integer(job),
        info=integer(1),
        DUP=FALSE)[c("d","u","v","info")]
    if(z$info)
        stop(paste("error ",z$info," in dsvdc"))
    z$d <- z$d[1:mn]
    if(nv && nv < p) z$v <- z$v[, 1:nv]
    z[c("d", if(nu) "u" else NULL, if(nv) "v" else NULL)]
}
sweep <-
function(x, MARGIN, STATS, FUN = "-", ...)
{
    if(is.character(FUN))
        FUN <- get(FUN)
    dims <- dim(x)
    perm <- c(MARGIN, (1:length(dims))[ - MARGIN])
    FUN(x, aperm(array(STATS, dims[perm]), order(perm)), ...)
}
switch <- function(EXPR,...)
    .Internal(switch(EXPR,...))
sys.call <-function(which=0)
 .Internal(sys.call(which))
sys.calls <-function()
 .Internal(sys.calls())
sys.frame <-function(which=0)
 .Internal(sys.frame(which))
sys.function <-function(which=0)
 .Internal(sys.function(which))
sys.frames <-function()
 .Internal(sys.frames())
sys.nframe <- function()
 .Internal(sys.nframe())
sys.parent <- function(n = 1)
 .Internal(sys.parent(n))
sys.parents <- function()
 .Internal(sys.parents())
sys.status <- function()
 list(sys.calls=sys.calls(), sys.parents=sys.parents(), sys.frames=sys.frames())
sys.on.exit <- function()
 .Internal(sys.on.exit())
system <- function(call, intern=FALSE)
    .Internal(system(call,intern))
system.file <- function(dir, name)
{
    system(paste("$RHOME/cmd/filename", dir, name), intern=TRUE)
}
table <- function(x, ...)
{
    if (nargs() == 0)
        stop("no arguments")
    bin <- 0
    lens <- NULL
    dims <- integer(0)
    pd <- 1
    dn <- NULL
    if(nargs() == 1 && is.list(x))
        args <- x
    else
        args <- list(x, ...)
    for (a in args) {
        if (is.null(lens))
            lens <- length(a)
        if (length(a) != lens)
            stop("all arguments must have the same length")
        if (!is.factor(a))
            cat <- as.factor(a)
        else cat <- a
        l <- levels(cat)
        dims <- c(dims, nlevels(cat))
        dn <- c(dn, list(l))
        bin <- bin + pd * (as.integer(cat) - 1)
        pd <- pd * nlevels(cat)
    }
    bin <- bin[!is.na(bin)]
    array(tabulate(bin + 1, pd), dims, dimnames = dn)
}
tabulate <- function(bin, nbins = max(bin))
{
    if(!is.numeric(bin) && !is.factor(bin))
        stop("tabulate: bin must be numeric or a factor")
    if((n <- length(bin)) == 0) bin <- 1
    else bin <- as.integer(bin)
    .C("tabulate",
        ans = integer(nbins),
        bin,
        n)$ans
}
tapply <- function (x, INDEX, FUN, ...)
{
    if (is.character(FUN))
        FUN <- get(FUN, mode = "function")
    if (mode(FUN) != "function")
        stop(paste("\"", FUN, "\" is not a function"))
    if (!is.list(INDEX)) INDEX <- list(INDEX)
    namelist <- vector("list", length(INDEX))
    extent <- integer(length(INDEX))
    nx <- length(x)
    group <- rep(1, nx)
    ngroup <- 1
    for (i in seq(INDEX)) {
        index <- as.factor(INDEX[[i]])
        if (length(index) != nx)
            stop("arguments must have same length")
        namelist[[i]] <- levels(index)
        extent[[i]] <- nlevels(index)
        group <- group + ngroup * (codes(index) - 1)
        ngroup <- ngroup * nlevels(index)
    }
    if (missing(FUN)) return(group)
    ans <- lapply(split(x, group), FUN, ...)
    if(all(unlist(lapply(ans,length))==1))
        ans <- unlist(ans, recursive=FALSE)
    if(length(INDEX) == 1) {
        names(ans) <- namelist[[1]]
    }
    else {
        dim(ans) <- extent
        dimnames(ans) <- namelist
    }
    return(ans)
}
text <- function(x, y=NULL, text, ...)
    .Internal(text(xy.coords(x,y), as.character(text), ...))
title <- function(main=NULL, sub=NULL, xlab=NULL, ylab=NULL, ...)
    .Internal(title(
        as.character(main),
        as.character(sub),
        as.character(xlab),
        as.character(ylab),
        ...
    ))
traceback <- function() as.character(.Traceback)
trunc <- function(x) ifelse(x<0,ceiling(x),floor(x))
start <-
function(x, ...)
UseMethod("start")
end <-
function(x, ...)
UseMethod("end")
frequency <-
function(x, ...)
UseMethod("frequency")
time <-
function(x, ...)
UseMethod("time")
window <-
function(x, ...)
UseMethod("window")
ts <-
function(data=NA, start=1, end=numeric(0), frequency=1, deltat=1)
{
    ts.eps <- .Options$ts.eps
    if(is.null(ts.eps)) ts.eps <- 1.e-6
    if(is.matrix(data)) {
        nseries <- ncol(data)
        ndata <- nrow(data)
    }
    else {
        nseries <- 1
        ndata <- length(data)
    }
    if(missing(frequency)) frequency <- 1/deltat
    if(missing(deltat)) deltat <- 1/deltat
    if(frequency > 1 && abs(frequency - round(frequency)) < ts.eps)
        frequency <- round(frequency)
    if(length(start) > 1) {
        if(start[2] > frequency) stop("invalid start")
        start <- start[1] + (start[2] - 1)/frequency
    }
    if(length(end) > 1) {
        if(end[2] > frequency) stop("invalid end")
        end <- end[1] + (end[2] - 1)/frequency
    }
    if(missing(end))
        end <- start + (ndata - 1)/frequency
    else if(missing(start))
        start <- end - (ndata - 1)/frequency
    nobs <- floor((end - start) * frequency + 1.01)
    if(nseries == 1) {
        if(ndata < nobs)
            data <- rep(data, length=nobs)
        else if(nobs > ndata)
            data <- data[1:nobs]
    }
    else {
        if(ndata < nobs)
            data <- data[rep(1:ndata, length=nobs)]
        else if(nobs > ndata)
            data <- data[1:nobs,]
    }
    attr(data, "tsp") <- c(start, end, frequency)
    attr(data, "class") <- "ts"
    data
}
tsp <-
function(x)
    attr(x, "tsp")
"tsp<-" <-
function(x, tsp)
{
    attr(x,"tsp") <- tsp
    class(x) <- "ts"
    x
}
is.ts <-
function (x)
inherits(x, "ts")
as.ts <-
function (x)
if (is.ts(x)) x else ts(x)
start.ts <-
function(x)
{
    ts.eps <- .Options$ts.eps
    if(is.null(ts.eps)) ts.eps <- 1.e-6
    tsp <- attr(as.ts(x), "tsp")
    is <- tsp[1]*tsp[3]
    if(abs(is-round(is)) < ts.eps) {
        is <- floor(tsp[1])
        fs <- floor(tsp[3]*(tsp[1] - is)+0.001)
        c(is,fs+1)
    }
    else ts[1]
}
end.ts <-
function(x)
{
    ts.eps <- .Options$ts.eps
    if(is.null(ts.eps)) ts.eps <- 1.e-6
    tsp <- attr(as.ts(x), "tsp")
    is <- tsp[2]*tsp[3]
    if(abs(is-round(is)) < ts.eps) {
        is <- floor(tsp[2])
        fs <- floor(tsp[3]*(tsp[2] - is)+0.001)
        c(is, fs+1)
    }
    else ts[2]
}
frequency.ts <-
function(x)
{
    tsp <- attr(as.ts(x), "tsp")
    tsp[3]
}
time.ts <-
function (x)
{
    x <- as.ts(x)
    if(is.matrix(x)) n <- nrow(x)
    else n <- length(x)
    xtsp <- attr(x, "tsp")
    ts(seq(xtsp[1], xtsp[2], length=n),
        start=start(x), end=end(x), frequency=frequency(x))
}
print.ts <- function(x, calender, ...)
{
    if(missing(calender))
        calender <- any(frequency(x)==c(4,12))
    if(all(frequency(x) != c(1,4,12)))
        calender <- FALSE
    if(!is.matrix(x) && calender) {
        if(frequency(x) == 12) {
            start.pad <- start(x)[2] - 1
            end.pad <- 12 - end(x)[2]
            data <- matrix(c(rep(NA, start.pad), x,
                rep(NA, end.pad)), nc=12, byrow=T)
            dimnames(data) <- list(
                as.character(start(x)[1]:end(x)[1]),
                month.abb)
        }
        else if(frequency(x) == 4) {
            start.pad <- start(x)[2] - 1
            end.pad <- 4 - end(x)[2]
            data <- matrix(c(rep(NA, start.pad), x,
                rep(NA, end.pad)), nc=4, byrow=T)
            dimnames(data) <- list(
                    paste(start(x)[1]:end(x)[1], ":" , sep=""),
                    c("Qtr1", "Qtr2", "Qtr3", "Qtr4"))
        }
        else if(frequency(x) == 1) {
            data <- x
            attributes(data) <- NULL
            names(data) <- time(x)
        }
    }
    else  {
        cat("Time-Series:\nStart =", deparse(start(x)),
            "\nEnd =", deparse(end(x)),
            "\nFrequency =", deparse(frequency(x)), "\n")
        data <- x
        attr(data, "tsp") <- NULL
        attr(data, "class") <- NULL
        # something like this is needed here
        # if(is.matrix(x)) rownames(data) <- time(x)
    }
    print(data, ...)
    invisible(x)
}
plot.ts <-
function (x, type="l", xlim, ylim, xlab, ylab, log="",
    col=par("col"), bg=NA, pch=par("pch"), lty=par("lty"), ...)
{
    time.x <- time(x)
    if(missing(xlim)) xlim <- range(time.x)
    if(missing(ylim)) ylim <- range(x, na.rm=TRUE)
    if(missing(xlab)) xlab <- "Time"
    if(missing(ylab)) ylab <- deparse(substitute(x))
    plot.new()
    plot.window(xlim, ylim, log)
    if(is.matrix(x)) {
        for(i in 1:ncol(x))
            lines.default(time.x, x[,i],
                col=col[(i-1)%%length(col) + 1],
                lty=lty[(i-1)%%length(lty) + 1],
                bg=bg[(i-1)%%length(bg) + 1],
                pch=pch[(i-1)%%length(pch) + 1],
                type=type)
    }
    else {
        lines.default(time.x, x, col=col[1], bg=bg, lty=lty[1],
            pch=pch[1], type=type)
    }
    title(xlab=xlab, ylab=ylab, ...)
    axis(1, ...)
    axis(2, ...)
    box(...)
}
window.ts <-
function(x, start, end)
{
    x <- as.ts(x)
    xtsp <- tsp(x)
    freq <- xtsp[3]
    xtime <- time(x)
    ts.eps <- .Options$ts.eps
        if (is.null(ts.eps))
                ts.eps <- 1e-06
    if(missing(start))
        start <- xtsp[1]
    else start <- switch(length(start),
            start,
            start[1] + (start[2] - 1)/freq,
            stop("Bad value for start"))
    if(start < xtsp[1]) {
        start <- xtsp[1]
        warning("start value not changed")
    }
    if(missing(end))
        end <- xtsp[2]
    else end <- switch(length(end),
            end,
            end[1] + (end[2] - 1)/freq,
            stop("Bad value for end"))
    if(end > xtsp[2]) {
        end <- xtsp[2]
        warning("end value not changed")
    }
    if(start > end)
        stop("start cannot be after end")
    if(all(abs(start - xtime) > abs(start) * ts.eps)) {
        start <- xtime[(xtime > start) & ((start + 1/freq) > xtime)]
    }
    if(all(abs(end - xtime) > abs(end) * ts.eps)) {
        end <- xtime[(xtime < end) & ((end - 1/freq) < xtime)]
    }
    if(is.matrix(x))
        x <- x[(trunc((start - xtsp[1]) * freq + 1.5):trunc((end -
            xtsp[1]) * freq + 1.5)), , drop = F]
    else x <- x[trunc((start - xtsp[1]) * freq + 1.5):trunc((end - xtsp[1])
            * freq + 1.5)]
    tsp(x) <- c(start, end, freq)
    x
}
"[.ts" <-
function(x, i, j)
{
    y <- NextMethod("[")
    if(is.matrix(x) & missing(i))
        ts(y, start=start(x), freq=frequency(x))
    else y
}
t.test <- function(x, y=NULL, alternative="two.sided",mu=0, paired = FALSE, var.equal = FALSE,  conf.level = 0.95) {
    choices<-c("two.sided","greater","less")
    alt<- pmatch(alternative,choices)
    alternative<-choices[alt]
    if( length(alternative)>1 || is.na(alternative) )
        stop("alternative must be one \"greater\", \"less\", \"two.sided\"")
    if( !missing(mu) )
        if( length(mu) != 1  || is.na(mu) )
            stop("mu must be a single number")
    if( !missing(conf.level) )
        if( length(conf.level) !=1 || is.na(conf.level) || conf.level<0 || conf.level > 1)
            stop("conf.level must be a number between 0 and 1")
    if( !is.null(y) ) {
        dname<-paste(deparse(substitute(x)),"and",paste(deparse(substitute(y))))
        if(paired)
            xok<-yok<-complete.cases(x,y)
        else {
            yok<-!is.na(y)
            xok<-!is.na(x)
        }
        y<-y[yok]
    }
    else {
        dname<-deparse(substitute(x))
        if( paired ) stop("y is missing for paired test")
        xok<-!is.na(x)
        yok<-NULL
    }
    x<-x[xok]
    if( paired ) {
        x<- x-y
        y<- NULL
    }
    nx <- length(x)
    if(nx <= 2) stop("not enough x observations")
    mx <- mean(x)
    vx <- var(x)
    estimate<-mx
    if(is.null(y)) {
        df <- length(x)-1
        stderr<-sqrt(vx/nx)
        tstat <- (mx-mu)/stderr
        method<-ifelse(paired,"Paired t-test","One Sample t-test")
        names(estimate)<-ifelse(paired,"mean of the differences","mean of x")
    } else {
        ny <- length(y)
        if(ny <= 2) stop("not enough y observations")
        my <- mean(y)
        vy <- var(y)
        method<-ifelse(var.equal,"Two Sample t-test","Welch Two Sample t-test")
        estimate<-c(mx,my)
        names(estimate)<-c("mean of x","mean of y")
        if(var.equal) {
            df <- nx+ny-2
            v <- ((nx-1)*vx + (ny-1)*vy)/df
            stderr <- sqrt(v*(1/nx+1/ny))
            tstat <- (mx-my-mu)/stderr
        } else {
            stderrx <-sqrt(vx/nx)
            stderry <-sqrt(vy/ny)
            stderr <- sqrt(stderrx^2 + stderry^2)
            df <- stderr^4/(stderrx^4/(nx-1) + stderry^4/(ny-1))
            tstat <- (mx - my - mu)/stderr
        }
    }
    if (alternative == "less") {
        pval <- pt(tstat, df)
        cint <- c(NA, tstat * stderr + qt(conf.level, df) * stderr)
    }
    else if (alternative == "greater") {
        pval <- 1 - pt(tstat, df)
        cint <- c(tstat * stderr - qt(conf.level, df) * stderr, NA)
    }
    else {
        pval <- 2 * pt(-abs(tstat), df)
        alpha <- 1 - conf.level
        cint <- c(tstat * stderr - qt((1 - alpha/2), df) * stderr,
            tstat * stderr + qt((1 - alpha/2), df) * stderr)
    }
    names(tstat)<-"t"
    names(df)<-"df"
    if(paired || !is.null(y) )
        names(mu)<-"difference in means"
    else
        names(mu)<- "mean"
    attr(cint,"conf.level")<-conf.level
    rval<-list(statistic = tstat, parameter = df, p.value = pval,
conf.int=cint, estimate=estimate, null.value = mu, alternative=alternative,
method=method, data.name=dname)
    attr(rval,"class")<-"htest"
    return(rval)
}
cm <- function(x) 2.54*x
xinch <- function(x=1)
    x * diff(par("usr")[1:2])/par("pin")[1]
yinch <- function(x=1)
    x * diff(par("usr")[3:4])/par("pin")[2]
unix.time <- function(expr)
{
  ## Purpose: Return CPU (and other) times that 'expr' used ..
  ##    Modelled after S`s  "unix.time"; to be used with R (rel. 0.4 & later)
  ## -------------------------------------------------------------------------
  ## Arguments: expr: 'any' valid R expression
  ## -------------------------------------------------------------------------
  if(!exists("proc.time", mode = "function", inherits=TRUE))
    stop(paste("proc.time  must be enabled at configuration / compile time\n",
           "    [add '-DProctime' to SYSTEM in src/Systems/<YOURSYS>]"))
  loc.frame <- sys.parent(1)
##-S   if(loc.frame == 1)
##-S    loc.frame <- F
  on.exit(cat("Timing stopped at:", proc.time() - time, "\n"))
  expr <- substitute(expr)
  time <- proc.time()
  eval(expr, envir = loc.frame) #<-- 'R'
  new.time <- proc.time()
  on.exit()
  if(length(new.time) == 3)    new.time <- c(new.time, 0, 0)
  if(length(time) == 3)        time     <- c(time, 0, 0)
  new.time - time
}
upper.tri <- function(x, diag = FALSE)
{
    x <- as.matrix(x)
        if(diag) row(x) <= col(x)
        else row(x) < col(x)
}
mat.or.vec <- function(nr,nc)
        if(nc==1) numeric(nr) else matrix(0,nr,nc)
var <- function(x, y=x, na.rm = FALSE) {
    if(is.matrix(x) || !missing(y)) cov(x,y)
    else {
        if (na.rm) x <- x[!is.na(x)]
        sum((x - mean(x))^2)/(length(x) - 1)
        }
}
logical <- function(n=0) vector("logical",n)
integer <- function(n=0) vector("integer",n)
real <- function(n=0) vector("real", n)
double <- function(n=0) vector("real", n)
numeric <- double
complex <- function(n=0, real=numeric(), imag=numeric())
.Internal(complex(n, real, imag))
character <- function(n=0) vector("character",n)
which <- function(x) {
    if(is.logical(x)) seq(x)[x]
    else stop("argument to \"which\" is not logical")
}
windows<- function() .Internal(device("Win32","",c(0,0,0)))
write <- function(x, file="data",ncolumns=if(is.character(x)) 1 else 5, append=FALSE)
        cat(x, file=file, sep=c(rep(" ",ncolumns-1), "\n"), append=append)
x11 <-
function(display="", width=7, height=7, ps=12,
    printcmd=options("printcmd")$printcmd,
    paper=options("papersize")$papersize,
    orientation="flexible")
{
    if(is.na(match(paper, c("none", "a4", "letter"))))
        stop("unsupported paper size in x11")
    orientation <- match(orientation,c("portrait", "landscape", "flexible"))
    if(is.na(orientation))
        stop("unknown page orientation in x11")
    .Internal(
        device("X11",
        as.character(c(display[1], paper)),
        as.double(c(width, height, ps, orientation))))
}