线性回归Statsmodels模型报告(重在Statsmodels.OLS.summary())

写作原因:最近看了一下Statsmodels.OLS,即用Statsmodels使用最小二乘法获得线性回归的系数、截距,主要有一个model.summary(),其中有一些参数想深入弄明白,将学习结果分享:

如果用python,有很多种方法实现线性回归(带不带常数项截距都无所谓):

从计算原理上来分:一般经常使用正规方程(最小二乘法)和梯度下降。

比如博客园这篇文章:机器学习算法的Python实现(一):线性回归 - Geek-Thomas - 博客园 (cnblogs.com)

说了八种方法,其本质还是上述两种方法的不同库的实现;

本文内容:

这里只用最小二乘法,毕竟咱们目的是研究模型报告,但是过程还是要走一遍做做样子,写三种方法的计算实现,带上代码和对比思考:

1.sklearn.linear_model ,里面各种线性回归的库,具体计算有最小二乘、梯度下降、批量、随机梯度下降等;

2.numpy.linalg,也是最小二乘法,相当于借助numpy逐步实现(也算半个手动挡)

3.Statsmodels.OLS,也是最小二乘法,其中的summary()模型报告是本文的重点;

当然还有纯手动挡自己写方法计算的。

一、回顾最小二乘法

原理是通过矩阵形式,求MSE关于系数矩阵的偏导,令其为0,求解:

推导公式比较简单,但是敲latex太长嫌麻烦,此处从上链接那里截个图,该公式随便都能搜到。

梯度下降原理略,本文不讲这个。

二、构造数据与使用四种方法具体过程

我们考虑一个带常数项、二个系数、噪音的数据,因为三个维度可以画出来。

2.1生成数据及展示

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import time
import random
from collections import Counter
import itertools

np.set_printoptions(suppress=True)
plt.rcParams['font.sans-serif'] = ['SimHei']	# 显示中文
plt.rcParams['axes.unicode_minus'] = False		# 显示负号
np.random.seed(1999)
# 生成数据
x1 = np.linspace(1,30,num=100)
x2 = np.linspace(-10,10,num=100)
w1 = 2
w2 = 3
noise = np.random.normal(0,1,size=100)
y= x1*w1 + x2*w2 + noise +5
y_dream = x1*w1 + x2*w2+5
# print(y)
# 先画个3D图
ax = plt.subplot(111,projection='3d')
ax.scatter(x1,x2,y,color='navy',s=10,label='生成的数据点')
ax.plot(x1,x2,y_dream,color='red',lw=2,label='理想的拟合曲线(非最佳)')
ax.set_xlabel('x1',fontsize=20)
ax.set_ylabel('x2',fontsize=20)
ax.set_zlabel('y',fontsize=20)
plt.title('原数据3D展示图')
plt.legend(loc='upper left')
plt.show()

差不多数据就是这样,噪音稍微弄小点,一开始设置标准差设得太大,发现对初始的点影响极大;由于存在噪音,所以我们认为的这个红色的线,极大可能并不是最好的。

2.2sklearn.linear_model

from sklearn.linear_model import LinearRegression
modle_1 = LinearRegression(fit_intercept=True)
X=np.c_[x1,x2]
modle_1.fit(X,y)
coef_1 = modle_1.coef_
intercept_1 = modle_1.intercept_
print('系数是:',coef_1)
print('截距是:',intercept_1)

#系数是: [2.75724514 1.90154838]
#截距是: -6.760579051034362

2.3 手算

此时我们要手动添加一个截距项,系数全部为1,可知我们把系数放在了后面,注意要用np.linalg.pinv,如果用np.linalg.inv可能在矩阵奇异的情况下,解不出来或者与不太准确!

X_intercept = np.c_[x1,x2,np.ones(100)]
# 对应了最小二乘法的几个部分
theta = np.linalg.pinv(X_intercept.T.dot(X_intercept)).dot(X_intercept.T).dot(y)
print(theta)

# array([ 2.58969124,  2.89297533, 17.43027338])

2.4 statsmodels OLS方法

OLS means:(ordinary least square)

import statsmodels.api as sm
import statsmodels.formula.api as smf
# 我们使用自带的添加截距系数方法add_constant
X_new= sm.add_constant(X)
X_new[:10]

看来这个方法,把截距项的系数放在了前面 

model_3 = sm.OLS(y,X_new)
result = model_3.fit()
print(result.params)
# 前面是截距,后面两个是系数
# [-0.08619059  2.32663944  2.52592665]

2.5 三个方法比较

我们构建的公式y=2*x_1 + 3*x_2 + 5 +noise

但由于噪音干扰,实际上最佳的曲线(MSE最小)已经不是原本的曲线了。

首先看系数和截距:

# 三组系数和截距
print('Sklearn前面系数,最后截距',coef_1,intercept_1)
print('Numpy.linalg前面系数,最后截距',theta)
print('statsmodels.OLS','切片调下位置同上',result.params[1:],result.params[0])
Sklearn前面系数,最后截距 [2.75724514 1.90154838] -6.760579051034362
Numpy.linalg前面系数,最后截距 [ 2.32663944  2.52592665 -0.08619059]
statsmodels.OLS 切片调下位置同上 [2.32663944 2.52592665] -0.08619059391555109

可以看到numpy.linalg和statsmodels.OLS方法结果一样,说明都是矩阵求逆再乘两下得出的结果,sklearn的默认纯线性回归(不带正则项的),里面没有solver可选(岭回归Ridge就有),看来十有八九计算实现并不是用的正规方程,毕竟正规方程在系数多的情况下,很可能求不出来。

那么sklearn会不会结果比他们两个差呢?

residual_df = pd.DataFrame()
residual_df['y_real'] = y

