可以
文成网站制作公司哪家好，找成都创新互联公司！从网页设计、网站建设、微信开发、APP开发、成都响应式网站建设公司等网站项目制作，到程序开发，运营维护。成都创新互联公司于2013年创立到现在10年的时间，我们拥有了丰富的建站经验和运维经验，来保证我们的工作的顺利进行。专注于网站建设就选成都创新互联公司。
最近项目中涉及基于Gradient Boosting Regression 算法拟合时间序列曲线的内容，利用python机器学习包 scikit-learn 中的GradientBoostingRegressor完成
因此就学习了下Gradient Boosting算法，在这里分享下我的理解
Boosting 算法简介
Boosting算法，我理解的就是两个思想：
1）“三个臭皮匠顶个诸葛亮”，一堆弱分类器的组合就可以成为一个强分类器；
2）“知错能改，善莫大焉”，不断地在错误中学习，迭代来降低犯错概率
当然，要理解好Boosting的思想，首先还是从弱学习算法和强学习算法来引入：
1）强学习算法：存在一个多项式时间的学习算法以识别一组概念，且识别的正确率很高；
2）弱学习算法：识别一组概念的正确率仅比随机猜测略好；
Kearns Valiant证明了弱学习算法与强学习算法的等价问题，如果两者等价，只需找到一个比随机猜测略好的学习算法，就可以将其提升为强学习算法。
那么是怎么实现“知错就改”的呢？
Boosting算法，通过一系列的迭代来优化分类结果，每迭代一次引入一个弱分类器，来克服现在已经存在的弱分类器组合的shortcomings
在Adaboost算法中，这个shortcomings的表征就是权值高的样本点
而在Gradient Boosting算法中,这个shortcomings的表征就是梯度
无论是Adaboost还是Gradient Boosting，都是通过这个shortcomings来告诉学习器怎么去提升模型，也就是“Boosting”这个名字的由来吧
Adaboost算法
Adaboost是由Freund 和 Schapire在1997年提出的，在整个训练集上维护一个分布权值向量W,用赋予权重的训练集通过弱分类算法产生分类假设（基学习器）y(x),然后计算错误率,用得到的错误率去更新分布权值向量w,对错误分类的样本分配更大的权值,正确分类的样本赋予更小的权值。每次更新后用相同的弱分类算法产生新的分类假设,这些分类假设的序列构成多分类器。对这些多分类器用加权的方法进行联合,最后得到决策结果。
其结构如下图所示：
前一个学习器改变权重w，然后再经过下一个学习器，最终所有的学习器共同组成最后的学习器。
如果一个样本在前一个学习器中被误分，那么它所对应的权重会被加重，相应地，被正确分类的样本的权重会降低。
这里主要涉及到两个权重的计算问题：
1）样本的权值
1 没有先验知识的情况下,初始的分布应为等概分布,样本数目为n,权值为1/n
2 每一次的迭代更新权值，提高分错样本的权重
2）弱学习器的权值
1 最后的强学习器是通过多个基学习器通过权值组合得到的。
2 通过权值体现不同基学习器的影响,正确率高的基学习器权重高。实际上是分类误差的一个函数
Gradient Boosting
和Adaboost不同，Gradient Boosting 在迭代的时候选择梯度下降的方向来保证最后的结果最好。
损失函数用来描述模型的“靠谱”程度，假设模型没有过拟合，损失函数越大，模型的错误率越高
如果我们的模型能够让损失函数持续的下降，则说明我们的模型在不停的改进，而最好的方式就是让损失函数在其梯度方向上下降。
下面这个流程图是Gradient Boosting的经典图了，数学推导并不复杂，只要理解了Boosting的思想，不难看懂
这里是直接对模型的函数进行更新，利用了参数可加性推广到函数空间。
训练F0-Fm一共m个基学习器，沿着梯度下降的方向不断更新ρm和am
GradientBoostingRegressor实现
python中的scikit-learn包提供了很方便的GradientBoostingRegressor和GBDT的函数接口，可以很方便的调用函数就可以完成模型的训练和预测
GradientBoostingRegressor函数的参数如下：
class sklearn.ensemble.GradientBoostingRegressor(loss='ls', learning_rate=0.1, n_estimators=100, subsample=1.0, min_samples_split=2, min_samples_leaf=1, min_weight_fraction_leaf=0.0, max_depth=3, init=None, random_state=None, max_features=None, alpha=0.9, verbose=0, max_leaf_nodes=None, warm_start=False, presort='auto')[source]¶
loss: 选择损失函数，默认值为ls(least squres)
learning_rate: 学习率，模型是0.1
n_estimators: 弱学习器的数目，默认值100
max_depth: 每一个学习器的最大深度，限制回归树的节点数目，默认为3
min_samples_split: 可以划分为内部节点的最小样本数，默认为2
min_samples_leaf: 叶节点所需的最小样本数，默认为1
……
可以参考
官方文档里带了一个很好的例子，以500个弱学习器，最小平方误差的梯度提升模型，做波士顿房价预测，代码和结果如下：
1 import numpy as np 2 import matplotlib.pyplot as plt 3 4 from sklearn import ensemble 5 from sklearn import datasets 6 from sklearn.utils import shuffle 7 from sklearn.metrics import mean_squared_error 8 9 ###############################################################################10 # Load data11 boston = datasets.load_boston()12 X, y = shuffle(boston.data, boston.target, random_state=13)13 X = X.astype(np.float32)14 offset = int(X.shape[0] * 0.9)15 X_train, y_train = X[:offset], y[:offset]16 X_test, y_test = X[offset:], y[offset:]17 18 ###############################################################################19 # Fit regression model20 params = {'n_estimators': 500, 'max_depth': 4, 'min_samples_split': 1,21 'learning_rate': 0.01, 'loss': 'ls'}22 clf = ensemble.GradientBoostingRegressor(**params)23 24 clf.fit(X_train, y_train)25 mse = mean_squared_error(y_test, clf.predict(X_test))26 print("MSE: %.4f" % mse)27 28 ###############################################################################29 # Plot training deviance30 31 # compute test set deviance32 test_score = np.zeros((params['n_estimators'],), dtype=np.float64)33 34 for i, y_pred in enumerate(clf.staged_predict(X_test)):35 test_score[i] = clf.loss_(y_test, y_pred)36 37 plt.figure(figsize=(12, 6))38 plt.subplot(1, 2, 1)39 plt.title('Deviance')40 plt.plot(np.arange(params['n_estimators']) + 1, clf.train_score_, 'b-',41 label='Training Set Deviance')42 plt.plot(np.arange(params['n_estimators']) + 1, test_score, 'r-',43 label='Test Set Deviance')44 plt.legend(loc='upper right')45 plt.xlabel('Boosting Iterations')46 plt.ylabel('Deviance')47 48 ###############################################################################49 # Plot feature importance50 feature_importance = clf.feature_importances_51 # make importances relative to max importance52 feature_importance = 100.0 * (feature_importance / feature_importance.max())53 sorted_idx = np.argsort(feature_importance)54 pos = np.arange(sorted_idx.shape[0]) + .555 plt.subplot(1, 2, 2)56 plt.barh(pos, feature_importance[sorted_idx], align='center')57 plt.yticks(pos, boston.feature_names[sorted_idx])58 plt.xlabel('Relative Importance')59 plt.title('Variable Importance')60 plt.show()
可以发现，如果要用Gradient Boosting 算法的话，在sklearn包里调用还是非常方便的，几行代码即可完成，大部分的工作应该是在特征提取上。
感觉目前做数据挖掘的工作，特征设计是最重要的，据说现在kaggle竞赛基本是GBDT的天下，优劣其实还是特征上，感觉做项目也是，不断的在研究数据中培养对数据的敏感度。
1、什么是多元线性回归模型？
当y值的影响因素不唯一时,采用多元线性回归模型。
y =y=β0+β1x1+β2x2+...+βnxn
例如商品的销售额可能不电视广告投入,收音机广告投入,报纸广告投入有关系,可以有 sales =β0+β1*TV+β2* radio+β3*newspaper.
2、使用pandas来读取数据
pandas 是一个用于数据探索、数据分析和数据处理的python库
[python] view plain copy
import pandas as pd
[html] view plain copy
pre name="code" class="python"# read csv file directly from a URL and save the results
data = pd.read_csv('/home/lulei/Advertising.csv')
# display the first 5 rows
data.head()
上面代码的运行结果：
TV Radio Newspaper Sales
0 230.1 37.8 69.2 22.1
1 44.5 39.3 45.1 10.4
2 17.2 45.9 69.3 9.3
3 151.5 41.3 58.5 18.5
4 180.8 10.8 58.4 12.9
上面显示的结果类似一个电子表格，这个结构称为Pandas的数据帧(data frame)，类型全称：pandas.core.frame.DataFrame.
pandas的两个主要数据结构：Series和DataFrame：
Series类似于一维数组，它有一组数据以及一组与之相关的数据标签(即索引)组成。
DataFrame是一个表格型的数据结构，它含有一组有序的列，每列可以是不同的值类型。DataFrame既有行索引也有列索引，它可以被看做由Series组成的字典。
[python] view plain copy
# display the last 5 rows
data.tail()
只显示结果的末尾5行
TV Radio Newspaper Sales
195 38.2 3.7 13.8 7.6
196 94.2 4.9 8.1 9.7
197 177.0 9.3 6.4 12.8
198 283.6 42.0 66.2 25.5
199 232.1 8.6 8.7 13.4
[html] view plain copy
# check the shape of the DataFrame(rows, colums)
data.shape
查看DataFrame的形状,注意第一列的叫索引，和数据库某个表中的第一列类似。
(200,4)
3、分析数据
特征：
TV：对于一个给定市场中单一产品，用于电视上的广告费用（以千为单位）
Radio：在广播媒体上投资的广告费用
Newspaper：用于报纸媒体的广告费用
响应：
Sales：对应产品的销量
在这个案例中，我们通过不同的广告投入，预测产品销量。因为响应变量是一个连续的值，所以这个问题是一个回归问题。数据集一共有200个观测值，每一组观测对应一个市场的情况。
注意：这里推荐使用的是seaborn包。网上说这个包的数据可视化效果比较好看。其实seaborn也应该属于matplotlib的内部包。只是需要再次的单独安装。
[python] view plain copy
import seaborn as sns
import matplotlib.pyplot as plt
# visualize the relationship between the features and the response using scatterplots
sns.pairplot(data, x_vars=['TV','Radio','Newspaper'], y_vars='Sales', size=7, aspect=0.8)
plt.show()#注意必须加上这一句，否则无法显示。
[html] view plain copy
这里选择TV、Radio、Newspaper 作为特征，Sales作为观测值
[html] view plain copy
返回的结果：
seaborn的pairplot函数绘制X的每一维度和对应Y的散点图。通过设置size和aspect参数来调节显示的大小和比例。可以从图中看出，TV特征和销量是有比较强的线性关系的，而Radio和Sales线性关系弱一些，Newspaper和Sales线性关系更弱。通过加入一个参数kind='reg'，seaborn可以添加一条最佳拟合直线和95%的置信带。
[python] view plain copy
sns.pairplot(data, x_vars=['TV','Radio','Newspaper'], y_vars='Sales', size=7, aspect=0.8, kind='reg')
plt.show()
结果显示如下：
4、线性回归模型
优点：快速；没有调节参数；可轻易解释；可理解。
缺点：相比其他复杂一些的模型，其预测准确率不是太高，因为它假设特征和响应之间存在确定的线性关系，这种假设对于非线性的关系，线性回归模型显然不能很好的对这种数据建模。
线性模型表达式： y=β0+β1x1+β2x2+...+βnxn 其中
y是响应
β0是截距
β1是x1的系数，以此类推
在这个案例中： y=β0+β1∗TV+β2∗Radio+...+βn∗Newspaper
(1)、使用pandas来构建X(特征向量)和y(标签列)
scikit-learn要求X是一个特征矩阵，y是一个NumPy向量。
pandas构建在NumPy之上。
因此，X可以是pandas的DataFrame，y可以是pandas的Series，scikit-learn可以理解这种结构。
[python] view plain copy
#create a python list of feature names
feature_cols = ['TV', 'Radio', 'Newspaper']
# use the list to select a subset of the original DataFrame
X = data[feature_cols]
# equivalent command to do this in one line
X = data[['TV', 'Radio', 'Newspaper']]
# print the first 5 rows
print X.head()
# check the type and shape of X
print type(X)
print X.shape
输出结果如下：
TV Radio Newspaper
0 230.1 37.8 69.2
1 44.5 39.3 45.1
2 17.2 45.9 69.3
3 151.5 41.3 58.5
4 180.8 10.8 58.4
class 'pandas.core.frame.DataFrame'
(200, 3)
[python] view plain copy
# select a Series from the DataFrame
y = data['Sales']
# equivalent command that works if there are no spaces in the column name
y = data.Sales
# print the first 5 values
print y.head()
输出的结果如下：
0 22.1
1 10.4
2 9.3
3 18.5
4 12.9
Name: Sales
（2）、构建训练集与测试集
[html] view plain copy
pre name="code" class="python"span style="font-size:14px;"##构造训练集和测试集
from sklearn.cross_validation import train_test_split #这里是引用了交叉验证
X_train,X_test, y_train, y_test = train_test_split(X, y, random_state=1)
#default split is 75% for training and 25% for testing
[html] view plain copy
print X_train.shape
print y_train.shape
print X_test.shape
print y_test.shape
输出结果如下：
(150, 3)
(150,)
(50, 3)
(50,)
注：上面的结果是由train_test_spilit()得到的，但是我不知道为什么我的版本的sklearn包中居然报错：
ImportError Traceback (most recent call last)ipython-input-182-3eee51fcba5a in module() 1 ###构造训练集和测试集---- 2 from sklearn.cross_validation import train_test_split 3 #import sklearn.cross_validation 4 X_train,X_test, y_train, y_test = train_test_split(X, y, random_state=1) 5 # default split is 75% for training and 25% for testingImportError: cannot import name train_test_split
处理方法：1、我后来重新安装sklearn包。再一次调用时就没有错误了。
2、自己写函数来认为的随机构造训练集和测试集。(这个代码我会在最后附上。)
（3）sklearn的线性回归
[html] view plain copy
from sklearn.linear_model import LinearRegression
linreg = LinearRegression()
model=linreg.fit(X_train, y_train)
print model
print linreg.intercept_
print linreg.coef_
输出的结果如下：
LinearRegression(copy_X=True, fit_intercept=True, normalize=False)
2.66816623043
[ 0.04641001 0.19272538 -0.00349015]
[html] view plain copy
# pair the feature names with the coefficients
zip(feature_cols, linreg.coef_)
输出如下：
[('TV', 0.046410010869663267),
('Radio', 0.19272538367491721),
('Newspaper', -0.0034901506098328305)]
y=2.668+0.0464∗TV+0.192∗Radio-0.00349∗Newspaper
如何解释各个特征对应的系数的意义？
对于给定了Radio和Newspaper的广告投入，如果在TV广告上每多投入1个单位，对应销量将增加0.0466个单位。就是加入其它两个媒体投入固定，在TV广告上每增加1000美元（因为单位是1000美元），销量将增加46.6（因为单位是1000）。但是大家注意这里的newspaper的系数居然是负数，所以我们可以考虑不使用newspaper这个特征。这是后话，后面会提到的。
（4）、预测
[python] view plain copy
y_pred = linreg.predict(X_test)
print y_pred
[python] view plain copy
print type(y_pred)
输出结果如下：
[ 14.58678373 7.92397999 16.9497993 19.35791038 7.36360284
7.35359269 16.08342325 9.16533046 20.35507374 12.63160058
22.83356472 9.66291461 4.18055603 13.70368584 11.4533557
4.16940565 10.31271413 23.06786868 17.80464565 14.53070132
15.19656684 14.22969609 7.54691167 13.47210324 15.00625898
19.28532444 20.7319878 19.70408833 18.21640853 8.50112687
9.8493781 9.51425763 9.73270043 18.13782015 15.41731544
5.07416787 12.20575251 14.05507493 10.6699926 7.16006245
11.80728836 24.79748121 10.40809168 24.05228404 18.44737314
20.80572631 9.45424805 17.00481708 5.78634105 5.10594849]
type 'numpy.ndarray'
5、回归问题的评价测度
(1) 评价测度
对于分类问题，评价测度是准确率，但这种方法不适用于回归问题。我们使用针对连续数值的评价测度(evaluation metrics)。
这里介绍3种常用的针对线性回归的测度。
1)平均绝对误差(Mean Absolute Error, MAE)
(2)均方误差(Mean Squared Error, MSE)
(3)均方根误差(Root Mean Squared Error, RMSE)
这里我使用RMES。
[python] view plain copy
pre name="code" class="python"#计算Sales预测的RMSE
print type(y_pred),type(y_test)
print len(y_pred),len(y_test)
print y_pred.shape,y_test.shape
from sklearn import metrics
import numpy as np
sum_mean=0
for i in range(len(y_pred)):
sum_mean+=(y_pred[i]-y_test.values[i])**2
sum_erro=np.sqrt(sum_mean/50)
# calculate RMSE by hand
print "RMSE by hand:",sum_erro
最后的结果如下：
type 'numpy.ndarray' class 'pandas.core.series.Series'
50 50
(50,) (50,)
RMSE by hand: 1.42998147691
（2）做ROC曲线
[python] view plain copy
import matplotlib.pyplot as plt
plt.figure()
plt.plot(range(len(y_pred)),y_pred,'b',label="predict")
plt.plot(range(len(y_pred)),y_test,'r',label="test")
plt.legend(loc="upper right") #显示图中的标签
plt.xlabel("the number of sales")
plt.ylabel('value of sales')
plt.show()
显示结果如下：（红色的线是真实的值曲线，蓝色的是预测值曲线）
直到这里整个的一次多元线性回归的预测就结束了。
6、改进特征的选择
在之前展示的数据中，我们看到Newspaper和销量之间的线性关系竟是负关系（不用惊讶，这是随机特征抽样的结果。换一批抽样的数据就可能为正了），现在我们移除这个特征，看看线性回归预测的结果的RMSE如何？
依然使用我上面的代码，但只需修改下面代码中的一句即可：
[python] view plain copy
#create a python list of feature names
feature_cols = ['TV', 'Radio', 'Newspaper']
# use the list to select a subset of the original DataFrame
X = data[feature_cols]
# equivalent command to do this in one line
#X = data[['TV', 'Radio', 'Newspaper']]#只需修改这里即可pre name="code" class="python" style="font-size: 15px; line-height: 35px;"X = data[['TV', 'Radio']] #去掉newspaper其他的代码不变
# print the first 5 rowsprint X.head()# check the type and shape of Xprint type(X)print X.shape
最后的到的系数与测度如下：
LinearRegression(copy_X=True, fit_intercept=True, normalize=False)
2.81843904823
[ 0.04588771 0.18721008]
RMSE by hand: 1.28208957507
然后再次使用ROC曲线来观测曲线的整体情况。我们在将Newspaper这个特征移除之后，得到RMSE变小了，说明Newspaper特征可能不适合作为预测销量的特征，于是，我们得到了新的模型。我们还可以通过不同的特征组合得到新的模型，看看最终的误差是如何的。
备注：
之前我提到了这种错误：
注：上面的结果是由train_test_spilit()得到的，但是我不知道为什么我的版本的sklearn包中居然报错：
ImportError Traceback (most recent call last)ipython-input-182-3eee51fcba5a in module() 1 ###构造训练集和测试集---- 2 from sklearn.cross_validation import train_test_split 3 #import sklearn.cross_validation 4 X_train,X_test, y_train, y_test = train_test_split(X, y, random_state=1) 5 # default split is 75% for training and 25% for testingImportError: cannot import name train_test_split
处理方法：1、我后来重新安装sklearn包。再一次调用时就没有错误了。
2、自己写函数来认为的随机构造训练集和测试集。(这个代码我会在最后附上。)
这里我给出我自己写的函数：
利用python进行线性回归
理解什么是线性回归
线性回归也被称为最小二乘法回归（Linear Regression, also called Ordinary Least-Squares (OLS) Regression）。它的数学模型是这样的：
y = a+ b* x＋e
其中，a 被称为常数项或截距；b 被称为模型的回归系数或斜率；e 为误差项。a 和 b 是模型的参数。
当然，模型的参数只能从样本数据中估计出来：
y'= a' + b'* x
我们的目标是选择合适的参数，让这一线性模型最好地拟合观测值。拟合程度越高，模型越好。
那么，接下来的问题就是，我们如何判断拟合的质量呢？
这一线性模型可以用二维平面上的一条直线来表示，被称为回归线。
模型的拟合程度越高，也即意味着样本点围绕回归线越紧密。
如何计算样本点与回归线之间的紧密程度呢？
高斯和勒让德找到的方法是：被选择的参数，应该使算出来的回归线与观测值之差的平房和最小。用函数表示为：
这被称为最小二乘法。最小二乘法的原理是这样的：当预测值和实际值距离的平方和最小时，就选定模型中的两个参数（a 和 b）。这一模型并不一定反映解释变量和反应变量真实的关系。但它的计算成本低；相比复杂模型更容易解释。
模型估计出来后，我们要回答的问题是：
我们的模型拟合程度如何？或者说，这个模型对因变量的解释力如何？（R2）
整个模型是否能显著预测因变量的变化？（F 检验）
每个自变量是否能显著预测因变量的变化？（t 检验）
首先回答第一个问题。为了评估模型的拟合程度如何，我们必须有一个可以比较的基线模型。
如果让你预测一个人的体重是多少？在没有任何额外信息的情况下，你可能会用平均值来预测，尽管会存在一定误差，但总比瞎猜好。
现在，如果你知道他的身高信息，你的预测值肯定与平均值不一样。额外信息相比平均值更能准确地预测被预测的变量的能力，就代表模型的解释力大小。
上图中，SSA 代表由自变量 x 引起的 y 的离差平方和，即回归平方和，代表回归模型的解释力；SSE 代表由随机因素引起的 y 的离差平方和，即剩余平方和，代表回归模型未能解释的部分；SST 为总的离差平方和，即我们仅凭 y 的平均值去估计 y 时所产生的误差。
用模型能够解释的变异除以总的变异就是模型的拟合程度：
R2=SSA/SST=1-SSE
R2（R 的平方）也被称为决定系数或判定系数。
第二个问题，我们的模型是否显著预测了 y 的变化？
假设 y 与 x 的线性关系不明显，那么 SSA 相对 SSE 占有较大的比例的概率则越小。换句话说，在 y 与 x 无线性关系的前提下，SSA 相对 SSE 的占比越高的概率是越小的，这会呈现一定的概率分布。统计学家告诉我们它满足 F 分布，就像这样：
如果 SSA 相对 SSE 占比较大的情况出现了，比如根据 F 分布，这个值出现的概率小于 5%。那么，我们最好是拒绝 y 与 x 线性关系不显著的原始假设，认为二者存在显著的线性关系较为合适。
第三个问题，每个自变量是否能显著预测因变量的变化？换句话说，回归系数是否显著？
回归系数的显著性检验是围绕回归系数的抽样分布（t 分布）来进行的，推断过程类似于整个模型的检验过程，不赘言。
实际上，对于只有一个自变量的一元线性模型，模型的显著性检验和回归系数的检验是一致的，但对于多元线性模型来说，二者就不能等价了。
利用 statsmodels 进行最小二乘回归
＃导入相应模块
In [1]: import numpy as np
In [2]: import pandas as pd
In [3]: import statsmodels.api as sm
＃将数据导入 pandas 的 dataframe 对象，第一列（年份）作为行标签
In [4]: df=pd.read_csv('/Users/xiangzhendong/Downloads/vincentarelbundock-Rdatasets-1218370/csv/datasets/longley.csv', index_col=0)
＃查看头部数据
In [5]: df.head()
Out[5]:
GNP.deflator GNP Unemployed Armed.Forces Population Year \
1947 83.0 234.289 235.6 159.0 107.608 1947
1948 88.5 259.426 232.5 145.6 108.632 1948
1949 88.2 258.054 368.2 161.6 109.773 1949
1950 89.5 284.599 335.1 165.0 110.929 1950
1951 96.2 328.975 209.9 309.9 112.075 1951
Employed
1947 60.323
1948 61.122
1949 60.171
1950 61.187
1951 63.221
＃设置预测变量和结果变量，用 GNP 预测 Employed
In [6]: y=df.Employed ＃结果变量
In [7]: X=df.GNP ＃预测变量
＃为模型增加常数项，即回归线在 y 轴上的截距
In [8]: X=sm.add_constant(X)
＃执行最小二乘回归，X 可以是 numpy array 或 pandas dataframe（行数等于数据点个数，列数为预测变量个数），y 可以是一维数组（numpy array）或 pandas series
In [10]: est=sm.OLS(y,X)
使用 OLS 对象的 fit() 方法来进行模型拟合
In [11]: est=est.fit()
＃查看模型拟合的结果
In [12]: est.summary()
Out[12]:
＃查看最终模型的参数
In [13]: est.params
Out[13]:
const 51.843590
GNP 0.034752
dtype: float64
＃选择 100 个从最小值到最大值平均分布（equally spaced）的数据点
In [14]: X_prime=np.linspace(X.GNP.min(), X.GNP.max(),100)[:,np.newaxis]
In [15]: X_prime=sm.add_constant(X_prime)
＃计算预测值
In [16]: y_hat=est.predict(X_prime)
In [17]: plt.scatter(X.GNP, y, alpha=0.3) ＃画出原始数据
＃分别给 x 轴和 y 轴命名
In [18]: plt.xlabel("Gross National Product")
In [19]: plt.ylabel("Total Employment")
In [20]: plt.plot(X_prime[:,1], y_hat, 'r', alpha=0.9) ＃添加回归线，红色
多元线性回归（预测变量不止一个）
我们用一条直线来描述一元线性模型中预测变量和结果变量的关系，而在多元回归中，我们将用一个多维（p）空间来拟合多个预测变量。下面表现了两个预测变量的三维图形：商品的销量以及在电视和广播两种不同媒介的广告预算。
数学模型是：
Sales = beta_0 + beta_1＊TV + beta_2＊Radio
图中，白色的数据点是平面上的点，黑色的数据点事平面下的点。平面的颜色是由对应的商品销量的高低决定的，高是红色，低是蓝色。
利用 statsmodels 进行多元线性回归
In [1]: import pandas as pd
In [2]: import numpy as np
In [3]: import statsmodels.api as sm
In [4]: df_adv=pd.read_csv('g.csv',index_col=0)
In [6]: X=df_adv[['TV','Radio']]
In [7]: y=df_adv['Sales']
In [8]: df_adv.head()
Out[8]:
TV Radio Newspaper Sales
1 230.1 37.8 69.2 22.1
2 44.5 39.3 45.1 10.4
3 17.2 45.9 69.3 9.3
4 151.5 41.3 58.5 18.5
5 180.8 10.8 58.4 12.9
In [9]: X=sm.add_constant(X)
In [10]: est=sm.OLS(y,X).fit()
In [11]: est.summary()
Out[11]:
你也可以使用 statsmodels 的 formula 模块来建立多元回归模型
In [12]: import statsmodels.formula.api as smf
In [13]: est=smf.ols(formula='Sales ~ TV + Radio',data=df_adv).fit()
处理分类变量
性别或地域都属于分类变量。
In [15]: df= pd.read_csv('httd.edu/~tibs/ElemStatLearn/datasets/SAheart.data', index_col=0)
In [16]: X=df.copy()
利用 dataframe 的 pop 方法将 chd 列单独提取出来
In [17]: y=X.pop('chd')
In [18]: df.head()
Out[18]:
sbp tobacco ldl adiposity famhist typea obesity alcohol \
row.names
1 160 12.00 5.73 23.11 Present 49 25.30 97.20
2 144 0.01 4.41 28.61 Absent 55 28.87 2.06
3 118 0.08 3.48 32.28 Present 52 29.14 3.81
4 170 7.50 6.41 38.03 Present 51 31.99 24.26
5 134 13.60 3.50 27.78 Present 60 25.99 57.34
age chd
row.names
1 52 1
2 63 1
3 46 0
4 58 1
5 49 1
In [19]: y.groupby(X.famhist).mean()
Out[19]:
famhist
Absent 0.237037
Present 0.500000
Name: chd, dtype: float64
In [20]: import statsmodels.formula.api as smf
In [21]: df['famhist_ord']=pd.Categorical(df.famhist).labels
In [22]: est=smf.ols(formula="chd ~ famhist_ord", data=df).fit()
分类变量的编码方式有许多，其中一种编码方式是虚拟变量编码（dummy－encoding），就是把一个 k 个水平的分类变量编码成 k－1 个二分变量。在 statsmodels 中使用 C 函数实现。
In [24]: est=smf.ols(formula="chd ~ C(famhist)", data=df).fit()
In [26]: est.summary()
Out[26]:
处理交互作用
随着教育年限（education）的增长，薪酬 (wage) 会增加吗？这种影响对男性和女性而言是一样的吗？
这里的问题就涉及性别与教育年限的交互作用。
换言之，教育年限对薪酬的影响是男女有别的。
＃导入相关模块
In [1]: import pandas as pd
In [2]: import numpy as np
In [4]: import statsmodels.api as sm
＃导入数据，存入 dataframe 对象
In [5]: df=pd.read_csv('/Users/xiangzhendong/Downloads/pydatafromweb/wages.csv')
In [6]: df[['Wage','Education','Sex']].tail()
Out[6]:
Wage Education Sex
529 11.36 18 0
530 6.10 12 1
531 23.25 17 1
532 19.88 12 0
533 15.38 16 0
由于性别是一个二分变量，我们可以绘制两条回归线，一条是 sex＝0（男性），一条是 sex＝1（女性）
＃绘制散点图
In [7]: plt.scatter(df.Education,df.Wage, alpha=0.3)
In [9]: plt.xlabel('education')
In [10]: plt.ylabel('wage')
＃linspace 的作用是生成从最小到最大的均匀分布的 n 个数
In [17]: education_linspace=np.linspace(df.Education.min(), df.Education.max(),100)
In [12]: import statsmodels.formula.api as smf
In [13]: est=smf.ols(formula='Wage ~ Education + Sex', data=df).fit()
In [18]: plt.plot(education_linspace, est.params[0]+est.params[1]education_linspace+est.params[2]0, 'r')
In [19]: plt.plot(education_linspace, est.params[0]+est.params[1]education_linspace+est.params[2]1, 'g')
以上两条线是平行的。这是因为分类变量只影响回归线的截距，不影响斜率。
接下来我们可以为回归模型增加交互项来探索交互效应。也就是说，对于两个类别，回归线的斜率是不一样的。
In [32]: plt.scatter(df.Education,df.Wage, alpha=0.3)
In [33]: plt.xlabel('education')
In [34]: plt.ylabel('wage')
＃使用＊代表我们的回归模型中除了交互效应，也包括两个变量的主效应；如果只想看交互效应，可以用：代替，但通常不会只看交互效应
In [35]: est=smf.ols(formula='Wage ~ Sex*Education', data=df).fit()
In [36]: plt.plot(education_linspace, est.params[0]+est.params[1]0+est.params[2]education_linspace+est.params[3]0education_linspace, 'r')
In [37]: plt.plot(education_linspace, est.params[0]+est.params[1]1+est.params[2]education_linspace+est.params[3]1education_linspace, 'g')
参考资料：
DataRobot | Ordinary Least Squares in Python
DataRoboe | Multiple Regression using Statsmodels
AnalyticsVidhya | 7 Types of Regression Techniques you should know!
