python - round - media de los pandas y los diferentes números difieren




python truncate float (2)

Tengo una IMU de MEMS en la que he estado recopilando datos y estoy usando pandas para obtener algunos datos estadísticos. Hay 6 flotadores de 32 bits recogidos en cada ciclo. Las tasas de datos son fijas para una ejecución de recopilación determinada. Las velocidades de datos varían entre 100Hz y 1000Hz y los tiempos de recolección duran hasta 72 horas. Los datos se guardan en un archivo binario plano. Leí los datos de esta manera:

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 es el resultado correcto. Puedo crear un resultado similar que alguien debería poder repetir de esta manera:

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

He repetido esto en Linux y Windows, en procesadores AMD e Intel, en Python 2.7 y 3.5. Estoy perplejo. ¿Qué estoy haciendo mal? Y consigue esto:

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

Podría aceptar esta diferencia. Está en el límite de la precisión de los flotadores de 32 bits.

NO IMPORTA. Escribí esto el viernes y la solución me golpeó esta mañana. Es un problema de precisión de punto flotante exacerbado por la gran cantidad de datos. Necesitaba convertir los datos en una flotación de 64 bits en la creación del marco de datos de esta manera:

df=pd.DataFrame(np.fromfile('FILENAME',dataType),dtype='float64')

Dejaré la publicación si alguien más se encuentra con un problema similar.


Version corta:

La razón por la que es diferente es porque pandas usa un bottleneck (si está instalado) cuando llama a la operación mean , en lugar de simplemente confiar en el numpy . bottleneck se usa presumiblemente ya que parece ser más rápido que el numpy (al menos en mi máquina), pero a un costo de precisión. Suceden que coinciden con la versión de 64 bits, pero difieren en la tierra de 32 bits (que es la parte interesante).

Versión larga:

Es extremadamente difícil saber qué está pasando solo con inspeccionar el código fuente de estos módulos (son bastante complejos, incluso para cálculos simples como la mean , resulta que la computación numérica es difícil). Lo mejor es utilizar el depurador para evitar la compilación del cerebro y esos tipos de errores. El depurador no cometerá un error en la lógica, le dirá exactamente lo que está pasando.

Aquí hay algo de mi seguimiento de pila (los valores difieren ligeramente ya que no hay semilla para RNG):

Se puede reproducir (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

No hay nada extraordinario en la versión de numpy . Es la versión de los pandas que es un poco chiflada.

Echemos un vistazo dentro de 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)

Así que encontramos el punto problemático, pero ahora las cosas se ponen un poco raras:

(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

Tenga en cuenta que delegate.mean() y np.nanmean salida a -9.0000029 con el tipo float32 , no a -9.0 como lo hace pandas nanmean . Con un poco de curiosear, puedes encontrar la fuente de pandas nanmean en pandas.core.nanops . Curiosamente, en realidad parece que debería coincidir con numpy al principio. Echemos un vistazo a los 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)

Aquí hay una versión (corta) del decorador 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)

Esto se llama con alt como la función de pandas nanmean , por lo que bn_name es 'nanmean' , y este es el attr que se toma del módulo de 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

Imagina que el decorador bottleneck_switch() no existe por un segundo. Realmente podemos ver que llamar a ese paso manual a través de esta función (sin bottleneck ) le dará el mismo resultado que 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

Sin embargo, eso nunca se llama si tiene instalado un bottleneck . En su lugar, el decorador bottleneck_switch() explota la función nanmean con la versión del bottleneck de bottleneck . Aquí es donde radica la discrepancia (aunque es interesante que coincida con el caso de float64 ):

(Pdb) import bottleneck as bn
(Pdb) bn.nanmean(delegate)
-9.0
(Pdb) bn.nanmean(delegate.astype(np.float64))
-9.000003749978807

bottleneck se utiliza únicamente para la velocidad, por lo que puedo decir. Supongo que están tomando algún tipo de método abreviado con su función nanmean , pero no lo examiné mucho (consulte la respuesta de @ ead para obtener detalles sobre este tema). Puede ver que, por lo general, es un poco más rápido que sus numpy comparativos: https://github.com/kwgoodman/bottleneck . Claramente el precio a pagar por esta velocidad es la precisión.

