如何让我的循环在R中运行得更快?
问题描述:
我正在使用函数从多个HWE卡方检验中获取p值。我正在遍历一个名为geno.data
的大矩阵,它是(313行x 355232列)来执行此操作。我基本上是逐行循环两列矩阵。它运行非常缓慢。我怎样才能让它更快?由于如何让我的循环在R中运行得更快?
library(genetics)
geno.data<-matrix(c("a","c"), nrow=313,ncol=355232)
Num_of_SNPs<-ncol(geno.data) /2
alleles<- vector(length = nrow(geno.data))
HWE_pvalues<-vector(length = Num_of_SNPs)
j<- 1
for (count in 1:Num_of_SNPs){
for (i in 1:nrow(geno.data)){
alleles[i]<- levels(genotype(paste(geno.data[i,c(2*j -1, 2*j)], collapse = "/")))
}
g2 <- genotype(alleles)
HWE_pvalues[count]<-HWE.chisq(g2)[3]
j = j + 2
}
答
首先,注意张贴的代码将导致索引出界外的错误,因为毕竟主循环迭代Num_of_SNPs
您j
值将ncol(geno.data)-1
和你所访问的列2*j-1
和2*j
。我假设你想要删除列2*count-1
和2*count
和j
。
矢量化对于编写快速R代码非常重要。在你的代码中,你调用paste
函数313次,每次传递长度为1的向量。一旦传递了长度为313的向量,R中的paste
就会快得多。下面是main for循环的原始矢量化内部:
# Original
get.pval1 <- function(count) {
for (i in 1:nrow(geno.data)){
alleles[i]<- levels(genotype(paste(geno.data[i,c(2*count -1, 2*count)], collapse = "/")))
}
g2 <- genotype(alleles)
HWE.chisq(g2)[3]
}
# Vectorized
get.pval2 <- function(count) {
g2 <- genotype(paste0(geno.data[,2*count-1], "/", geno.data[,2*count]))
HWE.chisq(g2)[3]
}
我们得到关于从量化20倍速度提升:
library(microbenchmark)
all.equal(get.pval1(1), get.pval2(1))
# [1] TRUE
microbenchmark(get.pval1(1), get.pval2(1))
# Unit: milliseconds
# expr min lq mean median uq max neval
# get.pval1(1) 299.24079 304.37386 323.28321 307.78947 313.97311 482.32384 100
# get.pval2(1) 14.23288 14.64717 15.80856 15.11013 16.38012 36.04724 100
随着量化代码,你的代码应该在大约177616 * 0.01580856 =2807.853秒,或约45分钟完成(相比于16小时为原始代码)。如果对你来说这还不够快,那么我建议你看一下R中的parallel
包。mcmapply
应该给你一个很好的加速,因为外部for
循环的每次迭代都是独立的。
请参阅http://*.com/questions/5963269/how-to-make-a-great-r-reproducible-example – csgillespie 2014-11-24 15:47:37
所以你正在做'选择(355232,2)'chisq测试?你碰巧认识fortran吗? – rawr 2014-11-24 15:53:30
@rawr我不知道fortran。我正在使用R包中的函数来进行卡方检验。这是特定于我的问题。 – cooldood3490 2014-11-24 16:03:42