statsmodels库官方文档http://www.statsmodels.org/stable/,里面包含很多统计模型和相应计算结果;一些Linear Regression Models例子http://www.statsmodels.org/stable/examples/index.html#regression
下面主要陈述常用的回归分析中的OLS:Ordinary Least Squares。
给定kk组样本数据(yi,x(i)1,x(i)2,⋯,x(i)n),i=1,2,⋯,k(yi,x1(i),x2(i),⋯,xn(i)),i=1,2,⋯,k,
用n+1n+1维一次多项式回归模型
y(x)=α0+α1x1+α2x2+⋯+αnxn=(α0,α1,⋯,αn)⎛⎝⎜⎜⎜⎜⎜⎜1x1x2⋮xn⎞⎠⎟⎟⎟⎟⎟⎟≜αTx
y(x)=α0+α1x1+α2x2+⋯+αnxn=(α0,α1,⋯,αn)(1x1x2⋮xn)≜αTx
其中α=(α0,α1,⋯,αn),x=(1,x1,x2,⋯,xn)Tα=(α0,α1,⋯,αn),x=(1,x1,x2,⋯,xn)T。
这里将常数项并入α,xα,x以增广形式表出是为了和statsmodels.OLS模块编写源码对应,OLS里的多项式回归模型是没有常数项的,所以这里将常数项看作基为11的维度上的系数,OLS就是用样本数据拟合出最小二乘最小的系数组合,即求αα。
对上述kk组样本数据进行最小二乘拟合,即最小化
∑i=1k(yi−α0+α1x(i)1+⋯+αnx(i)n)2
∑i=1k(yi−α0+α1x1(i)+⋯+αnxn(i))2
statsmodels.OLS 的参数有endog, exog, missing, hasconst等 ,现在只考虑前两个。
第一个输入 endog 是回归模型中的因变量y(x)y(x), 输入是一个kk维向量(y1,y2,⋯,yk)T(y1,y2,⋯,yk)T。第二个输入 exog 是自变量,即kk个样本点构成的k×(n+1)k×(n+1)维数组
⎛⎝⎜⎜⎜⎜⎜11⋮1x(1)1x(2)1⋮x(k)1x(1)2x(2)2⋮x(k)2⋯⋯⋱⋯x(1)nx(2)n⋮x(k)n⎞⎠⎟⎟⎟⎟⎟
(1x1(1)x2(1)⋯xn(1)1x1(2)x2(2)⋯xn(2)⋮⋮⋮⋱⋮1x1(k)x2(k)⋯xn(k))
通常,我们使用的数据集的kk个样本点构成的数组第一列并不全是11,所以为了使用OLS模型函数,需要在数组左侧加上一列 1,就需要使用statmodels库的add_constant()函数,该函数的参数就是因变量数组(上述k×(n+1)k×(n+1)维数组去掉左侧一列11),也就是数据分析中用到的具有物理含义的list、pd.Series、pd.DataFrame;该函数的输出就如上述形式的k×(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