Эффективно вычислить гистограмму парных разностей в большом векторе в R?

Я работаю с большим вектором целых чисел в R (около 10 миллионов целых чисел), и мне нужно найти каждую отдельную пару целых чисел из этого вектора, которые отличаются на 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

Ответ 1