門限迴歸模型(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()

結果圖：

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