門(mén)限回歸模型(Threshold Regressive Model，簡(jiǎn)稱(chēng)TR模型或TRM)的基本思想是通過(guò)門(mén)限變量的控制作用，當(dāng)給出預(yù)報(bào)因子資料后，首先根據(jù)門(mén)限變量的門(mén)限閾值的判別控制作用，以決定不同情況下使用不同的預(yù)報(bào)方程，從而試圖解釋各種類(lèi)似于跳躍和突變的現(xiàn)象。其實(shí)質(zhì)上是把預(yù)報(bào)問(wèn)題按狀態(tài)空間的取值進(jìn)行分類(lèi)，用分段的線性回歸模式來(lái)描述總體非線性預(yù)報(bào)問(wèn)題。

多元門(mén)限回歸的建模步驟就是確實(shí)門(mén)限變量、率定門(mén)限數(shù)L、門(mén)限值及回歸系數(shù)的過(guò)程，為了計(jì)算方便，這里采用二分割（即L=2）說(shuō)明模型的建模步驟。

基本步驟如下(附代碼)：

1.讀取數(shù)據(jù)，計(jì)算預(yù)報(bào)對(duì)象與預(yù)報(bào)因子之間的互相關(guān)系數(shù)矩陣。

數(shù)據(jù)讀取#利用pandas讀取csv，讀取的數(shù)據(jù)為DataFrame對(duì)象data = pd.read_csv(’jl.csv’)# 將DataFrame對(duì)象轉(zhuǎn)化為數(shù)組,數(shù)組的第一列為數(shù)據(jù)序號(hào)，最后一列為預(yù)報(bào)對(duì)象，中間各列為預(yù)報(bào)因子data= data.values.copy()# print(data)# 計(jì)算互相關(guān)系數(shù)，參數(shù)為預(yù)報(bào)因子序列和滯時(shí)kdef get_regre_coef(X,Y,k): S_xy=0 S_xx=0 S_yy=0 # 計(jì)算預(yù)報(bào)因子和預(yù)報(bào)對(duì)象的均值 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)#計(jì)算相關(guān)系數(shù)矩陣def regre_coef_matrix(data): row=data.shape[1]#列數(shù) 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)#滯時(shí)為1 return r_matrixr_matrix=regre_coef_matrix(data)# print(r_matrix)###輸出####[[0.048979 0.07829989 0.19005705 0.27501209 0.28604638]]

2.對(duì)相關(guān)系數(shù)進(jìn)行排序，相關(guān)系數(shù)最大的因子作為門(mén)限元。

#對(duì)相關(guān)系數(shù)進(jìn)行排序找到相關(guān)系數(shù)最大者作為門(mén)限元def get_menxiannum(r_matrix): row=r_matrix.shape[1]#列數(shù) for i in range(row): if r_matrix.max()==r_matrix[0,i]: return i+1 return -1m=get_menxiannum(r_matrix)# print(m)##輸出##第五個(gè)因子的互相關(guān)系數(shù)最大#5

3.根據(jù)選取的門(mén)限元因子對(duì)數(shù)據(jù)進(jìn)行重新排序。

#根據(jù)門(mén)限元對(duì)因子序列進(jìn)行排序,m為門(mén)限變量的序號(hào)def resort_bymenxian(data,m): data=data.tolist()#轉(zhuǎn)化為列表 data.sort(key=lambda x: x[m])#列表按照m+1列進(jìn)行排序(升序) data=np.array(data) return datadata=resort_bymenxian(data,m)#得到排序后的序列數(shù)組

4.將排序后的序列按照門(mén)限元分割序列為兩段，第一分割第一段1個(gè)數(shù)據(jù)，第二段n-1（n為樣本容量）個(gè)數(shù)據(jù)；第二次分割第一段2個(gè)數(shù)據(jù)，第二段n-2個(gè)數(shù)據(jù)，一次類(lèi)推，分別計(jì)算出分割后的F統(tǒng)計(jì)量并選出最大統(tǒng)計(jì)量對(duì)應(yīng)的門(mén)限元的分割點(diǎn)作為門(mén)限值。

def get_var(x): return x.std() ** 2 * x.size # 計(jì)算總方差#統(tǒng)計(jì)量F的計(jì)算,輸入數(shù)據(jù)為按照門(mén)限元排序后的預(yù)報(bào)對(duì)象數(shù)據(jù)def get_F(Y): col=Y.shape[0]#行數(shù)，樣本容量 FF=np.ones((1,col-1))#存儲(chǔ)不同分割點(diǎn)的統(tǒng)計(jì)量 V=get_var(Y)#計(jì)算總方差 for i in range(1,col):#1到col-1 S=get_var(Y[0:i])+get_var(Y[i:col])#計(jì)算兩段的組內(nèi)方差和 F=(V-S)*(col-2)/S FF[0,i-1]=F#此步需要判斷是否通過(guò)F檢驗(yàn)，通過(guò)了才保留F統(tǒng)計(jì)量 return FFy=data[:,data.shape[1]-1]FF=get_F(y)def get_index(FF,element):#獲取element在一維數(shù)組FF中第一次出現(xiàn)的索引 i=-1 for item in FF.flat: i+=1 if item==element: return if_index=get_index(FF,np.max(FF))#獲取統(tǒng)計(jì)量F的最大索引# print(data[f_index,m-1])#門(mén)限元為第五個(gè)因子，代入索引得門(mén)限值 121

