R/bootfuns.q

Defines functions find_type ts.return scramble tsboot make.ends ts.array lines.saddle.distn print.saddle.distn saddle.distn saddle print.simplex simplex1 simplex zero iden inv.logit logit cum3 k3.linear var.linear control tilt.boot smooth.f imp.prob imp.quantile imp.reg imp.moments const imp.weights exp.tilt glm.diag.plots glm.diag envelope linear.approx positive.jack usual.jack empinf.reg inf.jack empinf cens.resamp cens.weird cens.case cens.return censboot abc.ci bca.ci perc.ci stud.ci basic.ci norm.inter norm.ci print.bootci boot.ci cv.glm permutation.array ordinary.array jack.after.boot index.array importance.array.bal importance.array freq.array extra.array corr print.boot plot.boot boot.array boot.return normalize boot balanced.array anti.arr antithetic.array rperm isMatrix bsample sample0

Documented in abc.ci boot boot.array boot.ci boot.return censboot cens.return control corr cum3 cv.glm empinf envelope exp.tilt freq.array glm.diag glm.diag.plots imp.moments imp.prob imp.quantile imp.reg imp.weights inv.logit jack.after.boot k3.linear linear.approx lines.saddle.distn logit norm.ci plot.boot print.boot print.bootci print.saddle.distn print.simplex saddle saddle.distn simplex smooth.f tilt.boot tsboot ts.return var.linear

# part of R package boot
# copyright (C) 1997-2001 Angelo J. Canty
# corrections (C) 1997-2014 B. D. Ripley
#
# Unlimited distribution is permitted

# safe version of sample
# needs R >= 2.9.0
# only works if size is not specified in R >= 2.11.0, but it always is in boot
sample0 <- function(x, ...) x[sample.int(length(x), ...)]
bsample <- function(x, ...) x[sample.int(length(x), replace = TRUE, ...)]

isMatrix <- function(x) length(dim(x)) == 2L

## random permutation of x.
rperm <- function(x) sample0(x, length(x))



antithetic.array <- function(n, R, L, strata)
#
#  Create an array of indices by antithetic resampling using the
#  empirical influence values in L.  This function just calls anti.arr
#  to do the sampling within strata.
#
{
    inds <- as.integer(names(table(strata)))
    out <- matrix(0L, R, n)
    for (s in inds) {
	gp <- seq_len(n)[strata == s]
        out[, gp] <- anti.arr(length(gp), R, L[gp], gp)
    }
    out
}

anti.arr <- function(n, R, L, inds=seq_len(n))
{
#  R x n array of bootstrap indices, generated antithetically
#  according to the empirical influence values in L.
    unique.rank <- function(x) {
# Assign unique ranks to a numeric vector
        ranks <- rank(x)
        if (any(duplicated(ranks))) {
            inds <- seq_along(x)
            uniq <- sort(unique(ranks))
            tab <- table(ranks)
            for (i in seq_along(uniq))
                if (tab[i] > 1L) {
                    gp <- inds[ranks == uniq[i]]
                    ranks[gp] <- rperm(inds[sort(ranks) == uniq[i]])
                }
        }
        ranks
    }
    R1 <- floor(R/2)
    mat1 <- matrix(bsample(inds, R1*n), R1, n)
    ranks <- unique.rank(L)
    rev <- inds
    for (i in seq_len(n)) rev[i] <- inds[ranks == (n+1-ranks[i])]
    mat1 <- rbind(mat1, matrix(rev[mat1], R1, n))
    if (R != 2*R1) mat1 <- rbind(mat1, bsample(inds, n))
    mat1
}




balanced.array <- function(n, R, strata)
{
#
# R x n array of bootstrap indices, sampled hypergeometrically
# within strata.
#
    output <- matrix(rep(seq_len(n), R), n, R)
    inds <- as.integer(names(table(strata)))
    for(is in inds) {
        group <- seq_len(n)[strata == is]
        if(length(group) > 1L) {
            g <- matrix(rperm(output[group,  ]), length(group), R)
            output[group,  ] <- g
        }
    }
    t(output)
}

boot <- function(data, statistic, R, sim = "ordinary",
                 stype = c("i", "f", "w"),
                 strata  =  rep(1, n), L = NULL, m = 0, weights = NULL,
		 ran.gen = function(d, p) d, mle = NULL, simple = FALSE, ...,
                 parallel = c("no", "multicore", "snow"),
                 ncpus = getOption("boot.ncpus", 1L), cl = NULL)
{
#
# R replicates of bootstrap applied to  statistic(data)
# Possible sim values are: "ordinary", "balanced", "antithetic",
#                     "permutation", "parametric"
# Various auxilliary functions find the indices to be used for the
# bootstrap replicates and then this function loops over those replicates.
#
    call <- match.call()
    stype <- match.arg(stype)
    if (missing(parallel)) parallel <- getOption("boot.parallel", "no")
    parallel <- match.arg(parallel)
    have_mc <- have_snow <- FALSE
    if (parallel != "no" && ncpus > 1L) {
        if (parallel == "multicore") have_mc <- .Platform$OS.type != "windows"
        else if (parallel == "snow") have_snow <- TRUE
        if (!have_mc && !have_snow) ncpus <- 1L
        loadNamespace("parallel") # get this out of the way before recording seed
    }
    if (simple && (sim != "ordinary" || stype != "i" || sum(m))) {
        warning("'simple=TRUE' is only valid for 'sim=\"ordinary\", stype=\"i\", n=0', so ignored")
        simple <- FALSE
    }
    if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) runif(1)
    seed <- get(".Random.seed", envir = .GlobalEnv, inherits = FALSE)
    n <- NROW(data)
    if ((n == 0) || is.null(n))
        stop("no data in call to 'boot'")
    temp.str <- strata
    strata <- tapply(seq_len(n),as.numeric(strata))
    t0 <- if (sim != "parametric") {
	if ((sim == "antithetic") && is.null(L))
            L <- empinf(data = data, statistic = statistic,
                        stype = stype, strata = strata, ...)
        if (sim != "ordinary") m <- 0
        else if (any(m < 0)) stop("negative value of 'm' supplied")
        if ((length(m) != 1L) && (length(m) != length(table(strata))))
            stop("length of 'm' incompatible with 'strata'")
        if ((sim == "ordinary") || (sim == "balanced")) {
            if (isMatrix(weights) && (nrow(weights) != length(R)))
                stop("dimensions of 'R' and 'weights' do not match")}
        else weights <- NULL
        if (!is.null(weights))
            weights <- t(apply(matrix(weights, n, length(R), byrow = TRUE),
                               2L, normalize, strata))
        if (!simple) i <- index.array(n, R, sim, strata, m, L, weights)

        original <- if (stype == "f") rep(1, n)
        else if (stype == "w") {
            ns <- tabulate(strata)[strata]
            1/ns
        } else seq_len(n)

        t0 <- if (sum(m) > 0L) statistic(data, original, rep(1, sum(m)), ...)
        else statistic(data, original, ...)
        rm(original)
        t0
    } else # "parametric"
	statistic(data, ...)

    pred.i <- NULL
    fn <- if (sim == "parametric") {
        ## force promises, so values get sent by parallel
        ran.gen; data; mle
        function(r) {
            dd <- ran.gen(data, mle)
            statistic(dd, ...)
        }
    } else {
        if (!simple && ncol(i) > n) {
            pred.i <- as.matrix(i[ , (n+1L):ncol(i)])
            i <- i[, seq_len(n)]
        }
        if (stype %in% c("f", "w")) {
            f <- freq.array(i)
            rm(i)
            if (stype == "w") f <- f/ns
            if (sum(m) == 0L) function(r) statistic(data, f[r,  ], ...)
            else function(r) statistic(data, f[r, ], pred.i[r, ], ...)
        } else if (sum(m) > 0L)
            function(r) statistic(data, i[r, ], pred.i[r,], ...)
        else if (simple)
            function(r)
                statistic(data,
                          index.array(n, 1, sim, strata, m, L, weights), ...)
        else function(r) statistic(data, i[r, ], ...)
    }
    RR <- sum(R)
    res <- if (ncpus > 1L && (have_mc || have_snow)) {
        if (have_mc) {
            parallel::mclapply(seq_len(RR), fn, mc.cores = ncpus)
        } else if (have_snow) {
            list(...) # evaluate any promises
            if (is.null(cl)) {
                cl <- parallel::makePSOCKcluster(rep("localhost", ncpus))
                if(RNGkind()[1L] == "L'Ecuyer-CMRG")
                    parallel::clusterSetRNGStream(cl)
                res <- parallel::parLapply(cl, seq_len(RR), fn)
                parallel::stopCluster(cl)
                res
            } else parallel::parLapply(cl, seq_len(RR), fn)
        }
    } else lapply(seq_len(RR), fn)
    t.star <- matrix(, RR, length(t0))
    for(r in seq_len(RR)) t.star[r, ] <- res[[r]]

    if (is.null(weights)) weights <- 1/tabulate(strata)[strata]
    boot0 <- boot.return(sim, t0, t.star, temp.str, R, data, statistic,
                         stype, call,
                         seed, L, m, pred.i, weights, ran.gen, mle)
    attr(boot0, "boot_type") <- "boot"
    boot0
}

normalize <- function(wts, strata)
{
#
# Normalize a vector of weights to sum to 1 within each strata.
#
    n <- length(strata)
    out <- wts
    inds <- as.integer(names(table(strata)))
    for (is in inds) {
        gp <- seq_len(n)[strata == is]
        out[gp] <- wts[gp]/sum(wts[gp]) }
    out
}

boot.return <- function(sim, t0, t, strata, R, data, stat, stype, call,
			seed, L, m, pred.i, weights, ran.gen, mle)
#
# Return the results of a bootstrap in the form of an object of class
# "boot".
#
{
    out <- list(t0=t0, t=t, R=R, data=data, seed=seed,
                statistic=stat, sim=sim, call=call)
    if (sim == "parametric")
        out <- c(out, list(ran.gen=ran.gen, mle=mle))
    else if (sim == "antithetic")
        out <- c(out, list(stype=stype, strata=strata, L=L))
    else if (sim == "ordinary") {
        if (sum(m) > 0)
            out <- c(out, list(stype=stype, strata=strata,
                               weights=weights, pred.i=pred.i))
        else 	out <- c(out, list(stype=stype, strata=strata,
                                   weights=weights))
    } else if (sim == "balanced")
        out <- c(out, list(stype=stype, strata=strata,
                           weights=weights ))
    else
        out <- c(out, list(stype=stype, strata=strata))
    class(out) <- "boot"
    out
}

c.boot <- function (..., recursive = TRUE)
{
    args <- list(...)
    nm <- lapply(args, names)
    if (!all(sapply(nm, function(x) identical(x, nm[[1]]))))
        stop("arguments are not all the same type of \"boot\" object")
    res <- args[[1]]
    res$R <- sum(sapply(args, "[[", "R"))
    res$t <- do.call(rbind, lapply(args, "[[", "t"))
    res
}

boot.array <- function(boot.out, indices=FALSE) {
#
#  Return the frequency or index array for the bootstrap resamples
#  used in boot.out
#  This function recreates such arrays from the information in boot.out
#
    if (exists(".Random.seed", envir=.GlobalEnv, inherits = FALSE))
        temp <- get(".Random.seed", envir=.GlobalEnv, inherits = FALSE)
    else temp<- NULL
    assign(".Random.seed",  boot.out$seed, envir=.GlobalEnv)
    n <- NROW(boot.out$data)
    R <- boot.out$R
    sim <- boot.out$sim
    type <- find_type(boot.out)
    if (type == "tsboot") {
#  Recreate the array for an object created by tsboot, The default for
#  such objects is to return the index array unless index is specifically
#  passed as F
        if (missing(indices)) indices <- TRUE
        if (sim == "model")
            stop("index array not defined for model-based resampling")
        n.sim <- boot.out$n.sim
        i.a <- ts.array(n, n.sim, R, boot.out$l,
			sim, boot.out$endcorr)
        out <- matrix(NA,R,n.sim)
        for(r in seq_len(R)) {
            if (sim == "geom")
                ends <- cbind(i.a$starts[r,  ],
                              i.a$lengths[r,  ])
            else
                ends <- cbind(i.a$starts[r,], i.a$lengths)
            inds <- apply(ends, 1L, make.ends, n)
            if (is.list(inds))
                inds <- unlist(inds)[seq_len(n.sim)]
            out[r,] <- inds
        }
    }
    else if (type == "censboot") {
#  Recreate the array for an object created by censboot as long
#  as censboot was called with sim = "ordinary"
        if (sim == "ordinary") {
            strata <- tapply(seq_len(n), as.numeric(boot.out$strata))
            out <- cens.case(n,strata,R)
        }
        else	stop("boot.array not implemented for this object")
    }
    else {
#  Recreate the array for objects created by boot or tilt.boot
        if (sim == "parametric")
            stop("array cannot be found for parametric bootstrap")
        strata <- tapply(seq_len(n),as.numeric(boot.out$strata))
        if (find_type(boot.out) == "tilt.boot")
            weights <- boot.out$weights
        else {
            weights <- boot.out$call$weights
            if (!is.null(weights))
                weights <- boot.out$weights
        }
        out <- index.array(n, R, sim, strata, 0, boot.out$L, weights)
    }
    if (!indices) out <- freq.array(out)
    if (!is.null(temp)) assign(".Random.seed", temp, envir=.GlobalEnv)
    else rm(.Random.seed, pos=1)
    out
}

