【问题标题】:Calculating Slopes in Numpy (or Scipy)在 Numpy(或 Scipy)中计算斜率
【发布时间】:2012-03-21 06:21:25
【问题描述】:

我正在尝试找到使用 Numpy 和 Scipy 计算斜率的最快和最有效的方法。我有一个包含三个 Y 变量和一个 X 变量的数据集,我需要计算它们各自的斜率。例如,我可以轻松地一次执行此操作,如下所示,但我希望有一种更有效的方法来执行此操作。我也不认为 linregress 是最好的方法,因为我的结果中不需要任何辅助变量,如截距、标准误差等。任何帮助是极大的赞赏。

    import numpy as np
    from scipy import stats

    Y = [[  2.62710000e+11   3.14454000e+11   3.63609000e+11   4.03196000e+11
        4.21725000e+11   2.86698000e+11   3.32909000e+11   4.01480000e+11
        4.21215000e+11   4.81202000e+11]
        [  3.11612352e+03   3.65968334e+03   4.15442691e+03   4.52470938e+03
        4.65011423e+03   3.10707392e+03   3.54692896e+03   4.20656404e+03
        4.34233412e+03   4.88462501e+03]
        [  2.21536396e+01   2.59098311e+01   2.97401268e+01   3.04784552e+01
        3.13667639e+01   2.76377113e+01   3.27846013e+01   3.73223417e+01
        3.51249997e+01   4.42563658e+01]]
    X = [ 1990.  1991.  1992.  1993.  1994.  1995.  1996.  1997.  1998.  1999.] 
    slope_0, intercept, r_value, p_value, std_err = stats.linregress(X, Y[0,:])
    slope_1, intercept, r_value, p_value, std_err = stats.linregress(X, Y[1,:])
    slope_2, intercept, r_value, p_value, std_err = stats.linregress(X, Y[2,:])
    slope_0 = slope/Y[0,:][0]
    slope_1 = slope/Y[1,:][0]
    slope_2 = slope/Y[2,:][0]
    b, a = polyfit(X, Y[1,:], 1)
    slope_1_a = b/Y[1,:][0]

