Niepoprawna wartość średnia?

Question

Zazwyczaj pracuję z ogromnymi symulacjami. Czasami muszę obliczyć środek masy zestawu cząsteczek. Zauważyłem, że w wielu sytuacjach średnia wartość zwracana przez numpy.mean() jest niepoprawna. Mogę zrozumieć, że jest to spowodowane nasyceniem akumulatora. Aby uniknąć tego problemu, mogę podzielić sumowanie na wszystkie cząstki w małym zestawie cząsteczek, ale jest niewygodne. Czy ktoś ma pomysł, jak rozwiązać ten problem w elegancki sposób?

Tylko dla piking się swoją ciekawość, następujący przykład produkować coś podobnego do tego, co obserwujemy w moich symulacji:

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

sprawdzając max i wartości min, otrzymasz:

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

jednak wartość średnia wynosi:

a.mean() 
30687.236328125 

można dowiedzieć się, że coś jest nie tak tutaj. Nie dzieje się tak podczas używania dtype = np.float64, więc powinno być fajnie rozwiązać problem z pojedynczą precyzją.

Answer 1

To nie jest problem NumPy, jest to problem zmiennoprzecinkowy. To samo ma miejsce w 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)

Problemem jest to, że zmiennoprzecinkowy nie ograniczają dokładność; gdy wartość akumulatora rośnie w stosunku do elementów dodawanych do niego, względna precyzja spada.

Jednym z rozwiązań jest ograniczenie względnego wzrostu poprzez utworzenie drzewa sumującego. Oto przykład w języku C (mój Python nie jest wystarczająco dobre ...):

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

Można zadzwonić np.mean z argumentem dtype słów kluczowych, które określa rodzaj akumulatora (która domyślnie jest tego samego typu, co tablica dla tablic zmiennoprzecinkowych).

Wywołanie numeru a.mean(dtype=np.float64) rozwiąże Twój przykład zabawek, a może Twój problem z większymi tablicami.

Answer 3

Szybki i brudny odpowiedź

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

To daje oczekiwanego rezultatu do 1024 * 1024 matrycy, ale oczywiście to nie będzie prawdziwa dla większych tablic ...

przypadku obliczania średniej będzie Nie wąskie gardło w twoim kodzie Wdałbym sobie algorytm ad-hoc w pythonie: szczegóły zależą jednak od twojej struktury danych.

Jeśli obliczanie średniej jest wąskim gardłem, wówczas może pomóc rozwiązać pewien wyspecjalizowany (równoległy) algorytm redukcji.

Edit

Takie podejście może wydawać się głupie, ale na pewno złagodzić ten problem i jest prawie tak samo skuteczne jak .mean().

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 

Udzielenie rozsądniejszej odpowiedzi wymaga więcej informacji o strukturach danych, rozmiarach i architekturze docelowej.

Answer 4

można częściowo zaradzić za pomocą wbudowanego math.fsum, który tropi sum częściowych (docs zawierać link do prototypu receptury AS):

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

O ile jestem świadomy , numpy nie ma analog.