Эффективно вычислить гистограмму парных разностей в большом векторе в 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
        }
    }
}

(Предположим, что каждое целое число уникально, так что разность> 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

10
задан Ryan C. Thompson 2 March 2012 в 19:41
поделиться