【问题讨论】:

    标签: python numpy scipy


    【解决方案1】:

    这是视觉示例,如何预测线性回归的系数。 如何为新手计算斜率和截距。快乐学习。

    【讨论】:

    • 您是在谈论“使用 Numpy 和 Scipy 计算斜率的最快和最有效的方法”吗?
    【解决方案2】:

    这取决于您拥有的点数。如果您有两点,请使用linregressstats 中的scipy。如果更多,请使用theilslope,因为它可以避免数据中多达 29% 的异常值并计算最佳斜率。前者只考虑所有样本,不担心异常值,计算出适合所有样本的最佳斜率。

    from scipy import stats
    
    slope1 = stats.linregress([2,4],[1,2])[0] # (ydata,xdata)
    
    slope2 = stats.theilslopes([0.2,0.5,0.9,0.4],[1,2,3,4],0.9) # (ydata,xdata,confidence)
    

    【讨论】:

      【解决方案3】:

      没有 scipy,这个清晰的单行应该足够高效:

      slope = np.polyfit(X,Y,1)[0]
      

      最后你应该得到

      import numpy as np
      
      Y = np.array([
          [  2.62710000e+11, 3.14454000e+11, 3.63609000e+11, 4.03196000e+11, 4.21725000e+11, 2.86698000e+11, 3.32909000e+11, 4.01480000e+11, 4.21215000e+11, 4.81202000e+11],
          [  3.11612352e+03, 3.65968334e+03, 4.15442691e+03, 4.52470938e+03, 4.65011423e+03, 3.10707392e+03, 3.54692896e+03, 4.20656404e+03, 4.34233412e+03, 4.88462501e+03],
          [  2.21536396e+01, 2.59098311e+01, 2.97401268e+01, 3.04784552e+01, 3.13667639e+01, 2.76377113e+01, 3.27846013e+01, 3.73223417e+01, 3.51249997e+01, 4.42563658e+01]]).T
      X = [ 1990,  1991,  1992,  1993,  1994,  1995,  1996,  1997,  1998,  1999] 
      
      print np.polyfit(X,Y,1)[0]
      

      输出为 [1.54983152e+10 9.98749876e+01 1.84564349e+00]

      【讨论】:

        【解决方案4】:

        我在其他答案和原始回归公式的基础上构建了一个适用于任何张量的函数。 它将计算数据沿给定轴的斜率。因此,如果您有任意张量 X[i,j,k,l], Y[i,j,k,l],并且您想知道沿第三轴数据的所有其他轴的斜率,您可以使用 calcSlopes( X, Y, axis = 2 ) 调用它。

        import numpy as np
        
        def calcSlopes( x = None, y = None, axis = -1 ):
            assert x is not None or y is not None
        
            # assume that the given single data argument are equally
            # spaced y-values (like in numpy plot command)
            if y is None:
                y = x
                x = None
        
            # move axis we wanna calc the slopes of to first
            # as is necessary for subtraction of the means
            # note that the axis 'vanishes' anyways, so we don't need to swap it back
            y = np.swapaxes( y, axis, 0 )
            if x is not None:
                x = np.swapaxes( x, axis, 0 )
        
            # https://en.wikipedia.org/wiki/Simple_linear_regression
            # beta = sum_i ( X_i - <X> ) ( Y_i - <Y> ) / ( sum_i ( X_i - <X> )^2 )
            if x is None:
                # axis with values to reduce must be trailing for broadcast_to,
                # therefore transpose
                x = np.broadcast_to( np.arange( y.shape[0] ), y.T.shape ).T
                x = x - ( x.shape[0] - 1 ) / 2. # mean of (0,1,...,n-1) is n*(n-1)/2/n
            else:
                x = x - np.mean( x, axis = 0 )
            y = y - np.mean( y, axis = 0 )
        
            # beta = sum_i x_i y_i / sum_i x_i*^2
            slopes = np.sum( np.multiply( x, y ), axis = 0 ) / np.sum( x**2, axis = 0 )
        
            return slopes
        

        它还具有仅使用等距 y 数据的噱头。比如:

        y = np.array( [
            [ 1, 2, 3, 4 ],
            [ 2, 4, 6, 8 ]
        ] )
        
        print( calcSlopes( y, axis = 0 ) )
        print( calcSlopes( y, axis = 1 ) )
        
        x = np.array( [
            [ 0, 2, 4, 6 ],
            [ 0, 4, 8, 12 ]
        ] )
        
        print( calcSlopes( x, y, axis = 1 ) )
        

        输出:

        [1. 2. 3. 4.]
        [1. 2.]
        [0.5 0.5]
        

        【讨论】:

        • 当你使用 numpy.. 它比 numpy.polyfit 好多少?
        • @kimstik 什么对你来说“更好”?我没有进行性能基准测试。我只是不想在循环中调用 polyfit,我不得不这样做,因为据我所知 polyfit 不能处理任意张量。因此,在这方面,使用纯 numpy 更好,因为它使我能够首先做我想做的事情。
        • “更好”是指“使用 Numpy 和 Scipy 计算斜率的最快和最有效的方法”。在这一点上,张量是题外话。 numpy.polyfit 仍然是纯 numpy。您的方法甚至不需要 numpy 并且可以是纯 python。
        • axis=0 上的除以 0 不受管理。您的示例的矩阵斜率计算可以简化为:“np.mean(np.diff(y,axis=1)/np.diff(x,axis=1),axis=1)”。顺便说一句.. 你有 4 轴张量的测试例子吗?
        【解决方案5】:

        如前所述,您可以使用 scipy 的 linregress。以下是如何让斜率消失:

            from scipy.stats import linregress
        
            x=[1,2,3,4,5]
            y=[2,3,8,9,22]
        
            slope, intercept, r_value, p_value, std_err = linregress(x, y)
            print(slope)
        

        请记住,由于您要计算 r_value 和 p_value 等额外值,因此这样做会比手动仅计算斜率花费更长的时间。但是,Linregress 非常快。

        来源:https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.linregress.html

        【讨论】:

          【解决方案6】:

          我的做法是使用 np.diff() 函数:

          dx = np.diff(xvals),

          dy = np.diff(yvals)

          斜率 = dy/dx

          【讨论】:

          • 你认为 dx 是非零的吗?
          【解决方案7】:

          比公认答案更简单的表示:

          x = np.linspace(0, 10, 11)
          y = np.linspace(0, 20, 11)
          y = np.c_[y, y,y]
          
          X = x - x.mean()
          Y = y - y.mean()
          
          slope = (X.dot(Y)) / (X.dot(X))
          

          斜率方程来自Vector notation for the slope of a line using simple regression

          【讨论】:

          • 我不认为这是正确的!这仅适用于您向 y 广播了 3 次相同的值。问题是默认情况下,平均值是在所有轴上计算的,而不仅仅是一个。所以你要减去所有不同 y 集的平均值,而不是每个 y 集的平均值。
          【解决方案8】:

          最快和最有效的方法是使用来自linregress 的本机 scipy 函数来计算所有内容:

          slope : 回归线的斜率

          intercept : 回归线的截距

          r 值:相关系数

          p 值:假设检验的双边 p 值,其原假设是斜率为零

          stderr : 估计的标准误差

          这是一个例子:

          a = [15, 12, 8, 8, 7, 7, 7, 6, 5, 3]
          b = [10, 25, 17, 11, 13, 17, 20, 13, 9, 15]
          from scipy.stats import linregress
          linregress(a, b)
          

          会还给你:

          LinregressResult(slope=0.20833333333333337, intercept=13.375, rvalue=0.14499815458068521, pvalue=0.68940144811669501, stderr=0.50261704627083648)
          

          附:只是斜率的数学公式:

          【讨论】:

            【解决方案9】:

            线性回归计算在一维上是vector calculation。这意味着我们可以将整个 Y 矩阵上的乘法相结合,然后使用 numpy 中的 axis 参数对拟合进行矢量化。在您的情况下,可以解决以下问题

            ((X*Y).mean(axis=1) - X.mean()*Y.mean(axis=1)) / ((X**2).mean() - (X.mean())**2)
            

            您对拟合质量参数不感兴趣,但大部分都可以通过类似方式获得。

            【讨论】:

              【解决方案10】:

              X 和 Y 的定义方式与您的问题相同,您可以使用:

              dY = (numpy.roll(Y, -1, axis=1) - Y)[:,:-1]
              dX = (numpy.roll(X, -1, axis=0) - X)[:-1]
              
              slopes = dY/dX
              

              numpy.roll() 帮助您将下一个观察结果与当前观察结果对齐,您只需删除最后一列,这是最后一个观察结果和第一个观察结果之间无用的区别。然后您可以一次计算所有斜率,无需 scipy。

              在您的示例中,dX 始终为 1,因此您可以通过计算 slopes = dY 来节省更多时间。

              【讨论】:

              • 我需要澄清一下,因为我只为所有点寻找一个斜率;当您在 X 上运行 Y 的线性回归时会得到什么。例如,slope, intercept = polyfit(X, Y[1,:], 1) 给出的斜率为 99.87。
              • 这为您提供了 Y (3) 中每组数据的斜率。
              猜你喜欢
              • 1970-01-01
              • 2014-01-03
              • 1970-01-01
              • 1970-01-01
              • 2021-03-09
              • 1970-01-01
              • 2012-02-23
              • 1970-01-01
              • 2018-05-11
              相关资源
              最近更新 更多