################################################################################ ## R Code für Promotionskolloquium vom 15.12.05 ################################################################################ ## load ROptEst - also loads: distr, evd, distrEx, RandVar require(ROptEst) ############### ## Optimally robust estimation ############### require(MASS) par(mfrow = c(2,1)) plot(chem, ylab = "parts per million") title("Copper in wholemeal flour") qqnorm(chem) qqline(chem) options(warn = -1) ## Step 1: ## initial estimate: Kolomogorov (-Smirnov) minimum distance estimator est0 <- ksEstimator(chem, Norm()) ## ideal model: normal location and scale N1 <- NormLocationScaleFamily(mean = est0$mean, sd = est0$sd) ## Step 2: ## radius-minimax IC: 1-2 observations contaminated ## takes about 130 sec. #IC1 <- radiusMinimaxIC(L2Fam=N1, neighbor=ContNeighborhood(),risk=asMSE(), # loRad=1/sqrt(24), upRad=2/sqrt(24)) ## Step 3: ## radius-minimax estimator via one-step method #est1 <- oneStepEstimator(chem, IC1, est0) ## special case: normal location and scale ## use package "RobLox"!!! - computations are clearly faster require(RobLox) ## Step 2 + 3: ## computes radius-minimax estimator via one-step method # takes about 6.5 sec. est2 <- roblox(chem, eps.lower = 1/24, eps.upper = 2/24, returnIC = TRUE) ## estimator comparison require(RColorBrewer) mypalette <- brewer.pal(4,"Dark2") palette(mypalette) par(mfrow=c(1,1)) plot(chem[-17]) abline(h = mean(chem), col = 1, lwd = 2) abline(h = mean(chem) - sd(chem), col = 1, lwd = 2) abline(h = mean(chem) + sd(chem), col = 1, lwd = 2) abline(h = median(chem), col = 2, lwd = 2) abline(h = median(chem) - mad(chem), col = 2, lty = 2, lwd = 2) abline(h = median(chem) + mad(chem), col = 2, lty = 2, lwd = 2) abline(h = est0$mean, col = 3, lwd = 2) abline(h = est0$mean - est0$sd, col = 3, lty = 2, lwd = 2) abline(h = est0$mean + est0$sd, col = 3, lty = 2, lwd = 2) abline(h = est2$mean, col = 4, lwd = 2) abline(h = est2$mean - est2$sd, col = 4, lty = 2, lwd = 2) abline(h = est2$mean + est2$sd, col = 4, lty = 2, lwd = 2) legend(1, max(chem[-17])-0.1, legend = c(expression("mean" %+-% "sd"), expression("median" %+-% "mad"), expression("ksMD$mean" %+-% "ksMD$sd"), expression("optRob$mean" %+-% "optRob$sd")), fill = 1:4) title("Optimally Robust Estimation")