Precisión, ¿por qué Matlab y Python numpy dan salidas tan diferentes?

Sé sobre los tipos de datos básicos y que los tipos flotantes (flotante, doble) no pueden contener algunos números exactamente.

Sin embargo, al portar algunos códigos de Matlab a Python (Numpy), encontré algunas diferencias significativas en los cálculos, y creo que está volviendo a la precisión.

Tome el siguiente código, z-normalizando un vector de 500 dimensiones con solo los dos primeros elementos que tienen un valor distinto de cero.

Matlab:

Z = repmat(0,500,1); Z(1)=3;Z(2)=1; Za = (Z-repmat(mean(Z),500,1)) ./ repmat(std(Z),500,1); Za(1) >>> 21.1694 

Pitón:

 from numpy import zeros,mean,std Z = zeros((500,)) Z[0] = 3 Z[1] = 1 Za = (Z - mean(Z)) / std(Z) print Za[0] >>> 21.1905669677 

Además de que el formato muestra un poco más de dígitos en Python, hay una gran diferencia (imho), más de 0.02

Tanto Python como Matlab están utilizando un tipo de datos de 64 bits (afaik). Python usa ‘numpy.float64’ y Matlab ‘double’.

¿Por qué la diferencia es tan grande? ¿Cuál es más correcto?

Tal vez la diferencia viene de las llamadas mean y std . Compara esos primero.

Hay varias definiciones para std , algunas usan la raíz sqaure de

 1 / n * sum((xi - mean(x)) ** 2) 

otros usan

 1 / (n - 1) * sum((xi - mean(x)) ** 2) 

en lugar.

Desde un punto matemático: estas fórmulas son estimadores de la varianza de una variable aleatoria distribuida normal. La distribución tiene dos parámetros sigma y mu . Si sabe mu exactamente, el estimador óptimo para sigma ** 2 es

 1 / n * sum((xi - mu) ** 2) 

Si tiene que estimar mu partir de los datos usando mu = mean(xi) , el estimador óptimo para sigma**2 es

 1 / (n - 1) * sum((xi- mean(x))**2) 

Para responder a su pregunta, no , este no es un problema de precisión. Como @rocksportrocker señala, hay dos estimadores populares para la desviación estándar . El estándar de MATLAB tiene ambos disponibles pero, como estándar, utiliza uno diferente al que usó en Python.

Pruebe std(Z,1) lugar de std(Z) :

 Za = (Z-repmat(mean(Z),500,1)) ./ repmat(std(Z,2),500,1);Za(1) sprintf('%1.10f', Za(1)) 

lleva a

Za (1) = 21.1905669677

en MATLAB. Lea la respuesta de rockspotrocker sobre cuál de los dos resultados es más apropiado para lo que quiere hacer ;-).

Según la documentación de std en SciPy , tiene un parámetro llamado ddof :

ddof : int, opcional
Significa Delta Grados de Libertad. El divisor utilizado en los cálculos es N – ddof, donde N representa el número de elementos. Por defecto, ddof es cero.

En numpy, ddof es cero por defecto, mientras que en MATLAB es uno. Entonces, creo que esto puede resolver el problema:

 std(Z,ddof=1)