Jason R 的答案有一个缺陷，这在 Knuth 的“计算机编程艺术”卷中进行了讨论。2. 如果您的标准偏差是平均值的一小部分，就会出现问题：E(x^2) - (E(x)^2) 的计算对浮点舍入误差非常敏感。

你甚至可以在 Python 脚本中自己尝试：

ofs = 1e9
A = [ofs+x for x in [1,-1,2,3,0,4.02,5]] 
A2 = [x*x for x in A]
(sum(A2)/len(A))-(sum(A)/len(A))**2

我得到 -128.0 作为答案，这显然在计算上无效，因为数学预测结果应该是非负的。

Knuth 引用了一种计算运行均值和标准差的方法（我不记得发明人的名字），如下所示：

 initialize:
    m = 0;
    S = 0;
    n = 0;

 for each incoming sample x:
    prev_mean = m;
    n = n + 1;
    m = m + (x-m)/n;
    S = S + (x-m)*(x-prev_mean);

然后在每一步之后， 的值m就是平均值，标准偏差可以根据您最喜欢的标准偏差定义sqrt(S/n)来计算。sqrt(S/n-1)

我在上面写的方程与 Knuth 中的方程略有不同，但它在计算上是等效的。

当我有更多的时间时，我将用 Python 编写上面的公式，并表明你会得到一个非否定的答案（希望接近正确的值）。

---

更新：在这里。

测试1.py：

import math

def stats(x):
  n = 0
  S = 0.0
  m = 0.0
  for x_i in x:
    n = n + 1
    m_prev = m
    m = m + (x_i - m) / n
    S = S + (x_i - m) * (x_i - m_prev)
  return {'mean': m, 'variance': S/n}

def naive_stats(x):
  S1 = sum(x)
  n = len(x)
  S2 = sum([x_i**2 for x_i in x])
  return {'mean': S1/n, 'variance': (S2/n - (S1/n)**2) }

x1 = [1,-1,2,3,0,4.02,5] 
x2 = [x+1e9 for x in x1]

print "naive_stats:"
print naive_stats(x1)
print naive_stats(x2)

print "stats:"
print stats(x1)
print stats(x2)

结果：

naive_stats:
{'variance': 4.0114775510204073, 'mean': 2.0028571428571427}
{'variance': -128.0, 'mean': 1000000002.0028572}
stats:
{'variance': 4.0114775510204073, 'mean': 2.0028571428571431}
{'variance': 4.0114775868357446, 'mean': 1000000002.0028571}

您会注意到仍然存在一些舍入误差，但还不错，而naive_stats只是呕吐。

---

编辑：刚刚注意到贝利撒留的评论引用了维基百科，其中确实提到了 Knuth 算法。

为实时应用找到信号的均值和标准偏差的理想方法是什么。当信号在一定时间内偏离平均值超过 3 个标准差时，我希望能够触发控制器。

在这种情况下，正确的方法通常是计算指数加权运行平均值和标准偏差。在指数加权平均值中，均值和方差的估计值偏向于最近的样本，从而为您提供最近的均值和方差的估计值$\tau$seconds，这可能是您想要的，而不是所有见过的样本的通常算术平均值。

在频域中，“指数加权运行平均值”只是一个真正的极点。在时域中实现很简单。

时域实现

设mean和meansq为信号平方的均值和均值的当前估计值。在每个周期中，使用新样本更新这些估计值x：

% update the estimate of the mean and the mean square:
mean = (1-a)*mean + a*x
meansq = (1-a)*meansq + a*(x^2)

% calculate the estimate of the variance:
var = meansq - mean^2;

% and, if you want standard deviation:
std = sqrt(var);

这里$0 < a < 1$是一个常数，它决定了移动平均线的有效长度。如何选择$a$在下面的“分析”中进行了描述。

上面表示为命令式程序的内容也可以描述为信号流程图：

分析

上述算法计算$y_i = a x_i + (1-a) y_{i-1}$在哪里$x_i$是样本的输入$i$， 和$y_i$是输出（即均值的估计）。这是一个简单的单极 IIR 滤波器。采取$z$变换，我们找到传递函数

$$H(z) = \frac{a}{1-(1-a)z^{-1}}$$ .

