掩码数组的Python IndexError

Question

我已将数据格式化为3维数组，并且我正在沿着一个轴对每个元素执行回归。以下代码按我的预期工作，并返回回归斜率的列表。

import numpy as np
from scipy import stats
import numpy.ma as ma

#make array
np.random.seed(0)
array = np.random.random((4,3,2))

def regress_slope(array):
  N=array.shape[0]
  alpha=0.9
  y = array[:,:,1]
  x = array[:,:,0] 
  result = [stats.mstats.theilslopes(y[i,...],x[i,...],alpha)[0] for i in range(0,N)]
  return result

result = regress_slope(array)
list(result)
print(result)

我的“真实”数据包括无效值，并且我定义了阈值（<0.1）并尝试从数组中屏蔽这些值。但是，当我使用蒙版数组时，它将引发此错误：

array2 = ma.masked_less_equal(array, 0.1)
result2 = regress_slope(array2)
list(result2)

IndexError: index 0 is out of bounds for axis 0 with size 0

我不确定此错误消息的含义，但我认为可能是因为没有足够的未屏蔽值来计算回归吗？如果是这种情况，如何调整代码以返回结果中的nan？

Answer 1

您是正确的，如果没有足够的未掩码值来计算回归，函数stats.mstats.theilslopes会失败并显示该错误消息。

一个最小的例子：

# this works
a = ma.masked_array([1, 2], mask=[0, 0])
b = ma.masked_array([1, 2], mask=[0, 0])
stats.mstats.theilslopes(a, b, 0.95)

# but this does not
b = ma.masked_array([1, 2], mask=[0, 1])
stats.mstats.theilslopes(a, b, 0.95)

错误消息表明，在计算过程中的某个位置，它试图访问没有元素的结果的第1轴上的第1个元素。

我对您试图做的事情的理论了解不足，以了解结果是否有用，但这将解决您的紧迫问题：

import numpy as np
from scipy import stats
import numpy.ma as ma

np.random.seed(0)
a = np.random.random((4, 3, 2))


def regress_slope(arr):
    def safe_first_theilslopes(arr1, arr2, a):
        try:
            return stats.mstats.theilslopes(arr1, arr2, a)[0]
        except IndexError:
            return np.NaN

    n = arr.shape[0]
    alpha = 0.9
    y = arr[:, :, 1]
    x = arr[:, :, 0]
    return [safe_first_theilslopes(y[i, ...], x[i, ...], alpha) for i in range(0, n)]


result = regress_slope(a)
print(result)

a2 = ma.masked_less_equal(a, 0.1)
result2 = regress_slope(a2)
print(result2)

请注意，我如何让函数返回函数结果的第一个元素（stats.mstats.theilslopes(arr1, arr2, a)[0]）或np.NaN，所以现在将其放入该函数中。

此代码有效，但是会引发一些警告，您可以取消警告，但应该首先考虑一下：

RuntimeWarning: Mean of empty slice.
  out=out, **kwargs)
\lib\site-packages\numpy\core\_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars
  ret = ret.dtype.type(ret / rcount)