relfun:

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

Usage

relfun(xv, yv, C = 36, epsilon = 1e-04, plotit = TRUE, pch = "*")

Arguments

xv
yv
C
epsilon
plotit
pch

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 (xv, yv, C = 36, epsilon = 1e-04, plotit = TRUE, pch = "*") 
{
    plotit <- as.logical(plotit)
    temp <- matrix(c(xv, yv), ncol = 2)
    temp <- elimna(temp)
    xv <- temp[, 1]
    yv <- temp[, 2]
    tx <- biloc(xv)
    ty <- biloc(yv)
    sx <- sqrt(bivar(xv))
    sy <- sqrt(bivar(yv))
    z1 <- (xv - tx)/sx + (yv - ty)/sy
    z2 <- (xv - tx)/sx - (yv - ty)/sy
    ee <- ((z1 - biloc(z1))/sqrt(bivar(z1)))^2 + ((z2 - biloc(z2))/sqrt(bivar(z2)))^2
    w <- (1 - ee/C)^2
    if (length(w[w == 0]) >= length(xv)/2) 
        warning("More than half of the w values equal zero")
    sumw <- sum(w[ee < C])
    tempx <- w * xv
    txb <- sum(tempx[ee < C])/sumw
    tempy <- w * yv
    tyb <- sum(tempy[ee < C])/sumw
    tempxy <- w * (xv - txb) * (yv - tyb)
    tempx <- w * (xv - txb)^2
    tempy <- w * (yv - tyb)^2
    sxb <- sum((tempx[ee < C]))/sumw
    syb <- sum((tempy[ee < C]))/sumw
    rb <- sum(tempxy[ee < C])/(sqrt(sxb * syb) * sumw)
    z1 <- ((xv - txb)/sqrt(sxb) + (yv - tyb)/sqrt(syb))/sqrt(2 * 
        (1 + rb))
    z2 <- ((xv - txb)/sqrt(sxb) - (yv - tyb)/sqrt(syb))/sqrt(2 * 
        (1 - rb))
    wo <- w
    ee <- z1^2 + z2^2
    w <- (1 - ee/C)^2
    sumw <- sum(w[ee < C])
    tempx <- w * xv
    txb <- sum(tempx[ee < C])/sumw
    tempy <- w * yv
    tyb <- sum(tempy[ee < C])/sumw
    tempxy <- w * (xv - txb) * (yv - tyb)
    tempx <- w * (xv - txb)^2
    tempy <- w * (yv - tyb)^2
    sxb <- sum((tempx[ee < C]))/sumw
    syb <- sum((tempy[ee < C]))/sumw
    rb <- sum(tempxy[ee < C])/(sqrt(sxb * syb) * sumw)
    z1 <- ((xv - txb)/sqrt(sxb) + (yv - tyb)/sqrt(syb))/sqrt(2 * 
        (1 + rb))
    z2 <- ((xv - txb)/sqrt(sxb) - (yv - tyb)/sqrt(syb))/sqrt(2 * 
        (1 - rb))
    iter <- 0
    while (iter <= 10) {
        iter <= iter + 1
        ee <- z1^2 + z2^2
        w <- (1 - ee/C)^2
        sumw <- sum(w[ee < C])
        tempx <- w * xv
        txb <- sum(tempx[ee < C])/sumw
        tempy <- w * yv
        tyb <- sum(tempy[ee < C])/sumw
        tempxy <- w * (xv - txb) * (yv - tyb)
        tempx <- w * (xv - txb)^2
        tempy <- w * (yv - tyb)^2
        sxb <- sum((tempx[ee < C]))/sumw
        syb <- sum((tempy[ee < C]))/sumw
        rb <- sum(tempxy[ee < C])/(sqrt(sxb * syb) * sumw)
        z1 <- ((xv - txb)/sqrt(sxb) + (yv - tyb)/sqrt(syb))/sqrt(2 * 
            (1 + rb))
        z2 <- ((xv - txb)/sqrt(sxb) - (yv - tyb)/sqrt(syb))/sqrt(2 * 
            (1 - rb))
        wo <- w
        ee <- z1^2 + z2^2
        w <- (1 - ee/C)^2
        dif <- w - wo
        crit <- sum(dif^2)/(mean(w))^2
        if (crit <= epsilon) 
            break
    }
    if (plotit) {
        em <- median(sqrt(ee))
        r1 <- em * sqrt((1 + rb)/2)
        r2 <- em * sqrt((1 - rb)/2)
        temp <- c(0:179)
        thet <- 2 * 3.141593 * temp/180
        theta1 <- r1 * cos(thet)
        theta2 <- r2 * sin(thet)
        xplot1 <- txb + (theta1 + theta2) * sqrt(sxb)
        yplot1 <- tyb + (theta1 - theta2) * sqrt(syb)
        emax <- max(sqrt(ee[ee < 7 * em^2]))
        r1 <- emax * sqrt((1 + rb)/2)
        r2 <- emax * sqrt((1 - rb)/2)
        theta1 <- r1 * cos(thet)
        theta2 <- r2 * sin(thet)
        xplot <- txb + (theta1 + theta2) * sqrt(sxb)
        yplot <- tyb + (theta1 - theta2) * sqrt(syb)
        totx <- c(xv, xplot, xplot1)
        toty <- c(yv, yplot, yplot1)
        plot(totx, toty, type = "n", xlab = "x", ylab = "y")
        points(xv, yv, pch = pch)
        points(xplot, yplot, pch = ".")
        points(xplot1, yplot1, pch = ".")
    }
    list(mest = c(txb, tyb), mvar = c(sxb, syb), mrho = rb)
  }

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