numpy 평균값이 잘못 되었습니까?

Question

저는 대개 거대한 시뮬레이션을합니다. 때로는 입자 집합의 질량 중심을 계산해야합니다. 나는 많은 경우에 numpy.mean()에 의해 반환 된 평균값이 잘못되었음을 지적했다. 나는 그것이 누산기의 포화로 인한 것임을 알 수있다. 문제를 피하기 위해 작은 입자 세트의 모든 입자에 대해 합계를 나눌 수 있지만 불편합니다. 누구든지이 문제를 어떻게 우아하게 해결할 것인가? 당신은 최대 및 최소 값을 확인하면, 당신이 얻을

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

:

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

그냥 호기심을 piking를 들어, 다음의 예는 내 시뮬레이션에서 관찰 것과 비슷한 생산

그러나, 평균 값은 다음과 같습니다

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의 예입니다 (내 파이썬은 ... 충분하지 않습니다) :

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

당신은 어큐뮬레이터의 유형 (지정하는 dtype 키워드 인수로 np.mean를 호출 할 수 있습니다 부동 소수점 배열의 배열과 동일한 유형으로 기본 설정됩니다).

그래서 a.mean(dtype=np.float64)을 호출하면 장난감 예제와 더 큰 배열 문제를 해결할 수 있습니다.

Answer 3

신속하고 더러운 답

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

이 1024 * 1024 행렬 예상 된 결과를 제공하지만, 물론이 큰 배열에 대한 사실이되지 않습니다 ...

것이다 평균을 계산하는 경우 귀하의 코드에서 병목 현상이되지 않습니다 나는 파이썬에서 애드혹 알고리즘을 구현할 것입니다. 그러나 세부 사항은 데이터 구조에 달려 있습니다.

평균 계산에 병목 현상이있는 경우 일부 특수 (병렬) 감소 알고리즘으로 문제를 해결할 수 있습니다. 확인 문제를 완화하고 거의 효율적 .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에는 아날로그가 없습니다.