python实现门限回归

2023-10-20 05:10
文章标签 python 实现 回归 门限

本文主要是介绍python实现门限回归,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

门限回归模型(Threshold Regressive Model,简称TR模型或TRM)的基本思想是通过门限变量的控制作用,当给出预报因子资料后,首先根据门限变量的门限阈值的判别控制作用,以决定不同情况下使用不同的预报方程,从而试图解释各种类似于跳跃和突变的现象。其实质上是把预报问题按状态空间的取值进行分类,用分段的线性回归模式来描述总体非线性预报问题。多元门限回归的建模步骤就是确实门限变量、率定门限数L、门限值及回归系数的过程,为了计算方便,这里采用二分割(即L=2)说明模型的建模步骤。

基本步骤如下(附代码):

1.读取数据,计算预报对象与预报因子之间的互相关系数矩阵。

数据读取
#利用pandas读取csv,读取的数据为DataFrame对象
data = pd.read_csv('jl.csv')
# 将DataFrame对象转化为数组,数组的第一列为数据序号,最后一列为预报对象,中间各列为预报因子
data= data.values.copy()
# print(data)
# 计算互相关系数,参数为预报因子序列和滞时k
def get_regre_coef(X,Y,k):S_xy=0S_xx=0S_yy=0# 计算预报因子和预报对象的均值X_mean = np.mean(X)Y_mean = np.mean(Y)for i in  range(len(X)-k):S_xy += (X[i] - X_mean) * (Y[i+k] - Y_mean)for i in range(len(X)):S_xx += pow(X[i] - X_mean, 2)S_yy += pow(Y[i] - Y_mean, 2)return S_xy/pow(S_xx*S_yy,0.5)
#计算相关系数矩阵
def regre_coef_matrix(data):row=data.shape[1]#列数r_matrix=np.ones((1,row-2))# print(row)for i in range(1,row-1):r_matrix[0,i-1]=get_regre_coef(data[:,i],data[:,row-1],1)#滞时为1return r_matrix
r_matrix=regre_coef_matrix(data)
# print(r_matrix)
###输出###
#[[0.048979   0.07829989 0.19005705 0.27501209 0.28604638]]

 2.对相关系数进行排序,相关系数最大的因子作为门限元。

#对相关系数进行排序找到相关系数最大者作为门限元
def get_menxiannum(r_matrix):row=r_matrix.shape[1]#列数for i in range(row):if r_matrix.max()==r_matrix[0,i]:return i+1return -1
m=get_menxiannum(r_matrix)
# print(m)
##输出##第五个因子的互相关系数最大
#5

3.根据选取的门限元因子对数据进行重新排序。

#根据门限元对因子序列进行排序,m为门限变量的序号
def resort_bymenxian(data,m):data=data.tolist()#转化为列表data.sort(key=lambda x: x[m])#列表按照m+1列进行排序(升序)data=np.array(data)return data
data=resort_bymenxian(data,m)#得到排序后的序列数组

4.将排序后的序列按照门限元分割序列为两段,第一分割第一段1个数据,第二段n-1(n为样本容量)个数据;第二次分割第一段2个数据,第二段n-2个数据,一次类推,分别计算出分割后的F统计量并选出最大统计量对应的门限元的分割点作为门限值。

def get_var(x):return x.std() ** 2 * x.size  # 计算总方差
#统计量F的计算,输入数据为按照门限元排序后的预报对象数据
def get_F(Y):col=Y.shape[0]#行数,样本容量FF=np.ones((1,col-1))#存储不同分割点的统计量V=get_var(Y)#计算总方差for i in range(1,col):#1到col-1S=get_var(Y[0:i])+get_var(Y[i:col])#计算两段的组内方差和F=(V-S)*(col-2)/SFF[0,i-1]=F#此步需要判断是否通过F检验,通过了才保留F统计量return FF
y=data[:,data.shape[1]-1]
FF=get_F(y)
def get_index(FF,element):#获取element在一维数组FF中第一次出现的索引i=-1for item in FF.flat:i+=1if item==element:return i
f_index=get_index(FF,np.max(FF))#获取统计量F的最大索引
# print(data[f_index,m-1])#门限元为第五个因子,代入索引得门限值 121

