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

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算法，在“naive methods”之后定义非常明确：

• 维基百科：Algorithms for calculating variance

它比其他响应中建议的两遍或在线简单平方和收集器在数值上更稳定。只有当您有很多彼此接近的值时，稳定性才真正重要，因为它们会导致浮点文献中所谓的“catastrophic cancellation”。

您可能还想复习除以样本数 (N) 和方差计算中的 N-1 之间的差异（平方偏差）。除以 N-1 会导致对样本方差的无偏估计，而除以 N 平均会低估方差（因为它没有考虑样本均值与真实均值之间的方差）。

我写了两篇关于该主题的博客文章，其中包含更多详细信息，包括如何在线删除以前的值：

• Computing Sample Mean and Variance Online in One Pass
• Deleting Values in Welford’s Algorithm for Online Mean and Variance

你也可以看看我的Java实现； javadoc、源代码和单元测试都在线：

• Javadoc: stats.OnlineNormalEstimator
• Source: stats.OnlineNormalEstimator.java
• JUnit Source: test.unit.stats.OnlineNormalEstimatorTest.java
• LingPipe Home Page

Answer 2

基本答案是累积x（称为'sum_x1'）和x²（称为'sum_x2 '） 当你去时。那么标准差的值为：

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

在哪里

mean = sum_x / n

这是样本标准差；您使用“n”而不是“n - 1”作为除数来获得总体标准差。

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

Answer 3

这里是来自http://www.johndcook.com/standard_deviation.html 的 Welford 算法实现的纯 Python 翻译：

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

import math

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()

print(f'Mean: {mean}, Variance: {variance}, Std. Dev.: {stdev}')

Answer 4

也许不是您要的，但是...如果您使用 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 ]

顺便说一下，这篇博文和 cmets 中有一些关于计算均值和方差的一次性方法的有趣讨论：

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

Answer 5

Python runstats Module 仅用于此类事情。 Install runstats 来自 PyPI：

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 卷，第 2 页中所述。 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]

查看PDL::Primitive 以获取有关 statsover 函数的更多信息。这似乎表明 ADEV 是“标准偏差”。

但是它可能是 PRMS（Sinan 的 Statistics::Descriptive 示例所示）或 RMS（ars 的 NumPy 示例所示）。我想这三个中的一个一定是对的；-)

有关更多 PDL 信息，请查看：

• pdl.perl.org（官方 PDL 页面）。
• PDL quick reference guide on PerlMonks
• Dr. Dobb's article on PDL
• PDL Wiki
• Wikipedia entry for PDL
• Sourceforge project page for PDL

Answer 8

您的阵列有多大？除非它是无数个元素，否则不要担心循环遍历它两次。代码简单，易于测试。

我的偏好是使用 numpy 数组数学扩展将您的数组数组转换为 numpy 二维数组并直接获取标准差：

>>> x = [ [ 1, 2, 4, 3, 4, 5 ], [ 3, 4, 5, 6, 7, 8 ] ] * 10
>>> import numpy
>>> a = numpy.array(x)
>>> a.std(axis=0) 
array([ 1. ,  1. ,  0.5,  1.5,  1.5,  1.5])
>>> a.mean(axis=0)
array([ 2. ,  3. ,  4.5,  4.5,  5.5,  6.5])

如果这不是一个选项，而您需要纯 Python 解决方案，请继续阅读...

如果你的数组是

x = [ 
      [ 1, 2, 4, 3, 4, 5 ],
      [ 3, 4, 5, 6, 7, 8 ],
      ....
]

那么标准差就是：

d = len(x[0])
n = len(x)
sum_x = [ sum(v[i] for v in x) for i in range(d) ]
sum_x2 = [ sum(v[i]**2 for v in x) for i in range(d) ]
std_dev = [ sqrt((sx2 - sx**2)/N)  for sx, sx2 in zip(sum_x, sum_x2) ]

如果您确定只循环一次数组，则可以合并运行总和。

