【问题标题】:mean from pandas and numpy differpandas 和 numpy 的意思不同
【发布时间】:2021-06-28 03:07:14
【问题描述】:

我有一个 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')

如果其他人遇到类似问题,我会离开帖子。

【问题讨论】:

  • 我无法重现您的第一次检查,我收到 float32 大小的错误。请注意,您的 x 包含双精度值,但您的 df 包含浮点数。这总会给你带来不同,但不如原来的那么大。您是否有可能缺少与平均值计算方式混淆的值?
  • 部分问题是 Pandas 使用了糟糕的算法来计算平均值;最终,随着总和的累积,接近-9.8 的值被重复添加到大于2**23 的值上,有限的float32 分辨率意味着对于大多数随机样本,实际总和恰好变化了-10.0。使用成对求和或 Kahan 求和而不是简单的累加求和会大大改善这里的结果。但是,是的,以双精度计算平均值是显而易见的快速解决方法。
  • @MarkDickinson,那么为什么df['x'].sum() / len(df.index) 不会出现问题,即使float32 也能给出正确的结果?
  • @jpp:好问题。我想你得问问 Pandas 的作者。 NumPy 确实 在某些(但不是全部)情况下对其sum 操作使用成对求和;无论出于何种原因,df['x'].sum() 的这种特殊用法都有可能最终出现在其中一种 NumPy 案例中。

标签: python pandas numpy floating-point floating-accuracy


【解决方案1】:

@Matt Messersmith 的回答是一个很好的调查,但我想补充一点我认为很重要的一点:两个结果(numpy 和 pandas)都是错误的。但是,numpy 比 panda 出错的概率更高。

使用float32float64 之间没有根本区别,但是对于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 使用稍微交替的pairwise summation,即(((x1+x2)+(x3+x4))+(...)),即使不完美,也可以说是相当不错的。另一方面,bottleneck 使用了朴素的求和 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,对于 float32,eps 约为 1e-7。这意味着经过 10^6 次求和后,我们可能会有大约 10% 的相对误差(eps*10^6,这是一个上限)。

对于 float64eps 大约是 1e-16 在 10^6 求和后的相对误差将仅约为 1e-10。这对我们来说似乎很精确,但以可能的精确度衡量它仍然是一场惨败!

另一方面,Numpy(至少对于手头的系列)将添加两个几乎相等的元素。在这种情况下,产生的相对误差的上限是eps*log_2(n),即

  • 2e-6 的最大 float32 和 10^6 个元素
  • float64 和 10^6 个元素的最大 2e-15

从以上内容,除其他外,还有以下值得注意的含义:

  • 如果分布的均值是0,那么pandas 和numpy 几乎同样精确——相加数字的大小大约是0.0,并且相加数之间没有太大差异,这会导致较大的舍入误差天真的总结。
  • 如果知道对均值的准确估计,计算x'i=xi-mean_estimate 的总和可能会更稳健,因为x'i 的均值是0.0
  • x=(.333*np.ones(1000000)).astype(np.float32) 之类的东西足以触发 pandas 版本的奇怪行为 - 不需要随机性,我们知道结果应该是什么,不是吗?重要的是,0.333 不能用浮点数精确表示。

注意:以上内容适用于一维 numpy 数组。对于多维 numpy 数组,沿轴求和的情况更为复杂,因为 numpy 有时会切换到朴素求和。如需更详细的调查,请参阅此SO-post,这也解释了@Mark Dickinson observation,即:

np.ones((2, 10**8), dtype=np.float32).mean(axis=1) 是准确的,但 np.ones((10**8, 2), dtype=np.float32).mean(axis=0) 不是

【讨论】:

    【解决方案2】:

    短版:

    之所以不同,是因为pandas 在调用mean 操作时使用bottleneck(如果已安装),而不是仅仅依赖numpybottleneck 可能被使用,因为它似乎比 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.core.nanops 中找到pandas nanmean 的来源。有趣的是,它实际上看起来应该首先匹配numpy。一起来看看pandasnanmean

    (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',这是从bottleneck 模块中获取的attr:

    (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() 装饰器使用bottleneck 的版本覆盖了nanmean 函数。这就是差异所在(有趣的是,它与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 在这里引入一个标志可能会很好(一个用于速度,另一个用于更好的精度,默认是速度,因为这是当前的实现)。有些用户更关心计算的准确性,而不是计算速度。

    HTH。

    【讨论】:

    • 您说“numpy 将其打入 float64 以提高精度”,但您显示的代码似乎不支持这一点。在numpy.core._methods._mean 中,总和(对umr_sum 的调用)最终以dtype=None 执行。
    • 啊,如果您正在查看x.mean(),那么x 首先具有dtype np.float64。这可以解释为什么您看到float64 结果在均值之内。
    • 如果你想让 NumPy 在执行求和之前自动从 float32 转换为 float64,请尝试执行 np.ones((10**8, 2), dtype=np.float32).mean(axis=0)。在 NumPy 的案例中,pairwise summation 的使用实际上对准确性产生了影响。 (至于 Pandas 在做什么:我不知道。)
    • 很好的答案+解释。我将给它一些播出时间,以便获得更多观看次数。我希望它能到达 Pandas 开发人员。似乎是一个意想不到的后果,可能会产生奇怪而重大的影响超越 float32 vs float64 精度,例如OP 的极端例子。
    • 嗯,NumPy 的行为也相当愚蠢。 np.ones((10**8, 1), dtype=np.float32).mean(axis=0)np.ones((2, 10**8), dtype=np.float32).mean(axis=1) 是准确的但 np.ones((10**8, 2), dtype=np.float32).mean(axis=0) 不准确的事实是愚蠢的。可以解释,当然,但仍然很愚蠢。
    猜你喜欢
    • 2022-11-07
    • 2020-10-14
    • 2017-11-01
    • 1970-01-01
    • 2014-09-18
    • 2020-05-22
    • 2016-11-21
    • 2021-10-08
    相关资源
    最近更新 更多