-----------------
#   Aufgabe 1   #
-----------------

# 1a)

x = seq(from=-5,to=5,by=0.5)
n = length(x)
n
noise = rnorm(n)
y = 3 - 0.5 * x + noise           
plot(x,y)


# 1b)

FL2 = function( beta0 , beta1 )
{
  res = sum( (y-beta0-beta1*x)^2 ) 
  return(res)
}
FL1 = function( beta0 , beta1 )
{
  res = sum( abs(y-beta0-beta1*x) ) 
  return(res)
}


# 1c)

beta0 = seq(from=0,to=6,by=0.01)
beta1 = seq(from=-3,to=3,by=0.01)

n0 = length(beta0)
n1 = length(beta1)

matFL2 = matrix(0,n0,n1)
matFL1 = matrix(0,n0,n1)

for(i in 1:n0)
{
  for(j in 1:n1)
  {
    matFL2[i,j] = FL2(beta0[i],beta1[j]) 
    matFL1[i,j] = FL1(beta0[i],beta1[j]) 
  }
}

contour(beta0,beta1,matFL2)
contour(beta0,beta1,matFL2,nlevels=50)

contour(beta0,beta1,matFL1)
contour(beta0,beta1,matFL1,nlevels=50)


# 1d)

minFL2 = min(matFL2)
minFL1 = min(matFL1)
minFL2
minFL1

# find position in matrix where
# the minimum is attained:
indicesL2 = which(matFL2 == minFL2, arr.ind = TRUE)
indicesL2
indicesL2[1]
indicesL2[2]

# also:
beta0_L2min = beta0[ indicesL2[1] ]
beta0_L2min                        # sieht ok aus
beta1_L2min = beta1[ indicesL2[2] ]
beta1_L2min                        # sieht ok aus

# L1-Minimum:
indicesL1 = which(matFL1 == minFL1, arr.ind = TRUE)
beta0_L1min = beta0[ indicesL1[1] ]
beta0_L1min                      
beta1_L1min = beta1[ indicesL1[2] ]
beta1_L1min                        


# 1e)

e = rep(1,n)
ee = sum(e*e)
ex = sum(e*x)
ey = sum(e*y)
xx = sum(x*x)
xy = sum(x*y)
b = c(ey,xy)
a1 = c(ee,ex)
a2 = c(ex,xx)
A = rbind(a1,a2)


# 1f)

betas = solve(A,b)
betas
beta0_L2min 
beta1_L2min 


# 1g)

# L^1-Regression: 

betas = c(0,0)
e = rep(1,n)
for(k in 1:100)
{
  eps = abs( y - betas[1] - betas[2]*x )
  eps = pmax( eps , 10^(-12) )
  ee = sum(e*e/eps)
  ex = sum(e*x/eps)
  xx = sum(x*x/eps)
  ey = sum(e*y/eps)
  xy = sum(x*y/eps)
  b_eps = c(ey,xy)
  a1 = c(ee,ex)
  a2 = c(ex,xx)
  A_eps = rbind(a1,a2)
  betas = solve(A_eps,b_eps)
  print(betas)
}
beta0 = betas[1]
beta1 = betas[2]
beta0
beta1
beta0_L1min 
beta1_L1min


# 1h)

Sind etwa x = (x1,x2,...,xn) und y = (y1,y2,...,yn) zwei 
Vektoren der Laenge n, dann liefert pmax(x,y) ebenfalls 
wieder einen Vektor der Laenge n, 

# pmax(x,y) = ( max(x1,y1) , max(x2,y2) , ... , max(xn,yn) ) 

Dabei steht das "p" etwa fuer "parallel", es wird parallel 
oder elementweise jeweils immer das Maximum genommen. Der 
Befehl max(x,y) liefert einfach eine Zahl, keinen Vektor, 
und zwar einfach die groesste unter allen xi und yj, 

# max(x,y) = max{x1,x2,...,xn,y1,y2,...,yn} . 


# 1i)

res = lm( y ~ x )
res
solve(A,b)          # passt

summary(res)
names(res)
res$coef
res$fit

mode(res)
str(res)


# 1j)

info = "L2 - Fit in red,   L1 - Fit in blue"
plot(x,y,main=info)
lines(x,res$fit,col="red") 
lines(x,beta0_L1min+beta1_L1min*x,col="blue")