library("matrixStats")
set.seed(1)
x <- rnorm(1e4)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Variance estimators
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
sigma2_a <- var(x)
cat(sprintf("var(x) = %g\n", sigma2_a))
sigma2_b <- varDiff(x)
cat(sprintf("varDiff(x) = %g\n", sigma2_b))
d <- abs(sigma2_b - sigma2_a)
cat(sprintf("Absolute difference = %g\n", d))
stopifnot(d < 0.02)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Standard deviation estimators
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
sigma_a <- sd(x)
cat(sprintf("sd(x) = %g\n", sigma_a))
sigma_b <- sdDiff(x)
cat(sprintf("sdDiff(x) = %g\n", sigma_b))
d <- abs(sigma_b - sigma_a)
cat(sprintf("Absolute difference = %g\n", d))
stopifnot(d < 0.01)
# Sanity checks
stopifnot(abs(sigma2_a - sigma_a ^ 2) < 1e-9)
stopifnot(abs(sigma2_b - sigma_b ^ 2) < 1e-9)
sigma_a2 <- mad(x)
cat(sprintf("mad(x) = %g\n", sigma_a2))
sigma_b2 <- madDiff(x)
cat(sprintf("madDiff(x) = %g\n", sigma_b2))
d <- abs(sigma_b2 - sigma_a2)
cat(sprintf("Absolute difference = %g\n", d))
stopifnot(d < 0.05)
sigma_a3 <- IQR(x)
cat(sprintf("IQR(x) = %g\n", sigma_a3))
sigma_b3 <- iqrDiff(x)
cat(sprintf("iqrDiff(x) = %g\n", sigma_b3))
d <- abs(sigma_b3 - sigma_a3)
cat(sprintf("Absolute difference = %g\n", d))
stopifnot(d < 0.05)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Trimmed estimators
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
y <- x
outliers <- sample(length(x), size = 0.1 * length(x))
y[outliers] <- 100 * y[outliers]
sigma_ao <- sd(y[-outliers])
cat(sprintf("sd(y) = %g\n", sigma_ao))
sigma_bo <- sdDiff(y[-outliers])
cat(sprintf("sdDiff(y) = %g\n", sigma_bo))
d <- abs(sigma_b - sigma_a)
cat(sprintf("Absolute difference = %g\n", d))
stopifnot(d < 0.01)
sigma_bot <- sdDiff(y, trim = 0.05)
cat(sprintf("sdDiff(y, trim = 0.05) = %g\n", sigma_bot))
d <- abs(sigma_bot - sigma_a)
cat(sprintf("Absolute difference = %g\n", d))
#stopifnot(d < 1e-3)
sigma_cot <- madDiff(y, trim = 0.05)
cat(sprintf("madDiff(y, trim = 0.05) = %g\n", sigma_cot))
sigma_dot <- iqrDiff(y, trim = 0.05)
cat(sprintf("iqrDiff(y, trim = 0.05) = %g\n", sigma_dot))
fcns <- list(
varDiff = varDiff,
sdDiff = sdDiff,
madDiff = madDiff,
iqrDiff = iqrDiff
)
for (name in names(fcns)) {
cat(sprintf("%s()...\n", name))
fcn <- fcns[[name]]
for (mode in c("integer", "double")) {
cat("mode: ", mode, "", sep = "")
for (n in 0:3) {
x <- runif(n, min = -5, max = 5)
storage.mode(x) <- mode
str(x)
y <- fcn(x)
yt <- fcn(x, trim = 0.1)
str(list("non-trimmed" = y, trimmed = yt))
} # for (mode ...)
}
cat(sprintf("%s()...DONE\n", name))
} # for (name ...)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.