plot.boot <- function(x,index=1, t0=NULL, t=NULL, jack=FALSE,
	qdist="norm",nclass=NULL,df, ...) {
#
#  A plot method for bootstrap output objects.  It produces a histogram
#  of the bootstrap replicates and a QQ plot of them.  Optionally it can
#  also produce a jackknife-after-bootstrap plot.
#
    boot.out <- x
    t.o <- t
    if (is.null(t)) {
        t <- boot.out$t[,index]
        if (is.null(t0)) t0 <- boot.out$t0[index]
    }
    t <- t[is.finite(t)]
    if (const(t, min(1e-8,mean(t, na.rm=TRUE)/1e6))) {
        print(paste("All values of t* are equal to ", mean(t, na.rm=TRUE)))
        return(invisible(boot.out))
    }
    if (is.null(nclass)) nclass <- min(max(ceiling(length(t)/25),10),100)
    if (!is.null(t0)) {
#  Calculate the breakpoints for the histogram so that one of them is
#  exactly t0.
        rg <- range(t)
        if (t0<rg[1L]) rg[1L] <- t0
        else if (t0 >rg[2L]) rg[2L] <- t0
        rg <- rg+0.05*c(-1,1)*diff(rg)
        lc <- diff(rg)/(nclass-2)
        n1 <- ceiling((t0-rg[1L])/lc)
        n2 <- ceiling((rg[2L]-t0)/lc)
        bks <- t0+(-n1:n2)*lc
    }
    R <- boot.out$R
    if (qdist == "chisq") {
        qq <- qchisq((seq_len(R))/(R+1),df=df)
        qlab <- paste("Quantiles of Chi-squared(",df,")",sep="")
    }
    else {
	if (qdist!="norm")
            warning(gettextf("%s distribution not supported: using normal instead", sQuote(qdist)), domain = NA)
        qq <- qnorm((seq_len(R))/(R+1))
        qlab <-"Quantiles of Standard Normal"
    }
    if (jack) {
        layout(mat = matrix(c(1,2,3,3), 2L, 2L, byrow=TRUE))
        if (is.null(t0))
            hist(t,nclass=nclass,probability=TRUE,xlab="t*")
        else	hist(t,breaks=bks,probability=TRUE,xlab="t*")
        if (!is.null(t0)) abline(v=t0,lty=2)
        qqplot(qq,t,xlab=qlab,ylab="t*")
        if (qdist == "norm") abline(mean(t),sqrt(var(t)),lty=2)
        else abline(0,1,lty=2)
        jack.after.boot(boot.out,index=index,t=t.o, ...)
    }
    else {
        par(mfrow=c(1,2))
        if (is.null(t0))
            hist(t,nclass=nclass,probability=TRUE,xlab="t*")
        else	hist(t,breaks=bks,probability=TRUE,xlab="t*")
        if (!is.null(t0)) abline(v=t0,lty=2)
        qqplot(qq,t,xlab=qlab,ylab="t*")
        if (qdist == "norm") abline(mean(t),sqrt(var(t)),lty=2)
        else abline(0,1,lty=2)
    }
    par(mfrow=c(1,1))
    invisible(boot.out)
}

print.boot <- function(x, digits = getOption("digits"),
                          index = 1L:ncol(boot.out$t), ...)
{
#
# Print the output of a bootstrap
#
    boot.out <- x
    sim <- boot.out$sim
    cl <- boot.out$call
    t <- matrix(boot.out$t[, index], nrow = nrow(boot.out$t))
    allNA <- apply(t,2L,function(t) all(is.na(t)))
    ind1 <- index[allNA]
    index <- index[!allNA]
    t <- matrix(t[, !allNA], nrow = nrow(t))
    rn <- paste("t",index,"*",sep="")
    if (length(index) == 0L)
        op <- NULL
    else if (is.null(t0 <- boot.out$t0)) {
        if (is.null(boot.out$call$weights))
            op <- cbind(apply(t,2L,mean,na.rm=TRUE),
                        sqrt(apply(t,2L,function(t.st) var(t.st[!is.na(t.st)]))))
        else {
            op <- NULL
            for (i in index)
                op <- rbind(op, imp.moments(boot.out,index=i)$rat)
            op[,2L] <- sqrt(op[,2])
        }
        dimnames(op) <- list(rn,c("mean", "std. error"))
    }
    else {
        t0 <- boot.out$t0[index]
        if (is.null(boot.out$call$weights)) {
            op <- cbind(t0,apply(t,2L,mean,na.rm=TRUE)-t0,
                        sqrt(apply(t,2L,function(t.st) var(t.st[!is.na(t.st)]))))
            dimnames(op) <- list(rn, c("original"," bias  "," std. error"))
        }
        else {
            op <- NULL
            for (i in index)
                op <- rbind(op, imp.moments(boot.out,index=i)$rat)
            op <- cbind(t0,op[,1L]-t0,sqrt(op[,2L]),
                        apply(t,2L,mean,na.rm=TRUE))
            dimnames(op) <- list(rn,c("original", " bias  ",
                                      " std. error", " mean(t*)"))
        }
    }
    type <- find_type(boot.out)
    if (type == "boot") {
        if (sim == "parametric")
            cat("\nPARAMETRIC BOOTSTRAP\n\n")
        else if (sim == "antithetic") {
            if (is.null(cl$strata))
                cat("\nANTITHETIC BOOTSTRAP\n\n")
            else
                cat("\nSTRATIFIED ANTITHETIC BOOTSTRAP\n\n")
        }
        else if (sim == "permutation") {
            if (is.null(cl$strata))
                cat("\nDATA PERMUTATION\n\n")
            else
                cat("\nSTRATIFIED DATA PERMUTATION\n\n")
        }
        else if (sim == "balanced") {
            if (is.null(cl$strata) && is.null(cl$weights))
                cat("\nBALANCED BOOTSTRAP\n\n")
            else if (is.null(cl$strata))
                cat("\nBALANCED WEIGHTED BOOTSTRAP\n\n")
            else if (is.null(cl$weights))
                cat("\nSTRATIFIED BALANCED BOOTSTRAP\n\n")
            else
                cat("\nSTRATIFIED WEIGHTED BALANCED BOOTSTRAP\n\n")
        }
        else {
            if (is.null(cl$strata) && is.null(cl$weights))
                cat("\nORDINARY NONPARAMETRIC BOOTSTRAP\n\n")
            else if (is.null(cl$strata))
                cat("\nWEIGHTED BOOTSTRAP\n\n")
            else if (is.null(cl$weights))
                cat("\nSTRATIFIED BOOTSTRAP\n\n")
            else
                cat("\nSTRATIFIED WEIGHTED BOOTSTRAP\n\n")
        }
    }
    else if (type == "tilt.boot") {
        R <- boot.out$R
        th <- boot.out$theta
        if (sim == "balanced")
            cat("\nBALANCED TILTED BOOTSTRAP\n\n")
        else	cat("\nTILTED BOOTSTRAP\n\n")
        if ((R[1L] == 0) || is.null(cl$tilt) || eval(cl$tilt))
            cat("Exponential tilting used\n")
        else	cat("Frequency Smoothing used\n")
        i1 <- 1
        if (boot.out$R[1L]>0)
            cat(paste("First",R[1L],"replicates untilted,\n"))
        else {
            cat(paste("First ",R[2L]," replicates tilted to ",
                      signif(th[1L],4),",\n",sep=""))
            i1 <- 2
        }
        if (i1 <= length(th)) {
            for (j in i1:length(th))
                cat(paste("Next ",R[j+1L]," replicates tilted to ",
                          signif(th[j],4L),
                          ifelse(j!=length(th),",\n",".\n"),sep=""))
        }
        op <- op[, 1L:3L]
    }
    else if (type == "tsboot") {
        if (!is.null(cl$indices))
            cat("\nTIME SERIES BOOTSTRAP USING SUPPLIED INDICES\n\n")
        else if (sim == "model")
            cat("\nMODEL BASED BOOTSTRAP FOR TIME SERIES\n\n")
        else if (sim == "scramble") {
            cat("\nPHASE SCRAMBLED BOOTSTRAP FOR TIME SERIES\n\n")
            if (boot.out$norm)
                cat("Normal margins used.\n")
            else
                cat("Observed margins used.\n")
        }
        else if (sim == "geom") {
            if (is.null(cl$ran.gen))
                cat("\nSTATIONARY BOOTSTRAP FOR TIME SERIES\n\n")
            else
                cat(paste("\nPOST-BLACKENED STATIONARY",
                          "BOOTSTRAP FOR TIME SERIES\n\n"))
            cat(paste("Average Block Length of",boot.out$l,"\n"))
        }
        else {
            if (is.null(cl$ran.gen))
                    cat("\nBLOCK BOOTSTRAP FOR TIME SERIES\n\n")
            else
                cat(paste("\nPOST-BLACKENED BLOCK",
                          "BOOTSTRAP FOR TIME SERIES\n\n"))
            cat(paste("Fixed Block Length of",boot.out$l,"\n"))
        }
    }
    else if (type == "censboot") {
        cat("\n")
        if (sim == "weird") {
            if (!is.null(cl$strata)) cat("STRATIFIED ")
            cat("WEIRD BOOTSTRAP FOR CENSORED DATA\n\n")
        }
        else if ((sim == "ordinary") ||
                 ((sim == "model") && is.null(boot.out$cox))) {
            if (!is.null(cl$strata)) cat("STRATIFIED ")
            cat("CASE RESAMPLING BOOTSTRAP FOR CENSORED DATA\n\n")
        }
        else if (sim == "model") {
            if (!is.null(cl$strata)) cat("STRATIFIED ")
            cat("MODEL BASED BOOTSTRAP FOR COX REGRESSION MODEL\n\n")
        }
        else if (sim == "cond") {
            if (!is.null(cl$strata)) cat("STRATIFIED ")
            cat("CONDITIONAL BOOTSTRAP ")
            if (is.null(boot.out$cox))
                cat("FOR CENSORED DATA\n\n")
            else
                cat("FOR COX REGRESSION MODEL\n\n")
        }
    } else warning('unknown type of "boot" object')
    cat("\nCall:\n")
    dput(cl, control=NULL)
    cat("\n\nBootstrap Statistics :\n")
    if (!is.null(op)) print(op,digits=digits)
    if (length(ind1) > 0L)
        for (j in ind1)
            cat(paste("WARNING: All values of t", j, "* are NA\n", sep=""))
    invisible(boot.out)
}




corr <- function(d, w=rep(1,nrow(d))/nrow(d))
{
#  The correlation coefficient in weighted form.
    s <- sum(w)
    m1 <- sum(d[, 1L] * w)/s
    m2 <- sum(d[, 2L] * w)/s
    (sum(d[, 1L] * d[, 2L] * w)/s - m1 * m2)/sqrt((sum(d[, 1L]^2 * w)/s - m1^2) * (sum(d[, 2L]^2 * w)/s - m2^2))
}


extra.array <- function(n, R, m, strata=rep(1,n))
{
#
# Extra indices for predictions.  Can only be used with
# types "ordinary" and "stratified".  For type "ordinary"
# m is a positive integer.  For type "stratified" m can
# be a positive integer or a vector of the same length as
# strata.
#
    if (length(m) == 1L)
        output <- matrix(sample.int(n, m*R, replace=TRUE), R, m)
    else {
        inds <- as.integer(names(table(strata)))
        output <- matrix(NA, R, sum(m))
        st <- 0
        for (i in inds) {
            if (m[i] > 0) {
                gp <- seq_len(n)[strata == i]
                inds1 <- (st+1):(st+m[i])
                output[,inds1] <- matrix(bsample(gp, R*m[i]), R, m[i])
                st <- st+m[i]
            }
        }
    }
    output
}

freq.array <- function(i.array)
{
#
# converts R x n array of bootstrap indices into
# R X n array of bootstrap frequencies
#
    result <- NULL
    n <- ncol(i.array)
    result <- t(apply(i.array, 1, tabulate, n))
    result
}



importance.array <- function(n, R, weights, strata){
#
#  Function to do importance resampling  within strata based
#  on the weights supplied.  If weights is a matrix with n columns
#  R must be a vector of length nrow(weights) otherwise weights
#  must be a vector of length n and R must be a scalar.
#
    imp.arr <- function(n, R, wts, inds=seq_len(n))
        matrix(bsample(inds, n*R, prob=wts), R, n)
    output <- NULL
    if (!isMatrix(weights))
        weights <- matrix(weights, nrow=1)
    inds <- as.integer(names(table(strata)))
    for (ir in seq_along(R)) {
        out <- matrix(rep(seq_len(n), R[ir]), R[ir], n, byrow=TRUE)
        for (is in inds) {
            gp <- seq_len(n)[strata == is]
            out[, gp] <- imp.arr(length(gp), R[ir],
                                 weights[ir,gp], gp)
        }
        output <- rbind(output, out)
    }
    output
}

importance.array.bal <- function(n, R, weights, strata) {
#
#  Function to do balanced importance resampling within strata
#  based on the supplied weights.  Balancing is achieved in such
#  a way that each index appears in the array approximately in
#  proportion to its weight.
#
    imp.arr.bal <- function(n, R, wts, inds=seq_len(n)) {
        if (sum (wts) != 1) wts <- wts / sum(wts)
        nRw1 <- floor(n*R*wts)
        nRw2 <- n*R*wts - nRw1
        output <- rep(inds, nRw1)
        if (any (nRw2 != 0))
            output <- c(output,
                        sample0(inds, round(sum(nRw2)), prob=nRw2))
        matrix(rperm(output), R, n)
    }
    output <- NULL
    if (!isMatrix(weights))
        weights <- matrix(weights, nrow = 1L)
    inds <- as.integer(names(table(strata)))
    for (ir in seq_along(R)) {
        out <- matrix(rep(seq_len(n), R[ir]), R[ir], n, byrow=TRUE)
        for (is in inds) {
            gp <- seq_len(n)[strata == is]
            out[,gp] <- imp.arr.bal(length(gp), R[ir], weights[ir,gp], gp)
        }
        output <- rbind(output, out)
    }
    output
}



index.array <- function(n, R, sim, strata=rep(1,n), m=0, L=NULL, weights=NULL)
{
#
#  Driver function for generating a bootstrap index array.  This function
#  simply determines the type of sampling required and calls the appropriate
#  function.
#
    indices <- NULL
    if (is.null (weights)) {
        if (sim == "ordinary") {
            indices <- ordinary.array(n, R, strata)
            if (sum(m) > 0)
                indices <- cbind(indices, extra.array(n, R, m, strata))
        }
    else if (sim == "balanced")
        indices <- balanced.array(n, R, strata)
    else if (sim == "antithetic")
        indices <- antithetic.array(n, R, L, strata)
    else if (sim == "permutation")
        indices <- permutation.array(n, R, strata)
    } else {
        if (sim == "ordinary")
            indices <- importance.array(n, R, weights, strata)
        else if (sim == "balanced")
            indices <- importance.array.bal(n, R, weights, strata)
    }
    indices
}

