按第一列定义的间隔有效地平均第二列答案

【问题标题】：Efficiently average the second column by intervals defined by the first column按第一列定义的间隔有效地平均第二列
【发布时间】：2023-04-06 01:17:02
【问题描述】：

数据文件中有两个数字列。我需要按第一列的间隔（例如100）计算第二列的平均值。

我可以在 R 中编写这个任务，但是对于相对较大的数据文件（数百万行，第一列的值在 1 到 33132539 之间变化），我的 R 代码确实很慢。

这里我展示了我的 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)
}

【问题讨论】：

请定义“第一列的间隔（如100）”。目前尚不清楚“100”是什么意思。这不是前 100 行。因此，也许您的意思是从初始值（5380 到 5480）具有 100 值的行。无论“第一列的间隔”的定义是什么，用文字写下该定义都会有很大帮助。不是每个人都会读 R。
同意 S.Lott。我不知道 R，也看不到给定输入和给定输出之间的关系。这是您期望从该输入获得的输出吗？
这是一个数值区间，例如第一列>= 1和=101 和 200 应该是 101 到 200 的区间，然后继续。并且，应该对归类到由第一列定义的不同间隔的行中的第二列中的所有值进行平均。我想要这些平均值对应于这样的间隔。感谢 S. Lott 和 Bill Ruppert 的回复。
将meanrho.chr1 <- NULL 更改为meanrho.chr1 <- numeric(nrow(interval.chr1)) 应该会显着加快速度。这称为预分配。

标签： python perl r awk data.table

【解决方案1】：

首先想到的是python生成器，它是内存高效的。

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 的所有数字。由于您的情况不遵循这种模式（来自 cmets 中的额外细节），这里是正确的函数：

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

【讨论】：

亲爱的 hzmloth，非常感谢您为我的问题提供的 python 解决方案。这对我很有启发性。

【解决方案2】：

根据您的代码，我猜这将适用于完整的数据集（取决于您系统的内存）：

chr1 <- 33132539 
window <- 100 

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

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

我认为您需要一个因子来定义第一列 (rho) 中每 100 个区间的组，然后您可以使用标准的 apply 系列函数来获取组内的平均值。

这是您以可复制形式发布的数据。

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)

（很多 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 了解生成分组因子的方法。

【讨论】：

尊敬的 mdsumner，非常感谢您的友好解释。我从你的代码中学到了很多，比如 cut()。

【解决方案3】：

您实际上不需要设置输出 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))

【讨论】：

【解决方案4】：

这是一个 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;

【讨论】：

【解决方案5】：

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

【讨论】：

【解决方案6】：

考虑到您的问题的规模，您需要使用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']

在配备 2.53Ghz 4GB RAM 的 Macbook Pro 上，这需要 20 秒。如果您的第二列中没有任何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 获得最佳效果。

【讨论】：

亲爱的 Ramnath，感谢您向我介绍 data.table。很有价值。

【解决方案7】：

Perl 中的 Oneliner 像往常一样简单高效：

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'

【讨论】：