Ускорение функций с участием (ы) применяются - PullRequest
Новая
       8

Ускорение функций с участием (ы) применяются

0 голосов
/ 04 ноября 2018

Я профилировал свой код с помощью пакета lineprof и определил узкие места в трех функциях perm.stat.list, G.hat и emp.FDR. Похоже, что общей темой является использование (s)apply, основанное на выводе профилировщика.

Ниже приведена упрощенная версия моих функций, а также код для генерации воспроизводимого примера, включающего три функции. Я добавил комментарии, чтобы лучше объяснить, что делает каждая функция, и необходимые входные данные.

Я бы хотел значительно ускорить мой код, потому что даже с B=10 этот процесс занимает почти полчаса вычислений. Ввод занимает большую матрицу (10000 x 10000), поэтому скорость важна. В идеале я хотел бы запустить B=5000 перестановок, что также увеличивает время вычислений.

Любые советы по улучшению моего кода приветствуются. 

### Functions ###
    perm.stat.list <- function(samp.dat,N1,N2,B){
      perm.list = NULL
      for (b in 1:B){
        #Permute the row "labels", preserving information across columns
        perm.dat.tmp = samp.dat[sample(nrow(samp.dat)),]

        #Compute the permutation-based test statistics
        #Need to save each (1 x M) permutation vector into a list
        perm.list[[b]] = apply(perm.dat.tmp,2,function(y) t.test(y[1:N1],y[(N1+1):(N1+N2)])$statistic)
      }
      return(perm.list)
    }

    G.hat = function(perm.mat,t){
      #Number of permutations
      B = nrow(perm.mat)
      #Compute an empirical distribution along each COLUMN of permutation matrix
      out = apply(perm.mat,2,function(x) sum(x>t,na.rm = TRUE))/B
      return(out)
    }

    emp.FDR <- function(t.vec,mat){
      #For each value in t.vec, apply G.hat function
      out = sapply(t.vec,function(i) sum(G.hat(mat,i),na.rm = TRUE)/max(sum(t.vec > i,na.rm = TRUE),1))
      return(out)
    }

. 

### Generate reproducible example ###

### Global variables ###
#Sample sizes (rows)
N1=3000
N2=7000
#Number of columns
M = 10000
#Number of permutations
B = 10

### Data ###
set.seed(1)
X1 = matrix(rnorm(N1*M),ncol=M)
X2 = matrix(rnorm(N2*M),ncol=M)

### Combine data in one large matrix of size (N1+N2) rows x M columns ###
samp.dat = rbind(X1,X2)

### Compute statistic for each column of samp.dat ###
t.stats = apply(samp.dat,2,
               function(x) t.test(x[1:N1],x[(N1+1):(N1+N2)])$statistic)

### Sort t.stats in decreasing order (not necessarily needed for computation) ###
t.vec = sort(t.stats,decreasing=TRUE)

### Permutation matrix based on the data ###
perm.mat = perm.stat.list(samp.dat=samp.dat,N1=N1,N2=N2,B=B)

eFDR = emp.FDR(t.vec=t.vec,mat=perm.mat)
...