请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 9306|回复: 10

[经验总结] 求回归系数

[复制链接]

新浪微博达人勋

发表于 2022-1-7 18:00:41 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 雨落森林 于 2022-4-22 00:30 编辑

关于做一元线性回归用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或者回归系数





Figure_1.png

20220107.docx

15.85 KB, 下载次数: 11, 下载积分: 金钱 -5

main20220107.py

1.81 KB, 下载次数: 12, 下载积分: 金钱 -5

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2022-1-10 11:16:52 | 显示全部楼层
楼主真给力!这样的原创多来点,,
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2022-1-11 16:21:21 | 显示全部楼层
puck66 发表于 2022-1-10 11:16
楼主真给力!这样的原创多来点,,

谢谢捧场
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-9-5 21:10:07 | 显示全部楼层
您好,我想问一下您在linear_regressor这个函数中是否要求y和x的形状是一样的?因为在计算相关系数中normx*normy这步是对normx和normy的形状有要求的
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-9-5 21:39:49 | 显示全部楼层
trojeee 发表于 2022-9-5 21:10
您好,我想问一下您在linear_regressor这个函数中是否要求y和x的形状是一样的?因为在计算相关系数中normx* ...

抱歉没看到,x和y是要求形状一样的
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2022-9-6 17:15:20 | 显示全部楼层
trojeee 发表于 2022-9-5 21:39
抱歉没看到,x和y是要求形状一样的

其实可以不一样,一维和三维数组可以计算出二维的相关系数场,但是需要改改代码
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-5-8 16:47:35 | 显示全部楼层
雨落森林 发表于 2022-9-6 17:15
其实可以不一样,一维和三维数组可以计算出二维的相关系数场,但是需要改改代码

你好 请问数组维度不一样时该怎么改代码,谢谢了~数组维度不一样的情况使用的比较多
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2023-5-9 11:11:18 | 显示全部楼层
紫藤萝 发表于 2023-5-8 16:47
你好 请问数组维度不一样时该怎么改代码,谢谢了~数组维度不一样的情况使用的比较多

什么意思,能举个例子吗?比如A数组(40, 12, 30, 20)和B数组(40,)做相关这种吗?还是别的什么?能说明一下数组的shape吗?
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-5-15 18:50:43 | 显示全部楼层
雨落森林 发表于 2023-5-9 11:11
什么意思,能举个例子吗?比如A数组(40, 12, 30, 20)和B数组(40,)做相关这种吗?还是别的什么?能说明一 ...

对 就是一个是一维的时间,另一个是三维的时间+区域。谢谢~~
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2023-5-16 09:33:52 | 显示全部楼层
紫藤萝 发表于 2023-5-15 18:50
对 就是一个是一维的时间,另一个是三维的时间+区域。谢谢~~

指定axis=0,然后那个检验形状是否一致,否则报错的那段删掉。直接用
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表