------------------
#   zu week4:    #
------------------


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

F2 = function(beta0,beta1)
{
  res = sum( (y - (beta0 + beta1*x))^2 )
  return(res)
}

F1 = function(beta0,beta1)
{
  res = sum( abs(y - (beta0 + beta1*x)) )
  return(res)
}

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

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

matF2 = matrix(0,nrow=n0,ncol=n1)
matF1 = matF2

for(i in 1:n0)
{
  for(j in 1:n1)
  {
    matF1[i,j] = F1(beta0[i],beta1[j])
    matF2[i,j] = F2(beta0[i],beta1[j])
  }
}

contour(beta0,beta1,matF2)

min1 = min(matF1)
min2 = min(matF2)

ind2 = which(matF2==min2,arr.ind=TRUE)

beta0[352] # 3.51
beta1[261] # -0.4

F2(beta0[274],beta1[254])

ind1 = which(matF1==min1,arr.ind=TRUE)

beta0[367] # 3.66
beta1[255] # -0.46

F1(beta0[269],beta1[261])
min1


e = rep(1,n)

ee = sum(e*e)
ex = sum(e*x)
xx = sum(x*x)
ey = sum(e*y)
xy = sum(x*y)

a1 = c(ee,ex)
a2 = c(ex,xx)
A = rbind( a1 , a2 )
b = c(ey,xy)

betaL2 = solve(A,b)
betaL2


beta0 = 0
beta1 = 0
res0 = rep(0,500)
res1 = rep(0,500)

for(k in 1:500)
{
  eps = abs( y - (beta0 + beta1*x) )
  eps = pmax( eps , 10^(-9) )

  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)

  a1 = c(ee,ex)
  a2 = c(ex,xx)
  A = rbind( a1 , a2 )
  b = c(ey,xy)

  betaL1 = solve(A,b)
  beta0 = betaL1[1]
  beta1 = betaL1[2]
  print( c(beta0,beta1) )
  res0[k] = beta0
  res1[k] = beta1
}

plot(res0)
plot(res1)


lm( y ~ x )