Эффективно усредните второй столбец по интервалам, определенным первым столбцом

Question

В файле данных есть два числовых столбца.Мне нужно вычислить среднее значение второго столбца по интервалам (например, 100) первого столбца.

Я могу запрограммировать эту задачу на R, но мой код R действительно медленный для относительно большого файла данных (миллионы строк, значение первого столбца меняется от 1 до 33132539).

Здесь я показываю свой R-код.Как я мог настроить это, чтобы быть быстрее?Другие решения на основе Perl, Python, AWK или оболочки приветствуются.

Заранее спасибо.

(1) мой файл данных (с разделителями табуляции, миллионы строк)

5380    30.07383\n
5390    30.87\n
5393    0.07383\n
5404    6\n
5428    30.07383\n
5437    1\n
5440    9\n
5443    30.07383\n
5459    6\n
5463    30.07383\n
5480    7\n
5521    30.07383\n
5538    0\n
5584    20\n
5673    30.07383\n
5720    30.07383\n
5841    3\n
5880    30.07383\n
5913    4\n
5958    30.07383\n

(2) что я хочу получить, здесь интервал = 100

intervals_of_first_columns, average_of_2nd column_by_the_interval
100, 0\n
200, 0\n
300, 20.34074\n
400, 14.90325\n
.....

(3) R код

chr1 <- 33132539 # set the limit for the interval
window <- 100 # set the size of interval

spe <- read.table("my_data_file", header=F) # read my data in
names(spe) <- c("pos", "rho") # name my data 

interval.chr1 <- data.frame(pos=seq(0, chr1, window)) # setup intervals
meanrho.chr1 <- NULL # object for the mean I want to get

# real calculation, really slow on my own data.
for(i in 1:nrow(interval.chr1)){
  count.sub<-subset(spe, chrom==1 & pos>=interval.chr1$pos[i] & pos<=interval.chr1$pos[i+1])
  meanrho.chr1[i]<-mean(count.sub$rho)
}

Answer 1

Вам не нужно настраивать выходной файл data.frame, но вы можете, если хотите.Вот как я бы его кодировал, и я гарантирую, что это будет быстро.

> dat$incrmt <- dat$V1 %/% 100
> dat
     V1       V2 incrmt
1  5380 30.07383     53
2  5390 30.87000     53
3  5393  0.07383     53
4  5404  6.00000     54
5  5428 30.07383     54
6  5437  1.00000     54
7  5440  9.00000     54
8  5443 30.07383     54
9  5459  6.00000     54
10 5463 30.07383     54
11 5480  7.00000     54
12 5521 30.07383     55
13 5538  0.00000     55
14 5584 20.00000     55
15 5673 30.07383     56
16 5720 30.07383     57
17 5841  3.00000     58
18 5880 30.07383     58
19 5913  4.00000     59
20 5958 30.07383     59

> with(dat, tapply(V2, incrmt, mean, na.rm=TRUE))
      53       54       55       56       57       58       59 
20.33922 14.90269 16.69128 30.07383 30.07383 16.53692 17.03692 

