一元线性回归-逐像素判断-多组同时运行-矩阵运算

目录

• 问题引入
• 线性回归函数
• 函数调用

问题引入

一元线性回归再简单不过了，实现的方式多种多样。调用scikit-learn linear_model.LinearRegression()、Scipy.polyfit( ) 或 numpy.polyfit( )、stats.linregress( )、optimize.curve_fit( )、numpy.linalg.lstsq、Statsmodels.OLS ( )使用矩阵求逆方法的解析解、高中数学讲的最小二乘法公式，详细请看博客：Python环境下的8种简单线性回归算法 。

但如果数据Y(y1,y2,y3,...yn)中出现随机位置与个数的无效值（0值），只有将对应的X(x1,x2,x3,...xn)与Y中的无效值（0值）剔除，得到的Y随X变化趋势才是准确的。对于较少组Y与X的回归，可以采用for循环一组一组的判断计算。但是当有数以亿计组Y与X，用for循环则显得效率低下。遇到这一问题，是由于我在计算多年NDVI序列趋势，每一年的NDVI会出现随机位置与个数的无效值（0值）。

经过探索，在健忘主义这篇博客启发下，找到了上述问题解决办法。下面分享编程实现思路

线性回归函数

通过改写健忘主义这篇博客，利用numpy 数组提供的求和平方等函数，设置操作在axis=0上对列进行求和等运算。通过数组运算提高操作效率。
函数输入数组X，Y需要提前进行处理，在下节函数调用中说明。

def lrg_ols(X,Y): 
    '''
    calculate slope and intercept by ols
    refer to https://blog.csdn.net/qq_35515661/article/details/84727471

    Parameters
    ----------
    X : 
        TYPE
            n*m ndarray
        DESCRIPTION.
            each column stands for a independent variable series, eg. x1,x2,x3,...,xn.
            which "x1,x2,x3,...,xn" is 1,2, 3 ...n insequence. X can be created by 
            X = np.arange(1,n+1).reshape(n,1)
            X = np.broadcast_to(X,(n,m))            
    Y : 
        TYPE
            n*m ndarray
        DESCRIPTION. 
            each column stands for a dependent variable series, eg. y1 y2 y3...yn.
            note the Y's dimensions must be the same as X's

    Returns
    -------
    k : TYPE
            1*m ndarray
        DESCRIPTION.
            corresponding slopes
        
    b : TYPE
            1*m ndarray
        DESCRIPTION.
            corresponding intercept
    '''
    
    X_size = np.count_nonzero(X, axis=0) # X_size should be equal Y_size
    # XYproduct = np.multiply(X,Y) #element-wise produt of X and Y
    # XXproduct = np.multiply(X,X)
    X_mean = np.sum(X,axis=0).reshape(1, -1) / X_size
    Y_mean = np.sum(Y,axis=0).reshape(1, -1) / X_size
    zi = np.sum(X*Y,axis=0).reshape(1, -1) - X_size * X_mean * Y_mean
    mu = np.sum(X*X,axis=0).reshape(1, -1) - X_size * X_mean * X_mean
    k = zi / mu
    b = Y_mean - k * X_mean
    
    # 计算决定系数
    Y_pred = X*k + b
    Y_pred[np.where(X==0)]=0   
    SSR = np.sum(np.square(Y_pred - Y_mean),axis=0).reshape(1, -1) # 回归平方和
    SSE = np.sum(np.square(Y - Y_pred),axis=0).reshape(1, -1) # 残差平方和
    SST = SSR + SSE # 总偏差平方和
    r2 = SSR / SST    
    return k,b,r2,X_size

函数调用

输入数组X，Y。因变量Y整理成二维数组，每一列是一个因变量时序；自变量整理成与X维度相同的数组，每一列是一个自变量数组。# 然后进行数据筛选，将因变量无效值处赋值为0，并将对应自变量处赋值为0。

shp = data.shape        
# reshape后成为每一列是一个时序
data = data.reshape(shp[0], -1)
# 构造自变量序列，1，2 ... n,
X = np.arange(1,shp[0]+1).reshape(shp[0],1)
X = np.repeat(X,shp[1]*shp[2],axis=1)
# 数据筛选，将负值赋值为0
cond = np.where(data<0)
data[cond]=0
X[cond]=0
k,b,r2,xsize = lrg_ols(X,data)