如何有效地计算运行标准偏差？

Question

我有一系列数字列表,例如:

[0] (0.01, 0.01, 0.02, 0.04, 0.03)
[1] (0.00, 0.02, 0.02, 0.03, 0.02)
[2] (0.01, 0.02, 0.02, 0.03, 0.02)
     ...
[n] (0.01, 0.00, 0.01, 0.05, 0.03)

我想要做的是有效地计算所有数组元素的列表的每个索引的平均值和标准差.

为了做到这一点,我一直循环遍历数组并将列表的给定索引处的值相加.最后,我将"平均值列表"中的每个值除以n.

为了做标准偏差,我再次循环,现在我已经计算了平均值.

我想避免两次通过阵列,一次为平均值,然后一次为SD(在我有一个意思之后).

是否有一种有效的方法来计算两个值,只通过一次数组？解释语言(例如Perl或Python)或伪代码中的任何代码都可以.

Answer 1

答案是使用Welford算法,该算法在"天真方法"之后非常明确地定义:

• 维基百科:计算方差的算法

它在数值上比在其他响应中建议的两遍或在线简单平方和收集器更稳定.只有当你有很多彼此接近的值时,稳定性才真正重要,因为它们会导致浮点文献中所谓的" 灾难性消除 ".

您可能还想要在方差计算中除以样本数(N)和N-1之间的差异(平方偏差).除以N-1导致对样本的方差的无偏估计,而除以N平均低估方差(因为它没有考虑样本均值和真实均值之间的方差).

我在这个主题上写了两篇关于更多细节的博客文章,包括如何在线删除以前的值:

• 一次通过计算在线样本均值和方差
• 删除Welford在线均值和方差算法中的值

您还可以查看我的Java工具; javadoc,源和单元测试都在线:

• 的Javadoc: stats.OnlineNormalEstimator
• 资源: stats.OnlineNormalEstimator.java
• JUnit来源: test.unit.stats.OnlineNormalEstimatorTest.java
• LingPipe主页

Answer 2

基本的答案是在你去的时候积累x(称之为'sum_x1')和x ²(称之为'sum_x2')之和.那么标准差的值是:

stdev = sqrt((sum_x2 / n) - (mean * mean)) 

哪里

mean = sum_x / n

这是样本标准差; 你使用'n'而不是'n - 1'作为除数得到总体标准差.

如果处理大样本,您可能需要担心两个大数之间取差异的数值稳定性.有关更多信息,请转到其他答案(维基百科等)中的外部参考.

Answer 3

也许不是你问的问题,但是......如果你使用一个numpy数组,它会为你做有效的工作:

from numpy import array

nums = array(((0.01, 0.01, 0.02, 0.04, 0.03),
              (0.00, 0.02, 0.02, 0.03, 0.02),
              (0.01, 0.02, 0.02, 0.03, 0.02),
              (0.01, 0.00, 0.01, 0.05, 0.03)))

print nums.std(axis=1)
# [ 0.0116619   0.00979796  0.00632456  0.01788854]

print nums.mean(axis=1)
# [ 0.022  0.018  0.02   0.02 ]

顺便说一句,在这篇博客文章中有一些有趣的讨论,以及对计算方法和差异的一次通过方法的评论:

• http://lingpipe-blog.com/2009/03/19/computing-sample-mean-variance-online-one-pass/

Answer 4

以下是来自http://www.johndcook.com/standard_deviation.html的Welford算法实现的文字纯Python翻译:

https://github.com/liyanage/python-modules/blob/master/running_stats.py

class RunningStats:

    def __init__(self):
        self.n = 0
        self.old_m = 0
        self.new_m = 0
        self.old_s = 0
        self.new_s = 0

    def clear(self):
        self.n = 0

    def push(self, x):
        self.n += 1

        if self.n == 1:
            self.old_m = self.new_m = x
            self.old_s = 0
        else:
            self.new_m = self.old_m + (x - self.old_m) / self.n
            self.new_s = self.old_s + (x - self.old_m) * (x - self.new_m)

            self.old_m = self.new_m
            self.old_s = self.new_s

    def mean(self):
        return self.new_m if self.n else 0.0

    def variance(self):
        return self.new_s / (self.n - 1) if self.n > 1 else 0.0

    def standard_deviation(self):
        return math.sqrt(self.variance())

