【问题标题】:What is the best way to take np.percentile along an axis ignoring nans?沿着忽略 nans 的轴获取 np.percentile 的最佳方法是什么?
【发布时间】:2017-08-30 10:21:44
【问题描述】:

对于包含 NaN 值的数据,是否有一种相当快速的方法来处理 np.percentile(ndarr, axis=0)

对于np.median,有对应的bottleneck.nanmedianhttps://pypi.python.org/pypi/Bottleneck)就不错了。

我想出的百分位数最好的是:

   from bottleneck import nanrankdata, nanmax, nanargmin
   def nanpercentile(x, q, axis):
       ranks = nanrankdata(x, axis=axis)
       peak = nanmax(ranks, axis=axis)
       pct = ranks/peak / 100. # to make a percentile
       wh = nanargmin(abs(pct-q),axis=axis)
       return x[wh]

这不起作用;真正需要的是某种方法来沿axis 获取第 n 个元素,但我还没有找到 numpy 切片技巧来做到这一点。

“相当快”意味着比循环遍历索引更好,例如:

q = 40
x = np.array([[[1,2,3],[6,np.nan,4]],[[0.5,2,1],[9,3,np.nan]]])
out = np.empty(x.shape[:-1])
for i in range(x.shape[0]):
   for j in range(x.shape[1]):
      d = x[i,j,:]
      out[i,j] = np.percentile(d[np.isfinite(d)], q)

print out

#array([[ 1.8,  4.8],
#       [ 0.9,  5.4]])

这可行,但可能非常慢。

np.ma 似乎没有按预期工作;它将nan 值视为inf

xm = np.ma.masked_where(np.isnan(x),x)
print np.percentile(xm,40,axis=2)

# array([[ 1.8,  5.6],
#        [ 0.9,  7.8]])

