按第一列定义的区间高效计算第二列的平均值

Question

在一个数据文件里，有两列数字。我需要通过第一列的区间（比如每100个）来计算第二列的平均值。

我可以用R语言来编写这个任务的程序，但我的R代码在处理相对较大的数据文件时（有几百万行，第一列的值在1到33132539之间变化）运行得非常慢。

这里我展示了我的R代码。请问我该如何优化它，让它运行得更快？如果有其他使用perl、python、awk或shell的解决方案也非常欢迎。

提前谢谢大家。

(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.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']

在我的Macbook Pro上，这个操作花了20秒，配置是2.53Ghz和4GB内存。如果你的第二列没有任何 NA 值，你可以通过把 mean 替换成 .Internal(mean) 来提高10倍的速度。

这里是使用 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  

---

来自Matthew的更新：

在v1.8.2版本中，这个优化（把 mean 替换成 .Internal(mean)）现在是自动进行的；也就是说，普通的 DT[,mean(somecol),by=] 现在运行速度是之前的10倍。我们会努力在未来做更多这样的便利改进，这样用户就不需要知道那么多技巧就能充分利用 data.table。

Answer 2

这段内容是关于编程问题的讨论，主要涉及一些技术细节和解决方案。它可能包含了代码示例和一些常见的错误提示。对于初学者来说，理解这些内容可能会有点挑战，但只要慢慢来，逐步掌握就能明白。

如果你在学习编程，遇到不懂的地方，可以尝试查阅相关资料，或者向更有经验的人请教。记住，编程是一个不断学习和实践的过程，遇到问题是很正常的。

总之，保持耐心，多做练习，你会逐渐变得更加熟练。

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 3

其实你并不一定要设置一个输出的数据框，但如果你想这样做也是可以的。下面是我会写的代码，我保证它会很快。

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