9

我有以下 Python 代码和输出：

>>> import numpy as np
>>> s = [12.40265325, -1.3362417499999921, 6.8768662500000062, 25.673127166666703, 19.733372250000002, 21.649556250000003, 7.1676752500000021, -0.85349583333329804, 23.130314250000012, 20.074925250000007, -0.29701574999999281, 17.078694250000012, 3.3652611666666985, 19.491246250000003, -0.76856974999999039, -1.8838917499999965, -6.8547018333333085, 4.5195052500000088, 5.9882702500000136, -9.5889237499999922, 13.98170916666669, -2.3662137499999929, 12.111165250000013, -6.8334957499999902, -21.379336749999993, 8.4651301666666967, 2.5094612500000082, -0.21386274999998989, 5.1226162500000072, 14.283680166666699, -4.3340977499999909, -2.7831607499999933, 8.2339832500000085, -12.841856749999991, -6.4984398333333075, -6.2773697499999912, -13.638411749999996, -15.90088974999999, -8.2505068333333043, -19.616496749999996, -4.4346607499999919, -10.056376749999991, -13.581729833333299, -8.2284047499999957, -4.5957137499999945, -5.3758427499999968, -12.254779833333302, 11.207287250000007, -12.848971749999997, -14.449801749999992, -17.247984749999993, -17.475253833333305]
>>> np.mean(s)
1.3664283380001927e-14
>>> np.std(s)
12.137473069268983
>>> (s - np.mean(s)) / np.std(s)
array([ 1.02184806, -0.11009225,  0.56658138,  2.1151954 , ...

当我在 R 中运行它时，结果不匹配：

> options(digits=16)
> s = c(12.40265325, -1.3362417499999921, 6.8768662500000062, 25.673127166666703, 19.733372250000002, 21.649556250000003, 7.1676752500000021, -0.85349583333329804, 23.130314250000012, 20.074925250000007, -0.29701574999999281, 17.078694250000012, 3.3652611666666985, 19.491246250000003, -0.76856974999999039, -1.8838917499999965, -6.8547018333333085, 4.5195052500000088, 5.9882702500000136, -9.5889237499999922, 13.98170916666669, -2.3662137499999929, 12.111165250000013, -6.8334957499999902, -21.379336749999993, 8.4651301666666967, 2.5094612500000082, -0.21386274999998989, 5.1226162500000072, 14.283680166666699, -4.3340977499999909, -2.7831607499999933, 8.2339832500000085, -12.841856749999991, -6.4984398333333075, -6.2773697499999912, -13.638411749999996, -15.90088974999999, -8.2505068333333043, -19.616496749999996, -4.4346607499999919, -10.056376749999991, -13.581729833333299, -8.2284047499999957, -4.5957137499999945, -5.3758427499999968, -12.254779833333302, 11.207287250000007, -12.848971749999997, -14.449801749999992, -17.247984749999993, -17.475253833333305)
> mean(s)
[1] 1.243449787580175e-14
> sd(s)
[1] 12.25589024484334
> (s - mean(s)) / sd(s)
 [1]  1.01197489551755737 -0.10902853430514588  2.09475824715945480  0.56110703609584245 ...

我知道差异很小，但这对我的应用程序来说有点问题。另外值得注意的是，R 结果也与 Stata 的结果相匹配。

注意：我正在使用 Python 2.7.2、NumpPy 1.6.1、R 2.15.2 GUI 1.53 Leopard 构建 64 位 (6335)

4

3 回答 3

14

对于std明显偏离了相当大数量的 in numpy，默认std返回sqrt(sum((x-x.mean())**2)) / (n-ddof)where ddof=0。我猜R假设ddof=1，因为：

In [7]: s.std()
Out[7]: 12.137473069268983

In [8]: s.std(ddof=1)
Out[8]: 12.255890244843339

和：

> sd(s)
[1] 12.25589

我无法解释mean，但由于在每种情况下它基本上为零，我将其称为精度问题。numpy会在默认容差下将它们报告为“足够接近”：

In [5]: np.isclose(s.mean(), 1.24345e-14)
Out[5]: True

其他答案比我更好地讨论了这个问题。

于 2013-10-04T22:41:23.227 回答
9

这揭示了其中的一些，使用纯 Python，s列表如原始帖子中给出：

>>> import math
>>> sum(s) / len(s)
1.3664283380001927e-14
>>> math.fsum(s) / len(s)
1.2434497875801753e-14

第一个输出重现np.mean()，第二个输出重现 R mean()（我敢肯定，如果使用 R 代码，options(digits=17)它们将是相同的）。

Python 中的不同之处在于，sum()在每次加法后添加“从左到右”会产生舍入误差，同时在math.fsum()概念上计算无限精度和，最后总共进行一次舍入，以用最接近的可表示形式替换无限精度和双精度数。

美元对甜甜圈说这也是 R 所做的。这可以解释为什么@John 报告 R 返回相同的平均值，而不管数字的顺序如何s（无限精度的和对被加数的顺序完全不敏感）。

不过，我不认为这就是结束。R 可能也在使用更好的数值方法来计算标准偏差——在较小的数值误差的意义上“更好”，但在需要更多时间来计算的意义上可能“更糟”。

请注意，PEP 450 - “将统计模块添加到标准库”最近被 Python 接受。这会将这些东西的一些高质量（数字）实现添加到标准库中。当然，这numpy取决于他们是否也想使用这些。

关于标准开发

因为无论如何计算均值都接近 0，而且其中的数字s根本不接近 0，所以计算出的均值的差异几乎无关紧要。为了证明这一点，这是一个进行无限精度计算的构建块（同样是普通的 Python）：

from fractions import Fraction
def sumsq(xs):
    fs = [Fraction(x) for x in xs]
    mean = sum(fs) / len(fs)
    return sum((f - mean)**2 for f in fs)

现在我们可以使用它来产生非常高质量（而且非常慢！）的总体和样本标准差估计：

>>> ss = sumsq(s)
>>> ss  # exact result:  no rounding errors so far!
Fraction(606931231449932225838747590566767, 79228162514264337593543950336)
>>> from math import sqrt
>>> sqrt(ss / len(s))  # population sdev with 2 roundings
12.137473069268983 
>>> sqrt(ss / (len(s) - 1))     # sample sdev with 2 roundings
12.255890244843338

所以 - 惊喜，惊喜 ;-) -np.std(s)计算了总体标准偏差的最佳可能双重近似值，并且 Rsd()计算了样本标准偏差的最佳可能双重近似值。

因此，在这种特定情况下，计算平均值之间的数值差异是一个红鲱鱼 -因为与原始数字相比，平均值很小，所以几乎任何计算标准偏差的方法都会给出良好的数值结果。

这里真正的区别只是默认np使用n分母（population sdev），而 R 默认使用n-1分母（sample sdev）。

于 2013-10-05T00:26:58.233 回答
5

请记住，64 位的精度仅为 2e-16 左右。如果将这些数字相加，您会发现总和与平均值一样非常接近 0。因此问题可能与该精度有关。您引用的每个函数都需要先对数字求和。于是，我又回到了起点。

在 RReduce('+', s)中产生与 python 函数相同的总和sum。在 R 和 Python 中，它们的总和实际上完全相同。但是，R 中的meanandsum函数使用更精确的方法来进行数学运算。当您以与在 numpy 中完成的相同方式在 R 中进行所有数学运算时，它是相同的。

有理由担心您正在使用的 python 计算。您使用的 R 代码实际上可以更好地处理事情。尝试：

# R
sum(s)
sum(s * 10000) / 10000
Reduce('+', s)
Reduce('+', s*10000)/10000

# python (numpy is the same here)
sum(s)
sum(s * 10000) / 10000

R 中的sumin 可以正确处理缩放，因为两个总和是相同的。但是，R 和 python 都无法使用顺序求和方法来处理它。您可以尝试的另一件事是打乱数字。我不会提供代码，但sum在 R 中始终给出相同的值，而Reduce在 R 中，sum在 python 中根据订单给出不同的值

所以你会怎么做？我建议您必须接受您的精度只有这么高，并将接近 0 的值视为 0。正如您所见，这会给您带来问题，因为函数在内部对这些数字求和，例如均值和标准偏差。当您开始然后进行方差时，源自总和的平均误差就会爆炸。也许更多关于为什么这些数字必须相同的信息将帮助您获得更准确的建议。

如果相同是最重要的，那么有一个替代方案应该有效。不要使用 R 的内置函数。它们质量太高，突出了 numpy 统计数据中的问题。如果您像我向您展示的那样滚动平均值和标准差，Reduce那么结果将是相同的。然而，你要做的是让 R 变得更慢，更不精确。如果您可以完全避免此选项，请这样做。例如：

npMean <- function(x) Reduce('+', x)/length(x)
npMean(s)
npSD <- function(x) {m <- npMean(x); sqrt( Reduce('+', (x - m)^2)/(length(x)) )}
npSD(s)

将准确给出 python 平均值和（不正确的）numpy SD。这些会起作用，但有时很难绕过 R 的胆量，让事情对你来说太精确了。当然，如果你能找到 python 函数来替换 numpy 函数并使你的 python 代码更准确，那就更好了。

于 2013-10-04T23:02:58.320 回答