5.以门限值为分割点将数据序列分割为两段,分别进行多元线性回归,此处利用sklearn.linear_model模块中的线性回归模块。再代入预报因子分别计算两段的预测值。

#以门限值为分割点将新data序列分为两部分,分别进行多元回归计算
def data_excision(data,f_index):f_index=f_index+1data1=data[0:f_index,:]data2=data[f_index:data.shape[0],:]return data1,data2
data1,data2=data_excision(data,f_index)
# 第一段
def get_XY(data):# 数组切片对变量进行赋值Y = data[:, data.shape[1] - 1]  # 预报对象位于最后一列X = data[:, 1:data.shape[1] - 1]#预报因子从第二列到倒数第二列return X, Y
X,Y=get_XY(data1)
regs=LinearRegression()
regs.fit(X,Y)
# print('第一段')
# print(regs.coef_)#输出回归系数
# print(regs.score(X,Y))#输出相关系数
#计算预测值
Y1=regs.predict(X)
# print('第二段')
X,Y=get_XY(data2)
regs.fit(X,Y)
# print(regs.coef_)#输出回归系数
# print(regs.score(X,Y))#输出相关系数
#计算预测值
Y2=regs.predict(X)
Y=np.column_stack((data[:,0],np.hstack((Y1,Y2)))).copy()
Y=np.column_stack((Y,data[:,data.shape[1]-1]))
Y=resort_bymenxian(Y,0)

6.将预测值和实际值按照年份序号从新排序,恢复其顺序,利用matplotlib模块做出预测值与实际值得对比图。

#恢复顺序
Y=resort_bymenxian(Y,0)
# print(Y.shape)
# 预测结果可视化
plt.plot(Y[:,0],Y[:,1],'b--',Y[:,0],Y[:,2],'g')
plt.title('Comparison of predicted and measured values',fontsize=20,fontname='Times New Roman')#添加标题
plt.xlabel('Years',color='gray')#添加x轴标签
plt.ylabel('Average traffic in December',color='gray')#添加y轴标签
plt.legend(['Predicted values','Measured values'])#添加图例
plt.show()

结果图:

所用数据:引自《现代中长期水文预报方法及其应用》汤成友 官学文 张世明 著

numx1x2x3x4x5y
196030830135231014980.5
19611821861651277042.9
1962195134134976143.9
196313637833430714887.4
196423063033216110066.6
196522533320936515282.9
1966296225317527228111
196732422917631715379.3
196827823035231714382
1969662442453381188103
197018713610312974.743
197128440460032716192.2
1972427430843448236144
197325840463927515698.9
197411316012817777.250.1
197514330033321410663
19761137419324110758.6
19772041401549055.140.2
197817444535126712070.3
1979939519721494.964.3
198021425035438517873
198123267648321811372.6
198226621614611282.861.4
1983210433803301166115
198426170251229115397.5
198519717823818094.258.9
198644225662331014684.3
19871369925323211462
198825622618532115180.1
198947340930029814179.6
199027729163930214984.6
199137218117410468.858.4
19922511421269559.451.4
199318112513024012164
199425327821618212482.4
199516821426517510168.1
199698.89792.78856.745.6
199725238531327011978.8
199824219813711471.951.8
199926817812710968.653.3
200086.228623313377.858.6
20011501681229362.842.9
200218015097.87848.241.9
20031662031661247053.7
200440020212615892.754.7
200579.882.612916076.653.7

 

这篇关于python实现门限回归的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/244888

相关文章

Python实现终端清屏的几种方式详解

