###
###  MODERN STATISTICAL METHODS (NMST424)
###

rm(list=ls());

# # # # # # # # # # # # # # # # # # # # # # # # # # # 
# # # 1. Robust estimators of location 
# # # # # # # # # # # # # # # # # # # # # # # # # # # 

# library(robustbase);
library(MASS)

mu <- log(45.523);
sigma <- sqrt(2*(log(52.283) - mu));

# Sample from a log-normal distribution 
n <- 10000;
x <- exp(rnorm(n, mean=mu, sd=sigma));
hist(x, breaks=50)

#  # # # # # # # # # # # # # # # # 
# Why do the following estimators differ? 
mean(x);
median(x);

# # Huber studentized estimator with 
# MAD taken as the estimator of the scale, k=1.5 
# psi <- function(x, k=1.345) x*(abs(x) <= k) + k*sign(x)*(abs(x) > k)
huber(x); 


# MAD
mad(x)
median(abs(x - median(x)))/qnorm(0.75)

# Other estimates of scale 
sd(x);

# Interquartile range 
IQR(x);
(quantile(x,0.75)-quantile(x,0.25));
# Standardized IQR
(quantile(x,0.75)-quantile(x,0.25))/(qnorm(0.75)-qnorm(0.25));


# # # # # # # # # # # # # # # # # # # # # # # # # # # 
# # # 2. Robust regression
# # # # # # # # # # # # # # # # # # # # # # # # # # # 

# # A. Illustrating robustness 

  n <- 100;
  x <- sort(rnorm(n));
  # x <- sort(rnorm(n));
  # u <- runif(n);

  Y <- 1 + 2*x + rnorm(n);
#   Y[1:5] <- Y[1:5] + 10
#   Y[(n-4):n] <- Y[(n-4):n] + 10
  Y[(n-9):n] <- Y[(n-9):n] + 10

  plot(x,Y)

  summary(fit.lm <- lm(Y~x))
  summary(fit.lad <- rq(Y~x))
  summary(fit.huber <- rlm(Y~x))

  plot(x, Y);
  abline(fit.lm, col="black", lwd=2);
  abline(fit.lad, col="blue", lwd=2);
  abline(fit.huber, col="red", lwd=2);
  abline(a=1, b=2, lty="dotted");
  legend("bottomright", col=c("black", "blue", "red"), lwd=rep(2,3), legend=c("LS", "LAD", "Huber"))


#  Illustrating the algorithm
beta.iter <- coef(fit.lm);
n.iter <- 10;
X <- cbind(1, x);

psi <- function(x, k=1.345) x*(abs(x) <= k) + k*sign(x)*(abs(x) > k)
 
#  Manual calculation just to illustrate the iterative algorithm 
 for(i in 1:n.iter){
  res <- Y - X%*%beta.iter;
  s <- mad(res,0);
  scaled.res <- res/s;
  W <- psi(scaled.res)/scaled.res;
  beta.iter <- coef(lm(Y~x, weights=W));
  print(c(beta.iter, s))
#   beta.iter <- solve(t(X)%*%W%*%X)%*%t(X)%*%W%*%Y; 
 }