5.以門(mén)限值為分割點(diǎn)將數(shù)據(jù)序列分割為兩段，分別進(jìn)行多元線性回歸，此處利用sklearn.linear_model模塊中的線性回歸模塊。再代入預(yù)報(bào)因子分別計(jì)算兩段的預(yù)測(cè)值。

#以門(mén)限值為分割點(diǎn)將新data序列分為兩部分，分別進(jìn)行多元回歸計(jì)算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,data2data1,data2=data_excision(data,f_index)# 第一段def get_XY(data): # 數(shù)組切片對(duì)變量進(jìn)行賦值 Y = data[:, data.shape[1] - 1] # 預(yù)報(bào)對(duì)象位于最后一列 X = data[:, 1:data.shape[1] - 1]#預(yù)報(bào)因子從第二列到倒數(shù)第二列 return X, YX,Y=get_XY(data1)regs=LinearRegression()regs.fit(X,Y)# print(’第一段’)# print(regs.coef_)#輸出回歸系數(shù)# print(regs.score(X,Y))#輸出相關(guān)系數(shù)#計(jì)算預(yù)測(cè)值Y1=regs.predict(X)# print(’第二段’)X,Y=get_XY(data2)regs.fit(X,Y)# print(regs.coef_)#輸出回歸系數(shù)# print(regs.score(X,Y))#輸出相關(guān)系數(shù)#計(jì)算預(yù)測(cè)值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.將預(yù)測(cè)值和實(shí)際值按照年份序號(hào)從新排序，恢復(fù)其順序，利用matplotlib模塊做出預(yù)測(cè)值與實(shí)際值得對(duì)比圖。

#恢復(fù)順序Y=resort_bymenxian(Y,0)# print(Y.shape)# 預(yù)測(cè)結(jié)果可視化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’)#添加標(biāo)題plt.xlabel(’Years’,color=’gray’)#添加x軸標(biāo)簽plt.ylabel(’Average traffic in December’,color=’gray’)#添加y軸標(biāo)簽plt.legend([’Predicted values’,’Measured values’])#添加圖例plt.show()

結(jié)果圖：

所用數(shù)據(jù)：引自《現(xiàn)代中長(zhǎng)期水文預(yù)報(bào)方法及其應(yīng)用》湯成友 官學(xué)文 張世明 著

num x1 x2 x3 x4 x5 y 1960 308 301 352 310 149 80.5 1961 182 186 165 127 70 42.9 1962 195 134 134 97 61 43.9 1963 136 378 334 307 148 87.4 1964 230 630 332 161 100 66.6 1965 225 333 209 365 152 82.9 1966 296 225 317 527 228 111 1967 324 229 176 317 153 79.3 1968 278 230 352 317 143 82 1969 662 442 453 381 188 103 1970 187 136 103 129 74.7 43 1971 284 404 600 327 161 92.2 1972 427 430 843 448 236 144 1973 258 404 639 275 156 98.9 1974 113 160 128 177 77.2 50.1 1975 143 300 333 214 106 63 1976 113 74 193 241 107 58.6 1977 204 140 154 90 55.1 40.2 1978 174 445 351 267 120 70.3 1979 93 95 197 214 94.9 64.3 1980 214 250 354 385 178 73 1981 232 676 483 218 113 72.6 1982 266 216 146 112 82.8 61.4 1983 210 433 803 301 166 115 1984 261 702 512 291 153 97.5 1985 197 178 238 180 94.2 58.9 1986 442 256 623 310 146 84.3 1987 136 99 253 232 114 62 1988 256 226 185 321 151 80.1 1989 473 409 300 298 141 79.6 1990 277 291 639 302 149 84.6 1991 372 181 174 104 68.8 58.4 1992 251 142 126 95 59.4 51.4 1993 181 125 130 240 121 64 1994 253 278 216 182 124 82.4 1995 168 214 265 175 101 68.1 1996 98.8 97 92.7 88 56.7 45.6 1997 252 385 313 270 119 78.8 1998 242 198 137 114 71.9 51.8 1999 268 178 127 109 68.6 53.3 2000 86.2 286 233 133 77.8 58.6 2001 150 168 122 93 62.8 42.9 2002 180 150 97.8 78 48.2 41.9 2003 166 203 166 124 70 53.7 2004 400 202 126 158 92.7 54.7 2005 79.8 82.6 129 160 76.6 53.7

以上這篇python實(shí)現(xiàn)門(mén)限回歸方式就是小編分享給大家的全部?jī)?nèi)容了，希望能給大家一個(gè)參考，也希望大家多多支持好吧啦網(wǎng)。