binary.v <- function(n) { x <- 1:(2^n) mx <- max(x) digits <- floor(log2(mx)) ans <- 0:(digits-1); lx <- length(x) x <- matrix(rep(x,rep(digits, lx)),ncol=lx) x <- (x %/% 2^ans) %% 2 return(x) } allperms <- function(n,k,index=1:n){ tmp <- binary.v(n) tmp <- tmp[,colSums(tmp)==k,drop=FALSE] apply(tmp,2,function(i) index[i==1]) } rowSD <- function(x) sqrt(ncol(x)/(ncol(x)-1)*(rowMeans(x^2)-rowMeans(x)^2)) myB <- function(x,y){ require(Biobase) require(limma) design <- model.matrix(~factor(c(rep(1,ncol(x)),rep(0,ncol(y))))) eset <- new("exprSet",exprs=cbind(x,y)) fit <- lmFit(eset, design) eBayes(fit)$t[,2] } myT <- function(x,y){ tmp1 <- rowMeans(x)-rowMeans(y) tmp2 <- rowSD(x)^2/ncol(x) + rowSD(y)^2/ncol(y) return(tmp1/sqrt(tmp2)) } rank.agreement <- function(x,y,max=200){ if(length(x)!=length(y)) stop("length of x and y must be the same.") if(max>length(x)) stop("max must be at most length(x)") tmp <- as.numeric(as.factor(rank(-abs(y)))) x <- tmp[order(-abs(x))] y <- sort(tmp) ret <- vector("numeric",length=max) for(i in 1:max) ret[i] <- sum(x[1:i]<=y[i]) return(ret) }