# R/showdown.R In RobLox: Optimally Robust Influence Curves and Estimators for Location and Scale

#### Documented in showdown

```###############################################################################
## Function to perform simulation study comparing some estimator with
## rmx estimators
###############################################################################
showdown <- function(n, M, eps, contD, seed = 123, estfun, estMean, estSd,
eps.lower = 0, eps.upper = 0.05, steps = 3L, fsCor = TRUE,
plot1 = FALSE, plot2 = FALSE, plot3 = FALSE){
stopifnot(n >= 3)
stopifnot(eps >= 0, eps <= 0.5)
if(plot1){
from <- min(-6, q.l(contD)(1e-15))
to <- max(6, q.l(contD)(1-1e-15))
curve(pnorm, from = from, to = to, lwd = 2, n = 201,
main = "Comparison: ideal vs. real", ylab = "cdf")
fun <- function(x) (1-eps)*pnorm(x) + eps*p(contD)(x)
curve(fun, from = from, to = to, add = TRUE, col = "orange",
lwd = 2, n = 201, ylab = "cdf")
legend("topleft", legend = c("ideal", "real"),
fill = c("black", "orange"))
}

set.seed(seed)
rad <- rbinom(n*M, prob = eps, size = 1)
Mid <- rnorm(n*M)
Mcont <- r(contD)(n*M)
ind <- rowSums(matrix(rad, ncol = n)) >= n/2
while(any(ind)){
M1 <- sum(ind)
cat("Samples to re-simulate:\t", M1, "\n")
rad <- rbinom(n*M1, prob = eps, size = 1)
Mid <- rnorm(n*M1)
Mcont <- r(contD)(n*M1)
ind[ind] <- rowSums(matrix(rad, ncol = n)) >= n/2
}

if(plot2){
ind <- if(M <= 20) 1:M else sample(1:M, 20)
if(plot1) dev.new()
M1 <- min(M, 20)
print(
stripplot(rep(1:M1, each = n) ~ as.vector(Mre[ind,]),
ylab = "samples", xlab = "x", pch = 20,
main = ifelse(M <= 20, "Samples", "20 randomly chosen samples"))
)
}

## ML-estimator: mean and sd
Mean <- rowMeans(Mre)
Sd <- sqrt(rowMeans((Mre-Mean)^2))
Median <- rowMedians(Mre)
## Competitor
if(missing(estfun)){
Comp <- apply(Mre, 1, estMean)
Comp <- cbind(Comp, apply(Mre, 1, estSd))
}else
Comp <- t(apply(Mre, 1, estfun))

RadMinmax <- estimate(rowRoblox(Mre, eps.lower = eps.lower,
eps.upper = eps.upper, k = steps,
fsCor = fsCor))

if(plot3){
Ergebnis1 <- list(Mean, Median, Comp[,1], RadMinmax[,1])
myCol <- brewer.pal(4, "Dark2")
if(plot1 || plot2) dev.new()
layout(matrix(c(1, 1, 1, 1, 3, 2, 2, 2, 2, 3), ncol = 2))
boxplot(Ergebnis1, col = myCol, pch = 20, main = "Location")
abline(h = 0)
boxplot(Ergebnis2, col = myCol, pch = 20, main = "Scale")
abline(h = 1)
op <- par(mar = rep(2, 4))
plot(c(0,1), c(1, 0), type = "n", axes = FALSE)
fill = myCol, ncol = 4, cex = 1.5)
on.exit(par(op))
}

## ML-estimator
MSE1.1 <- n*mean(Mean^2)
MSE2.1 <- n*mean(Median^2)
## Tukey
MSE3.1 <- n*mean(Comp[,1]^2)
empMSE <- data.frame(ML = MSE1.1, Med = MSE2.1, Competitor = MSE3.1, "rmx" = MSE4.1)
rownames(empMSE) <- "n x empMSE (loc)"
relMSE <- empMSE[1,]/empMSE[1,4]
empMSE <- rbind(empMSE, relMSE)
rownames(empMSE)[2] <- "relMSE (loc)"

## ML-estimator
MSE1.2 <- n*mean((Sd-1)^2)
## Tukey
MSE3.2 <- n*mean((Comp[,2]-1)^2)
empMSE <- rbind(empMSE, c(MSE1.2, MSE2.2, MSE3.2, MSE4.2))
rownames(empMSE)[3] <- "n x empMSE (scale)"
relMSE <- empMSE[3,]/empMSE[3,4]
empMSE <- rbind(empMSE, relMSE)
rownames(empMSE)[4] <- "relMSE (scale)"
empMSE <- rbind(empMSE, c(MSE1.1 + MSE1.2, MSE2.1 + MSE2.2, MSE3.1 + MSE3.2,
MSE4.1 + MSE4.2))
rownames(empMSE)[5] <- "n x empMSE (loc + scale)"
relMSE <- empMSE[5,]/empMSE[5,4]
empMSE <- rbind(empMSE, relMSE)
rownames(empMSE)[6] <- "relMSE (loc + scale)"

empMSE
}
```

## Try the RobLox package in your browser

Any scripts or data that you put into this service are public.

RobLox documentation built on April 6, 2019, 3:04 a.m.