jack.after.boot <- function(boot.out, index=1, t=NULL, L=NULL,
	useJ=TRUE, stinf = TRUE, alpha=NULL, main = "", ylab=NULL, ...)
{
# jackknife after bootstrap plot
    t.o <- t
    if (is.null(t)) {
        if (length(index) > 1L) {
            index <- index[1L]
            warning("only first element of 'index' used")
        }
        t <- boot.out$t[, index]
    }
    fins <- seq_along(t)[is.finite(t)]
    t <- t[fins]
    if (is.null(alpha)) {
        alpha <- c(0.05, 0.1, 0.16, 0.5, 0.84, 0.9, 0.95)
        if (is.null(ylab))
            ylab <- "5, 10, 16, 50, 84, 90, 95 %-iles of (T*-t)"
    }
    if (is.null(ylab)) ylab <- "Percentiles of (T*-t)"
    data <- boot.out$data
    n <- NROW(data)
    f <- boot.array(boot.out)[fins, , drop=TRUE]
    percentiles <- matrix(data = NA, length(alpha), n)
    J <- numeric(n)
    for(j in seq_len(n)) {
# Find the quantiles of the bootstrap distribution on omitting each point.
        values <- t[f[, j] == 0]
        J[j] <- mean(values)
        percentiles[, j] <- quantile(values, alpha) - J[j]
    }
# Now find the jackknife values to be plotted, and standardize them,
# if required.
    if (!useJ) {
        if (is.null(L))
            J <- empinf(boot.out, index=index, t=t.o, ...)
        else 	J <- L
    }
    else	J <- (n - 1) * (mean(J) - J)
    xtext <- "jackknife value"
    if (!useJ) {
        if (!is.null(L) || (is.null(t.o) && (boot.out$stype == "w")))
            xtext <- paste("infinitesimal", xtext)
        else	xtext <- paste("regression", xtext)
    }
    if (stinf) {
        J <- J/sqrt(var(J))
        xtext <- paste("standardized", xtext)
    }
    top <- max(percentiles)
    bot <- min(percentiles)
    ylts <- c(bot - 0.35 * (top - bot), top + 0.1 * (top - bot))
    percentiles <- percentiles[,order(J)]#
# Plot the overall quantiles and the delete-1 quantiles against the
# jackknife values.
    plot(sort(J), percentiles[1,  ], ylim = ylts, type = "n", xlab = xtext,
         ylab = ylab, main=main)
    for(j in seq_along(alpha))
        lines(sort(J), percentiles[j,  ], type = "b", pch = "*")
    percentiles <- quantile(t, alpha) - mean(t)
    for(j in seq_along(alpha))
        abline(h=percentiles[j], lty=2)
# Now print the observation numbers below the plotted lines.  They are printed
# in five rows so that all numbers can be read easily.
    text(sort(J), rep(c(bot - 0.08 * (top - bot), NA, NA, NA, NA), n, n),
         order(J), cex = 0.5)
    text(sort(J), rep(c(NA, bot - 0.14 * (top - bot), NA, NA, NA), n, n),
         order(J), cex = 0.5)
    text(sort(J), rep(c(NA, NA, bot - 0.2 * (top - bot), NA, NA), n, n),
         order(J), cex = 0.5)
    text(sort(J), rep(c(NA, NA, NA, bot - 0.26 * (top - bot), NA), n, n),
         order(J), cex = 0.5)
    text(sort(J), rep(c(NA, NA, NA, NA, bot - 0.32 * (top - bot)), n, n),
         order(J), cex = 0.5)
    invisible()
}


ordinary.array <- function(n, R, strata)
{
#
# R x n array of bootstrap indices, resampled within strata.
# This is the function which generates a regular bootstrap array
# using equal weights within each stratum.
#
    inds <- as.integer(names(table(strata)))
    if (length(inds) == 1L) {
        output <- sample.int(n, n*R, replace=TRUE)
        dim(output) <- c(R, n)
    } else {
        output <- matrix(as.integer(0L), R, n)
        for(is in inds) {
            gp <- seq_len(n)[strata == is]
            output[, gp] <- if (length(gp) == 1) rep(gp, R) else bsample(gp, R*length(gp))
        }
    }
    output
}

permutation.array <- function(n, R, strata)
{
#
# R x n array of bootstrap indices, permuted within strata.
# This is similar to ordinary array except that resampling is
# done without replacement in each row.
#
    output <- matrix(rep(seq_len(n), R), n, R)
    inds <- as.integer(names(table(strata)))
    for(is in inds) {
        group <- seq_len(n)[strata == is]
        if (length(group) > 1L) {
            g <- apply(output[group,  ], 2L, rperm)
            output[group,  ] <- g
        }
    }
    t(output)
}


cv.glm <- function(data, glmfit, cost=function(y,yhat) mean((y-yhat)^2),
		   K=n)
{
# cross-validation estimate of error for glm prediction with K groups.
# cost is a function of two arguments: the observed values and the
# the predicted values.
    call <- match.call()
    if (!exists(".Random.seed", envir=.GlobalEnv, inherits = FALSE)) runif(1)
    seed <- get(".Random.seed", envir=.GlobalEnv, inherits = FALSE)
    n <- nrow(data)
    if ((K > n) || (K <= 1))
        stop("'K' outside allowable range")
    K.o <- K
    K <- round(K)
    kvals <- unique(round(n/(1L:floor(n/2))))
    temp <- abs(kvals-K)
    if (!any(temp == 0))
        K <- kvals[temp == min(temp)][1L]
    if (K!=K.o) warning(gettextf("'K' has been set to %f", K), domain = NA)
    f <- ceiling(n/K)
    s <- sample0(rep(1L:K, f), n)
    n.s <- table(s)
#    glm.f <- formula(glmfit)
    glm.y <- glmfit$y
    cost.0 <- cost(glm.y, fitted(glmfit))
    ms <- max(s)
    CV <- 0
    Call <- glmfit$call
    for(i in seq_len(ms)) {
        j.out <- seq_len(n)[(s == i)]
        j.in <- seq_len(n)[(s != i)]
        ## we want data from here but formula from the parent.
        Call$data <- data[j.in, , drop=FALSE]
        d.glm <- eval.parent(Call)
        p.alpha <- n.s[i]/n
        cost.i <- cost(glm.y[j.out],
                       predict(d.glm, data[j.out, , drop=FALSE],
                               type = "response"))
        CV <- CV + p.alpha * cost.i
        cost.0 <- cost.0 - p.alpha *
            cost(glm.y, predict(d.glm, data, type = "response"))
    }
    list(call = call, K = K,
         delta = as.numeric(c(CV, CV + cost.0)),  # drop any names
         seed = seed)
}


boot.ci <- function(boot.out,conf = 0.95,type = "all",
		    index = 1L:min(2L, length(boot.out$t0)),
		    var.t0 = NULL ,var.t = NULL, t0 = NULL, t = NULL,
		    L = NULL, h = function(t) t,
                    hdot = function(t) rep(1, length(t)),
		    hinv = function(t) t, ...)
#
#  Main function to calculate bootstrap confidence intervals.
#  This function calls a number of auxilliary functions to do
#  the actual calculations depending on the type of interval(s)
#  requested.
#
{
    call <- match.call()
#  Get and transform the statistic values and their variances,
    if ((is.null(t) && !is.null(t0)) ||
        (!is.null(t) && is.null(t0)))
        stop("'t' and 't0' must be supplied together")
    t.o <- t; t0.o <- t0
#    vt.o <- var.t
    vt0.o <- var.t0
    if (is.null(t)) {
        if (length(index) == 1L) {
            t0 <- boot.out$t0[index]
            t <- boot.out$t[,index]
        }
        else if (ncol(boot.out$t)<max(index)) {
            warning("index out of bounds; minimum index only used.")
            index <- min(index)
            t0 <- boot.out$t0[index]
            t <- boot.out$t[,index]
        }
        else {
            t0 <- boot.out$t0[index[1L]]
            t <- boot.out$t[,index[1L]]
            if (is.null(var.t0)) var.t0 <- boot.out$t0[index[2L]]
            if (is.null(var.t)) var.t <- boot.out$t[,index[2L]]
        }
    }
    if (const(t, min(1e-8, mean(t, na.rm=TRUE)/1e6))) {
        print(paste("All values of t are equal to ", mean(t, na.rm=TRUE),
                    "\n Cannot calculate confidence intervals"))
        return(NULL)
    }
    if (length(t) != boot.out$R)
        stop(gettextf("'t' must of length %d", boot.out$R), domain = NA)
    if (is.null(var.t))
        fins <- seq_along(t)[is.finite(t)]
    else {
        fins <- seq_along(t)[is.finite(t) & is.finite(var.t)]
        var.t <- var.t[fins]
    }
    t <- t[fins]
    R <- length(t)
    if (!is.null(var.t0)) var.t0 <- var.t0*hdot(t0)^2
    if (!is.null(var.t))  var.t <- var.t*hdot(t)^2
    t0 <- h(t0); t <- h(t)
    if (missing(L)) L <- boot.out$L
    output <- list(R = R, t0 = hinv(t0), call = call)
    #  Now find the actual intervals using the methods listed in type
    if (any(type == "all" | type == "norm"))
        output <- c(output,
                    list(normal = norm.ci(boot.out, conf,
                         index[1L], var.t0=vt0.o, t0=t0.o, t=t.o,
                         L=L, h=h, hdot=hdot, hinv=hinv)))
    if (any(type == "all" | type == "basic"))
        output <- c(output, list(basic=basic.ci(t0,t,conf,hinv=hinv)))
    if (any(type == "all" | type == "stud")) {
        if (length(index)==1L)
            warning("bootstrap variances needed for studentized intervals")
        else
            output <- c(output, list(student=stud.ci(c(t0,var.t0),
                                     cbind(t ,var.t), conf, hinv=hinv)))
    }
    if (any(type == "all" | type == "perc"))
        output <- c(output, list(percent=perc.ci(t,conf,hinv=hinv)))
    if (any(type == "all" | type == "bca")) {
        if (find_type(boot.out) == "tsboot")
            warning("BCa intervals not defined for time series bootstraps")
        else
            output <- c(output, list(bca=bca.ci(boot.out,conf,
                                     index[1L],L=L,t=t.o, t0=t0.o,
                                     h=h,hdot=hdot, hinv=hinv, ...)))
    }
    class(output) <- "bootci"
    output
}

print.bootci <- function(x, hinv = NULL, ...) {
#
#  Print the output from boot.ci
#
    ci.out <- x
    cl <- ci.out$call
    ntypes <- length(ci.out)-3L
    nints <- nrow(ci.out[[4L]])
    t0 <- ci.out$t0
    if (!is.null(hinv)) t0 <- hinv(t0)  #
#  Find the number of decimal places which should be used
    digs <- ceiling(log10(abs(t0)))
    if (digs <= 0) digs <- 4
    else if (digs >= 4) digs <- 0
    else digs <- 4-digs
    intlabs <- NULL
    basrg <- strg <- perg <- bcarg <- NULL
    if (!is.null(ci.out$normal))
        intlabs <- c(intlabs,"     Normal        ")
    if (!is.null(ci.out$basic)) {
        intlabs <- c(intlabs,"     Basic         ")
        basrg <- range(ci.out$basic[,2:3]) }
    if (!is.null(ci.out$student)) {
        intlabs <- c(intlabs,"   Studentized     ")
        strg <- range(ci.out$student[,2:3]) }
    if (!is.null(ci.out$percent)) {
        intlabs <- c(intlabs,"    Percentile     ")
        perg <- range(ci.out$percent[,2:3]) }
    if (!is.null(ci.out$bca)) {
        intlabs <- c(intlabs,"      BCa          ")
        bcarg <- range(ci.out$bca[,2:3]) }
    level <- 100*ci.out[[4L]][, 1L]
    if (ntypes == 4L) n1 <- n2 <- 2L
    else if (ntypes == 5L) {n1 <- 3L; n2 <- 2L}
    else {n1 <- ntypes; n2 <- 0L}
    ints1 <- matrix(NA,nints,2L*n1+1L)
    ints1[,1L] <- level
    n0 <- 4L
#  Re-organize the intervals and coerce them into character data
    for (i in n0:(n0+n1-1)) {
        j <- c(2L*i-6L,2L*i-5L)
        nc <- ncol(ci.out[[i]])
        nc <- c(nc-1L,nc)
        if (is.null(hinv))
            ints1[,j] <- ci.out[[i]][,nc]
        else	ints1[,j] <- hinv(ci.out[[i]][,nc])
    }
    n0 <- 4L+n1
    ints1 <- format(round(ints1,digs))
    ints1[,1L] <- paste("\n",level,"%  ",sep="")
    ints1[,2*(1L:n1)] <- paste("(",ints1[,2*(1L:n1)],",",sep="")
    ints1[,2*(1L:n1)+1L] <- paste(ints1[,2*(1L:n1)+1L],")  ")
    if (n2 > 0) {
        ints2 <- matrix(NA,nints,2L*n2+1L)
        ints2[,1L] <- level
        j <- c(2L,3L)
        for (i in n0:(n0+n2-1L)) {
            if (is.null(hinv))
                ints2[,j] <- ci.out[[i]][,c(4L,5L)]
            else	ints2[,j] <- hinv(ci.out[[i]][,c(4L,5L)])
            j <- j+2L
        }
        ints2 <- format(round(ints2,digs))
        ints2[,1L] <- paste("\n",level,"%  ",sep="")
        ints2[,2*(1L:n2)] <- paste("(",ints2[,2*(1L:n2)],",",sep="")
        ints2[,2*(1L:n2)+1L] <- paste(ints2[,2*(1L:n2)+1L],")  ")
    }
    R <- ci.out$R                       #
#  Print the intervals
    cat("BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS\n")
    cat(paste("Based on",R,"bootstrap replicates\n\n"))
    cat("CALL : \n")
    dput(cl, control=NULL)
    cat("\nIntervals : ")
    cat("\nLevel",intlabs[1L:n1])
    cat(t(ints1))
    if (n2 > 0) {
        cat("\n\nLevel",intlabs[(n1+1):(n1+n2)])
        cat(t(ints2))
    }
    if (!is.null(cl$h)) {
        if (is.null(cl$hinv) && is.null(hinv))
            cat("\nCalculations and Intervals on ",
                "Transformed Scale\n")
        else	cat("\nCalculations on Transformed Scale;",
                    " Intervals on Original Scale\n")
    }
    else if (is.null(cl$hinv) && is.null(hinv))
        cat("\nCalculations and Intervals on Original Scale\n")
    else 	cat("\nCalculations on Original Scale",
                    " but Intervals Transformed\n")#
#  Print any warnings about extreme values.
    if (!is.null(basrg)) {
        if ((basrg[1L] <= 1) || (basrg[2L] >= R))
            cat("Warning : Basic Intervals used Extreme Quantiles\n")
        if ((basrg[1L] <= 10) || (basrg[2L] >= R-9))
            cat("Some basic intervals may be unstable\n")
    }
    if (!is.null(strg)) {
        if ((strg[1L] <= 1) || (strg[2L] >= R))
            cat("Warning : Studentized Intervals used Extreme Quantiles\n")
        if ((strg[1L] <= 10) || (strg[2L] >= R-9))
            cat("Some studentized intervals may be unstable\n")
    }
    if (!is.null(perg)) {
        if ((perg[1L] <= 1) || (perg[2L] >= R))
            cat("Warning : Percentile Intervals used Extreme Quantiles\n")
        if ((perg[1L] <= 10) || (perg[2L] >= R-9))
            cat("Some percentile intervals may be unstable\n")
    }
    if (!is.null(bcarg)) {
        if ((bcarg[1L] <= 1) || (bcarg[2L] >= R))
            cat("Warning : BCa Intervals used Extreme Quantiles\n")
        if ((bcarg[1L] <= 10) || (bcarg[2L] >= R-9))
            cat("Some BCa intervals may be unstable\n")
    }
    invisible(ci.out)
}

