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实现一个图片拆分工具,可以根据需要的行数和列数进行拆分,感兴趣的小伙伴可以跟随小编一起学习一下... 简单介绍先自己选择输入的图片,默认是输出到项目文件夹中,可以自己选择其他的文件夹,选择需要拆分的行数和列数,可以通过

Python中反转字符串的常见方法小结

《Python中反转字符串的常见方法小结》在Python中,字符串对象没有内置的反转方法,然而,在实际开发中,我们经常会遇到需要反转字符串的场景,比如处理回文字符串、文本加密等,因此,掌握如何在Pyt... 目录python中反转字符串的方法技术背景实现步骤1. 使用切片2. 使用 reversed() 函

Python中将嵌套列表扁平化的多种实现方法

《Python中将嵌套列表扁平化的多种实现方法》在Python编程中,我们常常会遇到需要将嵌套列表(即列表中包含列表)转换为一个一维的扁平列表的需求,本文将给大家介绍了多种实现这一目标的方法,需要的朋... 目录python中将嵌套列表扁平化的方法技术背景实现步骤1. 使用嵌套列表推导式2. 使用itert

使用Docker构建Python Flask程序的详细教程

《使用Docker构建PythonFlask程序的详细教程》在当今的软件开发领域,容器化技术正变得越来越流行,而Docker无疑是其中的佼佼者,本文我们就来聊聊如何使用Docker构建一个简单的Py... 目录引言一、准备工作二、创建 Flask 应用程序三、创建 dockerfile四、构建 Docker

Python使用vllm处理多模态数据的预处理技巧

《Python使用vllm处理多模态数据的预处理技巧》本文深入探讨了在Python环境下使用vLLM处理多模态数据的预处理技巧,我们将从基础概念出发,详细讲解文本、图像、音频等多模态数据的预处理方法,... 目录1. 背景介绍1.1 目的和范围1.2 预期读者1.3 文档结构概述1.4 术语表1.4.1 核

Python使用pip工具实现包自动更新的多种方法

《Python使用pip工具实现包自动更新的多种方法》本文深入探讨了使用Python的pip工具实现包自动更新的各种方法和技术,我们将从基础概念开始,逐步介绍手动更新方法、自动化脚本编写、结合CI/C... 目录1. 背景介绍1.1 目的和范围1.2 预期读者1.3 文档结构概述1.4 术语表1.4.1 核

在Linux中改变echo输出颜色的实现方法

《在Linux中改变echo输出颜色的实现方法》在Linux系统的命令行环境下,为了使输出信息更加清晰、突出,便于用户快速识别和区分不同类型的信息,常常需要改变echo命令的输出颜色,所以本文给大家介... 目python录在linux中改变echo输出颜色的方法技术背景实现步骤使用ANSI转义码使用tpu

Conda与Python venv虚拟环境的区别与使用方法详解

《Conda与Pythonvenv虚拟环境的区别与使用方法详解》随着Python社区的成长,虚拟环境的概念和技术也在不断发展,:本文主要介绍Conda与Pythonvenv虚拟环境的区别与使用... 目录前言一、Conda 与 python venv 的核心区别1. Conda 的特点2. Python v

Python使用python-can实现合并BLF文件

《Python使用python-can实现合并BLF文件》python-can库是Python生态中专注于CAN总线通信与数据处理的强大工具,本文将使用python-can为BLF文件合并提供高效灵活... 目录一、python-can 库:CAN 数据处理的利器二、BLF 文件合并核心代码解析1. 基础合

Python使用OpenCV实现获取视频时长的小工具

《Python使用OpenCV实现获取视频时长的小工具》在处理视频数据时,获取视频的时长是一项常见且基础的需求,本文将详细介绍如何使用Python和OpenCV获取视频时长,并对每一行代码进行深入解析... 目录一、代码实现二、代码解析1. 导入 OpenCV 库2. 定义获取视频时长的函数3. 打开视频文