pisim:

In musto101/wilcox_R: What the package does (short line)

Usage

pisim(n = 100, q = 7, nruns = 100, alpha = 0.05, eps = 0.1, shift = 9, type = 1)

Arguments

n
q
nruns
alpha
eps
shift
type

Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (n = 100, q = 7, nruns = 100, alpha = 0.05, eps = 0.1, 
    shift = 9, type = 1) 
{
    b <- 0 * 1:q + 1
    cpicov <- 0
    npicov <- 0
    acpicov <- 0
    opicov <- 0
    val3 <- 1:nruns
    val4 <- val3
    val5 <- val3
    pilen <- matrix(0, nrow = nruns, ncol = 4)
    coef <- matrix(0, nrow = nruns, ncol = q + 1)
    corfac <- (1 + 15/n) * sqrt(n/(n - q - 1))
    corfac2 <- sqrt(n/(n - q - 1))
    for (i in 1:nruns) {
        x <- matrix(rnorm(n * q), nrow = n, ncol = q)
        if (type == 1) {
            y <- 1 + x %*% b + rnorm(n)
            xf <- rnorm(q)
            yf <- 1 + xf %*% b + rnorm(1)
        }
        if (type == 2) {
            y <- 1 + x %*% b + rt(n, df = 3)
            xf <- rnorm(q)
            yf <- 1 + xf %*% b + rt(1, df = 3)
        }
        if (type == 3) {
            y <- 1 + x %*% b + rexp(n) - 1
            xf <- rnorm(q)
            yf <- 1 + xf %*% b + rexp(1) - 1
        }
        if (type == 4) {
            y <- 1 + x %*% b + runif(n, min = -1, max = 1)
            xf <- rnorm(q)
            yf <- 1 + xf %*% b + runif(1, min = -1, max = 1)
        }
        if (type == 5) {
            err <- rnorm(n, sd = 1 + rbinom(n, 1, eps) * shift)
            y <- 1 + x %*% b + err
            xf <- rnorm(q)
            yf <- 1 + xf %*% b + rnorm(1, sd = 1 + rbinom(1, 
                1, eps) * shift)
        }
        out <- lsfit(x, y)
        fres <- out$resid
        coef[i, ] <- out$coef
        yfhat <- out$coef[1] + xf %*% out$coef[-1]
        w <- cbind(1, x)
        xtxinv <- solve(t(w) %*% w)
        xf <- c(1, xf)
        hf <- xf %*% xtxinv
        hf <- hf %*% xf
        val <- sqrt(1 + hf)
        mse <- sum(fres^2)/(n - q - 1)
        val2 <- qt(1 - alpha/2, n - q - 1) * sqrt(mse) * val
        up <- yfhat + val2
        low <- yfhat - val2
        pilen[i, 1] <- up - low
        if (low < yf && up > yf) 
            cpicov <- cpicov + 1
        val2 <- quantile(fres, c(alpha/2, 1 - alpha/2))
        val3[i] <- as.single(corfac * val2[1] * val)
        val4[i] <- as.single(corfac * val2[2] * val)
        up <- yfhat + val4[i]
        low <- yfhat + val3[i]
        pilen[i, 2] <- up - low
        if (low < yf && up > yf) 
            npicov <- npicov + 1
        val6 <- corfac2 * max(abs(val2))
        val5[i] <- val6 * val
        up <- yfhat + val5[i]
        low <- yfhat - val5[i]
        pilen[i, 3] <- up - low
        if (low < yf && up > yf) 
            acpicov <- acpicov + 1
        sres <- sort(fres)
        cc <- ceiling(n * (1 - alpha))
        rup <- sres[cc]
        rlow <- sres[1]
        olen <- rup - rlow
        if (cc < n) {
            for (j in (cc + 1):n) {
                zlen <- sres[j] - sres[j - cc + 1]
                if (zlen < olen) {
                  olen <- zlen
                  rup <- sres[j]
                  rlow <- sres[j - cc + 1]
                }
            }
        }
        up <- yfhat + corfac * val * rup
        low <- yfhat + corfac * val * rlow
        pilen[i, 4] <- up - low
        if (low < yf && up > yf) 
            opicov <- opicov + 1
    }
    pimnlen <- apply(pilen, 2, mean)
    mnbhat <- apply(coef, 2, mean)
    lcut <- mean(val3)
    hcut <- mean(val4)
    accut <- mean(val5)
    cpicov <- cpicov/nruns
    npicov <- npicov/nruns
    acpicov <- acpicov/nruns
    opicov <- opicov/nruns
    list(mnbhat = mnbhat, pimenlen = pimnlen, cpicov = cpicov, 
        npicov = npicov, acpicov = acpicov, opicov = opicov, 
        lcut = lcut, hcut = hcut, accut = accut)
  }

musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.