我在R中处理一个包含大量整数的向量(约1000万个整数),我需要找到每对差值小于等于500的不同整数,并绘制它们之间差异的直方图(即对于每一对,第二个数减去第一个数的差值)。
以下是非矢量化代码,可以极其缓慢地实现我的目标:
# Generate some random example data
V <- round(rnorm(100) * 1000)
# Prepare the histogram
my.hist <- rep(0, 500)
names(my.hist) <- as.character(seq(1,500))
for (x1 in V) {
for (x2 in V) {
difference = x2 - x1
if (difference > 0 && difference <= 500) {
my.hist[difference] = my.hist[difference] + 1
}
}
}
假设每个整数都是唯一的,因此difference > 0
位是可以的。这是允许的,因为我实际上不关心差值为零的情况。
以下是向量化内部循环的代码:
my.hist2 <- rep(0, 500)
names(my.hist2) <- as.character(seq(1,500))
for (x1 in V) {
differences <- V[V > x1 & V <= x1+500] - x1
difftable <- table(differences)
my.hist2[names(difftable)] = my.hist2[names(difftable)] + difftable
}
这肯定比第一个版本快。然而,即使是这个变体当V
只包含500000个元素(半百万)时也已经太慢了。
我可以按照以下方式无需任何显式循环来完成:
X <- combn(V, 2)
# X is a matrix with two rows where each column represents a pair
diffs <- abs(X[2,] - X[1,])
my.hist3 <- table(diffs[diffs <= 500])
然而,矩阵 X 将包含 10e6 * (10e6 - 1) / 2
,即约为 50,000,000,000,000 列,这可能是个问题。
那么有没有一种方法可以在不使用显式循环(太慢)或展开所有对的矩阵(太大)的情况下完成此操作?
如果你想知道我为什么需要这样做,我正在实现这个: http://biowhat.ucsd.edu/homer/chipseq/qc.html#Sequencing_Fragment_Length_Estimation
combn
将返回一个2 x 50E12的对矩阵(所有可能的组合)。 - jbaums