计算大数据的分位数

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. 分块读入矩阵。对于每个块，计算高于范围和低于范围的观察值数量。然后保留该范围内的所有观察结果。
   3. 当您读完最后一个块时，您现在拥有三条信息（上面的计数、下面的计数、内部的观察向量）。获取分位数的一种方法是对整个向量进行排序并找到第 n 个观测值，您可以使用上述信息来实现这一点：对范围内的观测值进行排序，并找到第 (n-count_below) 个观测值。

编辑：(3) 的示例。

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

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 )

最后，该值与实际版本略有偏差。我怀疑我被一个或一些同样愚蠢的解释所改变，但也许我错过了一些基本的东西。