估计统计中位数,模式,偏度,峰度的“在线”(迭代器)algorithm?

是否有algorithm来估计一组值的中位数,模式,偏度和/或峰度,但不要求将所有的值一次存储在内存中?

我想计算一下基本的统计数据:

  • 意思是:算术平均值
  • 方差:平均偏差的平均值
  • 标准偏差:方差的平方根
  • 中值:将较大一半的数字与较小的一半分开的值
  • 模式:在集合中find的最频繁的值
  • 偏度:tl; 博士
  • 峰度:tl; 博士

计算任何这些的基本公式是小学算术,我知道他们。 有很多统计库也可以实现它们。

我的问题是我正在处理的集合中有大量数值(数十亿):使用Python,我不能仅仅制作一个包含数十亿个元素的列表或哈希值。 即使我用C语言写这个,十亿个元素的数组也不太实际。

数据没有sorting。 它是由其他进程随机,随机产生的。 每一套的大小是非常可变的,大小不会事先知道。

我已经知道如何很好地处理均值和方差,以任何顺序遍历集合中的每个值。 (实际上,就我而言,我将它们按照生成顺序排列)。下面是我正在使用的algorithm, http : //en.wikipedia.org/wiki/Algorithms_for_calculating_variance#On-line_algorithm :

  • 初始化三个variables:count,sum和sum_of_squares
  • 对于每个值:
    • 增量计数。
    • 将该值添加到总和。
    • 将值的平方加到sum_of_squares。
  • 除数和存储为variables的意思。
  • 按count计算sum_of_squares,作为variablesmean_of_squares存储。
  • 正方形的意思是,存储为square_of_mean。
  • 从mean_of_squares减去square_of_mean,作为方差存储。
  • 产出均值和方差。

这种“在线”algorithm存在缺陷(例如,由于sum_of_squares快速增长大于整数范围或浮点精度),但基本上给了我所需要的,而不必在每个集合中存储每个值。

但我不知道是否有类似的技术来估计额外的统计数据(中位数,模式,偏度,峰度)。 只要处理N值所需的内存大大小于O(N),我就可以忍受一个有偏差的估计器,甚至是一个在一定程度上降低精度的方法。

