我有一系列数字列表,例如:
[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)或伪代码中的任何代码都可以。
答案是使用Welford算法,该算法在“天真方法”之后非常明确地定义:
它在数值上比在其他响应中建议的两遍或在线简单平方和收集器更稳定。当你有很多彼此接近的值时,稳定性才真正重要,因为它们会导致浮点文献中所谓的“catastrophic cancellation”。
您可能还想要在方差计算中除以样本数(N)和N-1之间的差异(平方偏差)。除以N-1导致对样本的方差的无偏估计,而将N除以平均低估方差(因为它没有考虑样本均值和真实均值之间的方差)。
我在这个主题上写了两篇关于更多细节的博客文章,包括如何在线删除以前的值:
您还可以查看我的Java工具; javadoc,源和单元测试都在线:
您可以查看关于Standard Deviation的维基百科文章,特别是关于快速计算方法的部分。
还有一篇我发现使用Python的文章,您应该能够使用其中的代码而不需要太多更改:Subliminal Messages - Running Standard Deviations。
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
如下面的答案描述:Does pandas/scipy/numpy provide a cumulative standard deviation function? Python Pandas模块包含一个计算运行或cumulative standard deviation的方法。为此,您必须将数据转换为pandas数据帧(如果是1D则转换为系列),但有一些功能。
这是一个“单行”,分布在多行,功能编程风格:
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)))
我喜欢这样表达更新:
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)
基本答案是在你去的时候积累x(称之为'sum_x1')和x2(称之为'sum_x2')之和。那么标准差的值是:
stdev = sqrt((sum_x2 / n) - (mean * mean))
哪里
mean = sum_x / n
这是样本标准差;你使用'n'而不是'n - 1'作为除数得到总体标准差。
如果处理大样本,您可能需要担心两个大数之间取差异的数值稳定性。有关更多信息,请转到其他答案(维基百科等)中的外部参考。
也许不是你问的问题,但是......如果你使用一个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://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();
Python runstats Module就是出于这种情况。来自PyPI的Install 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模块的作者。
看看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信息,请查看:
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
你的阵列有多大?除非它是数以万计的元素,否则不要担心循环两次。代码简单,易于测试。
我的偏好是使用numpy数组数学扩展来将您的数组数组转换为numpy 2D数组并直接获得标准偏差:
>>> 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
这并不像上面的列表理解解决方案那样优雅。
我认为这个问题会对你有所帮助。 Standard deviation