# slr2.r
 # first probe on SimStudy problem 2
 # normality of L1 regression
 #
  y <- c(4.5,5.5,6.5,8,10,12)
  x <- (1:6)
 #
  g <- function(b) {
         res <- y - b*x
         med <- median(res)
         g <- sum(abs(res-med)) }
 # 
 # simple linear regression
  slr <- function(y,x) {
    xb <- mean(x)
    yb <- mean(y)
    sxx <- sum((x-xb)*(x-xb))
    b1 <- sum((x-xb)*(y-yb))/sxx
    b0 <- yb - b1*xb
    sse <- sum((y-b0-b1*x)*(y-b0-b1*x))
    slr <- c(b0,b1,sqrt(sse/((length(x)-2)*sxx)))   }
 #
    rslt <- slr(y,x)
    rslt
    int <- c(rslt[2]-4*rslt[3],rslt[2]+4*rslt[3]) # +/- 4 se's
    int
 # use optimize on smoother function 
    this <- optimize(g,int)
    this
 # different approach
   c2 <- .1
   h <- function(u){sqrt(c2+u*u)}
   hp <- function(u){u/sqrt(c2+u*u)}
   hpp <- function(u){c2/((c2+u*u)*sqrt(c2+u*u))}
    b <- c(1,1)
   for ( i in(1:5) ) {        # iterate on i
    e <- y - b[1] - x*b[2]
    print(e)
    print(sapply(e,h))
    print(sum(sapply(e,h)))
    hpatb <- sapply(e,hp)
    hppatb <- sapply(e,hpp)
    print(hpatb)
    print(hppatb)
    grad <- -c(sum(hpatb),sum(x*hpatb))
    print(grad)
    hess <- matrix(c(sum(hppatb),sum(x*hppatb),sum(x*hppatb),
                  sum(x*x*hppatb)),2,2)
    print(hess)
    b <- b - solve(hess,grad)
    print(c(i,b)) }
 # done
  rm(list=ls())
 q()