Rob*_*Rob 29 python floating-point numpy floating-accuracy pandas
我有一个MEMS IMU,我一直在收集数据,我正在使用pandas从中获取一些统计数据.每个循环收集6个32位浮点数.对于给定的集合运行,数据速率是固定的.数据速率在100Hz和1000Hz之间变化,收集时间长达72小时.数据保存在平面二进制文件中.我这样读了数据:
import numpy as np
import pandas as pd
dataType=np.dtype([('a','<f4'),('b','<f4'),('c','<f4'),('d','<f4'),('e','<f4'),('e','<f4')])
df=pd.DataFrame(np.fromfile('FILENAME',dataType))
df['c'].mean()
-9.880581855773926
x=df['c'].values
x.mean()
-9.8332081
-9.833是正确的结果.我可以创建一个类似的结果,有人应该能够这样重复:
import numpy as np
import pandas as pd
x=np.random.normal(-9.8,.05,size=900000)
df=pd.DataFrame(x,dtype='float32',columns=['x'])
df['x'].mean()
-9.859579086303711
x.mean()
-9.8000648778888628
我在Linux和Windows上,在AMD和Intel处理器上,在Python 2.7和3.5中重复了这一点.我很难过.我究竟做错了什么？得到这个:
x=np.random.normal(-9.,.005,size=900000)
df=pd.DataFrame(x,dtype='float32',columns=['x'])
df['x'].mean()
-8.999998092651367
x.mean()
-9.0000075889406528
我可以接受这种差异.它处于32位浮点数精度的极限.
没关系.我星期五写了这个,解决方案今天早上给了我.这是由大量数据加剧的浮点精度问题.我需要在创建数据帧时将数据转换为64位浮点数:
df=pd.DataFrame(np.fromfile('FILENAME',dataType),dtype='float64')
如果其他人遇到类似问题,我会留下帖子.
Mat*_*ith 19
之所以不同,是因为在调用操作时pandas使用bottleneck(如果已安装)mean,而不是仅仅依赖numpy.bottleneck可能是因为它似乎比numpy(至少在我的机器上)更快,但是以精确度为代价.它们碰巧匹配64位版本,但32位不同(这是有趣的部分).
通过检查这些模块的源代码来判断发生了什么是非常困难的(它们非常复杂,即使对于简单的计算,例如mean数值计算也很难).最好使用调试器来避免大脑编译和那些类型的错误.调试器不会在逻辑上出错,它会告诉你究竟发生了什么.
这是我的一些堆栈跟踪(由于没有RNG的种子,值略有不同):
可以重现(Windows):
>>> import numpy as np; import pandas as pd
>>> x=np.random.normal(-9.,.005,size=900000)
>>> df=pd.DataFrame(x,dtype='float32',columns=['x'])
>>> df['x'].mean()
-9.0
>>> x.mean()
-9.0000037501099754
>>> x.astype(np.float32).mean()
-9.0000029
没有什么特别numpy的版本.这个pandas版本有点古怪.
我们来看看里面df['x'].mean():
>>> def test_it_2():
... import pdb; pdb.set_trace()
... df['x'].mean()
>>> test_it_2()
... # Some stepping/poking around that isn't important
(Pdb) l
2307
2308 if we have an ndarray as a value, then simply perform the operation,
2309 otherwise delegate to the object
2310
2311 """
2312 -> delegate = self._values
2313 if isinstance(delegate, np.ndarray):
2314 # Validate that 'axis' is consistent with Series's single axis.
2315 self._get_axis_number(axis)
2316 if numeric_only:
2317 raise NotImplementedError('Series.{0} does not implement '
(Pdb) delegate.dtype
dtype('float32')
(Pdb) l
2315 self._get_axis_number(axis)
2316 if numeric_only:
2317 raise NotImplementedError('Series.{0} does not implement '
2318 'numeric_only.'.format(name))
2319 with np.errstate(all='ignore'):
2320 -> return op(delegate, skipna=skipna, **kwds)
2321
2322 return delegate._reduce(op=op, name=name, axis=axis, skipna=skipna,
2323 numeric_only=numeric_only,
2324 filter_type=filter_type, **kwds)
所以我们找到了麻烦点,但现在事情变得有些奇怪了:
(Pdb) op
<function nanmean at 0x000002CD8ACD4488>
(Pdb) op(delegate)
-9.0
(Pdb) delegate_64 = delegate.astype(np.float64)
(Pdb) op(delegate_64)
-9.000003749978807
(Pdb) delegate.mean()
-9.0000029
(Pdb) delegate_64.mean()
-9.0000037499788075
(Pdb) np.nanmean(delegate, dtype=np.float64)
-9.0000037499788075
(Pdb) np.nanmean(delegate, dtype=np.float32)
-9.0000029
请注意,delegate.mean()并np.nanmean输出-9.0000029型float32,不 -9.0为pandas nanmean呢.只要一点闲逛的,你可以找到源头pandas nanmean在pandas.core.nanops.有趣的是,它实际上似乎应该首先匹配numpy.我们来看看pandas nanmean:
(Pdb) import inspect
(Pdb) src = inspect.getsource(op).split("\n")
(Pdb) for line in src: print(line)
@disallow('M8')
@bottleneck_switch()
def nanmean(values, axis=None, skipna=True):
values, mask, dtype, dtype_max = _get_values(values, skipna, 0)
dtype_sum = dtype_max
dtype_count = np.float64
if is_integer_dtype(dtype) or is_timedelta64_dtype(dtype):
dtype_sum = np.float64
elif is_float_dtype(dtype):
dtype_sum = dtype
dtype_count = dtype
count = _get_counts(mask, axis, dtype=dtype_count)
the_sum = _ensure_numeric(values.sum(axis, dtype=dtype_sum))
if axis is not None and getattr(the_sum, 'ndim', False):
the_mean = the_sum / count
ct_mask = count == 0
if ct_mask.any():
the_mean[ct_mask] = np.nan
else:
the_mean = the_sum / count if count > 0 else np.nan
return _wrap_results(the_mean, dtype)
这是bottleneck_switch装饰器的(短)版本:
import bottleneck as bn
...
class bottleneck_switch(object):
def __init__(self, **kwargs):
self.kwargs = kwargs
def __call__(self, alt):
bn_name = alt.__name__
try:
bn_func = getattr(bn, bn_name)
except (AttributeError, NameError): # pragma: no cover
bn_func = None
...
if (_USE_BOTTLENECK and skipna and
_bn_ok_dtype(values.dtype, bn_name)):
result = bn_func(values, axis=axis, **kwds)
这就是所谓的有alt作为pandas nanmean的功能,所以bn_name就是'nanmean',这是是一个从抢过ATTR bottleneck模块:
(Pdb) l
93 result = np.empty(result_shape)
94 result.fill(0)
95 return result
96
97 if (_USE_BOTTLENECK and skipna and
98 -> _bn_ok_dtype(values.dtype, bn_name)):
99 result = bn_func(values, axis=axis, **kwds)
100
101 # prefer to treat inf/-inf as NA, but must compute the fun
102 # twice :(
103 if _has_infs(result):
(Pdb) n
> d:\anaconda3\lib\site-packages\pandas\core\nanops.py(99)f()
-> result = bn_func(values, axis=axis, **kwds)
(Pdb) alt
<function nanmean at 0x000001D2C8C04378>
(Pdb) alt.__name__
'nanmean'
(Pdb) bn_func
<built-in function nanmean>
(Pdb) bn_name
'nanmean'
(Pdb) bn_func(values, axis=axis, **kwds)
-9.0
假装bottleneck_switch()装饰器暂时不存在.我们实际上可以看到调用手动单步执行此函数(不带bottleneck)将获得与以下相同的结果numpy:
(Pdb) from pandas.core.nanops import _get_counts
(Pdb) from pandas.core.nanops import _get_values
(Pdb) from pandas.core.nanops import _ensure_numeric
(Pdb) values, mask, dtype, dtype_max = _get_values(delegate, skipna=skipna)
(Pdb) count = _get_counts(mask, axis=None, dtype=dtype)
(Pdb) count
900000.0
(Pdb) values.sum(axis=None, dtype=dtype) / count
-9.0000029
但是,如果你已bottleneck安装,那就永远不会被调用.相反,bottleneck_switch()装饰器反而nanmean使用bottleneck版本来破坏函数.这就是差异所在(有趣的是它与float64案例相符):
(Pdb) import bottleneck as bn
(Pdb) bn.nanmean(delegate)
-9.0
(Pdb) bn.nanmean(delegate.astype(np.float64))
-9.000003749978807
bottleneck据我所知,它仅用于速度.我假设他们正在使用某种类型的快捷方式nanmean,但是我没有对它进行过多的研究(有关此主题的详细信息,请参阅@ ead的答案).您可以看到它通常比numpy基准测试快一点:https://github.com/kwgoodman/bottleneck.显然,为这个速度付出的代价是精确的.
瓶颈实际上更快吗？
当然看起来像(至少在我的机器上).
In [1]: import numpy as np; import pandas as pd
In [2]: x=np.random.normal(-9.8,.05,size=900000)
In [3]: y_32 = x.astype(np.float32)
In [13]: %timeit np.nanmean(y_32)
100 loops, best of 3: 5.72 ms per loop
In [14]: %timeit bn.nanmean(y_32)
1000 loops, best of 3: 854 µs per loop
pandas在这里引入一个标志可能会很好(一个用于速度,另一个用于更好的精度,默认用于速度,因为那是当前的impl).一些用户更关心计算的准确性而不是它发生的速度.
HTH.
ead*_*ead 16
@Matt Messersmith的答案是一个很好的调查,但我想在我看来重要的一点:两个结果(numpy和pandas')是错误的.然而,numpy比熊猫更容易出错.
存在使用之间没有根本的区别float32和float64,但是,对于float32可以用于较小的数据集比被观察的问题float64.
它没有真正定义,mean应该如何计算 - 给定的数学定义对于无限精确的数字只是明确的,而不是我们的PC正在使用的浮点运算.
那么"正确"的公式是什么？
mean = (x0+..xn)/n
or
mean = [(x0+x1)+(x2+x3)+..]/n
or
mean = 1.0/n*(x0+..xn)
and so on...
显然,当在现代硬件上计算时,它们都将给出不同的结果 - 理想情况下,可以查看一个公式,该公式使得与理论正确值(以无限精度计算)相比的最小误差.
Numpy使用略微交替的成对求和,即(((x1+x2)+(x3+x4))+(...)),即使不完美,也已知非常好.另一方面,瓶颈使用天真求和x1+x2+x3+...:
REDUCE_ALL(nanmean, DTYPE0)
{
...
WHILE {
FOR {
ai = AI(DTYPE0);
if (ai == ai) {
asum += ai; <---- HERE WE GO
count += 1;
}
}
NEXT
}
...
}
我们可以很容易地看到发生的事情:经过一些步骤,bottleneck总结一个大的(以前的所有元素的总和,比例-9.8*number_of_steps)和一个小元(约-9.8),这导致约相当一个舍入误差big_number*eps,与EPS被周围1e-7的float32.这意味着在10 ^ 6个求和之后,我们可能会有大约10%的相对误差(eps*10^6这是一个上限).
为float64和eps为约1e-16相对误差将是只有大约1e-1010 ^ 6个求和之后.它对我们来说似乎是精确的,但是根据可能的精度来衡量它仍然是一场惨败!
另一方面,Numpy(至少对于手头的系列)将添加两个几乎相等的元素.在这种情况下,得到的相对误差的上限是eps*log_2(n),即
2e-6为float3210 ^ 6个元素2e-15为float6410 ^ 6个元素.除上述之外,还有以下值得注意的含义:
0,那么大熊猫和numpy几乎同样精确 - 总0.0和数的大小是大约的,并且在求和之间没有大的差异,这将导致天真求和的大舍入误差.x'i=xi-mean_estimate,因为它x'i具有平均值0.0.x=(.333*np.ones(1000000)).astype(np.float32)足以触发大熊猫版本的奇怪行为 - 不需要随机性,我们知道结果应该是什么,不是吗？重要的是,0.333不能用浮点精确呈现.