任何想法为什么R和Python的NumPy向量缩放不匹配？

Question

我有以下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 build 64-bit(6335)

Answer 1

对于std,这显然是关闭的一些实质性的量,在numpy,std返回sqrt(sum((x-x.mean())**2)) / (n-ddof)其中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

其他答案比我能更好地讨论这个问题.

Answer 2

这使用纯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可能也使用更好的数值方法来计算std dev - 从更小的数值误差的意义上说"更好",但在花费更多时间计算的意义上可能"更差".

请注意,PEP 450 - "将统计模块添加到标准库"最近被Python接受.这将为标准库添加一些高质量(数字)的实现.当然,numpy还要决定是否要使用这些.

关于std dev

因为无论如何计算均值接近于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)计算出最佳的双逼近总体标准差,和R的sd()计算最佳的双重逼近样本标准差.

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

这里真正的区别仅仅是np使用n默认的分母(人口发展局局长),而R使用n-1分母(样品SDEV)默认情况下.

Answer 3

请记住,64位的精度仅约为2e-16.如果你对这些数字求和,你会发现总和,就像平均值一样,非常接近于0.所以这个问题很可能与那个精度有关.您引用的每个函数都需要先对数字求和.所以,我回到了开始.

在R中Reduce('+', s)产生与python函数相同的总和sum.在R和Python中,它们实际上总结完全相同.但是,R中的函数mean和sum函数使用更准确的方法来进行数学运算.当你在R中完成所有的数学运算时,它就像在numpy中完成的那样,那么它就是相同的.

有理由担心你正在使用的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

在sumR中处理的比例适当为两个和是相同的.但是,R和python都无法使用顺序求和方法来处理它.你可以尝试的另一件事是争抢数字.我不会提供代码但sum在R中始终给出相同的值,而Reduce在R和sumpython中根据订单给出不同的值

所以你会怎么做？我建议你必须接受你的精度只有这么高,并将你的值接近0视为0.这会给你带来问题,正如你所见,这些函数在内部将这些数字相加为均值和标准偏差.当你开始做差异时,从总和得到的平均误差就会爆炸.也许有关这些数字必须相同的确切原因的更多信息可以帮助您获得更精确的建议.

有一个替代方案应该有效,如果相同是重要的.不要使用R的内置功能.他们的质量太高,突出了numpy统计数据中的问题.如果你滚动一个平均值和sd,就像我给你看的那样,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代码更准确,那就更好了.