statsmodels库官方文档http://www.statsmodels.org/stable/,里面包含很多统计模型和相应计算结果;一些Linear Regression Models例子http://www.statsmodels.org/stable/examples/index.html#regression
下面主要陈述常用的回归分析中的OLS:Ordinary Least Squares。
给定
k
k
k组样本数据
(
y
i
,
x
1
(
i
)
,
x
2
(
i
)
,
⋯
,
x
n
(
i
)
)
,
i
=
1
,
2
,
⋯
,
k
(y_{i},x_{1}^{(i)},x_{2}^{(i)},\cdots,x_{n}^{(i)}),i=1,2,\cdots,k
(yi,x1(i),x2(i),⋯,xn(i)),i=1,2,⋯,k,
用
n
+
1
n+1
n+1维一次多项式回归模型
y
(
x
)
=
α
0
+
α
1
x
1
+
α
2
x
2
+
⋯
+
α
n
x
n
=
(
α
0
,
α
1
,
⋯
,
α
n
)
(
1
x
1
x
2
⋮
x
n
)
≜
α
T
x
y(x)=\alpha_{0}+\alpha_{1}x_{1}+\alpha_{2}x_{2}+\cdots+\alpha_{n}x_{n}\\ =(\alpha_{0},\alpha_{1},\cdots,\alpha_{n})\left( \begin{array}{ccc} 1\\ x_{1}\\ x_{2}\\ \vdots\\ x_{n} \end{array} \right )\triangleq\alpha^{T}x
y(x)=α0+α1x1+α2x2+⋯+αnxn=(α0,α1,⋯,αn)
1x1x2⋮xn
≜αTx
其中
α
=
(
α
0
,
α
1
,
⋯
,
α
n
)
,
x
=
(
1
,
x
1
,
x
2
,
⋯
,
x
n
)
T
\alpha=(\alpha_{0},\alpha_{1},\cdots,\alpha_{n}),x=(1,x_{1},x_{2},\cdots,x_{n})^{T}
α=(α0,α1,⋯,αn),x=(1,x1,x2,⋯,xn)T。
这里将常数项并入 α , x \alpha,x α,x以增广形式表出是为了和statsmodels.OLS模块编写源码对应,OLS里的多项式回归模型是没有常数项的,所以这里将常数项看作基为 1 1 1的维度上的系数,OLS就是用样本数据拟合出最小二乘最小的系数组合,即求 α \alpha α。
对上述
k
k
k组样本数据进行最小二乘拟合,即最小化
∑
i
=
1
k
(
y
i
−
α
0
+
α
1
x
1
(
i
)
+
⋯
+
α
n
x
n
(
i
)
)
2
\sum_{i=1}^{k}(y_{i}-\alpha_{0}+\alpha_{1}x_{1}^{(i)}+\cdots+\alpha_{n}x_{n}^{(i)})^2
i=1∑k(yi−α0+α1x1(i)+⋯+αnxn(i))2
statsmodels.OLS 的参数有endog, exog, missing, hasconst等 ,现在只考虑前两个。
第一个输入 endog 是回归模型中的因变量
y
(
x
)
y(x)
y(x), 输入是一个
k
k
k维向量
(
y
1
,
y
2
,
⋯
,
y
k
)
T
(y_{1},y_{2},\cdots,y_{k})^{T}
(y1,y2,⋯,yk)T。第二个输入 exog 是自变量,即
k
k
k个样本点构成的
k
×
(
n
+
1
)
k\times(n+1)
k×(n+1)维数组
(
1
x
1
(
1
)
x
2
(
1
)
⋯
x
n
(
1
)
1
x
1
(
2
)
x
2
(
2
)
⋯
x
n
(
2
)
⋮
⋮
⋮
⋱
⋮
1
x
1
(
k
)
x
2
(
k
)
⋯
x
n
(
k
)
)
\left( \begin{array}{ccc} 1&x_{1}^{(1)}&x_{2}^{(1)}&\cdots&x_{n}^{(1)}\\ 1&x_{1}^{(2)}&x_{2}^{(2)}&\cdots&x_{n}^{(2)}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&x_{1}^{(k)}&x_{2}^{(k)}&\cdots&x_{n}^{(k)} \end{array} \right )
11⋮1x1(1)x1(2)⋮x1(k)x2(1)x2(2)⋮x2(k)⋯⋯⋱⋯xn(1)xn(2)⋮xn(k)
通常,我们使用的数据集的
k
k
k个样本点构成的数组第一列并不全是
1
1
1,所以为了使用OLS模型函数,需要在数组左侧加上一列 1,就需要使用statmodels库的add_constant()函数,该函数的参数就是因变量数组(上述
k
×
(
n
+
1
)
k\times(n+1)
k×(n+1)维数组去掉左侧一列
1
1
1),也就是数据分析中用到的具有物理含义的list、pd.Series、pd.DataFrame;该函数的输出就如上述形式的
k
×
(
n
+
1
)
k\times(n+1)
k×(n+1)维数组。
import pandas as pd # 读取数据到DataFrame
import urllib # 获取网络数据
import shutil # 文件操作
import zipfile # 压缩解压
import os
# 建立临时目录
try:
os.system('mkdir bike_data')
except:
os.system('rm -rf bike_data; mkdir bike_data')
data_source = 'http://archive.ics.uci.edu/ml/machine-learning-databases/00275/Bike-Sharing-Dataset.zip' # 网络数据地址
zipname = 'bike_data/Bike-Sharing-Dataset.zip' # 拼接文件和路径
urllib.request.urlretrieve(data_source, zipname) # 获得数据
zip_ref = zipfile.ZipFile(zipname, 'r') # 创建一个ZipFile对象处理压缩文件
#zip_ref.extractall(temp_dir) # 解压
zip_ref.extractall('bike_data')
zip_ref.close()
daily_path = 'bike_data/day.csv'
daily_data = pd.read_csv(daily_path) # 读取csv文件
daily_data['dteday'] = pd.to_datetime(daily_data['dteday']) # 把字符串数据传换成日期数据
drop_list = ['instant', 'season', 'yr', 'mnth', 'holiday', 'workingday', 'weathersit', 'atemp', 'hum'] # 不关注的列
daily_data.drop(drop_list, inplace = True, axis = 1) # inplace=true在对象上直接操作
daily_data.head() # 看一看数据~
import statsmodels.api as sm
#最小二乘
from statsmodels.stats.outliers_influence import summary_table
#获得汇总信息
x=sm.add_constant(daily_data['temp'])
#线性回归增加常数项 y=kx+b
y=daily_data['cnt']
regr=sm.OLS(y,x)
res=regr.fit()
st, data, ss2 = summary_table(res, alpha=0.05)
#置信水平alpha=5%,st数据汇总,data数据详情,ss2数据列名
fitted_values = data[:,2]
#等价于res.fittedvalues
res.model.endog ==y.values
#拟合回归模型的endog值就是因变量y
res.fittedvalues #获取拟合y值
res.params #拟合回归模型参数
res.params[0]+res.params[1]*daily_data['temp']==res.fittedvalues #验证二维回归模型的拟合y值计算原理
总结下,常用的OLS模型模板
import statsmodels.api as sm # 最小二乘
from statsmodels.stats.outliers_influence import summary_table # 获得汇总信息
x = sm.add_constant(daily_data['temp']) # 线性回归增加常数项 y=kx+b
y = daily_data['cnt']
regr = sm.OLS(y, x) # 普通最小二乘模型,ordinary least square model
res = regr.fit() #res.model.endog
# 从模型获得拟合数据
st, data, ss2 = summary_table(res, alpha=0.05) # 置信水平alpha=5%,st数据汇总,data数据详情,ss2数据列名
fitted_values = data[:,2] #等价于res.fittedvalues
以上内容参考statsmodels官网和pip install的源码
欢迎关注公众号:瑞行AI,欢迎交流AI算法、数据分析、leetcode刷题等技术,提供技术方案咨询和就业指导服务!