# (1)

set.seed(1234)
n=30
X=rnorm(n)
#X=rcauchy(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)
}

curve(psi,-5,5) #graf funkce 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)
}

# curve(rho,-5,5) #graf funkce 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)

##########################################################################
# skorova funkce Cauchyho rozdeleni
####################################################################
# nejprve definujeme psi funkci
psi=function(x){
psi=rep(0,length(x))
for (i in 1:length(x)){
psi[i]=2*x[i]/(1+x[i]^2)
}
return(psi)
}

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

# prislusny M-odhad
theta=uniroot(cauchy,c(-5,5))$root
theta

ind=seq(-5,5, by=0.01)
plot(Vectorize(cauchy)(ind)~ind, type="l")
abline(h=0,col=2)
abline(v=theta,col=2)

##########################################################################
# Tukeyho biweight
####################################################################
k=10

# nejprve definujeme 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]*(1-(x[i]/k)^2)
else psi[i]=0
}
return(psi)
}

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

# prislusny M-odhad
theta=uniroot(Tukey,c(-5,5))$root
theta

ind=seq(-15,15, by=0.01)
plot(Vectorize(Tukey)(ind)~ind, type="l")
abline(h=0,col=2)
abline(v=theta,col=2)


# (2)
library(MASS)

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,k=1,s=1)$mu, hubers(chem2,k=1,s=1)$mu, hubers(chem3,k=1,s=1)$mu)
c(hubers(chem1,s=1)$mu, hubers(chem2,s=1)$mu, hubers(chem3,s=1)$mu)
c(hubers(chem1,k=2,s=1)$mu, hubers(chem2,k=2,s=1)$mu, hubers(chem3,k=2,s=1)$mu)