这不是一个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）

您可以使用dtype关键字参数调用np.mean，该参数指定累加器的类型（对于浮点数组，其默认类型与array相同）。
因此，调用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矩阵的预期结果，但当然对于较大的阵列这将不是真的...
如果计算平均值不会成为代码的瓶颈，我会在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

给出一个更合理的答案需要更多关于数据结构、大小和目标体系结构的信息。