将 IIR 滤波器压缩到它们自己的块中，该图现在看起来像这样：

要进入连续域，我们进行替换$z = e^{s T}$在哪里$T$是采样时间和$f_s = 1/T$是采样率。求解$1-(1-a)e^{-sT}=0$，我们发现连续系统在$s = \frac{1}{T} \log (1-a)$.

选择$a$：

$$a = 1 - \exp \left\{2\pi\frac{T}{\tau}\right\}$$

参考

• Simulink 图源可以从https://gist.github.com/1942771下载

我以前在嵌入式处理应用程序中使用过的一种方法是维护感兴趣信号的总和和平方和的累加器：

$$A_{x,i} = \sum_{k=0}^{i}x[k] = A_{x,i-1} + x[i], A_{x,-1} = 0$$

$$A_{x^2,i} = \sum_{k=0}^{i}x^2[k] = A_{x^2,i-1} + x^2[i], A_{x^2,-1} = 0$$

此外，跟踪当前时刻$i$在上述等式中（即，注意您添加到累加器中的样本数量）。然后，时间的样本均值和标准差$i$是：

$$\tilde\mu = \frac{A_{x_i}}{i+1}$$

$$\tilde\sigma = \sqrt{\frac{A_{x^2_i}}{i+1} - \tilde\mu^2}$$

或者您可以使用：

$$\tilde\sigma = \sqrt{\frac{A_{x^2_i}}{i} - \tilde\mu^2}$$

取决于您喜欢哪种标准偏差估计方法。这些方程基于方差的定义：

$$\sigma^2 = \operatorname{E}(X^2) - (\operatorname{E}(X))^2$$

我过去曾成功使用过这些（尽管我只关心方差估计，而不是标准差），尽管如果你要求和，你必须小心你用来保存累加器的数字类型很长一段时间；你不想溢出。

编辑：除了上述关于溢出的评论外，应该注意的是，当在浮点算术中实现时，这不是一个数值稳健的算法，可能会导致估计统计数据中的大错误。在这种情况下，请查看 Jason S 的答案以获得更好的方法。

Jason 和 Nibot 的答案在一个重要方面有所不同：Jason 的方法计算整个信号的标准偏差和均值（因为 y = 0），而 Nibot 的方法是“运行”计算，即它比来自遥远的过去。

由于该应用程序似乎需要标准差和均值作为时间的函数，因此 Nibot 的方法可能是更合适的方法（对于这个特定的应用程序）。然而，真正棘手的部分将是正确设置时间加权部分。Nibot 的示例使用了一个简单的单极滤波器。

描述这一点的正确方法是我们得到一个估计$E[x]$通过过滤$x[n]$和估计$E[x^2]$通过过滤$x[n]^2$. 估计滤波器通常是低通滤波器。这些滤波器应缩放为在 0 Hz 时具有 0dB 增益。否则会有一个恒定的增益误差。

可以根据您对信号的了解以及估计所需的时间分辨率来选择低通滤波器。较低的截止频率和较高的阶数将导致更好的精度但更慢的响应时间。

更复杂的是，一个滤波器应用于线性域，另一个滤波器应用于平方域。平方会显着改变信号的频谱内容，因此您可能希望在平方域中使用不同的滤波器。

这是一个关于如何估计平均值、rms 和 std dev 作为时间函数的示例。

%% example
fs = 44100; n = fs; % 44.1 kHz sample rate, 1 second
% signal: white noise plus a low frequency drift at 5 Hz)
x = randn(n,1) + sin(2*pi*(0:n-1)'*5/fs);
% mean estimation filter: since we are looking for effects in the 5 Hz range we use maybe a
% 25 Hz filter, 2nd order so it's not too sluggish
[b,a] = butter(2,25*2/fs);
xmeanEst = filter(b,a,x);
% now we estimate x^2, since most frequency double we use twice the bandwidth
[b,a] = butter(2,50*2/fs);
x2Est = filter(b,a,x.^2);
% std deviation estimate
xstd = sqrt(x2Est)-xmeanEst;
% and plot it
h = plot([x, xmeanEst sqrt(x2Est) xstd]);
grid on;
legend('x','E<x>','sqrt(E<x^2>)','Std dev');
set(h(2:4),'Linewidth',2);