# (1)

set.seed(1234)
n=30
X=rnorm(n)
k=1

# nejprve definujeme Huberovu psi funkci
psi=function(x){
psi=rep(0,length(x))
for (i in 1:length(x)){ 
if (abs(x[i])<k) psi[i]=x[i]
else psi[i]=k*sign(x[i])
}
return(psi)
}

# a pote rovnici, kterou budeme resit
huber1=function(x){
huber1=sum(psi(X-x))
return(huber1)}

# Huberuv M-odhad
theta=uniroot(huber1,c(-5,5))$root
theta

# vykreslime graf dane funkce
par(mfrow=c(1,2))
ind=seq(-5,5, by=0.01)
plot(Vectorize(huber1)(ind)~ind, type="l")
abline(h=0,col=2)
abline(v=theta,col=2)

# (b) nyni hledejme M-odhad jako reseni minimalizacni ulohy
# definujme funkci rho
rho=function(x){
rho=rep(0,length(x))
for (i in 1:length(x)){ 
if (abs(x[i])<k) rho[i]=x[i]^2/2
else rho[i]=k*abs(x[i])-k^2/2
}
return(rho)
}

# a prislusnou funkci, kterou budeme minimalizovat
huber2=function(x){
huber2=sum(rho(X-x))
return(huber2)
}

# Huberuv M-odhad podruhe
theta2=optimize(huber2,c(-5,5))$min
val2=optimize(huber2,c(-5,5))$obj
theta2



ind=seq(-5,5, by=0.01)
plot(Vectorize(huber2)(ind)~ind, type="l")
abline(h=val2,col=2)
abline(v=theta2,col=2)


# porovnani obou odhadu
theta
theta2

#pomoci funkci z R-ka
library(MASS)
hubers(X,k=1,s=1)

#Studentizovany M-odhad
huber(X,k=1/sd(X))

# (2)
library(MASS)
library(psych)
library(ICSNP) 

plot(chem)
boxplot(chem)

chem1=chem             # puvodni datovy soubor
chem2=chem[-17]        # vylouceni odlehleho pozorovani 
chem3=c(chem2,2.895)   # oprava odlehleho pozorovani

c(mean(chem1), mean(chem2), mean(chem3))
c(median(chem1), median(chem2), median(chem3))
c(hubers(chem1)$mu, hubers(chem2)$mu, hubers(chem3)$mu)
c(mean(chem1,trim=0.05), mean(chem2,trim=0.05), mean(chem3,trim=0.05))
c(mean(chem1,trim=0.1), mean(chem2,trim=0.1), mean(chem3,trim=0.1))
c(winsor.mean(chem1,0.05),winsor.mean(chem2,0.05),winsor.mean(chem3,0.05))
c(winsor.mean(chem1,0.1),winsor.mean(chem2,0.1),winsor.mean(chem3,0.1))
c(hl.loc(chem1),hl.loc(chem2),hl.loc(chem3))

# vlastni funkce pro vypocet winsorizovaneho prumeru
winsorizedMean<-function(x,alpha) {
 xq<-quantile(x,c(alpha,1-alpha))
 x[x < xq[1]]<-xq[1]
 x[x > xq[2]]<-xq[2]
 return(mean(x))
}

# vlastni funkce pro vypocet Hodgesova - Lehmannova odhadu
HL=function(x){
n=length(x)
Y=rep(0,n*(n+1)/2)
k=1
for (i in 1:n){
for (j in i: n){
Y[k]=(x[i]+x[j])/2
k=k+1
}}
return(median(Y))
}