Неверное среднее значение numpy?

Question

Я обычно работаю с огромными симуляторами. Иногда мне нужно вычислить центр масс набора частиц. Я заметил, что во многих ситуациях среднее значение, возвращаемое numpy.mean(), неверно. Я могу понять, что это связано с насыщением аккумулятора. Чтобы избежать проблемы, я могу разделить суммирование по всем частицам в небольшом наборе частиц, но это неудобно. У кого-нибудь есть и идея, как решить эту проблему элегантным способом?

Только для piking вашего любопытства, в следующем примере производить что-то подобное тому, что я наблюдаю в моей симуляции:

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

если вы проверить максимальные и минимальные значения, вы получите:

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

однако, среднее значение:

a.mean() 
30687.236328125 

Вы можете понять, что что-то не так Вот. Этого не происходит при использовании dtype = np.float64, поэтому должно быть хорошо решить проблему для одиночной точности.

Answer 1

Это не проблема NumPy, это проблема с плавающей точкой. То же самое происходит в 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)

Проблема заключается в том, что с плавающей точкой имеет ограниченную точность; так как значение аккумулятора растет относительно добавляемых к нему элементов, относительная точность падает.

Одним из решений является ограничение относительного роста путем создания дерева сумматора. Вот пример в C (мой Python не достаточно хорошо ...):

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

Вы можете вызвать np.mean с аргументом в dtype ключевого слова, который определяет тип аккумулятора (который по умолчанию имеет тот же тип, что и массив для массивов с плавающей запятой).

Так что вызов a.mean(dtype=np.float64) поможет решить вашу игрушку, и, возможно, ваша проблема с более крупными массивами.

Answer 3

Быстрый и грязный ответ

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

Это дает ожидаемый результат для * 1024 матрицы 1024, но, конечно, это не будет справедливо для больших массивов ...

Если вычисления среднего значения воли не будет узким местом в вашем коде, я бы реализовал себе ad-hoc-алгоритм в python: подробности, однако, зависят от вашей структуры данных.

Если вычислить среднее значение является узким местом, то некоторые специализированные (параллельные) алгоритмы сокращения могут решить проблему.

Редактировать

Такой подход может показаться глупым, но наверняка смягчить проблему и почти столь же эффективно, как .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 

Предоставление более разумный ответ требует больше информации о структурах данных, это размеры и цель архитектурных сооружений.

Answer 4

Вы можете частично исправить это с помощью встроенного math.fsum, который отслеживает частичные суммы (Документы содержат ссылку на рецепт прототип AS):

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

Насколько я знаю , numpy не имеет аналога.