數(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í)k
def 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_matrix
r_matrix=regre_coef_matrix(data)
# print(r_matrix)
###輸出###
#[[0.048979 0.07829989 0.19005705 0.27501209 0.28604638]]

#對(duì)相關(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 -1
m=get_menxiannum(r_matrix)
# print(m)
##輸出##第五個(gè)因子的互相關(guān)系數(shù)最大
#5

#根據(jù)門限元對(duì)因子序列進(jìn)行排序,m為門限變量的序號(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 data
data=resort_bymenxian(data,m)#得到排序后的序列數(shù)組

def get_var(x):
 return x.std() ** 2 * x.size # 計(jì)算總方差
#統(tǒng)計(jì)量F的計(jì)算,輸入數(shù)據(jù)為按照門限元排序后的預(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#此步需要判斷是否通過F檢驗(yàn)，通過了才保留F統(tǒng)計(jì)量
 return FF
y=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 i
f_index=get_index(FF,np.max(FF))#獲取統(tǒng)計(jì)量F的最大索引
# print(data[f_index,m-1])#門限元為第五個(gè)因子，代入索引得門限值 121

#以門限值為分割點(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,data2
data1,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, Y
X,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)

#恢復(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()