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

rm(list=ls());

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

# library(robustbase);
library(MASS)

# A. Basic illustration of mean, median and Huber's estimator 

set.seed(20220323)
# Sample size 
n <- 100;

# # A random sample from contaminated normal distribution
u <- runif(n);
x <- rnorm(n)*(u < 0.9) + 5*rnorm(n)*(u >= 0.9);
# hist(x)
hist(x, breaks=20) 
 
mean(x);
# [1] 0.04090808
median(x);
# [1] -0.08085571

# # 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); 
# $mu
# [1] -0.02661503
# 
# $s
# [1] 0.837996

# 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));

# B. Very heavily contaminated sample (breakdown point illustration)
x0 <- rnorm(51);
x0[52:100] <- 1e12;
hist(x0[1:51]);

mean(x0);
median(x0);
huber(x0);

mean(x0);
# [1] 4.9e+11
median(x0);
# [1] 1.970713
huber(x0);
# $mu
# [1] 8.635593
# 
# $s
# [1] 5.697667


# Adding one more contaminated point so that exactly half of the points are contaminated 
x1 <- x0;
x1[51] <- 1e12;
huber(x1);
median(x1);


# # C. Simulation study 
nopak <- 1000;

MEAN <- numeric(nopak);
MEDIAN <- numeric(nopak);
HUBER <- numeric(nopak);

for(i in 1:nopak){
	u <- runif(n);
# 	x <- rexp(n)*sign(u - 0.5)
	x <- rnorm(n)*(u < 0.9) + 5*rnorm(n)*(u >= 0.9);

# x <- rnorm(n);
	
	MEAN[i] <- mean(x);
	MEDIAN[i] <- median(x);
	HUBER[i] <- huber(x)$mu;
}

assess.theta <- function(theta, digits=3){
  ret <- c(min(theta), mean(theta), median(theta), max(theta), sd(theta), IQR(theta));
  names(ret) <- c("min", "mean", "median", "max", "sd", "IQR");
  return(round(ret, digits=3))
}

# print("Performance of mean:")
# assess.theta(MEAN);
# print("Performance of median:")
# assess.theta(MEDIAN);
# print("Performance of (studentized) Huber estimate:")
# assess.theta(HUBER);

res <- rbind(assess.theta(MEAN), assess.theta(MEDIAN), assess.theta(HUBER));
rownames(res) <- c("MEAN", "MEDIAN", "HUBER");
res



# # # # # # # # # # # # # # # # # 
# D. A sample from a not-symmetric distributions 
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);
# [1] 52.68367
median(x);
# [1] 45.62115
huber(x);
# $mu
# [1] 49.06906
# 
# $s
# [1] 23.02938