# mineral dataset from Maronna, Martin & Yohai
# Load the package robustbase
# help(package = ...) for descriptions

cu = c(102,  96, 265, 185, 229,  20,  49,  28, 128,  83, 126,  79, 
       116,  34, 633, 258, 264, 189,  70,  71, 121,  60,  37,  60,  
       23,  19,  35,  45,  52,  44,  24,  48,  42,  46,  99,  17,  
       33,  78,  201,  89,   4,  18,  43,  29,  26,  33,  24,  12,  
       14, 179,  68,  66, 102)

zn = c(4, 56, 2, 8, 26, 1, 9, 9, 28, 16, 8, 22, 12, 14, 140, 46, 32,
       19, 19, 19, 22, 19, 11, 17, 40, 17, 19, 27, 24 ,35, 24, 24, 
       27, 11, 10, 15, 33, 12, 6, 14, 4, 10, 13, 18, 12, 10, 10, 3, 
       10, 25, 17, 22, 19)
plot(cu,zn)
text(630,135,"15")

abline(lsfit(cu,zn))
text(600, 90, "LS", lty=1)

Mfit=mym(x=cu, y=zn, type="Huber", c=1.5, intercept=T) 
	#My own function to fit Huber M-estimates, as in Lecture 17
abline(Mfit$fit$coef)
text(600, 65, "Huber")

abline(lsfit(cu[-15],zn[-15]))
text(600, 30, "LS-15")

GMfit = myGM(x=cu, y=zn, type="Huber", c=1.5, intercept=T, K=3) 
  #My GM fitting function; on the course website
print(GMfit$wx)
abline(GMfit$coef, lty=2)
text(600, 44, "GM/Huber")

GMfit = myGM(x=cu, y=zn, type="bisquare", c=4.5, intercept=T, K=3) 
  #My GM fitting function; on the course website
abline(GMfit$coef, lty=6)

mmfit = lmrob(zn ~ cu)
abline(mmfit$coef, lty=4)
legend("topleft", legend=c("LS-15", "GM/bisquare", "MM"),lty=c(1,6,4))