# sklearn学习原数据
y_sklean = modle_1.predict(X)
residual_df['y_sklean'] =y_sklean

# 手算原数据
y_manual = x1*theta[0] + x2*theta[1] +  theta[2]
residual_df['y_manual'] =y_manual
# X_intercept.dot(theta) # 或者用矩阵乘法

# OLS模型计算结果
y_ols = X_new.dot(result.params)
residual_df['y_ols'] =y_ols
residual_df
# 或者用result.fittedvalues

乍一看,结果貌似完全一样,我们仔细看看:

residual_df.iloc[0,1],residual_df.iloc[0,2],residual_df.iloc[0,3]

# (-23.01881766462413, -23.018817664624127, -23.018817664624144)

结果的很多位后的不一样,还是稍微有一点点的差异;

from sklearn.metrics import mean_squared_error,mean_absolute_error
pd.set_option( 'display.precision',12) # 浮点数精度
error_df = pd.DataFrame()
# 计算绝对值误差
error_df.loc['Skleanr','MAE'] = mean_absolute_error(y,y_sklean)
error_df.loc['numpy','MAE'] = mean_absolute_error(y,y_manual)
error_df.loc['OLS','MAE'] = mean_absolute_error(y,y_ols)
# 计算均方误差
error_df.loc['Skleanr','MSE'] = mean_squared_error(y,y_sklean)
error_df.loc['numpy','MSE'] = mean_squared_error(y,y_manual)
error_df.loc['OLS','MSE'] = mean_squared_error(y,y_ols)
error_df

基本是一样的...,讲道理应该不是完全一样,不过差异肯定在小数点后面很多位。

三、model.summary()

终于要进入正题:

其实本文主要研究的是这个:

可以看到能分成三部分,上面是一些基础信息,中间是系数部分,下面是一些检验量。

其中一眼就能看出来的不讲:

(1)R-squared 即R方,越接近1证明拟合得越好

sum of squares of the regression,即预测数据与原始数据均值之差的平方和。

 SSR=\sum\nolimits_{i=1}\nolimits^{n}(\hat{y_i}-\bar{y})^2

Total sum of squares,即原始数据和均值之差的平方和

SST=\sum\nolimits_{i=1}\nolimits^{n}(y_i-\bar{y})^2

The sum of squares due to error,拟合数据和原始数据对应点的误差的平方和

SSE=\sum\nolimits_{i=1}\nolimits^{n}(y_i-\hat{y_i})^2

SST=SSR+SSE

R-squared=\frac{SSR}{SST}=\frac{SST-SSE}{SST}=1-\frac{SSE}{SST}

# 验证R方
ssr = ((y_ols - y.mean())**2).sum() # 118359
sst = ((y - y.mean())**2).sum() # 118438
sse = ((y - y_ols)**2).sum() # 78.57
print('R方')
print(ssr/sst,result.rsquared) # 模型属性可查看

R方
0.9993365473073768 0.9993365473073765

(2)Adj. R-squared 调整R方,根据自由度和样本数修正了一下

注意这里面减去模型中的变量,我们此处DF model =1

# 验证调整R方
1-(100-1)*(1-result.rsquared)/(100-1-1),result.rsquared_adj

(0.9993297773819416, 0.9993297773819416)

(3)F-statistic,可参考

统计|如何理解多元线性回归的F检验的作用与目的-CSDN博客

f_stats = (ssr/1)/sse*(100-1-1) # 注意k=1
print(f_stats) 
print(result.fvalue)
# 完美匹配,最后一点误差可忽略
147614.11435208278
147614.11435208275

(4)Prob (F-statistic)

自然是查表得到的,要手算也行

from scipy.stats import f
# 懂点统计学都知道,曲线下面积即两个面积的和=1
p_value = f.cdf(f_stats,1,98) # 累计分布
print(p_value)

p_value_2= f.sf(f_stats,1,98) # 剩下的面积
print(p_value_2)             
print(result.f_pvalue)

0.9999999999999999
1.4929610240226694e-157
1.492961024022682e-157

(5)Log-Likelihood,似然函数的LOG

可参考https://www.le.ac.uk/users/dsgp1/COURSES/MATHSTAT/13mlreg.pdf

其中T为样本数,本文为100;sigma方即方差,是除以样本数100的,最右边公式可以抵消,sigma方即是总体的误差SSE除了一个样本数。

# T为样本数
L = -(100/2)*np.log(2*np.pi)-(100/2)*np.log(sse/100) - 100/2
print(L),print(result.llf)

-129.8399875471549
-129.8399875471549

(6) AIC,BIC

可参考:

模型中AIC和BIC以及loglikelihood的关系-腾讯云开发者社区-腾讯云

aic= -2*result.llf + 2*2
bic = -2*result.llf + 2*np.log(100)
print(aic),print(result.aic)
print(bic),print(result.bic)

263.6799750943098
263.6799750943098
268.890315466286
268.890315466286

至此第一部分的手动核算已经完毕;

写在最后:

一开始,我查到一篇文章,讲了一部分,但是很多东西并没有完全详尽,想着补充写明白:

statsmodels中的summary解读(使用OLS)_sm.glm.summary 解读-CSDN博客

(公式不够全面,不想深入研究也差不多够了)

有时候就很奇怪,等我边写边查资料的时候,发现了一位大佬的文章,很详细地写明了全部公式并手算验证,顿时就感觉没有必要继续写下去了。

请看他的:

【超详细】多元线性回归模型statsmodels_ols_回归残差 omnibus检验-CSDN博客

其他参考内容:

Ordinary Least Squares Regression in Python | DataRobot Blog(这篇文章有模型报告各个内容的英文简单释义,不过没有数学公式)

  • 9
    点赞
  • 27
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值