如何利用python實(shí)現(xiàn)逐步回歸

小編給大家分享一下如何利用python實(shí)現(xiàn)逐步回歸，相信大部分人都還不怎么了解，因此分享這篇文章給大家參考一下，希望大家閱讀完這篇文章后大有收獲，下面讓我們一起去了解一下吧！

逐步回歸的基本思想是將變量逐個(gè)引入模型，每引入一個(gè)解釋變量后都要進(jìn)行F檢驗(yàn)，并對(duì)已經(jīng)選入的解釋變量逐個(gè)進(jìn)行t檢驗(yàn)，當(dāng)原來(lái)引入的解釋變量由于后面解釋變量的引入變得不再顯著時(shí)，則將其刪除。以確保每次引入新的變量之前回歸方程中只包含顯著性變量。這是一個(gè)反復(fù)的過(guò)程，直到既沒(méi)有顯著的解釋變量選入回歸方程，也沒(méi)有不顯著的解釋變量從回歸方程中剔除為止。以保證最后所得到的解釋變量集是最優(yōu)的。

本例的逐步回歸則有所變化，沒(méi)有對(duì)已經(jīng)引入的變量進(jìn)行t檢驗(yàn)，只判斷變量是否引入和變量是否剔除，“雙重檢驗(yàn)”逐步回歸，簡(jiǎn)稱(chēng)逐步回歸。例子的鏈接：（原鏈接已經(jīng)失效），4項(xiàng)自變量，1項(xiàng)因變量。下文不再進(jìn)行數(shù)學(xué)推理，進(jìn)對(duì)計(jì)算過(guò)程進(jìn)行說(shuō)明，對(duì)數(shù)學(xué)理論不明白的可以參考《現(xiàn)代中長(zhǎng)期水文預(yù)報(bào)方法及其應(yīng)用》湯成友，官學(xué)文，張世明著；論文《逐步回歸模型在大壩預(yù)測(cè)中的應(yīng)用》王曉蕾等；

逐步回歸的計(jì)算步驟：

1.計(jì)算第零步增廣矩陣。第零步增廣矩陣是由預(yù)測(cè)因子和預(yù)測(cè)對(duì)象兩兩之間的相關(guān)系數(shù)構(gòu)成的。

2.引進(jìn)因子。在增廣矩陣的基礎(chǔ)上，計(jì)算每個(gè)因子的方差貢獻(xiàn)，挑選出沒(méi)有進(jìn)入方程的因子中方差貢獻(xiàn)最大者對(duì)應(yīng)的因子，計(jì)算該因子的方差比，查F分布表確定該因子是否引入方程。

3.剔除因子。計(jì)算此時(shí)方程中已經(jīng)引入的因子的方差貢獻(xiàn)，挑選出方差貢獻(xiàn)最小的因子,計(jì)算該因子的方差比，查F分布表確定該因子是否從方程中剔除。

4.矩陣變換。將第零步矩陣按照引入方程的因子序號(hào)進(jìn)行矩陣變換，變換后的矩陣再次進(jìn)行引進(jìn)因子和剔除因子的步驟，直到無(wú)因子可以引進(jìn)，也無(wú)因子可以剔除為止，終止逐步回歸分析計(jì)算。

a.以下代碼實(shí)現(xiàn)了數(shù)據(jù)的讀取，相關(guān)系數(shù)的計(jì)算子程序和生成第零步增廣矩陣的子程序。

注意：pandas庫(kù)讀取csv的數(shù)據(jù)結(jié)構(gòu)為DataFrame結(jié)構(gòu)，此處轉(zhuǎn)化為numpy中的（n-dimension array，ndarray）數(shù)組進(jìn)行計(jì)算

import numpy as np
import pandas as pd
#數(shù)據(jù)讀取
#利用pandas讀取csv，讀取的數(shù)據(jù)為DataFrame對(duì)象
data = pd.read_csv('sn.csv')
# 將DataFrame對(duì)象轉(zhuǎn)化為數(shù)組,數(shù)組的最后一列為預(yù)報(bào)對(duì)象
data= data.values.copy()
# print(data)
 
# 計(jì)算回歸系數(shù)，參數(shù)
def get_regre_coef(X,Y):
  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)):
    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)
#構(gòu)建原始增廣矩陣
def get_original_matrix():
  # 創(chuàng)建一個(gè)數(shù)組存儲(chǔ)相關(guān)系數(shù),data.shape幾行(維)幾列，結(jié)果用一個(gè)tuple表示
  # print(data.shape[1])
  col=data.shape[1]
  # print(col)
  r=np.ones((col,col))#np.ones參數(shù)為一個(gè)元組(tuple)
  # print(np.ones((col,col)))
  # for row in data.T:#運(yùn)用數(shù)組的迭代，只能迭代行，迭代轉(zhuǎn)置后的數(shù)組，結(jié)果再進(jìn)行轉(zhuǎn)置就相當(dāng)于迭代了每一列
    # 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.第二部分主要是計(jì)算公差貢獻(xiàn)和方差比。

def get_vari_contri(r):
  col = data.shape[1]
   #創(chuàng)建一個(gè)矩陣來(lái)存儲(chǔ)方差貢獻(xiàn)值
  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
