python实现门限回归方式
门限回归模型(ThresholdRegressiveModel，简称TR模型或TRM)的基本思想是通过门限变量的控制作用，当给出预报因子资料后，首先根据门限变量的门限阈值的判别控制作用，以决定不同情况下使用不同的预报方程，从而试图解释各种类似于跳跃和突变的现象。其实质上是把预报问题按状态空间的取值进行分类，用分段的线性回归模式来描述总体非线性预报问题。
多元门限回归的建模步骤就是确实门限变量、率定门限数L、门限值及回归系数的过程，为了计算方便，这里采用二分割（即L=2）说明模型的建模步骤。
基本步骤如下(附代码)：
1.读取数据，计算预报对象与预报因子之间的互相关系数矩阵。
数据读取
#利用pandas读取csv，读取的数据为DataFrame对象
data=pd.read_csv('jl.csv')
#将DataFrame对象转化为数组,数组的第一列为数据序号，最后一列为预报对象，中间各列为预报因子
data=data.values.copy()
#print(data)
#计算互相关系数，参数为预报因子序列和滞时k
defget_regre_coef(X,Y,k):
S_xy=0
S_xx=0
S_yy=0
#计算预报因子和预报对象的均值
X_mean=np.mean(X)
Y_mean=np.mean(Y)
foriinrange(len(X)-k):
S_xy+=(X[i]-X_mean)*(Y[i+k]-Y_mean)
foriinrange(len(X)):
S_xx+=pow(X[i]-X_mean,2)
S_yy+=pow(Y[i]-Y_mean,2)
returnS_xy/pow(S_xx*S_yy,0.5)
#计算相关系数矩阵
defregre_coef_matrix(data):
row=data.shape[1]#列数
r_matrix=np.ones((1,row-2))
#print(row)
foriinrange(1,row-1):
r_matrix[0,i-1]=get_regre_coef(data[:,i],data[:,row-1],1)#滞时为1
returnr_matrix
r_matrix=regre_coef_matrix(data)
#print(r_matrix)
###输出###
#[[0.0489790.078299890.190057050.275012090.28604638]]
2.对相关系数进行排序，相关系数最大的因子作为门限元。
#对相关系数进行排序找到相关系数最大者作为门限元 defget_menxiannum(r_matrix): row=r_matrix.shape[1]#列数 foriinrange(row): ifr_matrix.max()==r_matrix[0,i]: returni+1 return-1 m=get_menxiannum(r_matrix) #print(m) ##输出##第五个因子的互相关系数最大 #5
3.根据选取的门限元因子对数据进行重新排序。
#根据门限元对因子序列进行排序,m为门限变量的序号 defresort_bymenxian(data,m): data=data.tolist()#转化为列表 data.sort(key=lambdax:x[m])#列表按照m+1列进行排序(升序) data=np.array(data) returndata data=resort_bymenxian(data,m)#得到排序后的序列数组
4.将排序后的序列按照门限元分割序列为两段，第一分割第一段1个数据，第二段n-1（n为样本容量）个数据；第二次分割第一段2个数据，第二段n-2个数据，一次类推，分别计算出分割后的F统计量并选出最大统计量对应的门限元的分割点作为门限值。
defget_var(x): returnx.std()**2*x.size#计算总方差 #统计量F的计算,输入数据为按照门限元排序后的预报对象数据 defget_F(Y): col=Y.shape[0]#行数，样本容量 FF=np.ones((1,col-1))#存储不同分割点的统计量 V=get_var(Y)#计算总方差 foriinrange(1,col):#1到col-1 S=get_var(Y[0:i])+get_var(Y[i:col])#计算两段的组内方差和 F=(V-S)*(col-2)/S FF[0,i-1]=F#此步需要判断是否通过F检验，通过了才保留F统计量 returnFF y=data[:,data.shape[1]-1] FF=get_F(y) defget_index(FF,element):#获取element在一维数组FF中第一次出现的索引 i=-1 foriteminFF.flat: i+=1 ifitem==element: returni f_index=get_index(FF,np.max(FF))#获取统计量F的最大索引 #print(data[f_index,m-1])#门限元为第五个因子，代入索引得门限值121
5.以门限值为分割点将数据序列分割为两段，分别进行多元线性回归，此处利用sklearn.linear_model模块中的线性回归模块。再代入预报因子分别计算两段的预测值。
#以门限值为分割点将新data序列分为两部分，分别进行多元回归计算
defdata_excision(data,f_index):
f_index=f_index+1
data1=data[0:f_index,:]
data2=data[f_index:data.shape[0],:]
returndata1,data2
data1,data2=data_excision(data,f_index)
#第一段
defget_XY(data):
#数组切片对变量进行赋值
Y=data[:,data.shape[1]-1]#预报对象位于最后一列
X=data[:,1:data.shape[1]-1]#预报因子从第二列到倒数第二列
returnX,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('Comparisonofpredictedandmeasuredvalues',fontsize=20,fontname='TimesNewRoman')#添加标题
plt.xlabel('Years',color='gray')#添加x轴标签
plt.ylabel('AveragetrafficinDecember',color='gray')#添加y轴标签
plt.legend(['Predictedvalues','Measuredvalues'])#添加图例
plt.show()
结果图：
所用数据：引自《现代中长期水文预报方法及其应用》汤成友官学文张世明著
以上这篇python实现门限回归方式就是小编分享给大家的全部内容了，希望能给大家一个参考，也希望大家多多支持毛票票。
声明：本文内容来源于网络，版权归原作者所有，内容由互联网用户自发贡献自行上传，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任。如果您发现有涉嫌版权的内容，欢迎发送邮件至：czq8825#qq.com（发邮件时，请将#更换为@）进行举报，并提供相关证据，一经查实，本站将立刻删除涉嫌侵权内容。