如何使用逐元素操作获取多个 numpy 保存的数组的均值和标准

Posted 2023-04-14

【中文标题】如何使用逐元素操作获取多个 numpy 保存的数组的均值和标准

【英文标题】：How to get Mean & Std of multiple numpy saved arrays using element-wise operation

【发布时间】：2019-01-03 18:21:43

【问题描述】：

我有一个包含 1000 个 numpy 压缩文件 (npz) 的文件夹，这些文件代表数据模拟的结果。每个文件有两个数组a 和b，具有相同的维度、形状、数据类型。我想要的最终输出是a、b 和c（我在下面的示例中创建）的元素均值和标准差数组，考虑到所有模拟，即：

mean_a = np.mean(a1,a2,a3,...a1000)

std_a = np.std(a1,a2,a3...a1000) 等。

我设法获得了平均值，但没有使用直接的元素操作。我最挣扎的是得到性病。我试图将所有数组附加到列表中，但我遇到了内存错误的问题。知道我该如何进行吗？请参阅下面我到目前为止所取得的成就。提前致谢！！

import glob
import numpy as np
import os 

simulation_runs = 10
simulation_range = np.arange(simulation_runs)

npFiles = [npFile for npFile in glob.iglob(os.path.join(outDir, "sc0*.npz"))]

a_accum = np.empty([885, 854], dtype=np.float32)
b_accum = np.empty([885, 854], dtype=np.float32)    
c_accum = np.empty([885, 854], dtype=np.float32)    

for run, i in enumerate(npFiles):
    npData = np.load(i)
    a = npData['scc'] 
    b = npData['bcc']
    c = a+b
    a_accum  = a + a_accum
    b_accum = b + b_accum   
    c_accum = c + b_accum   

aMean = a_accum/len(simulation_range)
bMean= b_accum/len(simulation_range)
cMean = c_accum/len(simulation_range)

【问题讨论】：

既然你用的是numpy，那你调查过np.mean()和np.std()吗？

是的，我想应该这样做，但是我遇到了如何将所有 1000 个数组放在一起以同时使用 np.std() 和 np.mean 的问题()。我试图将每个数组附加到一个列表中，但它出现了内存错误。因此，也许有一种方法可以在循环内甚至是另一种方式来做到这一点。我对编程很陌生，所以我无法想象其他选择 =/

【参考方案1】：

首先，如果您可以 (ssh) 访问具有更多内存的机器，那是最简单的。也许你甚至可以在没有一个的情况下进行管理。 885*854*(1000 次模拟)*(4 个字节/float32) = 2.8 GiB，所以如果你分别做 a、b 和 c，你应该在一台合理的机器上有足够的内存。在这种情况下，只需将它们放入一个数组中，然后使用 np.mean 和 np.std：

a = np.zeros((1000,885,854), dtype=np.float32)
for run, i in enumerate(npFiles):
    a[i]=np.load(run)['scc']
amean = a.mean(axis=0)
astd = a.std(axis=0)

对于 b 和 c 也是如此。

否则，最优雅的选择是以易于延迟加载的格式保存数据。 dask 是专门为此设计的，但可能需要一些时间来学习（但从长远来看可能是值得的）。您也可以将其存储在 netcat 文件中，并使用xarray 作为dask 的一种前端，甚至可能更方便。

如果你只需要平均值，标准，你可以手动完成。标准的公式是

std = sqrt(mean(abs(x - x.mean())**2))

因此，既然您已经有了方法，那么该过程将与您已经做过的非常相似：（未经测试）

import numpy as np
import os 

simulation_runs = 10
simulation_range = np.arange(simulation_runs)

npFiles = [npFile for npFile in glob.iglob(os.path.join(outDir, "sc0*.npz"))]

a_accum = np.empty([885, 854], dtype=np.float32)
b_accum = np.empty([885, 854], dtype=np.float32)    
c_accum = np.empty([885, 854], dtype=np.float32)    

for run, i in enumerate(npFiles):
    npData = np.load(i)
    a = npData['scc'] 
    b = npData['bcc']
    c = a+b
    a_accum  = a + a_accum
    b_accum = b + b_accum   
    c_accum = c + b_accum   

aMean = a_accum/len(simulation_range)
bMean= b_accum/len(simulation_range)
cMean = c_accum/len(simulation_range)


a_sumsq = np.empty([885, 854], dtype=np.float32)
b_sumsq = np.empty([885, 854], dtype=np.float32)    
c_sumsq = np.empty([885, 854], dtype=np.float32)    

for run, i in enumerate(npFiles):
    npData = np.load(i)
    a = npData['scc'] 
    b = npData['bcc']
    c = a+b
    a_sumsq += (a-aMean)**2
    b_sumsq += (b-bMean)**2
    c_sumsq += (c-cMean)**2

a_std = np.sqrt(a_sumsq/(len(npFiles)-1)) # The -1 is to get an unbiased estimator
b_std = np.sqrt(b_sumsq/(len(npFiles)-1))
c_std = np.sqrt(c_sumsq/(len(npFiles)-1))

【讨论】：

谢谢！我正在努力让它发挥作用，但我仍然只是平均水平。当解释器读取带有 a_sumsq += (a-aMean)**2... 的行时，它会返回：''FloatingPointError: 在添加中遇到无效值'''。知道如何解决吗？

您的数据中是否有任何 NaN？如果您尝试在错误之前打印所有变量怎么办，那里有什么问题吗？

当我开始计算每次运行的标准时，会显示 NaN 值。计算平均值时没有 Nan。关于第一种方法...我不得不将运行次数从 1000 更改为 10000，所以我尝试并想出了内存错误，不幸的是

在不查看数据的情况下，我无法帮助您弄清楚 NaN 的来源。但是，如果第一种方法有效，您可以一次取 1000 个文件中的 10 个块，将它们填充到某个数组中，然后计算这些块的标准差的均方根。这将与实际标准相同。

只是一个更新：我已经设法通过用 numpy.zeros 替换 numpy.empty 来摆脱 NaN（不知道为什么）......所以我成功地使用了你的第二个解决方案！再次感谢