"Uall.equal"<-
function(target, current, ...)
{
	if(is.all.equal(target, current, ...))
		TRUE
	else {
		newtar <- substitute(target)
		newcur <- substitute(current)
		if(is.name(newtar)) {
			if(!is.name(newcur))
				newobj <- deparse(newtar)
			else {
				warning("two names, don't know which to assign"
					)
				return(paste(newtar, newcur, sep = " ?? "))
			}
		}
		else if(is.name(newcur))
			newobj <- deparse(newcur)
		else {
			warning("No names to update")
			return("all.equal with no named objects")
		}
		assign(paste("N",newobj,sep=""), current, where = 1, immediate=T)
		newobj
	}
}
"Utestfit"<-
function(object, tolerance = .Machine$double.eps^0.5, data = list())
{
	fit <- object$call
	if(is.null(fit))
		fit <- attr(object, "call")
	if(!is.language(fit))
		stop("object must have the fitting expression as component \"call\""
			)	#nframe <- new.frame(data)
	nframe <- FALSE
	newfit <- eval(fit, nframe)
	if(is.all.equal(object, newfit, tolerance))
		TRUE
	else {
		newobj <- paste("N", deparse(substitute(object)), sep = "")
		assign(newobj, newfit, where = 1, immediate=T)
		newobj
	}
}
"bin.fam"<-
structure(.Data = list(family = structure(.Data = c("Binomial", 
	"Logit: log(mu/(1 - mu))", "Binomial: mu(1-mu)"), .Names = c("name", 
	"link", "variance")), names = "Logit: log(mu/(1 - mu))", link = 
	function(mu)
log(mu/(1 - mu)), inverse = function(eta)
{
	junk <- care.exp(eta)
	junk/(1 + junk)
}
, deriv = function(mu)
{
	d <- mu * (1 - mu)
	if(any(tiny <- (d < .Machine$double.eps))) {
		warning("Model unstable; fitted probabilities of 0 or 1")
		d[tiny] <- .Machine$double.eps
	}
	1/d
}
, initialize = expression({
	if(is.matrix(y)) {
		if(dim(y)[2] > 2)
			stop("only binomial response matrices (2 columns)")
		n <- drop(y %*% c(1, 1))
		y <- y[, 1]
	}
	else {
		if(is.category(y))
			y <- y != levels(y)[1]
		else y <- as.vector(y)
		n <- rep(1, length(y))
	}
	w <- w * n
	n[n == 0] <- 1
	y <- y/n
	mu <- y + (0.5 - y)/n
}
), variance = function(mu)
mu * (1 - mu), deviance = function(mu, y, w, residuals = F)
{
	devy <- y
	nz <- y != 0
	devy[nz] <- y[nz] * log(y[nz])
	nz <- (1 - y) != 0
	devy[nz] <- devy[nz] + (1 - y[nz]) * log(1 - y[nz])
	devmu <- y * log(mu) + (1 - y) * log(1 - mu)
	if(any(small <- mu * (1 - mu) < .Machine$double.eps)) {
		warning("fitted values close to 0 or 1")
		smu <- mu[small]
		sy <- y[small]
		smu <- ifelse(smu < .Machine$double.eps, .Machine$double.eps, 
			smu)
		onemsmu <- ifelse((1 - smu) < .Machine$double.eps, .Machine$
			double.eps, 1 - smu)
		devmu[small] <- sy * log(smu) + (1 - sy) * log(onemsmu)
	}
	devi <- 2 * (devy - devmu)
	if(residuals)
		sign(y - mu) * sqrt(abs(devi) * w)
	else sum(w * devi)
}
, weight = expression(w * mu * (1 - mu))), class = "family")
"coverage"<-
function(pos = 2, tests)
{
	if(is.numeric(pos))
		pos.names <- search()[pos]
	else pos.names <- as.character(pos)
	for(i in tests)
		do.test(i)
	pos <- match(pos.names, search())
	if(any(is.na(pos)))
		stop(paste(pos.names[is.na(pos)], "not in search"))
	if(length(pos) == 1) {
		zz <- database.dictionary(pos)
		structure(zz$used, names = zz$name)
	}
	else {
		value <- as.list(pos)
		names(value) <- pos.names
		for(i in seq(along = pos)) {
			zz <- database.dictionary(pos[i])
			value[[i]] <- structure(zz$used, names = zz$name)
		}
		value
	}
}
"do.trace.time"<-
function(expr, what)
{
	Tpdiff <- function(x)
	{
		y <- numeric()
		for(i in seq(1, length(x), by = 2))
			y[(i + 1)/2] <- x[i + 1] - x[i]
		round(y, 2)
	}
	if(!missing(what)) {
		trace(what, trace.time, exit = trace.time)
		on.exit(untrace(what))
	}
	time <- proc.time()[1]
	Time <<- numeric()
	Twhat <<- character()
	N <<- 0
	if(is.name(substitute(expr))) {
		update(expr)
		expr <- expr$call
	}
	else expr
	if(length(Time))
		time <- sapply(split(Time - time, Twhat), Tpdiff)
	attr(time, "expr") <- substitute(expr)
	time
}
"g.min"<-
function(x, A, B)
{
	z <- A %*% x - B
	structure(crossprod(z), gradient = crossprod(A, z) * 2)
}
"is.all.equal"<-
function(target, current, ...)
{
	TTT <- all.equal(target, current, ...)
	if(is.logical(TTT))
		TTT
	else FALSE
}
"lgt.inv"<-
function(x)
{
	y <- x
	pos <- x > 0
	neg <- !pos
	y[pos] <- 1/(1 + exp( - x[pos]))
	y[neg] <- exp(x[neg])
	y[neg] <- y[neg]/(1 + y[neg])
	y
}
"logexp"<-
function(x)
{
#  log( 1+exp(x) )
	y <- x
	pos <- x > 0
	y[pos] <- x[pos] + logexp.neg( - x[pos])
	y[!pos] <- logexp.neg(x[!pos])
	y
}
"logexp.neg"<-
function(x)
{
#log(1 + exp (x)) for negative arguments
	y <- exp(x)
	small <- y < 4 * .Machine$double.eps
	y[small] <- y[small] * (1 - 0.5 * y[small])
	y[!small] <- log(1 + y[!small])
	y
}
"lprob"<-
function(lp)
log(1 + exp(lp)) - lp
"lprob2"<-
function(lp, X)
{
	elp <- exp(lp)
	z <- 1 + elp
	value <- log(z) - lp
	attr(value, "gradient") <-  - X/z
	value
}
"mprompt"<-
function(name, filename = paste(name, ".d", sep = ""), detail = F)
{
	separate.dots <- function(name)
	{
		if(!exists(name, mode = "function"))
			stop(paste(name, "does not exist as a function"))
		string <- unix("sed 's/\\./\\n/g'", name)
		switch(length(string),
			stop(paste(name, "is not a method")),
			if(exists(string[1], mode = "function")) list(generic
					 = string[1], class = string[2]),
			for(i in 1:(length(string) - 1)) {
				j <- seq(i)
				TT <- paste(string[j], collapse = ".")
				if(exists(TT, mode = "function")) {
					return(generic = TT, class = paste(
						string[ - j], collapse = "."))
					
				}
			}
			)
		list(generic=string[1],class=paste(string[-1],collapse="."))
	}
	name <- substitute(name)
	if(is.language(name) && !is.name(name))
		name <- eval(name)
	name <- as.character(name)
	if(is.logical(filename) && filename)
		filename <- paste(".Data/.Help/", name, sep = "")
	TT <- separate.dots(name)
	cl <- TT$class
	gn <- TT$generic
	default <- paste(gn, ".default", sep = "")
	qcl <- paste("`", cl, "'", sep = "")
	gnp <- switch(gn,
		Math = "exp",
		Ops = "+",
		Summary = "sum",
		As = "as.vector",
		gn)
	generic <- get(gnp, mode = "function")
	is.group <- gnp!=gn
	fgn <- if(is.group) paste("the `", gn, "' group of functions",sep="")
	else paste("the function `", gn, "()'", sep = "")
	if(exists(default, mode = "function")) {
		default <- get(default, mode = "function")
		fdflt <- paste("`", gn, "' or `", gn, ".default'", sep = "")
	}
	else {
		default <- generic
		fdflt <- paste("`", gn, "'", sep = "")
	}
	fn <- get(name, mode = "function")
	n <- length(fn) - 1
	s <- 1:n
	args <- fn[s]
	arg.names <- names(args)
	ng <- length(generic) - 1
	generic.names <- names(generic[1:ng])
	add.names <- arg.names[!match(arg.names, generic.names, F)]
	common.names <- arg.names[match(names(default)[ - length(default)],
		arg.names, 0)]
	if(length(common.names)==0) common.names <- arg.names[1]
	else common.names <- common.names[common.names != "..."]
	file <- c(".BG", paste(".FN", name), ".TL", if(detail) 
			"~method to do ???" else paste("Use", fgn, "on a",
			qcl, "object"), ".CS")
	call <- paste(name, "(", sep = "")
	for(i in s) {
		if(!detail | mode(args[[i]]) == "missing")
			call <- paste(call, arg.names[i], sep = "")
		else call <- paste(call, arg.names[i], "=", deparse(args[[i]]),
				sep = "")
		if(i != n)
			call <- paste(call, ", ", sep = "")
	}
	file <- c(file, paste(call, ")", sep = ""))
	message <- c(".PP", paste("This is a method for ", fgn,
		" for objects inheriting from class ", qcl, ".", sep = ""),
		paste("See", fdflt, "for the general behavior of this function"
		), paste("and for the interpretation of ", paste("`", 
		common.names, "'", collapse = ", ", sep = ""), ".", sep = ""))
	file <- c(file, message)
	if(detail & length(add.names)) {
		file <- c(file, ".PP", 
			"The following arguments are particular to this method:"
			)
		for(i in add.names)
			file <- c(file, paste(".AG", i), paste("~Describe",
				i, "here"))
		file <- c(file, ".RT", "~Describe the value returned", ".EX",
			"# The function is currently defined as", deparse(
			fn), ".KW ~keyword", ".WR")
	}
	else file <- c(file, ".WR")
	cat(file, file = filename, sep = "\n")
}
"sprompt"<-
function(support, master, filename = paste(support, ".d", sep = ""), detail = F
	)
{
	support <- substitute(support)
	if(is.language(support) && !is.name(support))
		support <- eval(support)
	support <- as.character(support)
	master <- substitute(master)
	if(is.language(master) && !is.name(master))
		master <- eval(master)
	master <- as.character(master)
	master <- paste(master, "()", sep = "")
	if(is.logical(filename) && filename)
		filename <- paste(".Data/.Help/", support, sep = "")
	title <- if(detail) "~Support function to do ???" else if(length(master
		) == 1)
		paste("Support for Function", master)
	else paste("Support for Functions", paste(master, collapse = ", "))
	file <- c(".BG", paste(".FN", support), ".TL", title, ".CS")
	get.call <- function(args, name, detail)
	{
		arg.names <- names(args)
		call <- paste(name, "(", sep = "")
		n <- length(args) - 1
		for(i in seq(length = n)) {
			if(!detail | mode(args[[i]]) == "missing")
				call <- paste(call, arg.names[i], sep = "")
				
			else call <- paste(call, arg.names[i], "=", deparse(
					args[[i]]), sep = "")
			if(i != n)
				call <- paste(call, ", ", sep = "")
		}
		paste(call, ")", sep = "")
	}
	for(i in support)
		file <- c(file, get.call(get(i, mode = "function"), i, detail))
		
	message <- c(".PP", paste("This is support for the", if(length(master) ==
		1) paste("function ", master, ".", sep = "") else paste(
			"functions ", paste(master[ - length(master)], collapse
			 = ", "), " and ", master[length(master)], ".", sep = 
			"")), 
		"It is not intended to be called directly by users.")
	file <- c(file, message)
	if(detail) {
		ut <- get(support, mode = "function")
		arg.names <- names(ut)[-1]
		for(i in arg.names)
			file <- c(file, paste(".AG", i), paste("~Describe",
				i, "here"))
		file <- c(file, ".RT", "~Describe the value returned", ".EX",
			"# The function is currently defined as", deparse(
			args), ".KW ~keyword", ".WR")
	}
	else file <- c(file, ".WR")
	cat(file, file = filename, sep = "\n")
}
"testfit"<-
function(object, tolerance = .Machine$double.eps^0.5, data = list())
{
	fit <- object$call
	if(is.null(fit))
		fit <- attr(object, "call")
	if(!is.language(fit))
		stop(
			"object must have the fitting expression as component \"call\""
			)	#nframe <- new.frame(data)
	nframe <- FALSE
	newfit <- eval(fit, nframe)
	all.equal(object, newfit, tolerance)
}
"trace.time"<-
function()
{
	if(!exists("N"))
		return()
	N <<- N + 1
	Time[N] <<- proc.time()[1]
	Twhat[N] <<- deparse(sys.call(sys.parent())[[1]])[1]
}
"zipmle"<-
function(xp, xl, y, coefp, coefl)
{
#the xp and xl matrices need a column of ones if there's an intercept
	ncoefp <- ncol(xp)
	ncoefl <- ncol(xl)
	ncoefpl <- ncoefp + ncoefl
	info <- array(NA, c(ncoefpl, ncoefpl))
	pos <- y > 0.5
	zero <- y < 0.5
	Sp <- xp %*% coefp
	mup <- exp(Sp)
	Sl <- xl %*% coefl
	negllik <- logexp(Sl) + lgamma(y + 1)
	negllik[pos] <- negllik[pos] - y[pos] * Sp[pos] + mup[pos]
	negllik[zero] <- negllik[zero] - Sl[zero] - logexp( - Sl[zero] - mup[
		zero])
	negllik <- sum(negllik)
	p0 <- lgt.inv(Sl)
	p1 <- p0 * 0
	p1[zero] <- lgt.inv(Sl[zero] + mup[zero])
	v0 <- p0 * (1 - p0)
	v1 <- p1 * (1 - p1)
	score.neg <- c(t(xp) %*% (mup * (1 - p1) - y), t(xl) %*% (p0 - p1))
	w1 <- mup * (1 - p1) * (1 - mup * p1)
	w2 <- ( - mup) * v1
	w3 <- v0 - v1
	info[1:ncoefp, 1:ncoefp] <- t(xp) %*% (as.vector(w1) * xp)
	info[(ncoefp + 1):(ncoefl + ncoefp), 1:ncoefp] <- t(xp) %*% (as.vector(
		w2) * xl)
	info[1:ncoefp, (ncoefp + 1):(ncoefpl)] <- t(info[(ncoefp + 1):(ncoefpl),
		1:ncoefp])
	info[(ncoefp + 1):(ncoefpl), (ncoefp + 1):(ncoefpl)] <- t(xl) %*% (
		as.vector(w3) * xl)
	attr(negllik, "gradient") <- score.neg
	attr(negllik, "hessian") <- info
	negllik
}