今天小编就为大家分享一篇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=0
 S_xx=0
 S_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)#滞时为1
 return 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+1
 return -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-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统计量
 return FF
y=data[:,data.shape[1]-1]
FF=get_F(y)
def get_index(FF,element):#获取element在一维数组FF中第一次出现的索引
 i=-1
 for item in FF.flat:
  i+=1
  if 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+1
 data1=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实现门限回归方式就是小编分享给大家的全部内容了,希望能给大家一个参考,也希望大家多多支持爱安网。

最新资讯
微软承认Win10存在严重NTFS漏洞:解压缩包等会损坏硬盘

微软承认Win10存在严

现在微软证实了这个消息,并且表示,会在后续的更新中修复
滴滴在北京已完成46787名司机接种疫苗 预约接种超10万人

滴滴在北京已完成4678

截止1月16日14点,在防疫部门的指导部署下,已完成46787名
高德打车北京升级疫情防控:要求合作平台完成驾驶员全员接种疫苗

高德打车北京升级疫情

为严格执行北京新冠肺炎疫情防控要求,保障广大驾驶员和
Uber计划分拆Postmates旗下快递机器人部门

Uber计划分拆Postmate

据Tech Crunch报道,Uber计划分拆旗下的快递机器人部门P
华为将开设海外最大旗舰店:选址沙特

华为将开设海外最大旗

中国科技公司华为计划在利雅得开设一家旗舰店,这是华为
巴菲特获权证可折价买入2307.69万股传媒公司股票

巴菲特获权证可折价买

伯克希尔·哈撒韦向SEC提交的信息披露报告,巴菲特获权
最新文章
在pycharm中为项目导入anacodna环境的操作方法

在pycharm中为项目导

这篇文章主要介绍了在pycharm中为项目导入anacodna环
tensorflow的ckpt及pb模型持久化方式及转化详解

tensorflow的ckpt及pb

今天小编就为大家分享一篇tensorflow的ckpt及pb模型持
PyTorch笔记之scatter()函数的使用

PyTorch笔记之scatter

这篇文章主要介绍了PyTorch笔记之scatter()函数的使用
python3实现网页版raspberry pi(树莓派)小车控制

python3实现网页版ras

这篇文章主要为大家详细介绍了python3实现网页版raspb
完美解决pycharm导入自己写的py文件爆红问题

完美解决pycharm导入

今天小编就为大家分享一篇完美解决pycharm导入自己写
pycharm内无法import已安装的模块问题解决

pycharm内无法import

今天小编就为大家分享一篇pycharm内无法import已安装