[R] Re: Oja median

[Roger Koenker]

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