oja <- function(x1, x2) 
{
# Pocita dvojrozmerny Ojov median pre body v rovine pomocou linearneho programu.
# x1 ... vektor prvych zloziek bodov
# x2 ... vektor druhych zloziek bodov

  n <- length(x1)
  m <- choose(n, 2)

  # Posun a preusporiadanie bodov

  shift.x1 <- 1 - min(x1)
  shift.x2 <- 1 - min(x2)
  x1 <- x1 + shift.x1
  x2 <- x2 + shift.x2
  rat <- x1 / x2
  x1 <- x1[order(rat)]
  x2 <- x2[order(rat)]

  # Priprava a vyriesenie prislusneho problemu LP

  h1 <- h2 <- b <- rep(0, m)
  k <- 0

  for(i in 1:(n - 1)) {
    for(j in (i + 1):n) {
      k <- k + 1
      h1[k] <- x2[i] - x2[j]
      h2[k] <- x1[j] - x1[i]
      b[k] <- x1[j] * x2[i] - x1[i] * x2[j]
      }
    }

  A1 <- cbind(h1, h2, -diag(m))
  A2 <- cbind(h1, h2, diag(m))
  a <- c(0, 0, rep(1, m))
  res <- simplex(a, A1=A1, b1=b, A2=A2, b2=b)

  list(x1=res$soln[1] - shift.x1, x2=res$soln[2] - shift.x2)
}


# Demonstracia programu:

x1 <- rnorm(20)
x2 <- rnorm(20)
z <- rnorm(20)
x1 <- x1 / z
x2 <- x2 / z 

plot(x1, x2, pch=19)
points(mean(x1), mean(x2), col="blue", pch=19)        # Tazisko bodov
points(median(x1), median(x2), col="green", pch=19)   # Pozlozkovy median

library(boot)
res <- oja(x1, x2)
points(res$x1, res$x2, col="red", pch=19)             # Ojov median