norm.ci <-
    function(boot.out = NULL,conf = 0.95,index = 1,var.t0 = NULL, t0 = NULL,
             t = NULL, L = NULL, h = function(t) t, hdot = function(t) 1,
             hinv = function(t) t)
#
#  Normal approximation method for confidence intervals.  This can be
#  used with or without a bootstrap object.  If a bootstrap object is
#  given then the intervals are bias corrected and the bootstrap variance
#  estimate can be used if none is supplied.
#
{
    if (is.null(t0))  {
        if (!is.null(boot.out)) t0 <-boot.out$t0[index]
        else stop("bootstrap output object or 't0' required")
    }
    if (!is.null(boot.out) && is.null(t))
        t <- boot.out$t[,index]
    if (!is.null(t)) {
        fins <- seq_along(t)[is.finite(t)]
        t <- h(t[fins])
    }
    if (is.null(var.t0)) {
        if (is.null(t)) {
            if (is.null(L))
                stop("unable to calculate 'var.t0'")
            else	var.t0 <- sum((hdot(t0)*L/length(L))^2)
        }
        else	var.t0 <- var(t)
    }
    else	var.t0 <- var.t0*hdot(t0)^2
    t0 <- h(t0)
    if (!is.null(t))
        bias <- mean(t)-t0
    else	bias <- 0
    merr <- sqrt(var.t0)*qnorm((1+conf)/2)
    out <- cbind(conf,hinv(t0-bias-merr),hinv(t0-bias+merr))
    out
}

norm.inter <- function(t,alpha)
#
#  Interpolation on the normal quantile scale.  For a non-integer
#  order statistic this function interpolates between the surrounding
#  order statistics using the normal quantile scale.  See equation
#  5.8 of Davison and Hinkley (1997)
#
{
    t <- t[is.finite(t)]
    R <- length(t)
    rk <- (R+1)*alpha
    if (!all(rk>1 & rk<R))
        warning("extreme order statistics used as endpoints")
    k <- trunc(rk)
    inds <- seq_along(k)
    out <- inds
    kvs <- k[k>0 & k<R]
    tstar <- sort(t, partial = sort(union(c(1, R), c(kvs, kvs+1))))
    ints <- (k == rk)
    if (any(ints)) out[inds[ints]] <- tstar[k[inds[ints]]]
    out[k == 0] <- tstar[1L]
    out[k == R] <- tstar[R]
    not <- function(v) xor(rep(TRUE,length(v)),v)
    temp <- inds[not(ints) & k != 0 & k != R]
    temp1 <- qnorm(alpha[temp])
    temp2 <- qnorm(k[temp]/(R+1))
    temp3 <- qnorm((k[temp]+1)/(R+1))
    tk <- tstar[k[temp]]
    tk1 <- tstar[k[temp]+1L]
    out[temp] <- tk + (temp1-temp2)/(temp3-temp2)*(tk1 - tk)
    cbind(round(rk, 2), out)
}

basic.ci <- function(t0, t, conf = 0.95, hinv = function(t) t)
#
#  Basic bootstrap confidence method
#
{
    qq <- norm.inter(t,(1+c(conf,-conf))/2)
    cbind(conf, matrix(qq[,1L],ncol=2L), matrix(hinv(2*t0-qq[,2L]),ncol=2L))
}

stud.ci <- function(tv0, tv, conf = 0.95, hinv=function(t) t)
#
#  Studentized version of the basic bootstrap confidence interval
#
{
    if ((length(tv0) < 2) || (ncol(tv) < 2)) {
        warning("variance required for studentized intervals")
        NA
    } else {
        z <- (tv[,1L]-tv0[1L])/sqrt(tv[,2L])
        qq <- norm.inter(z, (1+c(conf,-conf))/2)
        cbind(conf, matrix(qq[,1L],ncol=2L),
              matrix(hinv(tv0[1L]-sqrt(tv0[2L])*qq[,2L]),ncol=2L))
    }
}

perc.ci <- function(t, conf = 0.95, hinv = function(t) t)
#
#  Bootstrap Percentile Confidence Interval Method
#
{
    alpha <- (1+c(-conf,conf))/2
    qq <- norm.inter(t,alpha)
    cbind(conf,matrix(qq[,1L],ncol=2L),matrix(hinv(qq[,2]),ncol=2L))
}

bca.ci <-
    function(boot.out,conf = 0.95,index = 1,t0 = NULL,t = NULL, L = NULL,
             h = function(t) t, hdot = function(t) 1, hinv = function(t) t,
             ...)
#
#  Adjusted Percentile (BCa) Confidence interval method.  This method
#  uses quantities calculated from the empirical influence values to
#  improve on the precentile interval.  Usually the required order
#  statistics for this method will not be integers and so norm.inter
#  is used to find them.
#
{
    t.o <- t
    if (is.null(t) || is.null(t0)) {
        t <- boot.out$t[,index]
        t0 <- boot.out$t0[index]
    }
    t <- t[is.finite(t)]
    w <- qnorm(sum(t < t0)/length(t))
    if (!is.finite(w)) stop("estimated adjustment 'w' is infinite")
    alpha <- (1+c(-conf,conf))/2
    zalpha <- qnorm(alpha)
    if (is.null(L))
        L <- empinf(boot.out, index=index, t=t.o, ...)
    a <- sum(L^3)/(6*sum(L^2)^1.5)
    if (!is.finite(a)) stop("estimated adjustment 'a' is NA")
    adj.alpha <- pnorm(w + (w+zalpha)/(1-a*(w+zalpha)))
    qq <- norm.inter(t,adj.alpha)
    cbind(conf, matrix(qq[,1L],ncol=2L), matrix(hinv(h(qq[,2L])),ncol=2L))
}



abc.ci <- function(data, statistic, index = 1, strata = rep(1, n), conf = 0.95,
                   eps = 0.001/n, ...)
#
#   Non-parametric ABC method for constructing confidence intervals.
#
{
    y <- data
    n <- NROW(y)
    strata1 <- tapply(strata,as.numeric(strata))
    if (length(index) != 1L) {
	warning("only first element of 'index' used in 'abc.ci'")
        index <- index[1L]
    }
    S <- length(table(strata1))
    mat <- matrix(0,n,S)
    for (s in 1L:S) {
        gp <- seq_len(n)[strata1 == s]
        mat[gp,s] <- 1
    }
#  Calculate the observed value of the statistic
    w.orig <- rep(1/n,n)
    t0 <- statistic(y,w.orig/(w.orig%*%mat)[strata1], ...)[index]#
#  Now find the linear and quadratic empirical influence values through
#  numerical differentiation
    L <- L2 <- numeric(n)
    for (i in seq_len(n)) {
        w1 <- (1-eps)*w.orig
        w1[i] <- w1[i]+eps
        w2 <- (1+eps)*w.orig
        w2[i] <- w2[i] - eps
        t1 <- statistic(y,w1/(w1%*%mat)[strata1], ...)[index]
        t2 <- statistic(y,w2/(w2%*%mat)[strata1], ...)[index]
        L[i] <- (t1-t2)/(2*eps)
        L2[i] <- (t1-2*t0+t2)/eps^2
    }
#  Calculate the required quantities for the intervals
    temp1 <- sum(L*L)
    sigmahat <- sqrt(temp1)/n
    ahat <- sum(L^3)/(6*temp1^1.5)      # called a in the text
    bhat <- sum(L2)/(2*n*n)             # called b in the text
    dhat <- L/(n*n*sigmahat)            # called k in the text
    w3 <- w.orig+eps*dhat
    w4 <- w.orig-eps*dhat
    chat <- (statistic(y,w3/(w3%*%mat)[strata1], ...)[index]-2*t0 +
             statistic(y,w4/(w4%*%mat)[strata1], ...)[index]) /
                 (2*eps*eps*sigmahat)   # called c in the text
    bprime <- ahat-(bhat/sigmahat-chat) # called w in the text
    alpha <- (1+as.vector(rbind(-conf,conf)))/2
    zalpha <- qnorm(alpha)
    lalpha <- (bprime+zalpha)/(1-ahat*(bprime+zalpha))^2#
#  Finally calculate the interval endpoints by calling the statistic with
#  various weight vectors.
    out <- seq(alpha)
    for (i in seq_along(alpha)) {
        w.fin <- w.orig+lalpha[i]*dhat
        out[i] <- statistic(y,w.fin/(w.fin%*%mat)[strata1], ...)[index]
    }
    out <- cbind(conf,matrix(out,ncol=2L,byrow=TRUE))
    if (length(conf) == 1L) out <- as.vector(out)
    out
}

censboot <-
    function(data, statistic, R, F.surv, G.surv, strata = matrix(1, n, 2),
             sim = "ordinary", cox = NULL, index = c(1, 2), ...,
             parallel = c("no", "multicore", "snow"),
             ncpus = getOption("boot.ncpus", 1L), cl = NULL)
{
#
#  Bootstrap replication for survival data.  Possible resampling
#  schemes are case, model-based, conditional bootstrap (with or without
#  a model) and the weird bootstrap.
#
    mstrata <- missing(strata)
    if (any(is.na(data)))
        stop("missing values not allowed in 'data'")
    if ((sim != "ordinary") && (sim != "model") && (sim != "cond")
        && (sim != "weird")) stop("unknown value of 'sim'")
    if ((sim == "model") && (is.null(cox))) sim <- "ordinary"
    if (missing(parallel)) parallel <- getOption("boot.parallel", "no")
    parallel <- match.arg(parallel)
    have_mc <- have_snow <- FALSE
    if (parallel != "no" && ncpus > 1L) {
        if (parallel == "multicore") have_mc <- .Platform$OS.type != "windows"
        else if (parallel == "snow") have_snow <- TRUE
        if (!have_mc && !have_snow) ncpus <- 1L
        loadNamespace("parallel") # get this out of the way before recording seed
    }
    if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) runif(1)
    seed <- get(".Random.seed", envir = .GlobalEnv, inherits = FALSE)
    call <- match.call()
    if (isMatrix(data)) n <- nrow(data)
    else stop("'data' must be a matrix with at least 2 columns")
    if (ncol(data) < 2L)
        stop("'data' must be a matrix with at least 2 columns")
    if (length(index) < 2L)
        stop("'index' must contain 2 elements")
    if (length(index) > 2L) {
        warning("only first 2 elements of 'index' used")
        index <- index[1L:2L]
    }
    if (ncol(data) < max(index))
        stop("indices are incompatible with 'ncol(data)'")
    if (sim == "weird") {
        if (!is.null(cox))
            stop("sim = \"weird\" cannot be used with a \"coxph\" object")
        if (ncol(data) > 2L)
            warning(gettextf("only columns %s and %s of 'data' used",
                             index[1L], index[2L]), domain = NA)
        data <- data[,index]
    }
    if (!is.null(cox) && is.null(cox$coefficients) &&
        ((sim == "cond") || (sim == "model"))) {
        warning("no coefficients in Cox model -- model ignored")
        cox <- NULL
    }
    if ((sim != "ordinary")  && missing(F.surv))
        stop("'F.surv' is required but missing")
    if (missing(G.surv) && ((sim == "cond") || (sim == "model")))
        stop("'G.surv' is required but missing")
    if (NROW(strata) != n) stop("'strata' of wrong length")
    if (!isMatrix(strata)) {
        if (!((sim == "weird") || (sim == "ordinary")))
            strata <- cbind(strata, 1)
    } else {
        if ((sim == "weird") || (sim == "ordinary")) strata <- strata[, 1L]
        else  strata <- strata[, 1L:2L]
    }
    temp.str <- strata
    strata <- if (isMatrix(strata))
        apply(strata, 2L, function(s, n) tapply(seq_len(n), as.numeric(s)), n)
    else  tapply(seq_len(n), as.numeric(strata))
    t0 <- if ((sim == "weird") && !mstrata) statistic(data, temp.str, ...)
    else  statistic(data, ...)
    ## Calculate the resampled data sets.  For ordinary resampling this
    ## involves finding the matrix of indices of the case to be resampled.
    ## For the conditional bootstrap or model-based we must find an array
    ## consisting of R matrices containing the resampled times and their
    ## censoring indicators.  The data sets for the weird bootstrap must be
    ## calculated individually.
    fn <- if (sim == "ordinary") {
        bt <- cens.case(n, strata, R)
        function(r) statistic(data[sort(bt[r, ]), ], ...)
    } else if (sim == "weird") {
        ## force promises
        data; F.surv
        if (!mstrata) {
            function(r) {
                bootdata <- cens.weird(data, F.surv, strata)
                statistic(bootdata[, 1:2], bootdata[, 3L], ...)
            }
        } else  {
            function(r) {
                bootdata <- cens.weird(data, F.surv, strata)
                statistic(bootdata[, 1:2], ...)
            }
        }
    } else {
        bt <- cens.resamp(data, R, F.surv, G.surv, strata, index, cox, sim)
        function(r) {
            bootdata <- data
            bootdata[, index] <- bt[r, , ]
            oi <- order(bt[r, , 1L], 1-bt[r, , 2L])
            statistic(bootdata[oi, ], ...)
        }
    }
    rm(mstrata)

    res <- if (ncpus > 1L && (have_mc || have_snow)) {
        if (have_mc) {
            ## ... omitted as from 1.3-27
            parallel::mclapply(seq_len(R), fn, mc.cores = ncpus)
        } else if (have_snow) {
            list(...) # evaluate any promises
            if (is.null(cl)) {
                cl <- parallel::makePSOCKcluster(rep("localhost", ncpus))
                if(RNGkind()[1L] == "L'Ecuyer-CMRG")
                    parallel::clusterSetRNGStream(cl)
                parallel::clusterEvalQ(cl, library(survival))
                res <- parallel::parLapply(cl, seq_len(R), fn)
                parallel::stopCluster(cl)
                res
            } else {
                parallel::clusterEvalQ(cl, library(survival))
                parallel::parLapply(cl, seq_len(R), fn)
            }
       }
    } else lapply(seq_len(R), fn)

    t <- matrix(, R, length(t0))
    for(r in seq_len(R)) t[r, ] <- res[[r]]

    cens.return(sim, t0, t, temp.str, R, data, statistic, call, seed)
}

