logo标签列表

通过第一列定义的间隔有效地平均第二列。

pythonperlrawkdata.table
7

数据文件中有两个数字列。我需要按第一列的间隔（例如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)
}
3
请定义“第一列的间隔（例如100）”。“100”是什么意思不清楚。它并不代表前100行。因此，也许你指的是从初始值（5380到5480）开始具有值为100的行。无论“第一列的间隔”是什么定义，用文字写下这个定义会非常有帮助。并非每个人都阅读R。 - S.Lott
2
同意S.Lott的观点。我不了解R，也看不出给定的输入和输出之间的关系。那是您从该输入中期望得到的输出吗？ - Bill Ruppert
这是一个数字区间，例如，第一列中所有值>=1且<=100的应该是1到100的区间；而值>=101和200的应该是101到200的区间，以此类推。并且，所有分类到由第一列定义的不同区间的第二列中的值应该被平均。我想要这些平均值对应到这些区间。感谢S. Lott和Bill Ruppert的回复。 - jianfeng.mao
将“meanrho.chr1 <- NULL”更改为“meanrho.chr1 <- numeric(nrow(interval.chr1))”应该显著加快速度。这被称为预分配。 - Marek
7个回答
7

您不一定需要设置一个输出数据框，但如果需要的话可以设置。以下是我编码的方式，保证速度快。

> 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变量，从而减少设置步骤（保留HTML标记）：

您可以使用以下代码来跳过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))
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";
}
3

鉴于您面临的问题规模，您需要使用 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次复制的速度比较。请注意，使用.Internal(mean)的data.table速度快了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中，这种优化（用.Internal(mean)替换mean）现在会自动进行；也就是说，常规的DT[,mean(somecol),by=]现在以10倍的速度运行。我们将来会尝试进行更多这样的便利性改变，让用户不需要了解太多技巧就能从data.table中获得最佳效果。

亲爱的Ramnath，感谢您向我介绍data.table。这非常有价值。 - jianfeng.mao
2

首先想到的是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的所有数字。由于您的情况不遵循此模式（根据评论中的额外细节），因此这是正确的函数：

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解决方案。这对我来说很有启发性。 - jianfeng.mao
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)

由于我们没有从0开始，因此出现了很多NA值。

(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()函数。 - jianfeng.mao
2

这里有一个 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;
2
在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'
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接