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

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

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

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

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

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

#對相關(guān)系數(shù)進(jìn)行排序找到相關(guān)系數(shù)最大者作為門限元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)##輸出##第五個因子的互相關(guān)系數(shù)最大#5

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

#根據(jù)門限元對因子序列進(jìn)行排序,m為門限變量的序號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.將排序后的序列按照門限元分割序列為兩段，第一分割第一段1個數(shù)據(jù)，第二段n-1（n為樣本容量）個數(shù)據(jù)；第二次分割第一段2個數(shù)據(jù)，第二段n-2個數(shù)據(jù)，一次類推，分別計(jì)算出分割后的F統(tǒng)計(jì)量并選出最大統(tǒng)計(jì)量對應(yīng)的門限元的分割點(diǎn)作為門限值。

def get_var(x): return x.std() ** 2 * x.size # 計(jì)算總方差#統(tǒng)計(jì)量F的計(jì)算,輸入數(shù)據(jù)為按照門限元排序后的預(yù)報對象數(shù)據(jù)def get_F(Y): col=Y.shape[0]#行數(shù)，樣本容量 FF=np.ones((1,col-1))#存儲不同分割點(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#此步需要判斷是否通過F檢驗(yàn)，通過了才保留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])#門限元為第五個因子，代入索引得門限值 121

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

#以門限值為分割點(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ù)組切片對變量進(jìn)行賦值 Y = data[:, data.shape[1] - 1] # 預(yù)報對象位于最后一列 X = data[:, 1:data.shape[1] - 1]#預(yù)報因子從第二列到倒數(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ù)測值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ù)測值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ù)測值和實(shí)際值按照年份序號從新排序，恢復(fù)其順序，利用matplotlib模塊做出預(yù)測值與實(shí)際值得對比圖。

#恢復(fù)順序Y=resort_bymenxian(Y,0)# print(Y.shape)# 預(yù)測結(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)代中長期水文預(yù)報方法及其應(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)門限回歸方式就是小編分享給大家的全部內(nèi)容了，希望能給大家一個參考，也希望大家多多支持好吧啦網(wǎng)。