cens.return <- function(sim, t0, t, strata, R, data, statistic, call, seed) {
#
#  Create an object of class "boot" from the output of a censored bootstrap.
#
    out <- list(t0 = t0, t = t, R = R, sim = sim, data = data, seed = seed,
                statistic = statistic, strata = strata, call = call)
    class(out) <- "boot"
    attr(boot, "boot_type") <- "censboot"
    out
}

cens.case <- function(n, strata, R) {
#
#  Simple case resampling.
#
    out <- matrix(NA, nrow = R, ncol = n)
    for (s in seq_along(table(strata))) {
        inds <- seq_len(n)[strata == s]
        ns <- length(inds)
        out[, inds] <- bsample(inds,  ns*R)
    }
    out
}


cens.weird <- function(data, surv, strata) {
#
#  The weird bootstrap.  Censoring times are fixed and the number of
#  failures at each failure time are sampled from a binomial
#  distribution.  See Chapter 3 of Davison and Hinkley (1997).
#
#  data is a two column matrix containing the times and censoring
#    indicator.
#  surv is a survival object giving the failure time distribution.
#  strata is a the strata vector used in surv or a vector of 1's if no
#    strata were used.
#
    m <- length(surv$time)
    if (is.null(surv$strata)) {
        nstr <- 1
        str <- rep(1, m)
    } else {
        nstr <- length(surv$strata)
        str <- rep(1L:nstr, surv$strata)
    }
    n.ev <- rbinom(m, surv$n.risk, surv$n.event/surv$n.risk)
    while (any(tapply(n.ev, str, sum) == 0))
        n.ev <- rbinom(m, surv$n.risk, surv$n.event/surv$n.risk)
    times <- rep(surv$time, n.ev)
    str <- rep(str, n.ev)
    out <- NULL
    for (s in 1L:nstr) {
        temp <- cbind(times[str == s], 1)
        temp <- rbind(temp,
                      as.matrix(data[(strata == s&data[, 2L] == 0), , drop=FALSE]))
        temp <- cbind(temp, s)
        oi <- order(temp[, 1L], 1-temp[, 2L])
        out <- rbind(out, temp[oi, ])
    }
    if (is.data.frame(data)) out <- as.data.frame(out)
    out
}



cens.resamp <- function(data, R, F.surv, G.surv, strata, index = c(1,2),
                        cox = NULL, sim = "model")
{
#
#  Other types of resampling for the censored bootstrap.  This function
#  uses some local functions to implement the conditional bootstrap for
#  censored data and resampling based on a Cox regression model.  This
#  latter method of sampling can also use conditional sampling to get the
#  censoring times.
#
#  data is the data set
#  R is the number of replicates
#  F.surv is a survfit object for the failure time distribution
#  G.surv is a survfit object for the censoring time distribution
#  strata is a two column matrix, the first column gives the strata
#     gives the strata for the failure times and the second for the
#     censoring times.
#  index is a vector with two integer components giving the position
#     of the times and censoring indicators in data
#  cox is an object returned by the coxph function to give the Cox
#     regression model for the failure times.
#  sim is the simulation type which will always be "model" or "cond"
#
    gety1 <- function(n, R, surv, inds) {
# Sample failure times from the product limit estimate of the failure
# time distribution.
        survival <- surv$surv[inds]
        time <- surv$time[inds]
        n1 <- length(time)
        if (survival[n1] > 0L) {
            survival <- c(survival, 0)
            time <- c(time, Inf)
        }
        probs <- diff(-c(1, survival))
        matrix(bsample(time, n*R, prob = probs), R, n)
    }
    gety2 <- function(n, R, surv, eta, inds) {
# Sample failure times from the Cox regression model.
        F0 <- surv$surv[inds]
        time <- surv$time[inds]
        n1 <- length(time)
        if (F0[n1] > 0) {
            F0 <- c(F0, 0)
            time <- c(time, Inf)
        }
        ex <- exp(eta)
        Fh <- 1 - outer(F0, ex, "^")
        apply(rbind(0, Fh), 2L,
              function(p, y, R) bsample(y, R, prob = diff(p)), time, R)
    }
    getc1 <- function(n, R, surv, inds) {
# Sample censoring times from the product-limit estimate of the
# censoring distribution.
        cens <- surv$surv[inds]
        time <- surv$time[inds]
        n1 <- length(time)
        if (cens[n1] > 0) {
            cens <- c(cens, 0)
            time <- c(time, Inf)
        }
        probs <- diff(-c(1, cens))
        matrix(bsample(time, n*R, prob = probs), nrow = R)
    }
    getc2 <- function(n, R, surv, inds, data, index) {
# Sample censoring times form the conditional distribution.  If a failure
# was observed then sample from the product-limit estimate of the censoring
# distribution conditional on the time being greater than the observed
# failure time.  If the observation is censored then resampled time is the
# observed censoring time.
        cens <- surv$surv[inds]
        time <- surv$time[inds]
        n1 <- length(time)
        if (cens[n1] > 0) {
            cens <- c(cens, 0)
            time <- c(time, Inf)
        }
        probs <- diff(-c(1, cens))
        cout <- matrix(NA, R, n)
        for (i in seq_len(n)) {
            if (data[i, 2] == 0) cout[, i] <- data[i, 1L]
            else {
                pri <- probs[time > data[i, 1L]]
                ti <- time[time > data[i, 1L]]
                if (length(ti) == 1L) cout[, i] <- ti
                else cout[, i] <- bsample(ti, R, prob = pri)
            }
        }
        cout
    }
    n <- nrow(data)
    Fstart <- 1
    Fstr <- F.surv$strata
    if (is.null(Fstr)) Fstr <- length(F.surv$time)
    Gstart <- 1
    Gstr <- G.surv$strata
    if (is.null(Gstr)) Gstr <- length(G.surv$time)
    y0 <- matrix(NA, R, n)
    for (s in seq_along(table(strata[, 1L]))) {
# Find the resampled failure times within strata for failures
        ns <- sum(strata[, 1L] == s)
        inds <- Fstart:(Fstr[s]+Fstart-1)
        y0[, strata[, 1L] == s] <- if (is.null(cox)) gety1(ns, R, F.surv, inds)
        else  gety2(ns, R, F.surv, cox$linear.predictors[strata[, 1L] == s], inds)
        Fstart <- Fstr[s]+Fstart
    }
    c0 <- matrix(NA, R, n)
    for (s in seq_along(table(strata[, 2L]))) {
# Find the resampled censoring times within strata for censoring times
        ns <- sum(strata[, 2] == s)
        inds <- Gstart:(Gstr[s]+Gstart-1)
        c0[, strata[, 2] == s] <- if (sim != "cond") getc1(ns, R, G.surv, inds)
        else  getc2(ns, R, G.surv, inds, data[strata[,2] == s, index])
        Gstart <- Gstr[s]+Gstart
    }
    infs <- (is.infinite(y0) & is.infinite(c0))
    if (sum(infs) > 0) {
# If both the resampled failure time and the resampled censoring time
# are infinite then set the resampled time to be a failure at the largest
# failure time in the failure time stratum containing the observation.
        evs <- seq_len(n)[data[, index[2L]] == 1]
        maxf <- tapply(data[evs, index[1L]], strata[evs, 1L], max)
        maxf <- matrix(maxf[strata[, 1L]], nrow = R, ncol = n, byrow = TRUE)
        y0[infs] <- maxf[infs]
    }
    array(c(pmin(y0, c0), 1*(y0 <= c0)), c(dim(y0), 2))
}

empinf <- function(boot.out = NULL, data = NULL, statistic = NULL,
                   type = NULL, stype = NULL ,index = 1, t = NULL,
                   strata = rep(1, n), eps = 0.001, ...)
{
#
#   Calculation of empirical influence values.  Possible types are
#   "inf" = infinitesimal jackknife (numerical differentiation)
#   "reg" = regression based estimation
#   "jack" = usual jackknife estimates
#   "pos" = positive jackknife estimates
#
    if (!is.null(boot.out))
    {
	if (boot.out$sim == "parametric")
            stop("influence values cannot be found from a parametric bootstrap")
        data <- boot.out$data
        if (is.null(statistic))
            statistic <- boot.out$statistic
        if (is.null(stype))
            stype <- boot.out$stype
        if (!is.null(boot.out$strata))
            strata <- boot.out$strata
    }
    else
    {
	if (is.null(data))
            stop("neither 'data' nor bootstrap object specified")
        if (is.null(statistic))
            stop("neither 'statistic' nor bootstrap object specified")
        if (is.null(stype)) stype <- "w"
    }
    n <- NROW(data)
    if (is.null(type)) {
        if (!is.null(t)) type <- "reg"
        else if (stype == "w") type <- "inf"
        else if (!is.null(boot.out) &&
                 (boot.out$sim != "parametric") &&
                 (boot.out$sim != "permutation")) type <- "reg"
        else type <- "jack"
    }

    if (type == "inf") {
# calculate the infinitesimal jackknife values by numerical differentiation
    	if (stype !="w") stop("'stype' must be \"w\" for type=\"inf\"")
        if (length(index) != 1L) {
            warning("only first element of 'index' used")
            index <- index[1L]
        }
        if (!is.null(t))
            warning("input 't' ignored; type=\"inf\"")
        L <- inf.jack(data, statistic, index, strata, eps, ...)
    } else if (type == "reg") {
# calculate the regression estimates of the influence values
        if (is.null(boot.out))
            stop("bootstrap object needed for type=\"reg\"")
        if (is.null(t)) {
            if (length(index) != 1L) {
                warning("only first element of 'index' used")
                index <- index[1L]
            }
            t <- boot.out$t[,index]
        }
        L <- empinf.reg(boot.out, t)
    } else if (type == "jack") {
        if (!is.null(t))
            warning("input 't' ignored; type=\"jack\"")
        if (length(index) != 1L) {
            warning("only first element of 'index' used")
            index <- index[1L]
        }
        L <- usual.jack(data, statistic, stype, index, strata, ...)
    } else if (type == "pos") {
        if (!is.null(t))
            warning("input 't' ignored; type=\"pos\"")
        if (length(index) != 1L) {
            warning("only first element of 'index' used")
            index <- index[1L]
        }
        L <- positive.jack(data, statistic, stype, index, strata, ...)
    }
    L
}

inf.jack <-
    function(data, stat, index = 1, strata  =  rep(1, n), eps  =  0.001, ...)
{
#
#   Numerical differentiation to get infinitesimal jackknife estimates
#   of the empirical influence values.
#
    n <- NROW(data)
    L <- seq_len(n)
    eps <- eps/n
    strata <- tapply(strata, as.numeric(strata))
    w.orig <- 1/table(strata)[strata]
    tobs <- stat(data, w.orig, ...)[index]
    for(i in seq_len(n)) {
        group <- seq_len(n)[strata == strata[i]]
        w <- w.orig
        w[group] <- (1 - eps)*w[group]
        w[i] <- w[i] + eps
        L[i] <- (stat(data, w, ...)[index] - tobs)/eps
    }
    L
}

empinf.reg <- function(boot.out, t = boot.out$t[,1L])
#
#  Function to estimate empirical influence values using regression.
#  This method regresses the observed bootstrap values on the bootstrap
#  frequencies to estimate the empirical influence values
#
{
    fins <- seq_along(t)[is.finite(t)]
    t <- t[fins]
    R <- length(t)
    n <- NROW(boot.out$data)
    strata <- boot.out$strata
    if (is.null(strata))
        strata <- rep(1,n)
    else 	strata <- tapply(strata,as.numeric(strata))
    ns <- table(strata)
#    S <- length(ns)
    f <- boot.array(boot.out)[fins,]
    X <- f/matrix(ns[strata], R, n ,byrow=TRUE)
    out <- tapply(seq_len(n), strata, min)
    inc <- seq_len(n)[-out]
    X <- X[,inc]
    beta <- coefficients(glm(t ~ X))[-1L]
    l <- rep(0, n)
    l[inc] <- beta
    l <- l - tapply(l,strata,mean)[strata]
    l
}

usual.jack <- function(data, stat, stype = "w", index = 1,
                       strata = rep(1, n), ...)
#
#  Function to use the normal (delete 1) jackknife method to estimate the
#  empirical influence values
#
{
    n <- NROW(data)
    l <- rep(0,n)
    strata <- tapply(strata,as.numeric(strata))
    if (stype == "w") {
        w0 <- rep(1, n)/table(strata)[strata]
        tobs <- stat(data, w0, ...)[index]
        for (i in seq_len(n)) {
            w1 <- w0
            w1[i] <- 0
            gp <- strata == strata[i]
            w1[gp] <- w1[gp]/sum(w1[gp])
            l[i] <- (sum(gp)-1)*(tobs - stat(data,w1, ...)[index])
        }
    } else if (stype == "f") {
        f0 <- rep(1,n)
        tobs <- stat(data, f0, ...)[index]
        for (i in seq_len(n)) {
            f1 <- f0
            f1[i] <- 0
            gp <- strata == strata[i]
            l[i] <- (sum(gp)-1)*(tobs - stat(data, f1, ...)[index])
        }
    } else {
        i0 <- seq_len(n)
        tobs <- stat(data, i0, ...)[index]
        for (i in seq_len(n)) {
            i1 <- i0[-i]
            gp <- strata == strata[i]
            l[i] <- (sum(gp)-1)*(tobs - stat(data, i1, ...)[index])
        }
    }
    l
}