用法:

rs = RunningStats()
rs.push(17.0);
rs.push(19.0);
rs.push(24.0);

mean = rs.mean();
variance = rs.variance();
stdev = rs.standard_deviation();

Answer 5

在Python的RUNSTATS模块是只是这样的事情.从PyPI 安装runstats:

pip install runstats

Runstats摘要可以在一次数据传递中产生均值,方差,标准差,偏度和峰度.我们可以用它来创建你的"运行"版本.

from runstats import Statistics

stats = [Statistics() for num in range(len(data[0]))]

for row in data:

    for index, val in enumerate(row):
        stats[index].push(val)

    for index, stat in enumerate(stats):
        print 'Index', index, 'mean:', stat.mean()
        print 'Index', index, 'standard deviation:', stat.stddev()

统计概要基于用于计算一次通过中的标准偏差的Knuth和Welford方法,如计算机程序设计,第2卷,第7页中所述.232,第3版.这样做的好处是数值稳定和准确的结果.

免责声明:我是Python runstats模块的作者.

Answer 6

Statistics :: Descriptive是一个非常适合这些类型计算的Perl模块:

#!/usr/bin/perl

use strict; use warnings;

use Statistics::Descriptive qw( :all );

my $data = [
    [ 0.01, 0.01, 0.02, 0.04, 0.03 ],
    [ 0.00, 0.02, 0.02, 0.03, 0.02 ],
    [ 0.01, 0.02, 0.02, 0.03, 0.02 ],
    [ 0.01, 0.00, 0.01, 0.05, 0.03 ],
];

my $stat = Statistics::Descriptive::Full->new;
# You also have the option of using sparse data structures

for my $ref ( @$data ) {
    $stat->add_data( @$ref );
    printf "Running mean: %f\n", $stat->mean;
    printf "Running stdev: %f\n", $stat->standard_deviation;
}
__END__

输出:

C:\Temp> g
Running mean: 0.022000
Running stdev: 0.013038
Running mean: 0.020000
Running stdev: 0.011547
Running mean: 0.020000
Running stdev: 0.010000
Running mean: 0.020000
Running stdev: 0.012566

Answer 7

看看PDL(发音为"piddle!").

这是Perl数据语言,专为高精度数学和科学计算而设计.

这是一个使用你的数字的例子....

use strict;
use warnings;
use PDL;

my $figs = pdl [
    [0.01, 0.01, 0.02, 0.04, 0.03],
    [0.00, 0.02, 0.02, 0.03, 0.02],
    [0.01, 0.02, 0.02, 0.03, 0.02],
    [0.01, 0.00, 0.01, 0.05, 0.03],
];

my ( $mean, $prms, $median, $min, $max, $adev, $rms ) = statsover( $figs );

say "Mean scores:     ", $mean;
say "Std dev? (adev): ", $adev;
say "Std dev? (prms): ", $prms;
say "Std dev? (rms):  ", $rms;

哪个产生:

Mean scores:     [0.022 0.018 0.02 0.02]
Std dev? (adev): [0.0104 0.0072 0.004 0.016]
Std dev? (prms): [0.013038405 0.010954451 0.0070710678 0.02]
Std dev? (rms):  [0.011661904 0.009797959 0.0063245553 0.017888544]

有关statsover函数的更多信息, 请查看PDL :: Primitive.这似乎表明ADEV是"标准偏差".

然而,它可能是PRMS(Sinan的统计数据::描述性示例显示)或RMS(ars的NumPy示例显示).我猜这三个中的一个一定是对的;-)

有关更多PDL信息,请查看:

• pdl.perl.org (官方PDL页面).
• 关于PerlMonks的PDL快速参考指南
• Dobb博士关于PDL的文章
• PDL Wiki
• 维基百科的PDL条目
• PDfor的Sourceforge项目页面