門限回歸模型(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)# 計算互相關系數，參數為預報因子序列和滯時kdef 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_matrixr_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 -1m=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 datadata=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 FFy=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 if_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,data2data1,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, YX,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()

結果圖：

所用數據：引自《現代中長期水文預報方法及其應用》湯成友 官學文 張世明 著

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實現門限回歸方式就是小編分享給大家的全部內容了，希望能給大家一個參考，也希望大家多多支持好吧啦網。