¿Es realmente más rápido el cuello de botella?

Seguro que lo parece (al menos en mi máquina).

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

Puede ser bueno que los pandas introduzcan una bandera aquí (una para la velocidad, otra para una mejor precisión, la predeterminada es la velocidad ya que esa es la implícita actual). Algunos usuarios se preocupan mucho más por la precisión del cálculo que por la velocidad a la que ocurre.

HTH.


La respuesta de @Matt Messersmith es una gran investigación, pero me gustaría agregar un punto importante en mi opinión: ambos resultados (numpy y pandas) son incorrectos. Sin embargo, el número tiene mayor probabilidad de ser menos malo que el panda.

No hay una diferencia fundamental entre el uso de float32 y float64 , sin embargo, para float32 , se pueden observar problemas para conjuntos de datos más pequeños que para float64 .

No está realmente definido, cómo debe calcularse la mean : la definición matemática dada es solo inequívoca para números infinitamente precisos, pero no para las operaciones de punto flotante que están usando nuestras PC.

Entonces, ¿cuál es la fórmula "correcta"?

    mean = (x0+..xn)/n 
  or 
    mean = [(x0+x1)+(x2+x3)+..]/n
  or
    mean = 1.0/n*(x0+..xn)
  and so on...

Obviamente, cuando se calculan en el hardware moderno, todos darán resultados diferentes: lo ideal sería echar un vistazo a una fórmula que produce el error más pequeño en comparación con un valor teórico correcto (que se calcula con una precisión infinita).

Numpy usa una suma por pares ligeramente alternada, es decir (((x1+x2)+(x3+x4))+(...)) , que, aunque no sea perfecta, se sabe que es bastante buena. Por otro lado, el bottleneck utiliza la suma ingenua x1+x2+x3+... :

REDUCE_ALL(nanmean, DTYPE0)
{
    ...
    WHILE {
        FOR {
            ai = AI(DTYPE0);
            if (ai == ai) {
                asum += ai;   <---- HERE WE GO
                count += 1;
            }
        }
        NEXT
    }
    ...
}

y podemos ver fácilmente lo que sucede: después de algunos pasos, el bottleneck suma un gran (suma de todos los elementos anteriores, proporcional a -9.8*number_of_steps ) y un pequeño elemento (aproximadamente -9.8 ) que conduce a un error de redondeo de aproximadamente big_number*eps , con eps alrededor de 1e-7 para float32 . Esto significa que después de 10 ^ 6 sumas podríamos tener un error relativo de aproximadamente el 10% ( eps*10^6 , este es un límite superior).

Para float64 y eps siendo aproximadamente 1e-16 el error relativo sería solo alrededor de 1e-10 después de 10 ^ 6 sumas. Puede que nos parezca preciso, pero si lo comparamos con la posible precisión, ¡sigue siendo un fiasco!

Numpy por otro lado (al menos para la serie en cuestión) agregará dos elementos que son casi iguales. En este caso, el límite superior para el error relativo resultante es eps*log_2(n) , que es

  • float32 2e-6 para float32 y 10 ^ 6 elementos
  • float64 2e-15 para float64 y 10 ^ 6 elementos.

De lo anterior, entre otros, se destacan las siguientes implicaciones:

  • si la media de la distribución es 0 , entonces los pandas y los números son casi igualmente precisos: la magnitud de los números sumados es aproximadamente 0.0 y no existe una gran diferencia entre los sumandos, lo que llevaría a un gran error de redondeo para la suma ingenua.
  • si se conoce una buena estimación de la media, podría ser más robusto calcular la suma de x'i=xi-mean_estimate , porque x'i tendrá una media de 0.0 .
  • algo como x=(.333*np.ones(1000000)).astype(np.float32) es suficiente para desencadenar el extraño comportamiento de la versión de pandas: no hay necesidad de aleatoriedad y sabemos cuál debería ser el resultado, no nosotros Es importante, que 0.333 no se puede presentar precisamente con punto flotante.




floating-accuracy