贬低和标准化 3d 面具中的 4d 阵列?

Demean and standardise a 4d array within a 3d mask?

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),全部在掩码内。

我这样做正确吗?

有没有更快的方法?

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]

您有几个选项可供选择。

首先是使用masked arrays as you are doing, but provide a proper mask and use the masked functions. Right now, your code is computing all the means and standard deviations, and slapping a mask on the result. To skip masked elements, use np.ma.mean and 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是一个特殊的对象,表示"all the remaining dimensions"。它通常在切片符号中缩写为 ...np.newaxisNone 的别名。结合这些信息,你会发现 [: :, :, np.newaxis] 可以更干净优雅地写成 [..., None]。后者更通用,因为它适用于任意数量的维度。

Numpy 允许负轴索引。 "last axis" 更好的说法通常是 axis=-1.