逐步迴歸的基本思想是將變數逐個引入模型，每引入一個解釋變數後都要進行F檢驗，並對已經選入的解釋變數逐個進行t檢驗，當原來引入的解釋變數由於後面解釋變數的引入變得不再顯著時，則將其刪除。以確保每次引入新的變數之前回歸方程中只包含顯著性變數。這是一個反覆的過程，直到既沒有顯著的解釋變數選入迴歸方程，也沒有不顯著的解釋變數從迴歸方程中剔除為止。以保證最後所得到的解釋變數集是最優的。

本例的逐步迴歸則有所變化，沒有對已經引入的變數進行t檢驗，只判斷變數是否引入和變數是否剔除，“雙重檢驗”逐步迴歸，簡稱逐步迴歸。例子的連結：（原連結已經失效），4項自變數，1項因變數。下文不再進行數學推理，進對計算過程進行說明，對數學理論不明白的可以參考《現代中長期水文預報方法及其應用》湯成友，官學文，張世明著；論文《逐步迴歸模型在大壩預測中的應用》王曉蕾等；

逐步迴歸的計算步驟：

1.計算第零步增廣矩陣。第零步增廣矩陣是由預測因子和預測物件兩兩之間的相關係數構成的。

2.引進因子。在增廣矩陣的基礎上，計算每個因子的方差貢獻，挑選出沒有進入方程的因子中方差貢獻最大者對應的因子，計算該因子的方差比，查F分佈表確定該因子是否引入方程。

3.剔除因子。計算此時方程中已經引入的因子的方差貢獻，挑選出方差貢獻最小的因子,計算該因子的方差比，查F分佈表確定該因子是否從方程中剔除。

4.矩陣變換。將第零步矩陣按照引入方程的因子序號進行矩陣變換，變換後的矩陣再次進行引進因子和剔除因子的步驟，直到無因子可以引進，也無因子可以剔除為止，終止逐步迴歸分析計算。

a.以下程式碼實現了資料的讀取，相關係數的計運算元程式和生成第零步增廣矩陣的子程式。

注意：pandas庫讀取csv的資料結構為DataFrame結構，此處轉化為numpy中的（n-dimension array，ndarray）陣列進行計算

import numpy as np
import pandas as pd
#資料讀取
#利用pandas讀取csv，讀取的資料為DataFrame物件
data = pd.read_csv('sn.csv')
# 將DataFrame物件轉化為陣列,陣列的最後一列為預報物件
data= data.values.copy()
# print(data)
 
# 計算迴歸係數，引數
def get_regre_coef(X,Y):
  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)):
    S_xy += (X[i] - X_mean) * (Y[i] - Y_mean)
    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 get_original_matrix():
  # 建立一個數組儲存相關係數,data.shape幾行(維)幾列，結果用一個tuple表示
  # print(data.shape[1])
  col=data.shape[1]
  # print(col)
  r=np.ones((col,col))#np.ones引數為一個元組(tuple)
  # print(np.ones((col,col)))
  # for row in data.T:#運用陣列的迭代，只能迭代行，迭代轉置後的陣列，結果再進行轉置就相當於迭代了每一列
    # print(row.T)
  for i in range(col):
    for j in range(col):
      r[i,j]=get_regre_coef(data[:,i],data[:,j])
  return r

b.第二部分主要是計算公差貢獻和方差比。

def get_vari_contri(r):
  col = data.shape[1]
   #建立一個矩陣來儲存方差貢獻值
  v=np.ones((1,col-1))
  # print(v)
  for i in range(col-1):
    # v[0,i]=pow(r[i,col-1],2)/r[i,i]
    v[0,i] = pow(r[i,col - 1],2) / r[i,i]
  return v
#選擇因子是否進入方程，
#引數說明：r為增廣矩陣，v為方差貢獻值，k為方差貢獻值最大的因子下標,p為當前進入方程的因子數
def select_factor(r,v,k,p):
  row=data.shape[0]#樣本容量
  col=data.shape[1]-1#預報因子數
  #計算方差比
  f=(row-p-2)*v[0,k-1]/(r[col,col]-v[0,k-1])
  # print(calc_vari_contri(r))
  return f

c.第三部分呼叫定義的函式計算方差貢獻值

#計算第零步增廣矩陣
r=get_original_matrix()
# print(r)
#計算方差貢獻值
v=get_vari_contri(r)
print(v)
#計算方差比

計算結果：

此處沒有編寫判斷方差貢獻最大的子程式，因為在其他計算中我還需要變數的具體物理含義所以不能單純的由計算決定對變數的取捨，此處看出第四個變數的方查貢獻最大