sum_x  = [ 0 ] * d
sum_x2 = [ 0 ] * d
for v in x:
   for i, t in enumerate(v):
   sum_x[i] += t
   sum_x2[i] += t**2

这并不像上面的列表理解解决方案那么优雅。

Answer 9

我喜欢这样表达更新：

def running_update(x, N, mu, var):
    '''
        @arg x: the current data sample
        @arg N : the number of previous samples
        @arg mu: the mean of the previous samples
        @arg var : the variance over the previous samples
        @retval (N+1, mu', var') -- updated mean, variance and count
    '''
    N = N + 1
    rho = 1.0/N
    d = x - mu
    mu += rho*d
    var += rho*((1-rho)*d**2 - var)
    return (N, mu, var)

这样单程函数看起来像这样：

def one_pass(data):
    N = 0
    mu = 0.0
    var = 0.0
    for x in data:
        N = N + 1
        rho = 1.0/N
        d = x - mu
        mu += rho*d
        var += rho*((1-rho)*d**2 - var)
        # could yield here if you want partial results
   return (N, mu, var)

请注意，这是计算样本方差 (1/N)，而不是总体方差的无偏估计（使用 1/(N-1) 归一化因子）。与其他答案不同，跟踪运行方差的变量var 不会与样本数量成比例地增长。在任何时候，它只是迄今为止看到的一组样本的方差（在获得方差时没有最终的“除以 n”）。

在一个类中它看起来像这样：

class RunningMeanVar(object):
    def __init__(self):
        self.N = 0
        self.mu = 0.0
        self.var = 0.0
    def push(self, x):
        self.N = self.N + 1
        rho = 1.0/N
        d = x-self.mu
        self.mu += rho*d
        self.var += + rho*((1-rho)*d**2-self.var)
    # reset, accessors etc. can be setup as you see fit

这也适用于加权样本：

def running_update(w, x, N, mu, var):
    '''
        @arg w: the weight of the current sample
        @arg x: the current data sample
        @arg mu: the mean of the previous N sample
        @arg var : the variance over the previous N samples
        @arg N : the number of previous samples
        @retval (N+w, mu', var') -- updated mean, variance and count
    '''
    N = N + w
    rho = w/N
    d = x - mu
    mu += rho*d
    var += rho*((1-rho)*d**2 - var)
    return (N, mu, var)

Answer 10

这是一个“单行”，分布在多行，采用函数式编程风格：

def variance(data, opt=0):
    return (lambda (m2, i, _): m2 / (opt + i - 1))(
        reduce(
            lambda (m2, i, avg), x:
            (
                m2 + (x - avg) ** 2 * i / (i + 1),
                i + 1,
                avg + (x - avg) / (i + 1)
            ),
            data,
            (0, 0, 0)))

Answer 11

n=int(raw_input("Enter no. of terms:"))

L=[]

for i in range (1,n+1):

    x=float(raw_input("Enter term:"))

    L.append(x)

sum=0

for i in range(n):

    sum=sum+L[i]

avg=sum/n

sumdev=0

for j in range(n):

    sumdev=sumdev+(L[j]-avg)**2

dev=(sumdev/n)**0.5

print "Standard deviation is", dev

Answer 12

如以下答案所述： Does pandas/scipy/numpy provide a cumulative standard deviation function? Python Pandas 模块包含一个计算运行或cumulative standard deviation 的方法。 为此，您必须将数据转换为 pandas 数据框（如果是一维数据框，则为系列），但有一些功能。

Answer 13

这是一个实际示例，说明如何使用 python 和 numpy 实现运行标准差：

a = np.arange(1, 10)
s = 0
s2 = 0
for i in range(0, len(a)):
    s += a[i]
    s2 += a[i] ** 2 
    n = (i + 1)
    m = s / n
    std = np.sqrt((s2 / n) - (m * m))
    print(std, np.std(a[:i + 1]))

这将打印出计算出的标准偏差和用 numpy 计算的检查标准偏差：

