library(MASS)

#(1)
n=30
p=3
opak=10000

MSE1=rep(0,opak)
MSE2=rep(0,opak)

sigma=1
theta=rep(1,p)
Sigma=sigma^2* diag(p)

set.seed(1234)
for (i in 1:opak){
X=mvrnorm(n,theta,Sigma)

X_bar=apply(X, 2, mean)                            #vyberovy prumer
theta_hat=(1-(p-2)*sigma^2/n/sum(X_bar)^2)*X_bar   #J.-S. odhad

MSE1[i]=sum((X_bar-theta)^2)
MSE2[i]=sum((theta_hat-theta)^2)
}

sum(MSE1>MSE2)/opak # podil pripadu, kdy je J.-S. odhad lepsi
sum(MSE1<MSE2)/opak # podil pripadu, kdy je prumer lepsi

mean(MSE1)   # stredni ctvercova chyba prumeru
mean(MSE2)   # stredni ctvercova chyba J.-S. odhadu
sigma^2/n*p  # teoreticka stredni ctvercova chyba prumeru



library(depth)
library(HSAUR)


#(2) 
data=heptathlon
x=data$hurdles
y=data$shot
X=t(rbind(x,y))

#(a)
plot(X[,1],X[,2])
boxplot(X[,1])
boxplot(X[,2])

hloubka=rep(0,25)
for (i in 1:25){
print(depth(X[i,],X))
hloubka[i]=depth(X[i,],X)
}

#(b)
text(X[,1],X[,2],label=round(hloubka,2),col='blue', cex=0.7,pos=2)

#(c)
isodepth(X,dpth=c(1,2,3,4,5,6,7,8,9,10),colcontours=1:10)
med1=med(X,method="Tukey")$median
med1
#(d)
med2=med(X,method="Spatial")$median
med3=med(X,method="CWmed")$median
med2
med3

plot(X[,1],X[,2],pch=4)
points(med1[1],med1[2],col=2, pch=8)
points(med2[1],med2[2],col=3, pch=8)
points(med3[1],med3[2],col=4, pch=8)

legend("topright", legend=c("Tukey", "Spatial", "marginal"),pch=8,col=2:4)


#(3)
library(aplpack)

#(a)
par(mfrow=c(1,2))
boxplot(x)
boxplot(y)

#(b)
bagplot(X)