L1reg <- function(x, y) 
{
# Vypocet regresnej priamky minimalizaciou L1 normy odchyliek (L1 regresia)
# x ... vektor realnych hodnot vysvetlujucej premennej 
# y ... vektor realnych hodnot vysvetlovanej (pozorovanej) premennej 
#
# Poznamka: Tato ukazkova implementacia stopercentne funguje len ak y[i] su nezaporne
# pre vsetky i a ak koeficienty L1 priamky su tiez nezaporne.
# Pre seriozne vypocty odporucam funkciu rq v kniznici quantreg.

  n <- length(x)
  u <- rep(1, n)
  a <- c(0, 0, u)

  A1 <- cbind(x, u, -diag(n))
  A2 <- cbind(x, u, diag(n))
  res <- simplex(a, A1=A1, b1=y, A2=A2, b2=y)

  list(a=res$soln[1], b=res$soln[2], val=res$value)
}

# Demonstracia:

x <- 1:30
y <- 2 * x + 100 + rcauchy(30)  # Chyby maju Cauchyho rozdelenie s tazkymi chvostami

library(boot)
L1res <- L1reg(x,y)
print(L1res)

plot(x, y, pch=19)
abline(L1res$b, L1res$a, col="red")

library(quantreg)
summary(rq(y~x))
dat <- data.frame(x=x, y=y)
L2res <- lm(y~x, data=dat)

abline(L2res$coefficients[1], L2res$coefficients[2], col="blue")