【问题讨论】:

    标签: python numpy


    【解决方案1】:

    【讨论】:

      【解决方案2】:

      您可以使用在numpy.lib.stride_tricks 中找到的as_strided() 来操纵数组的步幅以更快地迭代它。

      您的计算可以被视为在阵列的 (1,1,3) 窗口上运行。我喜欢使用一个通用函数(sliding_window() 使用as_strided() 在 n 个窗口中创建 n。我在这里找到了它 - Efficient Overlapping Windows with Numpy;该函数的功劳显然归功于 johnvinyard。那个博客页面很好地描述了正在发生的事情。

      制作一些 1x1x3 的窗口

      import numpy as np
      x = np.array([[[1,2,3],[6,np.nan,4]],[[0.5,2,1],[9,3,np.nan]]])
      for thing in sliding_window(x, (1,1,3)):
          print thing
      
      # [ 1.  2.  3.]
      # [  6.  nan   4.]
      # [ 0.5  2.   1. ]
      # [  9.   3.  nan]
      

      应用 ```np.percentile()'' - 忽略 NaN

      for thing in sliding_window(x, (1,1,3)):
          print np.percentile(thing[np.isfinite(thing)], 40)
      
      # 1.8
      # 4.8
      # 0.9
      # 5.4
      

      制作一个结果数组:

      per_s = [np.percentile(thing[np.isfinite(thing)], 40)
               for thing in sliding_window(x, (1,1,3))]
      
      print per_s
      # [1.8, 4.8000000000000007, 0.90000000000000002, 5.4000000000000004]
      
      per_s = np.array(per_s)
      print per_s
      # array([ 1.8,  4.8,  0.9,  5.4])
      

      让它恢复到您期望的形状

      print per_s.reshape((2,2))
      # array([[ 1.8,  4.8],
      #        [ 0.9,  5.4]])
      
      print per_s.reshape(x.shape[:-1])
      # array([[ 1.8,  4.8],
      #        [ 0.9,  5.4]])
      

      这应该更快。我很好奇它是否会 - 我没有任何现实世界问题来测试它。

      对 numpy as_strided 的谷歌搜索出现了一些不错的结果:我有这个书签,http://scipy-lectures.github.io/advanced/advanced_numpy/

      sliding_window() 来自Efficient Overlapping Windows with Numpy

      from numpy.lib.stride_tricks import as_strided as ast
      from itertools import product
      
      def norm_shape(shape):
          '''
          Normalize numpy array shapes so they're always expressed as a tuple, 
          even for one-dimensional shapes.
      
          Parameters
              shape - an int, or a tuple of ints
      
          Returns
              a shape tuple
          '''
          try:
              i = int(shape)
              return (i,)
          except TypeError:
              # shape was not a number
              pass
      
          try:
              t = tuple(shape)
              return t
          except TypeError:
              # shape was not iterable
              pass
      
          raise TypeError('shape must be an int, or a tuple of ints')
      
      
      def sliding_window(a,ws,ss = None,flatten = True):
          '''
          Return a sliding window over a in any number of dimensions
      
          Parameters:
              a  - an n-dimensional numpy array
              ws - an int (a is 1D) or tuple (a is 2D or greater) representing the size 
                   of each dimension of the window
              ss - an int (a is 1D) or tuple (a is 2D or greater) representing the 
                   amount to slide the window in each dimension. If not specified, it
                   defaults to ws.
              flatten - if True, all slices are flattened, otherwise, there is an 
                        extra dimension for each dimension of the input.
      
          Returns
              an array containing each n-dimensional window from a
          '''
      
          if None is ss:
              # ss was not provided. the windows will not overlap in any direction.
              ss = ws
          ws = norm_shape(ws)
          ss = norm_shape(ss)
      
          # convert ws, ss, and a.shape to numpy arrays so that we can do math in every 
          # dimension at once.
          ws = np.array(ws)
          ss = np.array(ss)
          shape = np.array(a.shape)
      
      
          # ensure that ws, ss, and a.shape all have the same number of dimensions
          ls = [len(shape),len(ws),len(ss)]
          if 1 != len(set(ls)):
              raise ValueError(\
              'a.shape, ws and ss must all have the same length. They were %s' % str(ls))
      
          # ensure that ws is smaller than a in every dimension
          if np.any(ws > shape):
              raise ValueError('ws cannot be larger than a in any dimension. a.shape was %s and ws was %s' % (str(a.shape),str(ws)))
      
          # how many slices will there be in each dimension?
          newshape = norm_shape(((shape - ws) // ss) + 1)
          # the shape of the strided array will be the number of slices in each dimension
          # plus the shape of the window (tuple addition)
          newshape += norm_shape(ws)
          # the strides tuple will be the array's strides multiplied by step size, plus
          # the array's strides (tuple addition)
          newstrides = norm_shape(np.array(a.strides) * ss) + a.strides
          strided = ast(a,shape = newshape,strides = newstrides)
          if not flatten:
              return strided
      
          # Collapse strided so that it has one more dimension than the window.  I.e.,
          # the new array is a flat list of slices.
          meat = len(ws) if ws.shape else 0
          firstdim = (np.product(newshape[:-meat]),) if ws.shape else ()
          dim = firstdim + (newshape[-meat:])
          # remove any dimensions with size 1
          #dim = filter(lambda i : i != 1,dim)
          dim = tuple(thing for thing in dim if thing != 1)
          return strided.reshape(dim)
      

      【讨论】:

        【解决方案3】:

        如果您不需要超快速的解决方案,您可以先将您的数组传输到 pandas DataFrame 并进行分位数,然后再返回到 numpy 数组。

        df = pd.DataFrame(array.T).quantile()
        arr = np.array(df)
        

        【讨论】:

          【解决方案4】:

          您可以在 numpy 1.8 中使用 partition() 沿轴获取第 n 个元素,这里是获取沿最后一个轴的第二个元素的代码:

          x = np.array([[[1,2,3],[6,np.nan,4]],[[0.5,2,1],[9,3,np.nan]]])
          np.partition(x, 1)[..., 1]
          

          输出:

          array([[ 2.,  6.],
                 [ 1.,  9.]])
          

          【讨论】:

          • 有趣,我使用的是 1.7.1 - 返回什么?
          • 如果您希望n'th 元素沿着n 不恒定的轴工作,这是否可行?此外,partition 在内存中创建了一个副本,这是不可取的,但可能是不可避免的。
          猜你喜欢
          • 2015-05-10
          • 1970-01-01
          • 2020-10-08
          • 1970-01-01
          • 1970-01-01
          • 2010-10-11
          • 2016-11-20
          • 1970-01-01
          • 2019-02-22
          相关资源
          最近更新 更多