#選擇因子是否進(jìn)入方程，
#參數(shù)說(shuō)明：r為增廣矩陣，v為方差貢獻(xiàn)值，k為方差貢獻(xiàn)值最大的因子下標(biāo),p為當(dāng)前進(jìn)入方程的因子數(shù)
def select_factor(r,v,k,p):
  row=data.shape[0]#樣本容量
  col=data.shape[1]-1#預(yù)報(bào)因子數(shù)
  #計(jì)算方差比
  f=(row-p-2)*v[0,k-1]/(r[col,col]-v[0,k-1])
  # print(calc_vari_contri(r))
  return f

c.第三部分調(diào)用定義的函數(shù)計(jì)算方差貢獻(xiàn)值

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

計(jì)算結(jié)果：

此處沒(méi)有編寫(xiě)判斷方差貢獻(xiàn)最大的子程序，因?yàn)樵谄渌?jì)算中我還需要變量的具體物理含義所以不能單純的由計(jì)算決定對(duì)變量的取舍，此處看出第四個(gè)變量的方查貢獻(xiàn)最大

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

計(jì)算第四個(gè)預(yù)測(cè)因子的方差比(粘貼在了代碼中)，并查F分布表3.280進(jìn)行比對(duì),22.8>3.28,引入第四個(gè)預(yù)報(bào)因子。（前三次不進(jìn)行剔除椅子的計(jì)算）

d.第四部分進(jìn)行矩陣的變換。

#逐步回歸分析與計(jì)算
#通過(guò)矩陣轉(zhuǎn)換公式來(lái)計(jì)算各部分增廣矩陣的元素值
def convert_matrix(r,k):
  col=data.shape[1]
  k=k-1#從第零行開(kāi)始計(jì)數(shù)
  #第k行的元素單不屬于k列的元素
  r1 = np.ones((col, col)) # np.ones參數(shù)為一個(gè)元組(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,j]/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.進(jìn)行完矩陣變換就循環(huán)上面步驟進(jìn)行因子的引入和剔除

再次計(jì)算各因子的方差貢獻(xiàn)

前三個(gè)未引入方程的方差因子進(jìn)行排序，得到第一個(gè)因子的方差貢獻(xiàn)最大，計(jì)算第一個(gè)預(yù)報(bào)因子的F檢驗(yàn)值，大于臨界值引入第一個(gè)預(yù)報(bào)因子進(jìn)入方程。

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

進(jìn)行矩陣變換，計(jì)算方差貢獻(xiàn)

可以看出還沒(méi)有引入方程的因子2和3，方差貢獻(xiàn)較大的是因子2，計(jì)算因子2的f檢驗(yàn)值5.026>3.28,故引入預(yù)報(bào)因子2

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

繼續(xù)進(jìn)行矩陣轉(zhuǎn)換，計(jì)算方差貢獻(xiàn)

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

和因子4的剔除f檢驗(yàn)值1.863，均小于3.28（查F分布表）因子3不能引入，因子4需要剔除，此時(shí)方程中引入的因子數(shù)為2

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

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

#因子4的引進(jìn)檢驗(yàn)值1.8632624221880876，小于3.28不能引進(jìn)
f=select_factor(r,v,4,2)
print(f)
#因子1的剔除檢驗(yàn)值146.52265486251397,大于3.28不能剔除
f=delete_factor(r,v,1,2)
print(f)

不能剔除也不能引進(jìn)變量，此時(shí)停止逐步回歸的計(jì)算。引進(jìn)方程的因子為預(yù)報(bào)因子1和預(yù)報(bào)因子2，借助上一篇博客寫(xiě)的多元回歸。對(duì)進(jìn)入方程的預(yù)報(bào)因子和預(yù)報(bào)對(duì)象進(jìn)行多元回歸。輸出多元回歸的預(yù)測(cè)結(jié)果，一次為常數(shù)項(xiàng)，第一個(gè)因子的預(yù)測(cè)系數(shù)，第二個(gè)因子的預(yù)測(cè)系數(shù)。

#因子1和因子2進(jìn)入方程
#對(duì)進(jìn)入方程的預(yù)報(bào)因子進(jìn)行多元回歸
# regs=LinearRegression()
X=data[:,0:2]
Y=data[:,4]
X=np.mat(np.c_[np.ones(X.shape[0]),X])#為系數(shù)矩陣增加常數(shù)項(xiàng)系數(shù)
Y=np.mat(Y)#數(shù)組轉(zhuǎn)化為矩陣
#print(X)
B=np.linalg.inv(X.T*X)*(X.T)*(Y.T)
print(B.T)#輸出系數(shù),第一項(xiàng)為常數(shù)項(xiàng)，其他為回歸系數(shù)
###輸出##
#[[52.57734888 1.46830574 0.66225049]]

以上是“如何利用python實(shí)現(xiàn)逐步回歸”這篇文章的所有內(nèi)容，感謝各位的閱讀！相信大家都有了一定的了解，希望分享的內(nèi)容對(duì)大家有所幫助，如果還想學(xué)習(xí)更多知識(shí)，歡迎關(guān)注億速云行業(yè)資訊頻道！