akerdmul:

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

Usage

akerdmul(x, pts = NA, hval = NA, aval = 0.5, fr = 0.8, pr = FALSE, plotit = TRUE, theta = 50, phi = 25, expand = 0.5, scale = FALSE, xlab = "X", ylab = "Y", zlab = "", ticktype = "simple")

Arguments

x
pts
hval
aval
fr
pr
plotit
theta
phi
expand
scale
xlab
ylab
zlab
ticktype

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 (x, pts = NA, hval = NA, aval = 0.5, fr = 0.8, pr = FALSE, 
    plotit = TRUE, theta = 50, phi = 25, expand = 0.5, scale = FALSE, 
    xlab = "X", ylab = "Y", zlab = "", ticktype = "simple") 
{
    library(MASS)
    library(akima)
    if (is.na(pts[1])) 
        pts <- x
    if (ncol(x) != ncol(pts)) 
        stop("Number of columns for x and pts do not match")
    if (!is.matrix(x)) 
        stop("Data should be stored in a matrix")
    fhat <- rdplot(x, pyhat = TRUE, plotit = FALSE, fr = fr)
    n <- nrow(x)
    d <- ncol(x)
    pi <- gamma(0.5)^2
    cd <- c(2, pi)
    if (d == 2) 
        A <- 1.77
    if (d == 3) 
        A <- 2.78
    if (d > 2) {
        for (j in 3:d) cd[j] <- 2 * pi * cd[j - 2]/n
    }
    if (d > 3) 
        A <- (8 * d * (d + 2) * (d + 4) * (2 * sqrt(pi))^d)/((2 * 
            d + 1) * cd[d])
    if (is.na(hval)) 
        hval <- A * (1/n)^(1/(d + 4))
    svec <- NA
    for (j in 1:d) {
        sig <- sqrt(var(x[, j]))
        temp <- idealf(x[, j])
        iqr <- (temp$qu - temp$ql)/1.34
        A <- min(c(sig, iqr))
        x[, j] <- x[, j]/A
        svec[j] <- A
    }
    hval <- hval * sqrt(mean(svec^2))
    gm <- exp(mean(log(fhat[fhat > 0])))
    alam <- (fhat/gm)^(0 - aval)
    dhat <- NA
    nn <- nrow(pts)
    for (j in 1:nn) {
        temp1 <- t(t(x) - pts[j, ])/(hval * alam)
        temp1 <- temp1^2
        temp1 <- apply(temp1, 1, FUN = "sum")
        temp <- 0.5 * (d + 2) * (1 - temp1)/cd[d]
        epan <- ifelse(temp1 < 1, temp, 0)
        dhat[j] <- mean(epan/(alam * hval)^d)
    }
    if (plotit && d == 2) {
        fitr <- dhat
        iout <- c(1:length(fitr))
        nm1 <- length(fitr) - 1
        for (i in 1:nm1) {
            ip1 <- i + 1
            for (k in ip1:length(fitr)) if (sum(x[i, ] == x[k, 
                ]) == 2) 
                iout[k] <- 0
        }
        fitr <- fitr[iout >= 1]
        mkeep <- x[iout >= 1, ]
        fit <- interp(mkeep[, 1], mkeep[, 2], fitr)
        persp(fit, theta = theta, phi = phi, xlab = xlab, ylab = ylab, 
            zlab = zlab, expand = expand, scale = scale, ticktype = ticktype)
    }
    m <- "Done"
    if (pr) 
        m <- dhat
    m
  }

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