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

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

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

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

если вы проверите значения max и min, вы получите:

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

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

 a.mean() 30687.236328125

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

4 Solutions collect form web for “Неверное среднее значение numpy?”

Это не проблема 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 )

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

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

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

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

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

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

 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

Для получения более разумного ответа требуется дополнительная информация о структурах данных, их размерах и целевой архитектуре.