# #計算方差比
# print(data.shape[0])
f=select_factor(r,4,0)
print(f)
#######輸出##########
22.79852020138227

計算第四個預測因子的方差比(貼上在了程式碼中)，並查F分佈表3.280進行比對,22.8>3.28,引入第四個預報因子。（前三次不進行剔除椅子的計算）

d.第四部分進行矩陣的變換。

#逐步迴歸分析與計算
#通過矩陣轉換公式來計算各部分增廣矩陣的元素值
def convert_matrix(r,k):
  col=data.shape[1]
  k=k-1#從第零行開始計數
  #第k行的元素單不屬於k列的元素
  r1 = np.ones((col,col)) # np.ones引數為一個元組(tuple)
  for i in range(col):
    for j in range(col):
      if (i==k and j!=k):
        r1[i,j]=r[k,j]/r[k,k]
      elif (i!=k and j!=k):
        r1[i,j]=r[i,j]-r[i,k]*r[k,k]
      elif (i!= k and j== k):
        r1[i,j] = -r[i,k]/r[k,k]
      else:
        r1[i,j] = 1/r[k,k]
  return r1

e.進行完矩陣變換就迴圈上面步驟進行因子的引入和剔除

再次計算各因子的方差貢獻

前三個未引入方程的方差因子進行排序，得到第一個因子的方差貢獻最大，計算第一個預報因子的F檢驗值，大於臨界值引入第一個預報因子進入方程。

#矩陣轉換，計算第一步矩陣
r=convert_matrix(r,4)
# print(r)
#計算第一步方差貢獻值
v=get_vari_contri(r)
#print(v)
f=select_factor(r,1,1)
print(f)
#########輸出#####
108.22390933074443

進行矩陣變換，計算方差貢獻

可以看出還沒有引入方程的因子2和3，方差貢獻較大的是因子2，計算因子2的f檢驗值5.026>3.28,故引入預報因子2

f=select_factor(r,2,2)
print(f)
##########輸出#########
5.025864648951804

繼續進行矩陣轉換，計算方差貢獻

這一步需要考慮剔除因子了，有方差貢獻可以知道，已引入方程的因子中方差貢獻最小的是因子4，分別計算因子3的引進f檢驗值0.0183

和因子4的剔除f檢驗值1.863，均小於3.28（查F分佈表）因子3不能引入，因子4需要剔除，此時方程中引入的因子數為2

#選擇是否剔除因子，
#引數說明：r為增廣矩陣，v為方差貢獻值，k為方差貢獻值最大的因子下標,t為當前進入方程的因子數
def delete_factor(r,t):
  row = data.shape[0] # 樣本容量
  col = data.shape[1] - 1 # 預報因子數
  # 計算方差比
  f = (row - t - 1) * v[0,k - 1] / r[col,col]
  # print(calc_vari_contri(r))
  return f
#因子3的引進檢驗值0.018233473487350636
f=select_factor(r,3,3)
print(f)
#因子4的剔除檢驗值1.863262422188088
f=delete_factor(r,3)
print(f)

在此對矩陣進行變換，計算方差貢獻 ，已引入因子（因子1和2）方差貢獻最小的是因子1，為引入因子方差貢獻最大的是因子4，計算這兩者的引進f檢驗值和剔除f檢驗值

#因子4的引進檢驗值1.8632624221880876，小於3.28不能引進
f=select_factor(r,2)
print(f)
#因子1的剔除檢驗值146.52265486251397,大於3.28不能剔除
f=delete_factor(r,2)
print(f)

不能剔除也不能引進變數，此時停止逐步迴歸的計算。引進方程的因子為預報因子1和預報因子2，藉助上一篇部落格寫的多元迴歸。對進入方程的預報因子和預報物件進行多元迴歸。輸出多元迴歸的預測結果，一次為常數項，第一個因子的預測係數，第二個因子的預測係數。

#因子1和因子2進入方程
#對進入方程的預報因子進行多元迴歸
# regs=LinearRegression()
X=data[:,0:2]
Y=data[:,4]
X=np.mat(np.c_[np.ones(X.shape[0]),X])#為係數矩陣增加常數項係數
Y=np.mat(Y)#陣列轉化為矩陣
#print(X)
B=np.linalg.inv(X.T*X)*(X.T)*(Y.T)
print(B.T)#輸出係數,第一項為常數項，其他為迴歸係數
###輸出##
#[[52.57734888 1.46830574 0.66225049]]

以上這篇利用python實現逐步迴歸就是小編分享給大家的全部內容了，希望能給大家一個參考，也希望大家多多支援我們。