如果图书馆具有计算一个或多个“联机”操作的function,指向现有的统计图书馆也将有所帮助。

    偏度和峰度

    对于偏度和峰度的在线algorithm(沿方差线),请参见相同的wiki页面中的高时间统计的并行algorithm。

    中位数

    没有sorting的数据,中位数很难。 如果你知道,你有多less个数据点,理论上你只需要部分sorting,例如使用selectalgorithm 。 然而,这对数十亿的价值没有太大的帮助。 我build议使用频率计数,请参阅下一节。

    中位数和频率计数模式

    如果是整数,我会计算频率 ,可能会将最高和最低值切断,超出某个值,我相信这个值不再相关。 对于浮点数(或太多的整数),我可能会创build桶/间隔,然后使用与整数相同的方法。 (近似)模式和中位数计算比较容易,基于频率表。

    正态分布的随机variables

    如果它是正态分布的,我将使用总体样本均值 , 方差 , 偏度和峰度作为一个小子集的最大似然估计量。 (在线)algorithm来计算这些,你现在已经。 例如,读取几十万或数百万个数据点,直到您的估计误差变得足够小。 只要确保你从你的集合中随机select(例如,你不会通过挑选第一个100000的值来引入偏见)。 同样的方法也可以用于估计正常情况下的模式和中位数(因为样本均值是一个估计量)。

    进一步意见

    以上所有algorithm可以并行运行(包括许多sorting和selectalgorithm,例如QuickSort和QuickSelect),如果有帮助的话。

    我一直假设(除正态分布部分外)我们讨论样本矩,中位数和模式,而不是给定已知分布的理论时刻的估计量。

    一般来说,只要所有的观察值都是相同的随机variables(具有相同的分布)的实现以及矩,模和数的数据量,对数据进行抽样(即只查看一个子集)应该是相当成功的这个分布实际上存在中位数。 最后一个警告并不是无害的。 例如, Cauchy分布的均值(和所有更高的矩)不存在。 在这种情况下,“小”子集的样本均值可能大大偏离整个样本的样本均值。

    我使用这些增量/recursion均值和中值估计量,它们都使用不变的存储:

    mean += eta * (sample - mean) median += eta * sgn(sample - median) 

    其中eta是一个小的学习速率参数(例如0.001), sgn ()是返回{-1,0,1}之一的符号函数。 (如果数据是非平稳的,并且要跟踪随时间变化的数据,则使用常量eta ;否则,对于固定的数据源,您可以使用类似于eta = 1 / n的数据作为均值估计器,其中n是所见的样本数量很远…不幸的是,这似乎不适用于中位数估计器。)

    这种types的增量平均估计似乎被用在所有的地方,例如在无监督的neural network学习规则中,但是中值版本似乎远不常见,尽pipe它有好处(对exception值的鲁棒性)。 似乎在许多应用中,中值版本可以用作平均估计量的替代。

    我很想看到一个类似forms的增量模式估计…

    UPDATE

    我只是修改增量中值估计量来估计任意分位数。 一般来说,分位数函数( http://en.wikipedia.org/wiki/Quantile_function )会告诉您将数据分成两部分的值:p和1-p。 以下估计值递增:

     quantile += eta * (sgn(sample - quantile) + 2.0 * p - 1.0) 

    值p应该在[0,1]之内。 这实质上将sgn ()函数的对称输出{-1,0,1}偏移到一侧,将数据样本分割成两个不相等大小的区域(数据的小数p和1-p小于/大于分位数估计值)。 注意,对于p = 0.5,这减less到中值估计器。

    我实现了P-Squarealgorithm,用于dynamic计算分位数和直方图,而不需要在一个整洁的Python模块中存储观测数据 ,我称之为LiveStats 。 它应该很有效地解决你的问题。 该库支持你提到的每一个统计,除了模式。 我还没有find一个满意的模式估计的解决scheme。

    瑞恩,我恐怕你不是这样做的意思和方式…这几个星期前在这里出现了 。 而在线版本的其中一个优点(实际上是用了德福的方法)就是它特别准确和稳定的事实,请看这里的讨论。 其中一个优点是您不需要存储总数或总平方数

    我想不出任何联机模式和中位数的方法,这似乎需要一次考虑整个列表。 但是,也可能是一个类似的方法,而不是方差和平均值也适用于偏度和峰度…

    在这个问题中引用的维基百科文章包含了在线计算偏度和峰度的公式。

    对于模式 – 我相信 – 在网上无法做到这一点。 为什么? 假设你input的所有值除了最后一个重复前一个值之外都是不同的。 在这种情况下,您必须记住input中已经看到的所有值,以检测最后一个值是否与所看到的值重复,并使其成为最常见的值。

    对于中位数,它几乎是相同的 – 直到最后一个input,你不知道什么值将成为中位数,如果所有的input值是不同的,因为它可能在当前中位数之前或之后。 如果你知道input的长度,你可以find没有存储所有值的中间值,但是你仍然需要存储很多值(我想大概是一半),因为不好的input序列可能会在下半场可能会从中位数上半场做出任何价值。

    (请注意,我只是指确切的计算。)

    如果你有数十亿的数据点,那么你不可能需要确切的答案,而不是密切的答案。 一般来说,如果你有数十亿的数据点,产生它们的基本过程可能会服从某种统计平稳性/遍历性/混合性。 不pipe你是否期望分配合理地连续或不合理,这也许是重要的。

    在这些情况下,存在在线,低存储量,分位数估计 (中位数是0.5分位数的特殊情况)以及模式(如果不需要确切的答案)的algorithm。 这是一个积极的统计领域。

    分位数估计示例: http : //www.computer.org/portal/web/csdl/doi/10.1109/WSC.2006.323014

    模式估计的例子:Bickel DR。 连续数据模式和偏度的鲁棒估计。 计算统计和数据分析。 2002; 39:153-163。 doi:10.1016 / S0167-9473(01)00057-3。

    这些是计算统计的活跃领域。 你正在进入那些没有任何单一的最好的精确algorithm的领域,但是它们的多样性(实际上是统计估计器),它们具有不同的性质,假设和性能。 这是实验math。 关于这个问题可能有数百到数千篇论文。

    最后一个问题是你是否真的需要自己的偏度和峰度,或者更可能是其他一些可能更可靠地表征概率分布的参数(假设你有一个概率分布!)。 你期待高斯?

    你有清理/预处理数据的方法,使其大部分Gaussianish? (例如,取对数后金融交易金额往往是高斯的)。 你期望有限的标准偏差吗? 你期望胖尾巴? 你关心的是尾巴还是散装?

    大家一直在说,你不能用网上的方式来做这个模式,但事实并非如此。 这里有一篇文章描述了一个algorithm来解决耶鲁大学的Michael E. Fischer和Steven L. Salzberg于1982年发明的这个问题。 从文章:

    大多数发现algorithm使用其中一个寄存器临时存储stream中的单个项目; 这个项目是多数元素的当前候选人。 第二个寄存器是初始化为0的计数器。对于stream的每个元素,我们要求algorithm执行以下例程。 如果计数器读取0,则将当前stream元素安装为新的多数候选人(replace可能已经在寄存器中的任何其他元素)。 然后,如果当前元素与大多数候选人匹配,则增加计数器; 否则,递减计数器。 在这个循环的这一点上,如果到目前为止所看到的stream的部分具有多数元素,则该元素在候选寄存器中,并且计数器保持大于0的值。如果没有多数元素,该怎么办? 如果不在数据stream中进行第二遍处理(这在stream环境中是不可能的),在这种情况下,algorithm总是无法给出明确的答案。 它只是承诺正确识别多数元素,如果有的话。

    它也可以扩展到find更多内存的前N,但这应该解决它的模式。

    最终,如果你没有关于分布的先验参数知识,我认为你必须存储所有的值。

    这就是说,除非你正在处理某种病理情况,否则remedian(Rousseuw和Bassett 1990)可能足以满足你的目的。

    很简单,它涉及到计算中间批次的中位数。

    中位数和模式不能仅使用恒定的可用空间在线计算。 然而,因为中位数和模式无论如何比“定量”更具“描述性”,您可以通过对数据集进行抽样来估计它们。

    如果数据是长期正态分布的,那么你可以用你的平均数来估计中位数。

    您也可以使用以下技术估计中位数:为数据stream中每1,000,000个条目build立中值估计M [i],以便M [0]是前100万个条目的中值,M [1]第二百万条目的中位数等等。然后使用M [0] … M [k]的中间值作为中值估计量。 这当然可以节省空间,并且可以通过“调整”参数1,000,000来控制要使用空间的程度。 这也可以recursion地推广。

    好伙计试试这些:

    对于c ++:

     double skew(double* v, unsigned long n){ double sigma = pow(svar(v, n), 0.5); double mu = avg(v, n); double* t; t = new double[n]; for(unsigned long i = 0; i < n; ++i){ t[i] = pow((v[i] - mu)/sigma, 3); } double ret = avg(t, n); delete [] t; return ret; } double kurt(double* v, double n){ double sigma = pow(svar(v, n), 0.5); double mu = avg(v, n); double* t; t = new double[n]; for(unsigned long i = 0; i < n; ++i){ t[i] = pow( ((v[i] - mu[i]) / sigma) , 4) - 3; } double ret = avg(t, n); delete [] t; return ret; } 

    在那里你可以计算样本方差(svar)和平均值(avg),你可以将它们指向你的函数。

    另外,看看皮尔森的近似值。 在这样一个大的数据集上,它会非常相似。 3(平均数 – 中位数)/标准偏差您的中位数为max – min / 2

    对于浮动模式没有任何意义。 通常会把它们粘在一个尺寸很大的容器里(比如1/100 *(max – min))。

    我会倾向于使用桶,这可能是适应性的。 铲斗尺寸应该是您所需要的精度。 然后,当每个数据点进来,你加一个相关的桶计数。 这些应该给你简单的中位数和峰度的近似值,通过计算每个桶的值来计算它的值。

    一个问题可能是经过数十亿次操作后浮点分辨率的损失,即增加一个不会改变数值! 为了解决这个问题,如果最大桶大小超过了某个限制,那么你可以从所有计数中取一大堆数量。

     for j in range (1,M): y=np.zeros(M) # build the vector y y[0]=y0 #generate the white noise eps=npr.randn(M-1)*np.sqrt(var) #increment the y vector for k in range(1,T): y[k]=corr*y[k-1]+eps[k-1] yy[j]=y list.append(y)