positive.jack <- function(data, stat, stype = "w", index = 1,
                          strata = rep(1 ,n), ...)
{
#
#  Use the positive jackknife to estimate the empirical influence values.
#  The positive jackknife includes one observation twice to find its
#  influence.
#
    strata <- tapply(strata,as.numeric(strata))
    n <- NROW(data)
    L <- rep(0, n)
    if (stype == "w") {
        w0 <- rep(1, n)/table(strata)[strata]
        tobs <- stat(data, w0, ...)[index]
        for (i in seq_len(n)) {
            st1 <- c(strata,strata[i])
            w1 <- 1/table(st1)[strata]
            w1[i] <- 2*w1[i]
            gp <- strata == strata[i]
            w1[gp] <- w1[gp]/sum(w1[gp])
            L[i] <- (sum(gp)+1)*(stat(data, w1, ...)[index] - tobs)
        }
    } else if (stype == "f") {
        f0 <- rep(1,n)
        tobs <- stat(data, f0, ...)[index]
        for (i in seq_len(n)) {
            f1 <- f0
            f1[i] <- 2
            gp <- strata == strata[i]
            L[i] <- (sum(gp)+1)*(stat(data, f1, ...)[index] - tobs)
        }
    } else if (stype == "i") {
        i0 <- seq_len(n)
        tobs <- stat(data, i0, ...)[index]
        for (i in seq_len(n)) {
            i1 <- c(i0, i)
            gp <- strata == strata[i]
            L[i] <- (sum(gp)+1)*(stat(data, i1, ...)[index] - tobs)
        }
    }
    L
}

linear.approx <- function(boot.out, L = NULL, index = 1, type = NULL,
                          t0 = NULL, t = NULL, ...)
#
#  Find the linear approximation to the bootstrap replicates of a
#  statistic.  L should be the linear influence values which will
#  be found by empinf if they are not supplied.
#
{
    f <- boot.array(boot.out)
    n <- length(f[1,  ])
    if ((length(index) > 1L) && (is.null(t0) || is.null(t))) {
        warning("only first element of 'index' used")
        index <- index[1L]
    }
    if (is.null(t0)) {
        t0 <- boot.out$t0[index]
        if (is.null(L))
            L <- empinf(boot.out, index=index, type=type, ...)
    } else if (is.null(t) && is.null(L)) {
        warning("input 't0' ignored: neither 't' nor 'L' supplied")
        t0 <- t0[index]
        L <- empinf(boot.out, index=index, type=type, ...)
    }
    else if (is.null(L))
        L <- empinf(boot.out, type=type, t=t, ...)
    tL <- rep(t0, boot.out$R)
    strata <- boot.out$strata
    if (is.null(strata))
        strata <- rep(1, n)
    else 	strata <- tapply(strata,as.numeric(strata))
    S <- length(table(strata))
    for(s in 1L:S) {
        i.s <- seq_len(n)[strata == s]
        tL <- tL + f[, i.s] %*% L[i.s]/length(i.s)
    }
    as.vector(tL)
}

envelope <-
    function(boot.out = NULL, mat = NULL, level = 0.95, index = 1L:ncol(mat))
#
#  Function to estimate pointwise and overall confidence envelopes for
#  a function.
#
#  mat is a matrix of bootstrap values of the function at a number of
#     points.  The points at which they are evaluated are assumed to
#     be constant over the rows.
#
{
    emperr <- function(rmat, p = 0.05, k = NULL)
#  Local function to estimate the overall error rate of an envelope.
    {
        R <- nrow(rmat)
        if (is.null(k)) k <- p*(R+1)/2 else p <- 2*k/(R+1)
        kf <- function(x, k, R) 1*((min(x) <= k)|(max(x) >= R+1L-k))
        c(k, p, sum(apply(rmat, 1L, kf, k, R))/(R+1))
    }
    kfun <- function(x, k1, k2)
# Local function to find the cut-off points in each column of the matrix.
        sort(x ,partial = sort(c(k1, k2)))[c(k1, k2)]
    if (!is.null(boot.out) && isMatrix(boot.out$t)) mat <- boot.out$t
    if (!isMatrix(mat)) stop("bootstrap output matrix missing")
    if (length(index) < 2L) stop("use 'boot.ci' for scalar parameters")
    mat <- mat[,index]
    rmat <- apply(mat,2L,rank)
    R <- nrow(mat)
    if (length(level) == 1L) level <- rep(level,2L)
    k.pt <- floor((R+1)*(1-level[1L])/2+1e-10)
    k.pt <- c(k.pt, R+1-k.pt)
    err.pt <- emperr(rmat,k = k.pt[1L])
    ov <- emperr(rmat,k = 1)
    ee <- err.pt
    al <- 1-level[2L]
    if (ov[3L] > al)
        warning("unable to achieve requested overall error rate")
    else {
        continue <- !(ee[3L] < al)
        while(continue) {
#  If the observed error is greater than the level required for the overall
#  envelope then try another envelope.  This loop uses linear interpolation
#  on the integers between 1 and k.pt[1L] to find the required value.
            kk <- ov[1L]+round((ee[1L]-ov[1L])*(al-ov[3L])/ (ee[3L]-ov[3L]))
            if (kk == ov[1L]) kk <- kk+1
            else if (kk == ee[1L]) kk <- kk-1
            temp <- emperr(rmat, k = kk)
            if (temp[3L] > al) ee <- temp
            else ov <- temp
            continue <- !(ee[1L] == ov[1L]+1)
        }
    }
    k.ov <- c(ov[1L], R+1-ov[1L])
    err.ov <- ov[-1L]
    out <- apply(mat, 2L, kfun, k.pt, k.ov)
    list(point = out[2:1,], overall = out[4:3,], k.pt = k.pt,
         err.pt = err.pt[-1L], k.ov = k.ov, err.ov = err.ov, err.nom = 1-level)
}


glm.diag <- function(glmfit)
{
#
#  Calculate diagnostics for objects of class "glm".  The diagnostics
#  calculated are various types of residuals as well as the Cook statistics
#  and the leverages.
#
    w <- if (is.null(glmfit$prior.weights)) rep(1,length(glmfit$residuals))
         else glmfit$prior.weights
    sd <- switch(family(glmfit)$family[1L],
                 "gaussian" = sqrt(glmfit$deviance/glmfit$df.residual),
                 "Gamma" = sqrt(sum(w*(glmfit$y/fitted(glmfit) - 1)^2)/
                           glmfit$df.residual),
                 1)
##     sd <- ifelse(family(glmfit)$family[1L] == "gaussian",
##                  sqrt(glmfit$deviance/glmfit$df.residual), 1)
##     sd <- ifelse(family(glmfit)$family[1L] == "Gamma",
##                  sqrt(sum(w*(glmfit$y/fitted(glmfit) - 1)^2)/glmfit$df.residual), sd)
    dev <- residuals(glmfit, type = "deviance")/sd
    pear <- residuals(glmfit, type = "pearson")/sd
    ## R change: lm.influence drops 0-wt cases.
    h <- rep(0, length(w))
    h[w != 0] <- lm.influence(glmfit)$hat
    p <- glmfit$rank
    rp <- pear/sqrt(1 - h)
    rd <- dev/sqrt(1 - h)
    cook <- (h * rp^2)/((1 - h) * p)
    res <- sign(dev) * sqrt(dev^2 + h * rp^2)
    list(res = res, rd = rd, rp = rp, cook = cook, h = h, sd = sd)
}


glm.diag.plots <-
    function(glmfit, glmdiag = glm.diag(glmfit), subset  =  NULL,
             iden = FALSE, labels = NULL, ret = FALSE)
{
#  Diagnostic plots for objects of class "glm"
    if (is.null(glmdiag))
        glmdiag <- glm.diag(glmfit)
    if (is.null(subset))
        subset <- seq_along(glmdiag$h)
    else if (is.logical(subset))
        subset <- seq_along(subset)[subset]
    else if (is.numeric(subset) && all(subset<0))
        subset <- (1L:(length(subset)+length(glmdiag$h)))[subset]
    else if (is.character(subset)) {
        if (is.null(labels)) labels <- subset
        subset <- seq_along(subset)
    }
#	close.screen(all = T)
#	split.screen(c(2, 2))
#	screen(1) #
    par(mfrow = c(2,2))
#  Plot the deviance residuals against the fitted values
    x1 <- predict(glmfit)
    plot(x1, glmdiag$res, xlab = "Linear predictor", ylab = "Residuals")
    pars <- vector(4L, mode="list")
    pars[[1L]] <- par("usr")
#	screen(2) #
#  Plot a normal QQ plot of the standardized deviance residuals
    y2 <- glmdiag$rd
    x2 <- qnorm(ppoints(length(y2)))[rank(y2)]
    plot(x2, y2, ylab = "Quantiles of standard normal",
         xlab = "Ordered deviance residuals")
    abline(0, 1, lty = 2)
    pars[[2L]] <- par("usr")
#	screen(3) #
#  Plot the Cook statistics against h/(1-h) and draw line to highlight
#  possible influential and high leverage points.
    hh <- glmdiag$h/(1 - glmdiag$h)
    plot(hh, glmdiag$cook, xlab = "h/(1-h)", ylab = "Cook statistic")
    rx <- range(hh)
    ry <- range(glmdiag$cook)
    rank.fit <- glmfit$rank
    nobs <- rank.fit + glmfit$df.residual
    cooky <- 8/(nobs - 2 * rank.fit)
    hy <- (2 * rank.fit)/(nobs - 2 * rank.fit)
    if ((cooky >= ry[1L]) && (cooky <= ry[2L])) abline(h = cooky, lty = 2)
    if ((hy >= rx[1L]) && (hy <= rx[2L])) abline(v = hy, lty = 2)
    pars[[3L]] <- par("usr")
#	screen(4) #
#  Plot the Cook statistics against the observation number in the original
#  data set.
    plot(subset, glmdiag$cook, xlab = "Case", ylab = "Cook statistic")
    if ((cooky >= ry[1L]) && (cooky <= ry[2L])) abline(h = cooky, lty = 2)
    xx <- list(x1,x2,hh,subset)
    yy <- list(glmdiag$res, y2, glmdiag$cook, glmdiag$cook)
    pars[[4L]] <- par("usr")

    if (is.null(labels)) labels <- names(x1)
    while (iden) {
#  If interaction with the plots is required then ask the user which plot
#  they wish to interact with and then run identify() on that plot.
#  When the user terminates identify(), reprompt until no further interaction
#  is required and the user inputs a 0.
        cat("****************************************************\n")
        cat("Please Input a screen number (1,2,3 or 4)\n")
        cat("0 will terminate the function \n")
#		num <- scan(nmax=1)
        num <- as.numeric(readline())
        if ((length(num) > 0L) &&
            ((num == 1)||(num == 2)||(num == 3)||(num == 4))) {
            cat(paste("Interactive Identification for screen",
                      num,"\n"))
            cat("left button = Identify, center button = Exit\n")
#			screen(num, new=F)
            nm <- num+1
            par(mfg = c(trunc(nm/2),1 +nm%%2, 2, 2))
            par(usr = pars[[num]])
            identify(xx[[num]], yy[[num]], labels)
        }
        else 	iden <- FALSE
    }
#	close.screen(all=T)
    par(mfrow = c(1, 1))
    if (ret) glmdiag else invisible()
}

exp.tilt <- function(L, theta = NULL, t0 = 0, lambda = NULL,
                     strata = rep(1, length(L)) )
{
# exponential tilting of linear approximation to statistic
# to give mean theta.
#
    tilt.dis <- function(lambda)  {
#  Find the squared error in the mean using the multiplier lambda
#  This is then minimized to find the correct value of lambda
#  Note that the function should have minimum 0.
        L <- para[[2L]]
        theta <- para[[1L]]
        strata <- para[[3L]]
        ns <- table(strata)
        tilt <- rep(NA, length(L) )
        for (s in seq_along(ns)) {
            p <- exp(lambda*L[strata == s]/ns[s])
            tilt[strata == s] <- p/sum(p)
        }
        (sum(L*tilt) - theta)^2
    }
    tilted.prob <- function(lambda, L, strata)  {
#  Find the tilted probabilities for a given value of lambda
        ns <- table(strata)
        m <- length(lambda)
        tilt <- matrix(NA, m, length(L))
        for (i in 1L:m)
            for (s in seq_along(ns)) {
                p <- exp(lambda[i]*L[strata == s]/ns[s])
                tilt[i,strata == s] <- p/sum(p)
            }
        if (m == 1) tilt <- as.vector(tilt)
        tilt
    }
    strata <- tapply(strata, as.numeric(strata))
    if (!is.null(theta)) {
        theta <- theta-t0
        m <- length(theta)
        lambda <- rep(NA,m)
        for (i in 1L:m) {
            para <- list(theta[i],L,strata)
#			assign("para",para,frame=1)
#			lambda[i] <- nlmin(tilt.dis, 0 )$x
            lambda[i] <- optim(0, tilt.dis, method = "BFGS")$par
            msd <- tilt.dis(lambda[i])
            if (is.na(msd) || (abs(msd) > 1e-6))
                stop(gettextf("unable to find multiplier for %f", theta[i]),
                     domain = NA)
        }
    }
    else if (is.null(lambda))
        stop("'theta' or 'lambda' required")
    probs <- tilted.prob( lambda, L, strata )
    if (is.null(theta)) theta <- t0 + sum(probs * L)
    else theta <- theta+t0
    list(p = probs, theta = theta, lambda = lambda)
}


imp.weights <- function(boot.out, def = TRUE, q = NULL)
{
#
# Takes boot.out object and calculates importance weights
# for each element of boot.out$t, as if sampling from multinomial
# distribution with probabilities q.
# If q is NULL the weights are calculated as if
# sampling from a distribution with equal probabilities.
# If def=T calculates weights using defensive mixture
# distribution, if F uses weights knowing from which element of
# the mixture they come.
#
    R <- boot.out$R
    if (length(R) == 1L)
        def <- FALSE
    f <- boot.array(boot.out)
    n <- ncol(f)
    strata <- tapply(boot.out$strata,as.numeric(boot.out$strata))
#    ns <- table(strata)
    if (is.null(q))  q <- rep(1,ncol(f))
    if (any(q == 0)) stop("0 elements not allowed in 'q'")
    p <- boot.out$weights
    if ((length(R) == 1L) && all(abs(p - q)/p < 1e-10))
        return(rep(1, R))
    np <- length(R)
    q <- normalize(q, strata)
    lw.q <- as.vector(f %*% log(q))
    if (!isMatrix(p))
        p <- as.matrix(t(p))
    p <- t(apply(p, 1L, normalize, strata))
    lw.p <- matrix(NA, sum(R), np)
    for(i in 1L:np) {
        zz <- seq_len(n)[p[i,  ] > 0]
        lw.p[, i] <- f[, zz] %*% log(p[i, zz])
    }
    if (def)
        w <- 1/(exp(lw.p - lw.q) %*% R/sum(R))
    else {
        i <- cbind(seq_len(sum(R)), rep(seq_along(R), R))
        w <- exp(lw.q - lw.p[i])
    }
    as.vector(w)
}

