Неверное среднее значение 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 

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