降级和标准化一个4d阵列在一个三维掩膜？

Question

import numpy as np

ts = np.random.rand(40,45,40,1000)
mask = np.random.randint(2, size=(40,45,40),dtype=bool)

#creating a masked array
ts_m = np.ma.array(ts, mask=ts*~mask[:,:,:,np.newaxis])
#demeaning
ts_md = ts_m - ts_m.mean(axis=3)[:,:,:,np.newaxis]
#standardisation
ts_mds = ts_md / ts_md.std(ddof=1,axis=3)[:,:,:,np.newaxis]

我想降低t(沿轴3)，并除以它的标准差(沿轴3)，所有在掩码内。

我做得对吗？

有更快的方法吗？

Answer 1

你有几个选择可供选择。

第一种方法是像您正在做的那样使用蒙面阵列，但是提供一个适当的掩码并使用蒙面函数。现在，您的代码正在计算所有的均值和标准差，并对结果进行掩码。要跳过蒙面元素，请使用np.ma.mean和np.ma.std，从而避免做大量额外的工作。

正如您正确理解的，掩码的大小必须与数据的大小相匹配。当与数据相乘时，就会得到正确的大小，但是在一般情况下，它是昂贵的，并且给出了错误的结果，因为当数据或掩码为零时，掩码将为零。一种更好的方法是创建一个在上一个(新)维度重复的掩码视图。如果首先获得匹配的尾随尺寸，则可以使用np.broadcast_to：

ts = np.random.rand(40, 45, 40, 1000)
mask = np.random.randint(2, size=(40, 45, 40), dtype=np.bool)

#creating a masked array
ts_m = np.ma.array(ts, mask=np.broadcast_to(mask[..., None], ts.shape)
#demeaning
ts_md = ts_m - np.ma.mean(ts_m, axis=3)[..., None]
#standardisation
ts_mds = ts_md / np.ma.std(ts_m, ddof=1,axis=3)[..., None]

掩码是只读的，因为它很可能有一个零步长的维度，有时可以做一些意想不到的事情。这里的广播版本大致相当于

np.lib.stride_tricks.as_strided(mask, ts.shape, (*mask.strides, 0), writeable=False)

这两个版本都为原始数据创建视图，因此非常快速。它们只分配一个指向现有数据的新数组对象，这些数据不被复制。请记住，np.lib.stride_tricks.as_strided是一个大锤，应该非常小心地使用。如果你任由它的话，它随时都会毁了你的解释。

注意：蒙面数组中的掩码被解释为被蒙面的True，而布尔索引数组被解释为False掩码。根据它是如何获得的，以及它在您的真实代码中的意义，您可能想要反转掩码

mask=np.broadcast_to(~mask[..., None], ...)

另一种选择是自己实现掩蔽。有两种方法可以做到这一点。如果您预先这样做，则掩码将应用于数据的前导维度：

ts = np.random.rand(40, 45, 40, 1000)
mask = np.random.randint(2, size=(40, 45, 40), dtype=np.bool)

#creating a masked array
mask = ~mask  # optional, see note above
ts_m = ts[mask]
#demeaning
ts_md = ts_m - ts_m.mean(axis=-1)
#standardisation
ts_mds = ts_md / ts_md.std(ddof=1,axis=-1)
# reshaping
result = np.empty_like(ts)  # alternatively, np.zeros_like
result[mask] = ts_mds

此选项可能比蒙面数组便宜，因为初始掩蔽步骤创建了一个40*45*40-mask_size x 1000数组，并且只在完成后将其替换为结果的掩蔽区域，而不是对完整大小的数据进行操作并保持形状。

第三种选择只有在只有少量的元素被屏蔽的情况下才会真正有用。这基本上就是您的原始代码所做的工作:执行所有的交换，并将掩码应用于结果。

更多提示

Ellipsis是一个特殊的对象，意思是“所有剩余的维度”。在片表示法中，它通常是缩写为...。np.newaxis是None的别名。将这些信息组合在一起，您就可以得到[: :, :, np.newaxis]可以编写得更加简洁和优雅，就像[..., None]一样。后者更一般，因为它适用于任意数量的维度。

Numpy允许负轴指数。一个更好的方式说“最后一轴”一般是axis=-1。

Answer 2

import numpy as np

ts = np.random.rand(40,45,40,1000)
mask = np.random.randint(2, size=(40,45,40)).astype(bool)

#creating a masked array
ts_m = np.ma.array(ts, mask=np.broadcast_to(~mask.reshape(40,45,40,1),ts.shape))
#demeaning
ts_md = ts_m - ts_m.mean(axis=3)[:,:,:,np.newaxis]
#standardisation
ts_mds = ts_md / ts_md.std(ddof=1,axis=3)[:,:,:,np.newaxis]