Valore medio numpy errato?

Question

Di solito lavoro con enormi simulazioni. A volte, ho bisogno di calcolare il centro di massa dell'insieme di particelle. Ho notato che in molte situazioni, il valore medio restituito da numpy.mean() è sbagliato. Posso capire che è dovuto alla saturazione dell'accumulatore. Per evitare il problema, posso suddividere la somma su tutte le particelle in un piccolo insieme di particelle, ma è scomodo. Qualcuno ha idea di come risolvere questo problema in modo elegante?

Solo per piking la vostra curiosità, il seguente esempio produrre qualcosa di simile a quello che osservo nel mio simulazioni:

import numpy as np 
a = np.ones((1024,1024), dtype=np.float32)*30504.00005 

se si controlla l'valori massimi e minimi, si ottiene:

a.max() 
30504.0 
a.min() 
30504.0 

tuttavia, il valore medio è:

a.mean() 
30687.236328125 

si può capire che qualcosa non va Qui. Questo non sta accadendo quando si usa dtype = np.float64, quindi dovrebbe essere bello risolvere il problema per precisione singola.

Answer 1

Questo non è un problema di NumPy, è un problema in virgola mobile. Lo stesso accade in C:

float acc = 0; 
for (int i = 0; i < 1024*1024; i++) { 
    acc += 30504.00005f; 
} 
acc /= (1024*1024); 
printf("%f\n", acc); // 30687.304688 

(Live demo)

Il problema è che in virgola mobile è limitata precisione; man mano che il valore dell'accumulatore cresce in relazione agli elementi aggiunti ad esso, la precisione relativa diminuisce.

Una soluzione è limitare la crescita relativa, costruendo un albero ad albero. Ecco un esempio in C (il mio Python non è abbastanza buono ...):

float sum(float *p, int n) { 
    if (n == 1) return *p; 
    for (int i = 0; i < n/2; i++) { 
     p[i] += p[i+n/2]; 
    } 
    return sum(p, n/2); 
} 

float x[1024*1024]; 
for (int i = 0; i < 1024*1024; i++) { 
    x[i] = 30504.00005f; 
} 

float acc = sum(x, 1024*1024); 

acc /= (1024*1024); 
printf("%f\n", acc); // 30504.000000 

(Live demo)

Answer 2

È possibile chiamare np.mean con un argomento dtype parola chiave, che specifica il tipo di accumulatore (che per impostazione predefinita è dello stesso tipo dell'array per gli array a virgola mobile).

Quindi chiamando a.mean(dtype=np.float64) risolverà il vostro esempio di giocattoli, e, forse, il problema con gli array più grandi.

Answer 3

rapido e sporco risposta

assert a.ndim == 2 
a.mean(axis=-1).mean() 

Questo dà il risultato previsto per il * 1024 matrice di 1024, ma naturalmente questo non sarà vero per gli array più grandi ...

Se calcolando la media sarà non essere un collo di bottiglia nel codice Vorrei implementarmi un algoritmo ad hoc in python: i dettagli dipendono comunque dalla struttura dei dati.

Se calcolare la media è un collo di bottiglia, allora un algoritmo di riduzione (parallelo) specializzato potrebbe risolvere il problema.

Modifica

Questo approccio può sembrare sciocco, ma di sicuro a mitigare il problema ed è quasi efficiente come .mean() stessa.

In [65]: a = np.ones((1024,1024), dtype=np.float32)*30504.00005 

In [66]: a.mean() 
Out[66]: 30687.236328125 

In [67]: a.mean(axis=-1).mean() 
Out[67]: 30504.0 

In [68]: %timeit a.mean() 
1000 loops, best of 3: 894 us per loop 

In [69]: %timeit a.mean(axis=-1).mean() 
1000 loops, best of 3: 906 us per loop 

Dare una risposta più sensata richiede qualche informazione in più sulle strutture di dati, è formati, e architeture bersaglio.

Answer 4

Si può parzialmente ovviare a questo utilizzando un built-in math.fsum, che rintraccia le somme parziali (i documenti contengono un link ad una ricetta prototipo AS):

>>> fsum(a.ravel())/(1024*1024) 
30504.0 

Per quanto ne sono a conoscenza , numpy non ha un analogo.