const <- function(w, eps=1e-8) {
# Are all of the values of w equal to within the tolerance eps.
    all(abs(w-mean(w, na.rm=TRUE)) < eps)
}

imp.moments <- function(boot.out=NULL, index=1, t=boot.out$t[,index],
			w=NULL, def=TRUE, q=NULL )
{
# Calculates raw, ratio, and regression estimates of mean and
# variance of t using importance sampling weights in w.
    if (missing(t) && is.null(boot.out$t))
        stop("bootstrap replicates must be supplied")
    if (is.null(w))
        if (!is.null(boot.out))
            w <- imp.weights(boot.out, def, q)
        else	stop("either 'boot.out' or 'w' must be specified.")
    if ((length(index) > 1L) && missing(t)) {
        warning("only first element of 'index' used")
        t <- boot.out$t[,index[1L]]
    }
    fins <- seq_along(t)[is.finite(t)]
    t <- t[fins]
    w <- w[fins]
    if (!const(w)) {
        y <- t*w
        m.raw <- mean( y )
        m.rat <- sum( y )/sum( w )
        t.lm <- lm( y~w )
        m.reg <- mean( y ) - coefficients(t.lm)[2L]*(mean(w)-1)
        v.raw <- mean(w*(t-m.raw)^2)
        v.rat <- sum(w/sum(w)*(t-m.rat)^2)
        x <- w*(t-m.reg)^2
        t.lm2 <- lm( x~w )
        v.reg <- mean( x ) - coefficients(t.lm2)[2L]*(mean(w)-1)
    }
    else {	m.raw <- m.rat <- m.reg <- mean(t)
		v.raw <- v.rat <- v.reg <- var(t)
            }
    list( raw=c(m.raw,v.raw), rat = c(m.rat,v.rat),
         reg = as.vector(c(m.reg,v.reg)))
}


imp.reg <- function(w)
{
#  This function takes a vector of importance sampling weights and
#  returns the regression importance sampling weights.  The function
#  is called by imp.prob and imp.quantiles to enable those functions
#  to find regression estimates of tail probabilities and quantiles.
    if (!const(w)) {
        R <- length(w)
        mw <- mean(w)
        s2w <- (R-1)/R*var(w)
        b <- (1-mw)/s2w
        w <- w*(1+b*(w-mw))/R
    }
    cumsum(w)/sum(w)
}


imp.quantile <- function(boot.out=NULL, alpha=NULL, index=1,
			t=boot.out$t[,index], w=NULL, def=TRUE, q=NULL )
{
# Calculates raw, ratio, and regression estimates of alpha quantiles
#  of t using importance sampling weights in w.
    if (missing(t) && is.null(boot.out$t))
        stop("bootstrap replicates must be supplied")
    if (is.null(alpha)) alpha <- c(0.01,0.025,0.05,0.95,0.975,0.99)
    if (is.null(w))
        if (!is.null(boot.out))
            w <- imp.weights(boot.out, def, q)
        else	stop("either 'boot.out' or 'w' must be specified.")
    if ((length(index) > 1L) && missing(t)){
        warning("only first element of 'index' used")
        t <- boot.out$t[,index[1L]]
    }
    fins <- seq_along(t)[is.finite(t)]
    t <- t[fins]
    w <- w[fins]
    o <- order(t)
    t <- t[o]
    w <- w[o]
    cum <- cumsum(w)
    o <- rev(o)
    w.m <- w[o]
    t.m <- -rev(t)
    cum.m <- cumsum(w.m)
    cum.rat <- cum/mean(w)
    cum.reg <- imp.reg(w)
    R <- length(w)
    raw <- rat <- reg <- rep(NA,length(alpha))
    for (i in seq_along(alpha)) {
        if (alpha[i]<=0.5) raw[i] <-  max(t[cum<=(R+1)*alpha[i]])
        else raw[i] <- -max(t.m[cum.m<=(R+1)*(1-alpha[i])])
        rat[i] <- max(t[cum.rat <= (R+1)*alpha[i]])
        reg[i] <- max(t[cum.reg <= (R+1)*alpha[i]])
    }
    list(alpha=alpha, raw=raw, rat=rat, reg=reg)
}

imp.prob <- function(boot.out=NULL, index=1, t0=boot.out$t0[index],
			t=boot.out$t[,index], w=NULL,  def=TRUE, q=NULL)
{
# Calculates raw, ratio, and regression estimates of tail probability
#  pr( t <= t0 ) using importance sampling weights in w.
    is.missing <- function(x) length(x) == 0L || is.na(x)

    if (missing(t) && is.null(boot.out$t))
        stop("bootstrap replicates must be supplied")
    if (is.null(w))
        if (!is.null(boot.out))
            w <- imp.weights(boot.out, def, q)
        else	stop("either 'boot.out' or 'w' must be specified.")
    if ((length(index) > 1L) && (missing(t) || missing(t0))) {
        warning("only first element of 'index' used")
        index <- index[1L]
        if (is.missing(t)) t <- boot.out$t[,index]
        if (is.missing(t0)) t0 <- boot.out$t0[index]
    }
    fins <- seq_along(t)[is.finite(t)]
    t <- t[fins]
    w <- w[fins]
    o <- order(t)
    t <- t[o]
    w <- w[o]
    raw <- rat <- reg <- rep(NA,length(t0))
    cum <- cumsum(w)/sum(w)
    cum.r <- imp.reg(w)
    for (i in seq_along(t0)) {
        raw[i] <-sum(w[t<=t0[i]])/length(w)
        rat[i] <- max(cum[t<=t0[i]])
        reg[i] <- max(cum.r[t<=t0[i]])
    }
    list(t0=t0, raw=raw, rat=rat, reg=reg )
}

smooth.f <- function(theta, boot.out, index=1, t=boot.out$t[,index],
			width=0.5 )
{
# Does frequency smoothing of the frequency array for boot.out with
# bandwidth A to give frequencies for 'typical' distribution at theta
    if ((length(index) > 1L) && missing(t)) {
        warning("only first element of 'index' used")
        t <- boot.out$t[,index[1L]]
    }
    if (isMatrix(t)) {
        warning("only first column of 't' used")
        t <- t[,1L]
    }
    fins <- seq_along(t)[is.finite(t)]
    t <- t[fins]
    m <- length(theta)
    v <- imp.moments(boot.out, t=t)$reg[2L]
    eps <- width*sqrt(v)
    if (m  == 1)
        w <- dnorm((theta-t)/eps )/eps
    else {
        w <- matrix(0,length(t),m)
        for (i in 1L:m)
            w[,i] <- dnorm((theta[i]-t)/eps )/eps
    }
    f <- crossprod(boot.array(boot.out)[fins,] , w)
    strata <- boot.out$strata
    strata <- tapply(strata, as.numeric(strata))
    ns <- table(strata)
    out <- matrix(NA,ncol(f),nrow(f))
    for (s in seq_along(ns)) {
        ts <- matrix(f[strata == s,],m,ns[s],byrow=TRUE)
        ss <- apply(ts,1L,sum)
        out[,strata == s] <-  ts/matrix(ss,m,ns[s])
    }
    if (m == 1) out <- as.vector(out)
    out
}

tilt.boot <- function(data, statistic, R, sim="ordinary",
		stype="i", strata = rep(1, n), L = NULL, theta=NULL,
		alpha=c(0.025,0.975), tilt=TRUE, width=0.5, index=1, ... )
{
#  Does tilted bootstrap sampling of stat applied to data with strata strata
#  and simulation type sim.
#  The levels of R give the number of simulations at each level.  For example,
#  R=c(199,100,50) will give three separate bootstraps with 199, 100, 50
#  simulations.  If R[1L]>0 the first simulation is assumed to be untilted
#  and L can be estimated from it by regression, or it can be frequency
#  smoothed to give probabilities p.
#  If tilt=T use exponential tilting with empirical influence value L
#  given explicitly or estimated from boot0, but if tilt=F
#  (in which case R[1L] should be large) frequency smoothing of boot0 is used
#  with bandwidth A.
#  Tilting/frequency smoothing is to theta (so length(theta)=length(R)-1).
#  The function assumes at present that q=0 is the median of the distribution
#  of t*.
    if ((sim != "ordinary") && (sim != "balanced"))
        stop("invalid value of 'sim' supplied")
    if (!is.null(theta) && (length(R) != length(theta)+1))
        stop("'R' and 'theta' have incompatible lengths")
    if (!tilt && (R[1L] == 0))
        stop("R[1L] must be positive for frequency smoothing")
    call <- match.call()
    n <- NROW(data)
    if (R[1L]>0) {
# If required run an initial bootstrap with equal weights.
        if (is.null(theta) && (length(R) != length(alpha)+1))
            stop("'R' and 'alpha' have incompatible lengths")
        boot0 <- boot(data, statistic, R = R[1L], sim=sim, stype=stype,
                      strata = strata, ... )
        if (is.null(theta)) {
            if (any(c(alpha,1-alpha)*(R[1L]+1) <= 5))
                warning("extreme values used for quantiles")
            theta <- quantile(boot0$t[,index],alpha)
        }
    }
    else {
# If no initial bootstrap is run then exponential tilting must be
# used.  Also set up a dummy bootstrap object to hold the output.
        tilt <- TRUE
        if (is.null(theta))
            stop("'theta' must be supplied if R[1L] = 0")
        if (!missing(alpha))
            warning("'alpha' ignored; R[1L] = 0")
        if (stype == "i") orig <- seq_len(n)
        else if (stype == "f") orig <- rep(1,n)
        else orig <- rep(1,n)/n
        boot0 <- boot.return(sim=sim,t0=statistic(data,orig, ...),
                             t=NULL, strata=strata, R=0, data=data,
                             stat=statistic, stype=stype,call=NULL,
                             seed=get(".Random.seed", envir=.GlobalEnv, inherits = FALSE),
                             m=0,weights=NULL)
    }
# Calculate the weights for the subsequent bootstraps
    if (is.null(L) & tilt)
        if (R[1L] > 0) L <- empinf(boot0, index, ...)
        else L <- empinf(data=data, statistic=statistic, stype=stype,
                         index=index, ...)
    if (tilt) probs <- exp.tilt(L, theta, strata=strata, t0=boot0$t0[index])$p
    else probs <- smooth.f(theta, boot0, index, width=width)#
# Run the weighted bootstraps and collect the output.
    boot1 <- boot(data, statistic, R[-1L], sim=sim, stype=stype,
                  strata=strata, weights=probs, ...)
    boot0$t <- rbind(boot0$t, boot1$t)
    boot0$weights <- rbind(boot0$weights, boot1$weights)
    boot0$R <- c(boot0$R, boot1$R)
    boot0$call <- call
    boot0$theta <- theta
    attr(boot0, "boot_type") <- "tilt.boot"
    boot0
}


control <- function(boot.out, L=NULL, distn=NULL, index=1, t0=NULL, t=NULL,
                    bias.adj=FALSE, alpha=NULL, ... )
{
#
#  Control variate estimation.  Post-simulation balance can be used to
#  find the adjusted bias estimate.  Alternatively the linear approximation
#  to the statistic of interest can be used as a control variate and hence
#  moments and quantiles can be estimated.
#
    if (!is.null(boot.out$call$weights))
        stop("control methods undefined when 'boot.out' has weights")
    if (is.null(alpha))
        alpha <- c(1,2.5,5,10,20,50,80,90,95,97.5,99)/100
    tL <- dL <- bias <- bias.L <- var.L <- NULL
    k3.L <- q.out <- distn.L <- NULL
    stat <- boot.out$statistic
    data <- boot.out$data
    R <- boot.out$R
    f <- boot.array(boot.out)
    if (bias.adj) {
# Find the adjusted bias estimate using post-simulation balance.
        if (length(index) > 1L) {
            warning("only first element of 'index' used")
            index <- index[1L]
        }
        f.big <- apply(f, 2L, sum)
        if (boot.out$stype == "i")
        { 	n <- ncol(f)
                i.big <- rep(seq_len(n),f.big)
                t.big <- stat(data, i.big, ...)[index]
            }
        else if (boot.out$stype == "f")
            t.big <- stat(data, f.big, ...)[index]
        else if (boot.out$stype == "w")
            t.big <- stat(data, f.big/R, ...)[index]
        bias <- mean(boot.out$t[, index]) - t.big
        out <- bias
    }
    else {
# Using the linear approximation as a control variable, find estimates
# of the moments and quantiles of the statistic of interest.
        if (is.null(t) || is.null(t0)) {
            if (length(index) > 1L) {
                warning("only first element of 'index' used")
                index <- index[1L]
            }
            if (is.null(L))
                L <- empinf(boot.out, index=index, ...)
            tL <- linear.approx(boot.out, L, index, ...)
            t <- boot.out$t[,index]
            t0 <- boot.out$t0[index]
        }
        else {
            if (is.null(L))
                L <- empinf(boot.out, t=t, ...)
            tL <- linear.approx(boot.out, L, t0=t0, ...)
        }
        fins <- seq_along(t)[is.finite(t)]
        t <- t[fins]
        tL <- tL[fins]
        R <- length(t)
        dL <- t - tL                    #
# Find the moments of the statistic of interest.
        bias.L <- mean(dL)
        strata <- tapply(boot.out$strata, as.numeric(boot.out$strata))
        var.L <- var.linear(L, strata) + 2*var(tL, dL) + var(dL)
        k3.L <- k3.linear(L, strata) + 3 * cum3(tL, dL) +
            3 * cum3(dL, tL) + cum3(dL)
        if (is.null(distn)) {
# If distn is not supplied then calculate the saddlepoint approximation to
# the distribution of the linear approximation.
            distn <- saddle.distn((t0+L)/length(L),
                                  alpha = (1L:R)/(R + 1),
                                  t0=c(t0,sqrt(var.L)), strata=strata)
            dist.q <- distn$quantiles[,2]
            distn <- distn$distn
        }
        else	dist.q <- predict(distn, x=qnorm((1L:R)/(R+1)))$y#
# Use the quantiles of the distribution of the linear approximation and
# the control variates to estimate the quantiles of the statistic of interest.
        distn.L <- sort(dL[order(tL)] + dist.q)
        q.out <- distn.L[(R + 1) * alpha]
        out <- list(L=L, tL=tL, bias=bias.L, var=var.L, k3=k3.L,
                    quantiles=cbind(alpha,q.out), distn=distn)
    }
    out
}

