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

Вам действительно не нужно настраивать выходной файл 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)) 

 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"; } 

Учитывая размер вашей проблемы, вам нужно использовать data.table который быстро 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.53Ghz 4GB RAM. Если в вашем втором столбце нет NA , вы можете получить 10-кратное ускорение, заменив mean на .Internal(mean) .

Ниже приведено сравнение скорости с использованием rbenchmark и 5 повторений. Обратите внимание, что data.table с 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 

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

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

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

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

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

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

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

(Множество НС, так как мы не начинали с 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 для способов генерации факторов группировки.

Вот программа 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; 

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

 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() 

EDIT : вышеприведенное решение предполагало, что числа в первом столбце – это ВСЕ числа от 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() 

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'