如何有效地计算运行标准偏差？
Posted
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了如何有效地计算运行标准偏差？相关的知识，希望对你有一定的参考价值。
我有一系列数字列表，例如：
[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除以平均低估方差（因为它没有考虑样本均值和真实均值之间的方差）。
我在这个主题上写了两篇关于更多细节的博客文章，包括如何在线删除以前的值：
- 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
您可以查看关于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();