tests/huber2.R

In robsurvey: Robust Survey Statistics Estimation

library("robsurvey", quietly = TRUE)

# test on floating point equality
all_equal <- function(target, current, label,
    tolerance = sqrt(.Machine$double.eps), scale = NULL,
    check.attributes = FALSE)
{
    if (missing(label))
        stop("Argument 'label' is missing\n")
    res <- all.equal(target, current, tolerance, scale,
        check.attributes = check.attributes)
    if (is.character(res))
        cat(paste0(label, ": ", res, "\n"))
}

# test against MASS::hubers
if (requireNamespace("MASS", quietly = TRUE)) {
    library(MASS)

    # convergence tolerance
    TOLERANCE <- 1e-8

    # make a copy of the function MASS::hubers
    hubers_mod <- hubers
    # replace the mad by the (scaled) IQR as initial scale estimator
    body(hubers_mod)[[7]][[3]][[2]] <- substitute(s0 <- IQR(y, type = 2) *
        0.7413)

    # payroll data
    attach(workplace)
    all_equal(huber2(payroll, rep(1, length(payroll)), tol = TOLERANCE),
        hubers_mod(payroll, tol = 1e-8)$mu,
        label = "huber2: payroll data: test against MASS::hubers failed")

    # x11 data
    x11 <- c(1:10, 1e99); w11 <- rep(1, 11)
    all_equal(huber2(x11, w11, tol = TOLERANCE),
        hubers_mod(x11, tol = 1e-8)$mu,
        label = "huber2: x11 data: test against MASS::hubers failed")

    # x11 data (scaled)
    x11 <- c(1:10, 1e99) / 1e20; w11 <- rep(1, 11)
    all_equal(huber2(x11, w11, tol = TOLERANCE),
        hubers_mod(x11, tol = 1e-8)$mu,
        label = "huber2: x11 data (scale): test against MASS::hubers failed")
}