r - 计算大数据的分位数

Question

我有大约 300 个文件，每个文件包含 1000 个时间序列实现(每个文件约 76 MB)。

我想从全套 300000 个实现中计算每个时间步的分位数(0.05、0.50、0.95)。

我无法将 1 个文件中的实现合并在一起，因为它会变得太大。

最有效的方法是什么？

每个矩阵都是通过运行模型生成的，但是这里有一个包含随机数的示例:

x <- matrix(rexp(10000000, rate=.1), nrow=1000)

Answer 1

至少有三种选择:

1. 你确定它必须来自全套吗？ 10% 的样本在这里应该是一个非常非常好的近似值。
2. 300k 元素对于向量来说不算大，但 300k x 100+ 列矩阵就很大了。仅将您需要的列而不是整个矩阵拉入内存(如有必要，可以在每一列上重复)。
3. 按顺序进行，可能与较小的样本结合使用，让您在正确的范围内开始。对于第 5 个百分位数，您只需要知道有多少项目高于当前猜测，有多少项目低于当前猜测。所以像:
   1. 抽取 1% 的样本，找出它的第 5 个百分位数。上下跳跃一些公差，这样您就可以确定确切的第 5 个百分位数位于该范围内。
   2. 分 block 读入矩阵。对于每个组 block ，计算超出范围和低于范围的观察值的数量。然后保留范围内的所有观察结果。
   3. 当您读完最后一 block 时，您现在拥有三个信息(上面的计数、下面的计数、其中的观察向量)。获取分位数的一种方法是对整个向量进行排序并找到第 n 个观察值，您可以使用上述信息来做到这一点:对范围内的观察值进行排序，并找到第 (n-count_below) 个。

编辑:(3)的例子。

请注意，我不是冠军算法设计者，几乎可以肯定有人为此设计了更好的算法。此外，这种实现方式并不是特别有效。如果速度对您很重要，请考虑 Rcpp，或者为此更优化的 R。制作一堆列表然后从中提取值并不是那么聪明，但以这种方式制作原型(prototype)很容易，所以我采用了它。

library(plyr)

set.seed(1)

# -- Configuration -- #
desiredQuantile <- .25

# -- Generate sample data -- #

# Use some algorithm (sampling, iteration, or something else to come up with a range you're sure the true value lies within)
guessedrange <- c( .2, .3 )
# Group the observations to correspond to the OP's files
dat <- data.frame( group = rep( seq(100), each=100 ), value = runif(10000) )

# -- Apply the algorithm -- #

# Count the number above/below and return the values within the range, by group
res <- dlply( dat, .( group ), function( x, guessedrange ) {
  above <- x$value > guessedrange[2]
  below <- x$value < guessedrange[1]
  list(
    aboveCount  = sum( above ),
    belowCount = sum( below ),
    withinValues = x$value[ !above & !below ]
  )
}, guessedrange = guessedrange )
# Exract the count of values below and the values within the range
belowCount <- sum( sapply( res, function(x) x$belowCount ) )
belowCount
withinValues <- do.call( c, sapply( res, function(x) x$withinValues ) )
str(withinValues)
# Count up until we find the within value we want
desiredQuantileCount <- floor( desiredQuantile * nrow(dat) ) #! Should fix this so it averages when there's a tie
sort(withinValues)[ desiredQuantileCount - belowCount + 1 ]
# Compare to exact value
quantile( dat$value, desiredQuantile )

最终，该值与确切版本有一点偏差。我怀疑我被一个或一些同样愚蠢的解释转移了注意力，但也许我遗漏了一些基本的东西。