《Python实现终端清屏的几种方式详解》在使用Python进行终端交互式编程时,我们经常需要清空当前终端屏幕的内容,本文为大家整理了几种常见的实现方法,有需要的小伙伴可以参考下... 目录方法一:使用 `os` 模块调用系统命令方法二:使用 `subprocess` 模块执行命令方法三:打印多个换行符模拟

SpringBoot+EasyPOI轻松实现Excel和Word导出PDF

《SpringBoot+EasyPOI轻松实现Excel和Word导出PDF》在企业级开发中,将Excel和Word文档导出为PDF是常见需求,本文将结合​​EasyPOI和​​Aspose系列工具实... 目录一、环境准备与依赖配置1.1 方案选型1.2 依赖配置(商业库方案)二、Excel 导出 PDF

Python实现MQTT通信的示例代码

《Python实现MQTT通信的示例代码》本文主要介绍了Python实现MQTT通信的示例代码,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下面随着小编来一... 目录1. 安装paho-mqtt库‌2. 搭建MQTT代理服务器(Broker)‌‌3. pytho

基于Python开发一个图像水印批量添加工具

《基于Python开发一个图像水印批量添加工具》在当今数字化内容爆炸式增长的时代,图像版权保护已成为创作者和企业的核心需求,本方案将详细介绍一个基于PythonPIL库的工业级图像水印解决方案,有需要... 目录一、系统架构设计1.1 整体处理流程1.2 类结构设计(扩展版本)二、核心算法深入解析2.1 自

使用zip4j实现Java中的ZIP文件加密压缩的操作方法

《使用zip4j实现Java中的ZIP文件加密压缩的操作方法》本文介绍如何通过Maven集成zip4j1.3.2库创建带密码保护的ZIP文件,涵盖依赖配置、代码示例及加密原理,确保数据安全性,感兴趣的... 目录1. zip4j库介绍和版本1.1 zip4j库概述1.2 zip4j的版本演变1.3 zip4

从入门到进阶讲解Python自动化Playwright实战指南

《从入门到进阶讲解Python自动化Playwright实战指南》Playwright是针对Python语言的纯自动化工具,它可以通过单个API自动执行Chromium,Firefox和WebKit... 目录Playwright 简介核心优势安装步骤观点与案例结合Playwright 核心功能从零开始学习

Python 字典 (Dictionary)使用详解

《Python字典(Dictionary)使用详解》字典是python中最重要,最常用的数据结构之一,它提供了高效的键值对存储和查找能力,:本文主要介绍Python字典(Dictionary)... 目录字典1.基本特性2.创建字典3.访问元素4.修改字典5.删除元素6.字典遍历7.字典的高级特性默认字典

Python自动化批量重命名与整理文件系统

《Python自动化批量重命名与整理文件系统》这篇文章主要为大家详细介绍了如何使用Python实现一个强大的文件批量重命名与整理工具,帮助开发者自动化这一繁琐过程,有需要的小伙伴可以了解下... 目录简介环境准备项目功能概述代码详细解析1. 导入必要的库2. 配置参数设置3. 创建日志系统4. 安全文件名处

使用Python构建一个高效的日志处理系统

《使用Python构建一个高效的日志处理系统》这篇文章主要为大家详细讲解了如何使用Python开发一个专业的日志分析工具,能够自动化处理、分析和可视化各类日志文件,大幅提升运维效率,需要的可以了解下... 目录环境准备工具功能概述完整代码实现代码深度解析1. 类设计与初始化2. 日志解析核心逻辑3. 文件处

python生成随机唯一id的几种实现方法

《python生成随机唯一id的几种实现方法》在Python中生成随机唯一ID有多种方法,根据不同的需求场景可以选择最适合的方案,文中通过示例代码介绍的非常详细,需要的朋友们下面随着小编来一起学习学习... 目录方法 1:使用 UUID 模块(推荐)方法 2:使用 Secrets 模块(安全敏感场景)方法