首页 > 学术百科

python实现逐步回归

逐步回归的基本思想是将变量逐个引⼊模型，每引⼊⼀个解释变量后都要进⾏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= py()

# 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)

s((col,col))#np.ones参数为⼀个元组(tuple)

# s((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]

#创建⼀个矩阵来存储⽅差贡献值

s((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,v,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,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.进⾏完矩阵变换就循环上⾯步骤进⾏因⼦的引⼊和剔除

拟合优度检验

以科学发展观为统领再次计算各因⼦的⽅差贡献

前三个未引⼊⽅程的⽅差因⼦进⾏排序，得到第⼀个因⼦的⽅差贡献最⼤，计算第⼀个预报因⼦的F检验值，⼤于临界值引⼊第⼀个预报因⼦进⼊⽅程。

#矩阵转换，计算第⼀步矩阵

r=convert_matrix(r,4)

# print(r)

#计算第⼀步⽅差贡献值

v=get_vari_contri(r)

#print(v)

f=select_factor(r,v,1,1)

print(f)

>####输出>

108.22390933074443

进⾏矩阵变换，计算⽅差贡献

可以看出还没有引⼊⽅程的因⼦2和3，⽅差贡献较⼤的是因⼦2，计算因⼦2的f检验值5.026>3.28，故引⼊预报因⼦2

f=select_factor(r,v,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,v,k,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,v,3,3)

print(f)

#因⼦4的剔除检验值1.863262422188088

f=delete_factor(r,v,4,3)

print(f)

在此对矩阵进⾏变换，计算⽅差贡献，已引⼊因⼦（因⼦1和2）⽅差贡献最⼩的是因⼦1，为引⼊因⼦⽅差贡献最⼤的是因⼦4，计算这两者的引进f检验值和剔除f检验值

#因⼦4的引进检验值1.8632624221880876，⼩于3.28不能引进

f=select_factor(r,v,4,2)

print(f)

#因⼦1的剔除检验值146.52265486251397,⼤于3.28不能剔除

f=delete_factor(r,v,1,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]]

本文发布于:2024-09-22 01:39:59，感谢您对本站的认可！

本文链接：https://www.17tex.com/xueshu/433727.html

上一篇：Matlab-t检验之ttest,ttest2

下一篇：(完整word版)T检验、F检验和统计学意义

标签：计算变量矩阵剔除检验贡献数组

留言与评论（共有 0 条评论）