"predict.lm" <- 
function(object, newdata, type = c("response", "terms"), se.fit = F, terms = 
	labels(object))
{
	type <- match.arg(type)
	if(missing(newdata) && type != "terms" && !se.fit)
		return(object$fitted)
	Terms <- object$terms
	if(!inherits(Terms, "terms"))
		stop("invalid terms component of  object")
	offset <- attr(Terms, "offset")
	intercept <- attr(Terms, "intercept")
	xbar <- NULL
	if(missing(newdata)) {
		x <- model.matrix(object)
		#center x if terms are to be computed
		if(type == "terms" && intercept) {
			xbar <- apply(x, 2., mean)
			x <- sweep(x, 2., xbar)
		}
	}
	else if(!((is.atomic(newdata) && length(newdata) == 1. && length(object$
		coef) != 1. && newdata > 0. && (newdata - trunc(newdata) < 
		.Machine$single.eps)) | is.list(newdata))) {
		if(!is.null(offset)) {
			warning("Offset not included")
			#try and coerce newdata to look like the x matrix
			offset <- NULL
		}
		TT <- length(object$coef)
		if(is.matrix(newdata) && ncol(newdata) == TT)
			x <- newdata
		else if(length(newdata) == TT)
			x <- matrix(newdata, 1., TT)
		else stop("Argument \"newdata\" is not a data frame, and cannot be coerced to an appropriate model matrix"
				)
	}
	else {
		x <- model.matrix(delete.response(Terms), newdata, contrasts = 
			object$contrasts)
		#newdata is a list, data frame or frame number
		if(!is.null(offset)) offset <- eval(attr(Terms, "variables")[
				offset], newdata)
	}
	if(!missing(newdata) && type == "terms" && intercept) {
		xold <- model.matrix(object)
		#need to center x 
		xbar <- apply(xold, 2., mean)
		x <- sweep(x, 2., xbar)
	}
	coefs <- coef(object)
	asgn <- attr(coefs, "assign")
	if(type == "terms") {
		terms <- match.arg(terms, labels(object))
		asgn <- asgn[terms]
	}
	nac <- is.na(object$coef)
	if(any(nac)) {
		xbar <- xbar[!nac]
		x <- x[, !nac]
	}
	attr(x, "constant") <- xbar
	if(se.fit) {
		fit.summary <- summary.lm(object)
		pred <- Build.terms(x, coefs, fit.summary$cov * fit.summary$
			sigma^2., asgn, collapse = type != "terms")
		pred$residual.scale <- fit.summary$sigma
		pred$df <- object$df.resid
	}
	else pred <- Build.terms(x, coefs, NULL, assign = asgn, collapse = type !=
			"terms")
	if(!is.null(offset) && type != "terms") {
		if(missing(newdata))
			warning("Offset not included")
		else {
			if(se.fit)
				pred$fit <- pred$fit + offset
			else pred <- pred + offset
		}
	}
	pred
}

"print.message" <- 
function(obj)
{
	message(obj)
}

"ptri" <- 
function(x, a = 0., b = a + 2.)
{
	x <- (x - a)/(b - a)
	temp <- ifelse(x < 0.5, 2. * x^2., 1. - 2. * (1. - x)^2.)
	temp <- ifelse(x < 0., 0., temp)
	ifelse(x > 1., 1., temp)
}

"pweib" <- 
function(x, alpha = 1., beta = 1.)
{
	if(alpha < 0.5)
		stop(" Sorry you don't have a density if alpha is \nless than .5!"
			)
	x <- x/beta
	1. - exp( - x^alpha)
}

"qbinom" <- 
function(px, N, p)
{
	size <- N
	prob <- p
	q <- .Internal(qbinom(px, par1 = size, par2 = prob), "S_qfuns", T,
		14)
	if(!is.null(Names <- names(p)))
		names(q) <- rep(Names, length = length(q))
	q
}

