【问题标题】:Long (>20million element) array summation in python numpypython numpy中的长(> 2000万个元素)数组求和
【发布时间】:2011-12-22 04:11:02
【问题描述】:

我是 python 和 numpy 的新手,所以如果这个问题太初级,请原谅! 我有一个负值数组(已排序):

>>>neg
[ -1.53507843e+02  -1.53200012e+02  -1.43161987e+02 ...,  -6.37326136e-1 -3.97518490e-10  -3.73480691e-10]
>>>neg.shape
(12922508,)

我需要将此数组添加到其副本(但具有正值)以找到平均为零的分布的标准偏差。所以我做了以下事情:

>>>pos=-1*neg
>>>pos=pos[::-1] #Just to make it look symmetric for the display bellow!
>>>total=np.hstack((neg,pos))
>>>total
[-153.50784302 -153.20001221 -143.1619873  ...,  143.1619873   153.20001221  153.50784302]
>>>total.shape
(25845016,)

到目前为止一切都很好,但奇怪的是这个新数组的和不为零:

>>>numpy.sum(total)
11610.6

标准差也完全不符合我的预期,但我猜这个问题的根源与此相同:为什么总和结果不为零?

当我将此方法应用于一个小数组时;例如 [-5, -3, -2] 和变为零。所以我猜问题在于数组的长度(超过 2000 万个元素)。有什么办法可以解决这个问题吗?

如果有人可以帮助我,我将不胜感激。

【问题讨论】:

  • math.fsum(total) 是否返回 0
  • 是的!哇!!!你的意思是我根本不应该使用 numpy,我误以为 numpy 是处理数组的最佳工具!!!但是查看docs.python.org/py3k/library/math.html#module-math 我没有看到任何计算标准偏差的工具。你有什么建议?
  • 没有。 fsum() 只是为了检查您的代码除了在求和过程中失去精度之外没有其他错误。 numpy.std() 可用于标准偏差。试试np.std(total, dtype=np.float64)
  • 我猜你看到了溢出问题。注意sum([1e308, 1, -1e308]) == 0.0math.fsum([1e308, 1, -1e308]) == 1.0
  • 它有助于在求和之前按绝对值对数据进行排序(尤其是在这种情况下,您希望正负贡献相互抵消)。它还有助于首先以块(例如 100000 左右)进行部分总和,然后将部分总和相加。

标签: python numpy


【解决方案1】:

如 cmets 中所述,将数百万个等号数字相加会导致浮点舍入问题。解决此问题的一种可能方法是在组合数组中混合正数和负数,以便求和时的任何中间结果始终大致保持在同一数量级内:

neg = -100*numpy.random.rand(20e6)
pos = -neg
combined = numpy.zeros(len(neg)+len(pos))
combined[::2] = neg
combined[1::2] = pos

现在combined.sum() 应该非常接近于零。

也许这种方法也有助于提高计算标准差的精度。

【讨论】:

  • 总和现在正好为零,同样不为零的均值也是如此。但非常奇怪的是,标准差并没有改变以前的样子。除非numpy.std() 使用另一种方法来计算总和(例如在math.fsum() 模块中),否则此结果是不可接受的,因为numpy.sum()numpy.mean() 发生了变化,numpy.std() 没有!!!!
  • 您需要自己实现std(),使用我们描述的技术在公式中进行求和。
  • 我写了一个程序来计算std()我自己;在 10,000 个元素的包中找到每个值与平均值的差异,将结果求和,最后除以元素的数量并求平方根。计算所有超过 2500 万个元素的标准偏差大约需要 15 分钟,它找到了与 numpy.std() 完全相同的值(只花了几分之一秒!)。这个数据集是对我算法的测试:正如我所说,我从一开始就知道标准偏差。我会检查那个标准差的来源,看看它是否正确!
  • 我被告知的标准偏差值是在 3sigma 剪辑之后,我没有被告知,在应用 3sigma 剪辑之后,我得到了我想要的值。非常感谢大家。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2014-01-05
  • 1970-01-01
  • 2018-09-09
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多