Python统计分析库statsmodels的OLS

来源：互联网发布：流燕抢票软件下载编辑：程序博客网时间：2024/06/06 04:13

statsmodels库官方文档http://www.statsmodels.org/stable/，里面包含很多统计模型和相应计算结果；一些Linear Regression Models例子http://www.statsmodels.org/stable/examples/index.html#regression

下面主要陈述常用的回归分析中的OLS:Ordinary Least Squares。
给定k组样本数据(yi,x(i)1,x(i)2,⋯,x(i)n)，i=1,2,⋯,k，
用n+1维一次多项式回归模型

y (x) = α 0 + α 1 x 1 + α 2 x 2 + \dots + α n x n = (α 0, α 1, \dots, α n) ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 1 x 1 x 2 ⋮ x n ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ≜ α T x

其中

α=(α0,α1,⋯,αn),x=(1,x1,x2,⋯,xn)T。

这里将常数项并入α,x以增广形式表出是为了和statsmodels.OLS模块编写源码对应，OLS里的多项式回归模型是没有常数项的，所以这里将常数项看作基为1的维度上的系数，OLS就是用样本数据拟合出最小二乘最小的系数组合，即求α。

对上述k组样本数据进行最小二乘拟合，即最小化

\sum i = 1 k (y i - α 0 + α 1 x (i) 1 + \dots + α n x (i) n) 2

statsmodels.OLS 的参数有endog, exog, missing, hasconst等，现在只考虑前两个。
这里写图片描述
第一个输入 endog 是回归模型中的因变量y(x), 输入是一个k维向量(y1,y2,⋯,yk)T。第二个输入 exog 是自变量，即k个样本点构成的k×(n+1)维数组

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 11 ⋮ 1 x (1) 1 x (2) 1 ⋮ x (k) 1 x (1) 2 x (2) 2 ⋮ x (k) 2 \dots \dots ⋱ \dots x (1) n x (2) n ⋮ x (k) n ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

通常，我们使用的数据集的

k个样本点构成的数组第一列并不全是

1，所以为了使用OLS模型函数，需要在数组左侧加上一列 1，就需要使用statmodels库的add_constant()函数，该函数的参数就是因变量数组(上述

k×(n+1)维数组去掉左侧一列

1)，也就是数据分析中用到的具有物理含义的list、pd.Series、pd.DataFrame；该函数的输出就如上述形式的

k×(n+1)维数组。

import pandas as pd # 读取数据到DataFrameimport 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+by=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+by = daily_data['cnt']regr = sm.OLS(y, x) # 普通最小二乘模型，ordinary least square modelres = 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的源码

阅读全文

0 0