"qexp" <- 
function(p, gamma = 1., beta = 1.)
{
	if(missing(gamma)) {
		gamma <- beta
	}
	q <- .Internal(qexp(p, par1 = 1./gamma), "S_qfuns", T, 4)
	if(!is.null(Names <- names(p)))
		names(q) <- rep(Names, length = length(q))
	q
}

"qgam" <- 
function(q, alpha, beta = 1.)
{
	qgamma(q, alpha) * beta
}

"qhnorm" <- 
function(p)
{
	qnorm((p * 0.5) + 0.5)
}

"qtri" <- 
function(p, a = 0., b = a + 2.)
{
	ifelse(p < 0.5, sqrt(p/2.), 1. - sqrt((1. - p)/2.)) * (b - a) + a
}

"random.walk" <- 
function(n, p = 0.5)
{
	cumsum(rbern(n, p) * 2. - 1.)
}

"rbern" <- 
function(n, p = 0.5)
{
	as.integer(ifelse(runif(n) <= p, 1., 0.))
}

"rbinom" <- 
function(n, N, p)
{
	size <- N
	prob <- p
	.Internal(rbinom(n, par1 = size, par2 = prob), "S_rfuns", T, 14)
}

"rdexp" <- 
function(n, mu = 0., sigma = 1.)
{
	x <- rexp(n) * sigma
	ifelse(rbern(n) == 0., x + mu, mu - x)
}

"read.data" <- 
function(file = "", what = double(0), n = -1, sep = "", multi.line = F, flush
	 = F, append = F, skip = 0, text = F, mess = F)
{
	if(text) {
		what <- "a"
	}
	if(mess) {
		sep <- "\r"
		what <- "a"
	}
	if(file == "") {
		print.message(c("If you forgot to give an output data set.\n",
			"That is, if you only typed:", 
			"                read.data()", "instead of say:", 
			"                read.data()-> zork.", 
			"enter two returns and start again specifying where the data\nshould go!",
			" ", "If you want to enter text you should have used:",
			"                read.data(text=T)-> zork", 
			"so that S knows these will not be numbers", "\r",
			"To make small corrections in your data after", 
			" you are finished use the subscripts:", 
			"  for example 165.9-> zork[15]", 
			" will change the 15th value of zork to 165.9", "   ",
			"Just enter your data after the colons. ", 
			"Enter a blank line to stop reading in data"))
	}
	if(skip != 0) {
		tempname <- tempfile("S.")
		on.exit(unlink(tempname))
		if(unix(paste("tail +", skip + 1, " ", file, ">", tempname,
			sep = ""), output = F))
			stop(paste("Problem skipping", skip, "lines of file",
				file, "or writing temp file", tempname))
		file <- tempname
	}
	data <- .Internal(scan(file, what, n, sep, multi.line, flush, append),
		"S_scan", T, 0)
	if(mess)
		class(data) <- "message"
	data
}

"rexp" <- 
function(n, gamma = 1., beta = 1.)
{
	if(!missing(beta)) {
		gamma <- beta
	}
 - gamma * log(runif(n))
}

"rgam" <- 
function(n, alpha, beta = 1.)
{
	rgamma(n, alpha) * beta
}

"rnorm.2d" <- 
function(n, rho = 0., mean = c(0., 0.), sd = c(1., 1.))
{
	if(abs(rho) > 1.)
		stop(" rho must be between -1 and 1 !")
	u <- rnorm(n)
	v <- u * rho + (rnorm(n)) * sqrt(1. - rho^2.)
	list(x = u * sd[1.] + mean[1.], y = v * sd[2.] + mean[2.])
}

"rperm" <- 
function(x)
{
	x[order(runif(length(x)))]
}

"runs" <- 
function(x)
{
	ind <- diff(c(0., x))
	runs <- c(cumsum(x)[ind == 1.], (sum(x) + 1.))
	as.integer(diff(runs))
}