求回归系数

雨落森林

关于做一元线性回归用Python去实现，其实有很多方法，什么Scipy库啊、sklearn库都可以做，但是有个问题就是，它只能对两个一维的向量做回归。如果我有很多格点，每个格点在时间方向上有序列，我要做回归，于是就不得不用for循环，但是for循环很慢，所以我就想自己去写一个函数去求回归系数和截距。话不多说，咱们直接看文档和代码。


import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

def linear_regressor(y, x, axis=-1, skipna=False, stat=False):
    # 对于缺失值，选择合适的求平均方法
    mean = np.nanmean if skipna else np.mean
    # 参数axis这里默认最后一维是用来回归的
    # 求出格点数，用来计算自由度n
    ny, nx = y.shape[-1], x.shape[-1]
    n = nx
    # 如果形状不一，则报错
    if y.shape!=x.shape:
        raise ValueError("x and y must keep the same")   
    # 先求平均
    meanx, meany = mean(x, axis=axis,keepdims=True), mean(y, axis=axis,keepdims=True)
    # 求距平
    anox, anoy = x - meanx, y-meany
    # 求均方差(距平的平方求均开根)
    Sx = np.sqrt(mean(anox**2,axis=axis,keepdims=True))
    Sy = np.sqrt(mean(anoy**2,axis=axis,keepdims=True))
    # 标准化(距平除以均方差)
    normx, normy = anox/Sx, anoy/Sy
    # 求相关系数(标准化后的协方差,协方差就是相乘再求均)
    r = mean(normx * normy, axis=axis, keepdims=True)
    # 求回归系数a(相关系数乘以y的均方差除以x的均方差)
    a = r*Sy/Sx
    # 求截距(y平均减去a乘以x平均)
    b = meany - a*meanx
    # 如果需要检验的话
    if stat:
        # 求t统计量
        t = r*np.sqrt(n-2)/np.sqrt(1 - r**2)
        # 求p值
        p = stats.t.cdf(np.abs(t), n-2)
        return a, b, 2-2*p
    else:
        return a, b
if __name__ == "__main__":
    # 测试
    x = np.tile(np.linspace(-10,10,100),[2,1])
    y = 0.5*x + 4 + np.random.uniform(-3,3,[2,100])
    a, b = linear_regressor(y, x)
   
    plt.subplot(2,1,1)
    plt.plot(x[0],y[0],'k.')
    plt.plot(x[0],a[0]*x[0]+b[0],'b-')
    plt.subplot(2,1,2)
    plt.plot(x[1],y[1],'k.')
    plt.plot(x[1],a[1]*x[1]+b[1],'b-')
    plt.show()





##############             2022-3-16        ##############
发现有新需求：
似乎大家都是做空间场相关比较多一点，计算相关系数并且做检验。
所以，我把这个函数改改，返回一个相关系数的场和它的概率。


def correlation(y, x, axis=0, skipna=False, stat=False,keepdims=False):
    # 对于缺失值，选择合适的求平均方法
    mean = np.nanmean if skipna else np.mean
    # 参数axis这里默认第一维是用来回归的
    # 求出格点数，用来计算自由度n
    ny, nx = y.shape[0], x.shape[0]
    n = nx
    # 如果形状不一，则报错
    if y.shape!=x.shape:
        raise ValueError("x and y must keep the same")
    # 先求平均
    meanx, meany = mean(x, axis=axis,keepdims=keepdims), mean(y, axis=axis,keepdims=keepdims)
    # 求距平
    anox, anoy = x - meanx, y-meany
    # 求均方差(距平的平方求均开根)
    Sx = np.sqrt(mean(anox**2,axis=axis,keepdims=keepdims))
    Sy = np.sqrt(mean(anoy**2,axis=axis,keepdims=keepdims))
    # 标准化(距平除以均方差)
    normx, normy = anox/Sx, anoy/Sy
    # 求相关系数(标准化后的协方差,协方差就是相乘再求均)
    r = mean(normx * normy, axis=axis, keepdims=True)
    # 如果需要检验的话
    if stat:
        # 求t统计量
        t = r*np.sqrt(n-2)/np.sqrt(1 - r**2)
        # 求p值
        p = stats.t.cdf(np.abs(t), n-2)
        return r,  2-2*p
    else:
        return r





##############             2022-4-21        ##############
发现了一点问题，检验的时候，公式计算t的时候那个1-r2忘了加根号
另外，求p的时候，要用t的绝对值求p，求得的是累积的概率密度，就相当于单边的、原假设不成立的、概率
所以在返回原假设的概率的时候，是2(1-p)

注：
检验的原理，在于：
原假设：我们假设总体的相关系数ρ=0，就是不相关的意思
然后我们通过样本求得样本的相关系数r
我们用r求一下r在样本数n（也就是自由度为n-2无偏估计下）对应的t值
这个t值取绝对值的累计概率密度，可以理解成单边检验的、总体相关的、概率p
所以要想知道双边的原假设的概率，就是，1-p的两倍

使用这个函数算出来原假设的概率p值之后，如果这个p值小于（一般来说是0.05
那么就可以拒绝原假设，认为求到样本的r或者回归系数能够代表总体的r或者回归系数

puck66

楼主真给力！这样的原创多来点，，

雨落森林

puck66 发表于 2022-1-10 11:16
楼主真给力！这样的原创多来点，，

谢谢捧场

trojeee

您好，我想问一下您在linear_regressor这个函数中是否要求y和x的形状是一样的？因为在计算相关系数中normx*normy这步是对normx和normy的形状有要求的

trojeee

trojeee 发表于 2022-9-5 21:10
您好，我想问一下您在linear_regressor这个函数中是否要求y和x的形状是一样的？因为在计算相关系数中normx* ...

抱歉没看到，x和y是要求形状一样的

雨落森林

trojeee 发表于 2022-9-5 21:39
抱歉没看到，x和y是要求形状一样的

其实可以不一样，一维和三维数组可以计算出二维的相关系数场，但是需要改改代码

紫藤萝

雨落森林 发表于 2022-9-6 17:15
其实可以不一样，一维和三维数组可以计算出二维的相关系数场，但是需要改改代码

你好 请问数组维度不一样时该怎么改代码，谢谢了~数组维度不一样的情况使用的比较多

雨落森林

紫藤萝 发表于 2023-5-8 16:47
你好 请问数组维度不一样时该怎么改代码，谢谢了~数组维度不一样的情况使用的比较多

什么意思，能举个例子吗？比如A数组(40, 12, 30, 20)和B数组(40,)做相关这种吗？还是别的什么？能说明一下数组的shape吗？

紫藤萝

雨落森林 发表于 2023-5-9 11:11
什么意思，能举个例子吗？比如A数组(40, 12, 30, 20)和B数组(40,)做相关这种吗？还是别的什么？能说明一 ...

对 就是一个是一维的时间，另一个是三维的时间+区域。谢谢~~

雨落森林

紫藤萝 发表于 2023-5-15 18:50
对 就是一个是一维的时间，另一个是三维的时间+区域。谢谢~~

指定axis=0，然后那个检验形状是否一致，否则报错的那段删掉。直接用