var.linear <- function(L, strata = NULL)
{
#  estimate the variance of a statistic using its linear approximation
    vL <- 0
    n <- length(L)
    if (is.null(strata))
        strata <- rep(1, n)
    else 	strata <- tapply(seq_len(n),as.numeric(strata))
    S <- length(table(strata))
    for(s in 1L:S) {
        i.s <- seq_len(n)[strata == s]
        vL <- vL + sum(L[i.s]^2/length(i.s)^2)
    }
    vL
}

k3.linear <- function(L, strata = NULL)
{
#  estimate the skewness of a statistic using its linear approximation
    k3L <- 0
    n <- length(L)
    if (is.null(strata))
        strata <- rep(1, n)
    else	strata <- tapply(seq_len(n),as.numeric(strata))
    S <- length(table(strata))
    for(s in 1L:S) {
        i.s <- seq_len(n)[strata == s]
        k3L <- k3L + sum(L[i.s]^3/length(i.s)^3)
    }
    k3L
}

cum3 <- function(a, b=a, c=a, unbiased=TRUE)
# calculate third order cumulants.
{
    n <- length(a)
    if (unbiased) mult <- n/((n-1)*(n-2))
    else mult <- 1/n
    mult*sum((a - mean(a)) * (b - mean(b)) * (c - mean(c)))
}

logit <- function(p) qlogis(p)
#
#  Calculate the logit of a proportion in the range [0,1]
#
## {
##     out <- p
##     inds <- seq_along(p)[!is.na(p)]
##     if (any((p[inds] < 0) | (p[inds] > 1)))
##         stop("invalid proportions input")
##     out[inds] <- log(p[inds]/(1-p[inds]))
##     out[inds][p[inds] == 0] <- -Inf
##     out[inds][p[inds] == 1] <- Inf
##     out
## }

inv.logit <- function(x)
#
#  Calculate the inverse logit of a number
#
# {
#     out <- exp(x)/(1+exp(x))
#     out[x==-Inf] <- 0
#     out[x==Inf] <- 1
#     out
# }
plogis(x)

iden <- function(n)
#
#  Return the identity matrix of size n
#
    if (n > 0) diag(rep(1,n)) else NULL

zero <- function(n,m)
#
#  Return an n x m matrix of 0's
#
    if ((n > 0) & (m > 0)) matrix(0,n,m) else NULL


simplex <- function(a,A1=NULL,b1=NULL,A2=NULL,b2=NULL,A3=NULL,b3=NULL,
		maxi=FALSE, n.iter=n+2*m, eps=1e-10)
#
#   This function calculates the solution to a linear programming
#   problem using the tableau simplex method.  The constraints are
#   given by the matrices A1, A2, A3 and the vectors b1, b2 and b3
#   such that A1%*%x <= b1, A2%*%x >= b2 and A3%*%x = b3.  The 2-phase
#   Simplex method is used.
#
{
    call <- match.call()
    if (!is.null(A1))
        if (is.matrix(A1))
            m1 <- nrow(A1)
        else 	m1 <- 1
    else 	m1 <- 0
    if (!is.null(A2))
        if (is.matrix(A2))
            m2 <- nrow(A2)
        else 	m2 <- 1
    else 	m2 <- 0
    if (!is.null(A3))
        if (is.matrix(A3))
            m3 <- nrow(A3)
        else 	m3 <- 1
    else 	m3 <- 0
    m <- m1+m2+m3
    n <- length(a)
    a.o <- a
    if (maxi) a <- -a
    if (m2+m3 == 0)
# If there are no >= or = constraints then the origin is a feasible
# solution, and so only the second phase is required.
        out <- simplex1(c(a,rep(0,m1)), cbind(A1,iden(m1)), b1,
                        c(rep(0,m1),b1), n+(1L:m1), eps=eps)
    else {
        if (m2 > 0)
            out1 <- simplex1(c(a,rep(0,m1+2*m2+m3)),
                             cbind(rbind(A1,A2,A3),
                                   rbind(iden(m1),zero(m2+m3,m1)),
                                   rbind(zero(m1,m2),-iden(m2),
                                         zero(m3,m2)),
                                   rbind(zero(m1,m2+m3),
                                         iden(m2+m3))),
                             c(b1,b2,b3),
                             c(rep(0,n),b1,rep(0,m2),b2,b3),
                             c(n+(1L:m1),(n+m1+m2)+(1L:(m2+m3))),
                             stage=1, n1=n+m1+m2,
                             n.iter=n.iter, eps=eps)
        else
            out1 <- simplex1(c(a,rep(0,m1+m3)),
                             cbind(rbind(A1,A3),
                                   iden(m1+m3)),
                             c(b1,b3),
                             c(rep(0,n),b1,b3),
                             n+(1L:(m1+m3)), stage=1, n1=n+m1,
                             n.iter=n.iter, eps=eps)
#  In phase 1 use 1 artificial variable for each constraint and
#  minimize the sum of the artificial variables.  This gives a
#  feasible solution to the original problem as long as all
#  artificial variables are non-basic (and hence the value of the
#  new objective function is 0).  If this is true then optimize the
#  original problem using the result as the original feasible solution.
        if (out1$val.aux > eps)
            out <- out1
        else	out <- simplex1(out1$a[1L:(n+m1+m2)],
                                out1$A[,1L:(n+m1+m2)],
                                out1$soln[out1$basic],
                                out1$soln[1L:(n+m1+m2)],
                                out1$basic,
                                val=out1$value, n.iter=n.iter, eps=eps)
    }
    if (maxi)
        out$value <- -out$value
    out$maxi <- maxi
    if (m1 > 0L)
        out$slack <- out$soln[n+(1L:m1)]
    if (m2 > 0L)
        out$surplus <- out$soln[n+m1+(1L:m2)]
    if (out$solved == -1)
        out$artificial <- out$soln[-(1L:n+m1+m2)]
    out$obj <- a.o
    names(out$obj) <- paste("x",seq_len(n),sep="")
    out$soln <- out$soln[seq_len(n)]
    names(out$soln) <- paste("x",seq_len(n),sep="")
    out$call <- call
    class(out) <- "simplex"
    out
}




simplex1 <- function(a,A,b,init,basic,val=0,stage=2, n1=N, eps=1e-10,
                     n.iter=n1)
#
#  Tableau simplex function called by the simplex routine.  This does
#  the actual calculations required in each phase of the simplex method.
#
{
    pivot <- function(tab, pr, pc) {
#  Given the position of the pivot and the tableau, complete
#  the matrix operations to swap the variables.
        pv <- tab[pr,pc]
        pcv <- tab[,pc]
        tab[-pr,]<- tab[-pr,] - (tab[-pr,pc]/pv)%o%tab[pr,]
        tab[pr,] <- tab[pr,]/(-pv)
        tab[pr,pc] <- 1/pv
        tab[-pr,pc] <- pcv[-pr]/pv
        tab
    }
    N <- ncol(A)
    M <- nrow(A)
    nonbasic <- (1L:N)[-basic]
    tableau <- cbind(b,-A[,nonbasic,drop=FALSE])
#  If in the first stage then find the artifical objective function,
#  otherwise use the original objective function.
    if (stage == 2) {
        tableau <- rbind(tableau,c(val,a[nonbasic]))
        obfun <- a[nonbasic]
    }
    else {	obfun <- apply(tableau[(M+n1-N+1):M,,drop=FALSE],2L,sum)
		tableau <- rbind(c(val,a[nonbasic]),tableau,obfun)
		obfun <- obfun[-1L]
            }
    it <- 1
    while (!all(obfun> -eps) && (it <= n.iter))
#  While the objective function can be reduced
#	Find a pivot
#	complete the matrix operations required
#	update the lists of basic and non-basic variables
    {
        pcol <- 1+order(obfun)[1L]
        if (stage == 2)
            neg <- (1L:M)[tableau[1L:M,pcol]< -eps]
        else 	neg <- 1+ (1L:M)[tableau[2:(M+1),pcol] < -eps]
        ratios <- -tableau[neg,1L]/tableau[neg,pcol]
        prow <- neg[order(ratios)[1L]]
        tableau <- pivot(tableau,prow,pcol)
        if (stage == 1) {
            temp <- basic[prow-1L]
            basic[prow-1L] <- nonbasic[pcol-1L]
            nonbasic[pcol-1L] <- temp
            obfun <- tableau[M+2L,-1L]
        }
        else {	temp <- basic[prow]
                basic[prow] <- nonbasic[pcol-1L]
                nonbasic[pcol-1L] <- temp
                obfun <- tableau[M+1L,-1L]
            }
        it <- it+1
    }
#  END of while loop
    if (stage == 1) {
        val.aux <- tableau[M+2,1L]
# If the value of the auxilliary objective function is zero but some
# of the artificial variables are basic (with value 0) then switch
# them with some nonbasic variables (which are not artificial).
        if ((val.aux < eps) && any(basic>n1)) {
            ar <- (1L:M)[basic>n1]
            for (j in seq_along(temp)) {
                prow <- 1+ar[j]
                pcol <- 1 + order(
                                  nonbasic[abs(tableau[prow,-1L])>eps])[1L]
                tableau <- pivot(tableau,prow,pcol)
                temp1 <- basic[prow-1L]
                basic[prow-1L] <- nonbasic[pcol-1L]
                nonbasic[pcol-1L] <- temp1
            }
        }
        soln <- rep(0,N)
        soln[basic] <- tableau[2:(M+1L),1L]
        val.orig <- tableau[1L,1L]
        A.out <- matrix(0,M,N)
        A.out[,basic] <- iden(M)
        A.out[,nonbasic] <- -tableau[2L:(M+1L),-1L]
        a.orig <- rep(0,N)
        a.orig[nonbasic] <- tableau[1L,-1L]
        a.aux <- rep(0,N)
        a.aux[nonbasic] <- tableau[M+2,-1L]
        list(soln=soln, solved=-1, value=val.orig, val.aux=val.aux,
             A=A.out, a=a.orig, a.aux=a.aux, basic=basic)
    }
    else {
 	soln <- rep(0,N)
        soln[basic] <- tableau[1L:M,1L]
        val <- tableau[(M+1L),1L]
        A.out <- matrix(0,M,N)
        A.out[,basic] <- iden(M)
        A.out[,nonbasic] <- tableau[1L:M,-1L]
        a.out <- rep(0,N)
        a.out[nonbasic] <- tableau[M+1L,-1L]
        if (it <= n.iter) solved <- 1L
        else solved <- 0L
        list(soln=soln, solved=solved, value=val,  A=A.out,
             a=a.out, basic=basic)
    }
}

print.simplex <- function(x, ...) {
#
#  Print the output of a simplex solution to a linear programming problem.
#
    simp.out <- x
    cat("\nLinear Programming Results\n\n")
    cl <- simp.out$call
    cat("Call : ")
    dput(cl, control=NULL)
    if (simp.out$maxi) cat("\nMaximization ")
    else cat("\nMinimization ")
    cat("Problem with Objective Function Coefficients\n")
    print(simp.out$obj)
    if (simp.out$solved == 1) {
        cat("\n\nOptimal solution has the following values\n")
        print(simp.out$soln)
        cat(paste("The optimal value of the objective ",
                  " function is ",simp.out$value,".\n",sep=""))
    }
    else if (simp.out$solved == 0) {
        cat("\n\nIteration limit exceeded without finding solution\n")
        cat("The coefficient values at termination were\n")
        print(simp.out$soln)
        cat(paste("The objective function value was ",simp.out$value,
                  ".\n",sep=""))
    }
    else cat("\nNo feasible solution could be found\n")
    invisible(x)
}


saddle <-
    function(A = NULL, u = NULL, wdist = "m", type = "simp", d = NULL, d1 = 1,
             init = rep(0.1, d), mu = rep(0.5, n), LR = FALSE, strata = NULL,
             K.adj = NULL, K2 = NULL)
#
#  Saddle point function.  Standard multinomial saddlepoints are
#  computed using nlmin whereas the more complicated conditional
#  saddlepoints for Poisson and Binary cases are done by fitting
#  a GLM to a set of responses which, in turn, are derived from a
#  linear programming problem.
#
{
    det <- function(mat) {
#  absolute value of the determinant of a matrix.
        if (any(is.na(mat))) NA
        else if (!all(is.finite(mat))) Inf
        else  abs(prod(eigen(mat,only.values = TRUE)$values))
    }
    sgn <- function(x, eps = 1e-10)
#  sign of a real number.
        if (abs(x) < eps) 0 else 2*(x > 0) - 1

    if (!is.null(A)) {
        A <- as.matrix(A)
        d <- ncol(A)
        if (length(u) != d)
            stop(gettextf("number of columns of 'A' (%d) not equal to length of 'u' (%d)",
                          d, length(u)), domain = NA)
        n <- nrow(A)
    } else if (is.null(K.adj))
        stop("either 'A' and 'u' or 'K.adj' and 'K2' must be supplied")
    if (!is.null(K.adj)) {
#  If K.adj and K2 are supplied then calculate the simple saddlepoint.
        if (is.null(d)) d <- 1
        type <- "simp"
        wdist <- "o"
        speq <- suppressWarnings(optim(init, K.adj))
        if (speq$convergence == 0) {
            ahat <- speq$par
            Khat <- K.adj(ahat)
            K2hat <- det(K2(ahat))
            gs <- 1/sqrt((2*pi)^d*K2hat)*exp(Khat)
            if (d == 1) {
                r <- sgn(ahat)*sqrt(-2*Khat)
                v <- ahat*sqrt(K2hat)
                if (LR)	Gs <- pnorm(r)+dnorm(r)*(1/r + 1/v)
                else	Gs <- pnorm(r+log(v/r)/r)
            }
            else	Gs <- NA
        }
        else gs <- Gs <- ahat <- NA
    }
    else if (wdist == "m") {
#  Calculate the standard simple saddlepoint for the multinomial case.
        type <- "simp"
        if (is.null(strata)) {
            p <- mu/sum(mu)
            para <- list(p,A,u,n)
            K <- function(al) {
                w <- para[[1L]]*exp(al%*%t(para[[2L]]))
                para[[4L]]*log(sum(w))-sum(al*para[[3L]])
            }
            speq <- suppressWarnings(optim(init, K))
            ahat <- speq$par
            w <- as.vector(p*exp(ahat%*%t(A)))
            Khat <- n*log(sum(w))-sum(ahat*u)
            sw <- sum(w)
            if (d == 1)
                K2hat <- n*(sum(w*A*A)/sw-(sum(w*A)/sw)^2)
            else {
                saw <- w %*% A
                sa2w <-