[R] Re: Oja median

Roger Koenker roger at ysidro.econ.uiuc.edu
Fri Aug 22 14:27:48 CEST 2003

```A report on my Oja Median problem:

Gabor Grothendieck suggested replacing my loopy cofactor function with:

cofactors2 <- function(x) {
x <- t(rbind(1,cbind(x,1)))
p <- nrow(x)
solve( x, (seq(p)==p)+0) * det(x)
}

Which is much better, especially if the dimension p is large.  Unfortunately,
it is still quite slow inside the apply loop of the main function, and
this suggests that it is probably necessary to recode in C or Fortran
to get something that would work well in realistic problems.  In case
anyone else would like to play with this...here is the original version
again and a timing comparison:

"mean.wilks" <- function(x){
# Computes the column means of the matrix x -- very sloooowly.
n <- dim(x)[1]
p <- dim(x)[2]
A <- t(combn(1:n,p))
X <- NULL
for(i in 1:p)
X <- cbind(X,x[A[,i],])
oja.ize <- function(v)cofactors(matrix(v,sqrt(length(v))))
A <- t(apply(X,1,oja.ize))
coef(lm(-A[,1]~A[,-1]-1))
}
"cofactors" <- function(A){
B <- rbind(1,cbind(A,1))
p <- ncol(B)
x <- rep(0,p)
for(i in 1:p)
x[i] <- ((-1)^(i+p)) *det(B[-i,-p])
return(x)
}
______________________________________

X <- matrix(rnorm(150),50)
t1 <- unix.time(mean.wilks(X)->m1)
t2 <- unix.time(mean.wilks2(X)->m2) # this uses cofactors2
t3 <- unix.time(apply(X,2,mean)->m3)
> m1
A[, -1]1    A[, -1]2    A[, -1]3
0.09205914  0.05603060 -0.24290857
> m2
A[, -1]1    A[, -1]2    A[, -1]3
0.09205914  0.05603060 -0.24290857
> m3
[1]  0.09205914  0.05603060 -0.24290857
> t1
[1] 147.70  25.33 173.41   0.00   0.00
> t2
[1] 108.97  21.87 131.23   0.00   0.00
> t3
[1] 0 0 0 0 0

NB.  A reminder that this mean version is only a testing ground for
the median version of Oja which replaces the lm() call with a median
regression call.

url:	www.econ.uiuc.edu/~roger/my.html	Roger Koenker
email	rkoenker at uiuc.edu			Department of Economics
vox: 	217-333-4558				University of Illinois
fax:   	217-244-6678				Champaign, IL 61820

```