0.0 0.0
0.5 0.5
0.8164965809277263 0.816496580927726
1.118033988749895 1.118033988749895
1.4142135623730951 1.4142135623730951
1.707825127659933 1.707825127659933
2.0 2.0
2.29128784747792 2.29128784747792
2.5819888974716116 2.581988897471611

我只是使用这个帖子中描述的公式：

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

Answer 14

回应@Charlie Parker 的 2021 年问题：

我想要一个答案，我可以在 numpy.xml 中将粘贴复制到我的代码中。我的输入是一个大小为 [N, 1] 的矩阵，其中 N 是数据点的数量，我已经计算了运行平均值，并且我假设我们已经计算了运行标准/方差，如何更新我们的新一批数据。

这里我们有两个函数的实现，它采用原始均值、原始方差和原始大小以及新样本并返回组合的原始和新样本的总均值和总方差（要获得 std-dev 只需采用方差sqrt 使用**(1/2))。第一个使用 numpy，第二个 onde 使用 welford。您可以选择最适合您的情况的一种。

def mean_and_variance_update_numpy(previous_mean, previous_var, previous_size, sample_to_append):
    if type(sample_to_append) is np.matrix:
        sample_to_append = sample_to_append.A1
    else:
        sample_to_append = sample_to_append.flatten()
    sample_to_append_mean = np.mean(sample_to_append)
    sample_to_append_size = len(sample_to_append)
    total_size = previous_size+sample_to_append_size
    total_mean = (previous_mean*previous_size+sample_to_append_mean*sample_to_append_size)/total_size
    total_var = (((previous_var+(total_mean-previous_mean)**2)*previous_size)+((np.var(sample_to_append)+(sample_to_append_mean-tm)**2)*sample_to_append_size))/total_size
    return (total_mean, total_var)
    
def mean_and_variance_update_welford(previous_mean, previous_var, previous_size, sample_to_append):
    if type(sample_to_append) is np.matrix:
        sample_to_append = sample_to_append.A1
    else:
        sample_to_append = sample_to_append.flatten()
    pos = previous_size
    mean = previous_mean
    v = previous_var*previous_size
    for value in sample_to_append:
        pos += 1
        mean_next = mean + (value - mean) / pos
        v = v + (value - mean)*(value - mean_next)
        mean = mean_next
    return (mean, v/pos)

让我们检查它是否有效：

import numpy as np

def mean_and_variance_udpate_numpy:
    ...
def mean_and_variance_udpate_welford:
    ...

# making the samples and results deterministic
np.random.seed(0)

# our initial sample has 100 samples, we want to append 10
n0, n1 = 100, 10

# using np.matrix only because it was in the question, np.array is more common
s0 = np.matrix(1e3+np.random.random_sample(n0)*1e-3).T
s1 = np.matrix(1e3+np.random.random_sample(n1)*1e-3).T

# precalculating our mean and var for initial sample:
s0mean, s0var = np.mean(s0), np.var(s0)

# calculating mean and variance for s0+s1 using our numpy updater
mean_and_variance_update_numpy(s0mean, s0var, len(s0), s1)
# (1000.0004826329636, 8.24577589696613e-08)

# calculating mean and variance for s0+s1 using our welford updater
mean_and_variance_update_welford(s0mean, s0var, len(s0), s1)
# (1000.0004826329634, 8.245775896913623e-08)

# similar results, now checking with numpy's calculation over the concatenation of s0 and s1
s0s1 = np.concatenate([s0,s1])
(np.mean(s0s1), np.var(s0s1))
# (1000.0004826329638, 8.245775896917313e-08)

这里的 3 个结果更接近

# np(s0s1)        (1000.0004826329638, 8.245775896917313e-08)
# np(s0)updnp(s1) (1000.0004826329636, 8.245775896966130e-08)
# np(s0)updwf(s1) (1000.0004826329634, 8.245775896913623e-08)

可以看出结果非常相似。