Среднее от панд и numpy отличается

Укороченная версия:

Причина этого в том, что 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, как pandasnanmean. Немного покопавшись, вы можете найти источник pandasnanmean в pandas.core.nanops. Интересно, что сначала кажется, что должен соответствует 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 как функция pandasnanmean, поэтому 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 было бы неплохо ввести здесь флаг (один для скорости, другой для большей точности, по умолчанию - для скорости, поскольку это текущий имп). Некоторые пользователи гораздо больше заботятся о точности вычислений, чем о скорости, с которой они происходят.

HTH.

Ответ @Matt Messersmith - отличное расследование, но я хотел бы добавить, на мой взгляд, важный момент: оба результата (numpy и pandas) неверны. Однако numpy имеет более высокую вероятность ошибиться, чем panda.

Нет принципиальной разницы между использованием float32 и float64, однако для float32 проблемы могут наблюдаться для меньших наборов данных, чем для float64.

На самом деле не определено, как следует рассчитывать mean - данное математическое определение однозначно только для бесконечно точных чисел, но не для операций с плавающей запятой, которые используют наши ПК.

Итак, какова «правильная» формула?

    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-10 после суммирования 10 ^ 6. Нам это может показаться точным, но в сравнении с возможной точностью это все равно фиаско!

С другой стороны, Numpy (по крайней мере, для данной серии) добавит два почти равных элемента. В этом случае верхняя граница результирующей относительной ошибки равна eps*log_2(n), что составляет

• максимальный 2e-6 для float32 и 10 ^ 6 элементов
• максимальный 2e-15 для float64 и 10 ^ 6 элементов.

Из вышесказанного, среди прочего, есть следующие примечательные последствия:

• если среднее значение распределения равно 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 нельзя представить точно с плавающей точкой.

NB: вышесказанное справедливо для одномерных массивов numpy. Ситуация более сложная для суммирования по оси для многомерных массивов numpy, поскольку numpy иногда переключается на наивное суммирование. Для более подробного исследования см. Этот SO-сообщение, который также объясняет @Mark Dickinson наблюдение, то есть:

np.ones((2, 10**8), dtype=np.float32).mean(axis=1) are accurate but np.ones((10**8, 2), dtype=np.float32).mean(axis=0) aren't