0 
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]

我想贬低 ts（沿轴 3），并除以它的标准偏差（沿轴 3），所有这些都在掩码内。

我这样做正确吗？

有更快的方法吗？

4

2 回答 2

1

您有几个可用的选项。

第一种是像您一样使用掩码数组，但提供适当的掩码并使用掩码函数。现在，您的代码正在计算所有均值和标准差，并在结果上加上掩码。要跳过被屏蔽的元素，请使用np.ma.meanand 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.

于 2018-07-26T15:23:31.757 回答
0 
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]
于 2018-07-26T14:56:56.367 回答