9

简而言之：我如何执行a+b使得由于截断而导致的任何精度损失远离零而不是接近零？

长篇大论

我正在计算一长串浮点值的总和，以计算集合的样本均值和方差。由于Var(X) = E(X 2 ) - E(X) 2，保持所有数字的运行计数、到目前为止所有数字的总和以及迄今为止所有数字的平方和就足够了。

到现在为止还挺好。

但是，绝对需要E(X 2 ) > E(X) 2，由于浮点精度并非总是如此。在伪代码中，问题是这样的：

int count;
double sum, sumOfSquares;
...
double value = <current-value>;
double sqrVal = value*value; 

count++;
sum += value; //slightly rounded down since value is truncated to fit into sum
sumOfSquares += sqrVal; //rounded down MORE since the order-of-magnitude 
//difference between sqrVal and sumOfSquares is twice that between value and sum;

对于可变序列，这不是一个大问题——你最终会稍微低估方差，但这通常不是一个大问题。然而，对于具有非零均值的常数或几乎常数集，它可能意味着E(X 2 ) < E(X) 2，导致计算出的方差为负，这违反了使用代码的预期。

现在，我知道了 Kahan Summation，这不是一个有吸引力的解决方案。首先，它使代码容易受到优化变幻莫测的影响（取决于优化标志，代码可能会或可能不会出现这个问题），其次，问题并不是真正由于精度 - 这已经足够了 - 这是因为加法引入了系统误差趋于零。如果我能执行这条线

sumOfSquares += sqrVal;

为了确保 sqrVal 向上而不是向下舍入到 sumOfSquares 的精度，我会有一个数值上合理的解决方案。但我怎样才能做到这一点？

编辑：已完成的问题 - 为什么在标签字段的下拉列表中按 enter 无论如何都会提交问题？

4

3 回答 3

6

IEEE 提供四种舍入模式（朝向 -inf、朝向 +inf、朝向 0、最接近）。朝着 +inf 方向似乎是您想要的。C90 或 C++ 中没有标准控件。C99 添加了头文件，该头文件<fenv.h>在某些​​ C90 和 C++ 实现中也作为扩展存在。要遵守 C99 标准，您必须编写如下内容：

#include <fenv.h>
#pragma STDC FENV_ACCESS ON

int old_round_mode = fegetround();
int set_round_ok = fesetround(FE_UPWARD);
assert(set_round_ok == 0);
...
int set_round_ok = fesetround(old_round_mode);
assert(set_round_ok == 0);

众所周知，您使用的算法在数值上不稳定并且存在精度问题。最好对数据进行两次传递以提高精度。

于 2009-08-10T08:33:45.387 回答
6

还有另一种单程算法可以稍微重新安排计算。在伪代码中：

n = 0
mean = 0
M2 = 0

for x in data:
    n = n + 1
    delta = x - mean
    mean = mean + delta/n
    M2 = M2 + delta*(x - mean)  # This expression uses the new value of mean

variance_n = M2/n         # Sample variance
variance = M2/(n - 1)     # Unbiased estimate of population variance

（来源：http ://en.wikipedia.org/wiki/Algorithms_for_calculating_variance ）

对于您使用常用算法指出的问题，这似乎表现得更好。

于 2009-08-10T08:38:23.617 回答
2

如果您不担心精度，而只是担心负方差，您为什么不简单地做V(x) = Max(0, E(X^2) - E(X)^2)

于 2009-08-10T09:18:10.040 回答