Вы могли бы выполнить еще меньше настроек (пропустите переменную incrmt с этим кодом:

    > with(dat, tapply(V2, V1 %/% 100, mean, na.rm=TRUE))
      53       54       55       56       57       58       59 
20.33922 14.90269 16.69128 30.07383 30.07383 16.53692 17.03692 

И если вы хотите, чтобы результат был доступен для чего-либо:

by100MeanV2 <- with(dat, tapply(V2, V1 %/% 100, mean, na.rm=TRUE))

Answer 2

Учитывая размер вашей проблемы, вам нужно использовать data.table, который быстро высвечивается.

require(data.table)
N = 10^6; M = 33132539
mydt = data.table(V1 = runif(N, 1, M), V2 = rpois(N, lambda = 10))
ans  = mydt[,list(avg_V2 = mean(V2)),'V1 %/% 100']

Это заняло 20 секунд на моем Macbook Pro со спецификацией 2,53 ГГц 4 ГБ ОЗУ.Если у вас нет NA во втором столбце, вы можете получить ускорение в 10 раз, заменив mean на .Internal(mean).

Вот сравнение скорости с использованием rbenchmark и 5 повторений.Обратите внимание, что data.table с .Internal(mean) в 10 раз быстрее.

test        replications   elapsed   relative 
f_dt()            5         113.752   10.30736   
f_tapply()        5         147.664   13.38021   
f_dt_internal()   5          11.036    1.00000  

---

Обновление от Матфея:

Новое в v1.8.2, эта оптимизация (замена mean на .Internal(mean)) теперь производится автоматически;то есть обычный DT[,mean(somecol),by=] теперь работает на скорости, в 10 раз большей.Мы постараемся сделать больше удобных изменений, подобных этому, в будущем, чтобы пользователям не нужно было знать столько хитростей, чтобы получить максимум от data.table.

Answer 3

use strict;
use warnings;

my $BIN_SIZE = 100;
my %freq;

while (<>){
    my ($k, $v) = split;
    my $bin = $BIN_SIZE * int($k / $BIN_SIZE);
    $freq{$bin}{n} ++;
    $freq{$bin}{sum} += $v;
}

for my $bin (sort { $a <=> $b  } keys %freq){
    my ($n, $sum) = map $freq{$bin}{$_}, qw(n sum);
    print join("\t", $bin, $n, $sum, $sum / $n), "\n";
}

Answer 4

Oneliner в Perl прост и эффективен как обычно:

perl -F\\t -lane'BEGIN{$l=33132539;$i=100;$,=", "}sub p(){print$r*$i,$s/$n if$n;$r=int($F[0]/$i);$s=$n=0}last if$F[0]>$l;p if int($F[0]/$i)!=$r;$s+=$F[1];$n++}{p'

Answer 5

Вот программа на Perl, которая делает то, что я думаю, что вы хотите.Предполагается, что строки отсортированы по первому столбцу.

#!/usr/bin/perl
use strict;
use warnings;

my $input_name       = "t.dat";
my $output_name      = "t_out.dat";
my $initial_interval = 1;

my $interval_size    = 100;
my $start_interval   = $initial_interval;
my $end_interval     = $start_interval + $interval_size;

my $interval_total   = 0;
my $interval_count   = 0;

open my $DATA, "<", $input_name  or die "$input_name: $!";
open my $AVGS, ">", $output_name or die "$output_name: $!";

my $rows_in  = 0;
my $rows_out = 0;
$| = 1;

for (<$DATA>) {
    $rows_in++;

    # progress indicator, nice for big data
    print "*" unless $rows_in % 1000;
    print "\n" unless $rows_in % 50000;

    my ($key, $value) = split /\t/;

    # handle possible missing intervals
    while ($key >= $end_interval) {

        # put your value for an empty interval here...
        my $interval_avg = "empty";

        if ($interval_count) {
            $interval_avg = $interval_total/$interval_count;
        }
        print $AVGS $start_interval,"\t", $interval_avg, "\n";
        $rows_out++;

        $interval_count = 0;
        $interval_total = 0;

        $start_interval = $end_interval;
        $end_interval   += $interval_size;
    }

    $interval_count++;
    $interval_total += $value;
}

# handle the last interval
if ($interval_count) {
    my $interval_avg = $interval_total/$interval_count;
    print $AVGS $start_interval,"\t", $interval_avg, "\n";
    $rows_out++;
}

print "\n";
print "Rows in:  $rows_in\n";
print "Rows out: $rows_out\n";

exit 0;

Answer 6

Исходя из вашего кода, я бы предположил, что это сработает полный набор данных (в зависимости от памяти вашей системы):

chr1 <- 33132539 
window <- 100 

pos <- cut(1:chr1, seq(0, chr1, window))

meanrho.chr1 <- tapply(spe$rho, INDEX = pos, FUN = mean)

Я думаю, что вам нужен коэффициент, который определяет группы интервалов для каждых 100в первом столбце (rho), а затем вы можете использовать стандартное семейство применений для получения средств в группах.

Вот данные, которые вы разместили в воспроизводимой форме.

spe <- structure(list(pos = c(5380L, 5390L, 5393L, 5404L, 5428L, 5437L, 
5440L, 5443L, 5459L, 5463L, 5480L, 5521L, 5538L, 5584L, 5673L, 
5720L, 5841L, 5880L, 5913L, 5958L), rho = c(30.07383, 30.87, 0.07383, 
6, 30.07383, 1, 9, 30.07383, 6, 30.07383, 7, 30.07383, 0, 20, 
30.07383, 30.07383, 3, 30.07383, 4, 30.07383)), .Names = c("pos", 
"rho"), row.names = c(NA, -20L), class = "data.frame")

Определите интервалы с помощью cut, мы просто хотим каждое сотое значение (но вы можете захотеть настроить детали в соответствии с вашим кодом для вашего реального набора данных).

pos.index <- cut(spe$pos, seq(0, max(spe$pos), by = 100))

Теперь передайте желаемую функцию (mean) по каждой группе.

tapply(spe$rho, INDEX = pos.index, FUN = mean)

(много NA, так как мы не начинали с 0, тогда)

(5.2e+03,5.3e+03] (5.3e+03,5.4e+03] (5.4e+03,5.5e+03] (5.5e+03,5.6e+03] (5.6e+03,5.7e+03] (5.7e+03,5.8e+03] (5.8e+03,5.9e+03] 
   20.33922          14.90269          16.69128          30.07383          30.07383          16.53692 

(При необходимости добавьте другие аргументы в FUN, например, na.rm, например:)

## tapply(spe$rho, INDEX = pos.index, FUN = mean, na.rm = TRUE)

См. ?tapply применение к группам в векторе (рваный массив) и ?cut для способов генерации факторов группировки.

Answer 7

Первое, что приходит на ум, это генератор питона, который эффективно использует память.

def cat(data_file): # cat generator
    f = open(data_file, "r")
    for line in f:
        yield line

Затем поместите некоторую логику в другую функцию (и предположим, что вы сохранили результаты в файле)

def foo(data_file, output_file):
    f = open(output_file, "w")
    cnt = 0
    suma = 0
    for line in cat(data_file):
        suma += line.split()[-1]
        cnt += 1
        if cnt%100 == 0:
            f.write("%s\t%s\n" %( cnt, suma/100.0)
            suma = 0
    f.close()

РЕДАКТИРОВАТЬ : Приведенное выше решение предполагало, что числа в первом столбце - это ВСЕ числа от 1 до N. Поскольку ваш случай не следует этому шаблону (из дополнительных подробностей в комментариях), здесь правильная функция:

def foo_for_your_case(data_file, output_file):
    f = open(output_file, "w")
    interval = 100
    suma = 0.0
    cnt = 0 # keep track of number of elements in the interval

    for line in cat(data_file):
        spl = line.split()

        while int(spl[0]) > interval:
            if cnt > 0 : f.write("%s\t%s\n" %( interval, suma/cnt)
            else: f.write("%s\t0\n" %( interval )
            interval += 100   
            suma = 0.0
            cnt = 0

        suma += float(spl[-1])
        cnt += 1

    f.close()