#---------------------------#
#   week2: R-Beispiel zu    # 
#     Regression in 1D      #
#---------------------------#



# Generieren wir etwa die Daten 
# aus dem week1.pdf:

x = seq(from=-10,to=10,by=1)
nx = length(x)
y = 1/2 * x + 5 + rnorm(nx,mean=0,sd=1)

plot(x,y)

# die Regressionskoeffizienten sind also 
# beta1 = 1/2  und  beta0 = 5 , 
# die wollen wir jetzt aus den Daten x,y 
# zurueckgewinnen:

# Wir stellen das F(beta0,beta1) auf und 
# minimieren das dann. Dazu legen wir das 
# F als eine benutzerdefinierte Funktion 
# an. Wir nehmen ein kleines f, da das 
# grosse F schon fuer FALSE verbraucht ist

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

# Schauen wir mal, ob das funktioniert:

f(5,1/2)
f(7,1/2)
f(5,0)
f(7,0)

# Plotten wir mal die Hoehenlinien von f 
# in der (beta0,beta1)-Ebene, das geht mit 
# dem contour()-Befehl:

?contour()

# wir brauchen 'Test-betas', fuer die das f 
# berechnet werden soll, nehmen wir etwa

beta0 = seq(from=0,to=10,by=0.01)         # sind dann 1000 oder 1001 Stueck
beta1 = seq(from=-5,to=5,by=0.01)         # ebenfalls 1001 Stueck

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

# sind dann insgesamt etwa 1 Million Punkte 
# in der (beta0,beta1)-Ebene; wir wollen das 
# Maximum finden, indem wir das f einfach an 
# allen Punkten berechnen:

# wir packen das Ergebnis in eine n0 x n1 - 
# Matrix z, legen wir erstmal so ein z an: 

z = matrix(0,nrow=n0,ncol=n1)
z

# 1 Mio Punkte, das koennte etwas dauern, 
# messen wir die Zeit:

t0 = Sys.time()

for(i in 1:n0)
{
  for(j in 1:n1)
  {
    z[i,j] = f(beta0[i],beta1[j])
  }
}

tend = Sys.time()
tend - t0

summary(z)
max(z)
min(z)

# wo wird das Minimum angenommen?

which.min(z)
which.min(z,arr.ind=TRUE)      # Fehlermeldung
?which.min()

# ..ok, geht so:

which( z==min(z), arr.ind=TRUE )

b0 = beta0[489]           # die 489 haengt natuerlich von den Zufallszahlen ab, 
b1 = beta1[552]           # ist immer etwas anders..
b0
b1

f(b0,b1)
min(z)                    # ok, das passt..


# wir wollten uns ja die Hoehenlinien 
# anschauen, geht jetzt so:

contour(beta0,beta1,z)
contour(beta0,beta1,z,nlevels = 100)

# das f ist offensicht wesentlich sensitiver 
# bezueglich der beta1-Variablen, waehlen wir 
# mal ein schmaleres beta1-Intervall und machen 
# dasselbe nochmal:


beta0 = seq(from=0,to=10,by=0.01)              # sind dann 1000 oder 1001 Stueck
beta1 = seq(from=-0.5,to=1.5,by=0.002)         # ebenfalls 1001 Stueck

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

z = matrix(0,nrow=n0,ncol=n1)

t0 = Sys.time()

for(i in 1:n0)
{
  for(j in 1:n1)
  {
    z[i,j] = f(beta0[i],beta1[j])
  }
}

tend = Sys.time()
tend - t0


# jetzt nochmal die Hoehenlinien:

contour(beta0, beta1, z, nlevels = 100 )                  # ok, etwas huebscher..

# ignorieren wir alle Werte ueber 50:

contour(beta0, beta1, z, nlevels = 100, zlim=c(0,50) )    # 100 Levels sind vielleicht 
                                                          # ein bischen viel..
contour(beta0, beta1, z, zlim=c(0,50) )


# Wenn wir das f also brute-force-maessig 
# einfach fuer eine Menge von Test-betas 
# evaluieren und in eine Matrix z schreiben, 
# lautet der allgemeine Code zum Finden von 
# beta0 und beta1: 

indices = which(z==min(z),arr.ind=TRUE)
i = indices[1]
j = indices[2]
beta0[i]
beta1[j]



#-----------------------------------------------------#
#  Bestimmen wir jetzt die Regressionskoeffizienten   # 
#     beta0 und beta1 mit dem  lm() - Befehl:         #
#-----------------------------------------------------# 


lm( y ~ x )          # das geht also etwas schneller...

res = lm(y~x)
res                  # liefert dasselbe..

summary(res)         # da gibt es mehr Informationen


# Infos zur Datenstruktur von res:

mode(res)            # ist eine Liste

str(res)             # ist offensichtlich eine Liste 
                     # mit 12 Elementen 

names(res)
# coefficients sollten die beta's sein, also das erste 
# Listenelement:

res[[1]]             # ja, beta0 und beta1

res$coefficients     # ist dasselbe..
res$coef             # geht auch, kann man offensichtlich 
                     # abkuerzen
res$c                # nicht mehr eindeutig, es gibt noch res$call
res$co               # jetzt wieder eindeutig

# Schauen wir uns noch die "fitted.values" an:

res$fitted.values
res$fit              # wieder dasselbe

plot(x,y)
points(x,res$fit,col="red")              # machen wir vielleicht eher 
                                         # eine Linie draus:
plot(x,y)
points(x,res$fit,type="l",col="red")     # ok, das ist also die 
                                         # Regressionsgerade