回归分析--python应用篇(statsmodels)

实例一：用statsmodels库做一元回归分析

import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt

#构造变量
number=20 
x=np.linspace(0,10,number) #x值
X=sm.add_constant(x) #回归方程添加一列x0=1
bata=np.array([2,5]) #系数
e=np.random.normal(size=number) #误差
y=np.dot(X,bata)+e #观察值y

#建回归方程
#OLS（endog,exog=None,missing='none',hasconst=None) (endog:因变量，exog=自变量）
modle=sm.OLS(y,X) #最小二乘法
res=modle.fit() #拟合数据
Bata=res.params  #取系数
res.summary() #结果

#作图
Y=res.fittedvalues  #预测值
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(x, y, 'o', label='data')#原始数据
ax.plot(x, Y, 'r--.',label='test')#拟合数据
ax.legend(loc='best') #展示各点表示意思，即label
plt.show()

---

实例二：用statsmodels做多元回归分析

import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt

#构造变量
number=50
x=np.linspace(0,10,number)
X=np.column_stack((x,x**2))  #!!!
X=sm.add_constant(X)
bata=np.array([5,2,3])
e=np.random.normal(size=number)
y=np.dot(X,bata)+e

#建立方程
model=sm.OLS(y,X)
res=model.fit()
Bata=res.params
print(Bata)
res.summary()

#作图
y_=res.fittedvalues
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(x, y, 'o', label='data')
ax.plot(x, y_, 'r--.',label='OLS')
ax.legend(loc='best')
plt.show()

---

实例三：分类变量

在现实中，常出现分类问题，比如冰箱容量、耗能等级预测冰箱价格问题。

分类函数：sm.categorical(类别数组，drop=True)  【drop=True：去掉原类别数据，默认False，原数据会和分类变量合并】

例子：energy_level=np.array((2,1,3,2,3,1))  #由6台冰箱耗能分成3个等级

sm.categorical(energy_level)

返回：array([[2., 0., 1., 0.], [1., 1., 0., 0.], [3., 0., 0., 1.], [2., 0., 1., 0.], [3., 0., 0., 1.], [1., 1., 0., 0.]])

sm.categorical(energy_level,drop=True)

返回：array([[0., 1., 0.], [1., 0., 0.], [0., 0., 1.], [0., 1., 0.], [0., 0., 1.], [1., 0., 0.]])

如下图，categorical函数类似数据透视化，将等级化成0，1表示。

冰箱等级	1等级	2等级	3等级
2	0	1	0
1	1	0	0
3	0	0	1
2	0	1	0
3	0	0	1
1	1	0	0

import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt

#构造函数
number=50
group=np.zeros(number,int)
group[20:40]=1
group[40:]=2
category=sm.categorical(group,drop=True) #构造分类变量
x=np.linspace(0,10,number)
X=np.column_stack((x,category))
X=sm.add_constant(X)
bata=np.array([2,3,4,5,6])
e=np.random.normal(size=number)
y=np.dot(X,bata)+e

#建立方程
model=sm.OLS(y,X).fit()
model.summary()

#作图
Y=model.fittedvalues
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(x, y, 'o', label='data')
ax.plot(x, Y, 'r--.',label='OLS')
ax.legend(loc='best')
plt.show()