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 [EMAIL PROTECTED] Department of Economics
vox: 217-333-4558 University of Illinois
fax: 217-244-6678 Champaign, IL 61820
______________________________________________
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
