线性回归介绍和Python实现

一、线性回归介绍

线性回归是利用数理统计中回归分析来确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法。其表达形式为$y=w^{\prime} x+e$,e为误差(服从均值为0的正态分布),数学上e叫截距。

  • 当回归分析中只包括一个自变量和一个因变量,且二者的关系可用一条直线近似表示,称为一元线性回归分析
  • 当回归分析中包括两个或两个以上的自变量,且因变量和自变量之间是线性关系,则称为多元线性回归分析

1、回归分析流程

线性回归属于回归问题。回归问题的流程如下:

  1. 给定数据集中每个样本及其正确答案。
  2. 选择一个模型函数h,这里的h代表hypothesis(假设)。
  3. 为h找到此数据集的(未必是全局)最优解,即找出最优解下的h的参数。

这里给定的数据集是训练集(Training Set)。不能所有数据都拿来训练,要留一部分验证模型的准确率,这部分叫测试集。

2、典型的线性回归模型

下面是几个典型的线性回归模型:

  • 最基本的单变量线性回归:
    • $y=b x + e$
  • 多变量线性回归:
    • $\hat{h}=\theta_{0}+\theta_{1} x_{1}+\theta_{2} x_{2}+\cdots+\theta_{n} x_{n}+e$
  • 多项式回归(Polynomial Regression):
    • $\hat{h}=\theta_{0}+\theta_{1} x^{1}+\ldots+\theta_{n-1} x^{n-1}+\theta_{n} x^{n}$
    • 我们可以令$x_{2}=x^{2}$,$x_{3}=x^{3}$将其转化为了线性回归模型。

最终通用表达式就是:$\hat{h}=\theta^{T}X=\theta_{0}+\theta_{1} x_{1}+\theta_{2} x_{2}+\cdots+\theta_{n} x_{n}+e$

二、sklearn实现简单线性回归

在进行多元线性回归之间通过简单线性回归来展现线性回归的特性和结果,后面再延伸至多元线性回归。

1、包导入

这里导入进行线性回归的包,我们利用pandas和numpy对数据进行操作,使用matplotlib进行图像化,使用sklearn进行数据集训练与模型导入。

1
2
3
4
5
6
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression

2、数据集创建

这里我们创建一个数据集来描述学生学习时间与成绩的关系并且做简单的线性回归。

1
2
3
4
5
6
7
#创建数据集
examDict = {'学习时间':[0.50,0.75,1.00,1.25,1.50,1.75,1.75,2.00,2.25,2.50,2.75,3.00,3.25,3.50,4.00,4.25,4.50,4.75,5.00,5.50],
'分数':[10,22,13,43,20,22,33,50,62,48,55,75,62,73,81,76,64,82,90,93]}

#转换为DataFrame的数据格式
examDf = pd.DataFrame(examDict)
examDf.head()
分数学习时间
0100.50
1220.75
2131.00
3431.25
4201.50

从上面的数据可以看到数据的特征值标签,学生的学习时间就是所需要的特征值,而分数就是通过特征值所反应的标签。

3、数据集散点图

这里我们利用散点图来简单看下学习时间与成绩的情况。

1
2
3
4
5
6
7
#绘制散点图
plt.scatter(examDf['分数'],examDf['学习时间'],color = 'b',label = "Exam Data")
#添加图的标签(x轴、y轴)
plt.xlabel("Hours")
plt.ylabel("Score")
#显示图像
plt.show()

从上图可以看到,分数和时间存在一定的相关性。下面我们看下两个变量因素的相关密切程度。

4、相关系数

相关系数的公式如下:

$$
P_{x y}=\frac{\operatorname{cov}(X, Y)}{\sigma_{x} \sigma_{y}}
$$

我们可以简单的认为:

  • 0~0.3 弱相关
  • 0.3~0.6 中等程度相关
  • 0.6~1 强相关
1
2
rDf = examDf.corr()
print(rDf)
分数学习时间
分数1.0000000.923985
学习时间0.9239851.000000

pandas中的数学统计函数DataFrame.corr()可以反应数据间的相关性关系,可从表值中反应出学习时间与分数之间的相关性为强相关(0.6~1)。

5、最小二乘法

简单回归方程为:y = a + b*x,而这个最佳拟合线是通过最小二乘法来得到。

最小二乘法的几个概念:

  • 点误差点误差 = 实际值 - 预测值
  • 误差平方和(Sum of square error)SSE = Σ(实际值-预测值)^2

最小二乘法就是使得误差平方和最小来得到最佳拟合。

6、训练集、测试集划分

这里使用sklearn中的train_test_split函数来划分训练、测试集。

1
2
3
4
5
6
7
8
exam_X = examDf['学习时间']
exam_Y = examDf['分数']

# 原数据集拆分训练集和测试集
X_train,X_test,Y_train,Y_test = train_test_split(exam_X,exam_Y,train_size = 0.8,test_size = 0.2)

print("原始数据特征:",exam_X.shape,",训练数据特征:",X_train.shape,",测试数据特征:",X_test.shape)
print("原始数据标签:",exam_Y.shape,",训练数据标签:",Y_train.shape,",测试数据标签:",Y_test.shape)
1
2
原始数据特征: (20,) ,训练数据特征: (16,) ,测试数据特征: (4,)
原始数据标签: (20,) ,训练数据标签: (16,) ,测试数据标签: (4,)

下面是可视化。

1
2
3
4
5
6
7
8
9
10
#散点图
plt.scatter(X_train, Y_train, color="blue", label="train data")
plt.scatter(X_test, Y_test, color="red", label="test data")

#添加图标标签
plt.legend(loc=2)
plt.xlabel("Hours")
plt.ylabel("Pass")
#显示图像
plt.show()

注意:由于训练集随机分配,每一次运行的结果(点的分布情况,训练集内的情况,测试集内的情况)不都相同。

7、模型训练

在创建数据集之后我们需要将训练集放入skleran中的线性回归模型LinearRegression()进行训练,使用函数中的.fit方法进行模型的训练。

1
2
3
4
5
6
model = LinearRegression()

X_train = X_train.values.reshape(-1,1)
X_test = X_test.values.reshape(-1,1)

model.fit(X_train,Y_train)
1
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

注意:因为model需要二维的数组来进行拟合但是这里只有一个特征所以需要reshape来转换为二维数组。

在模型训练完成之后会得到所对应的方程式(线性回归方程式)需要利用函数中的intercept_coef_来得到。

1
2
3
a = model.intercept_  #截距
b = model.coef_ #回归系数
print("最佳拟合线:截距",a,",回归系数:",b)
1
最佳拟合线:截距 6.837689362410515 ,回归系数: [17.11669719]

由上述的最佳拟合线的截距和回归系数可以算出其线性回归线方程:y = 7.56 + 16.28*x

8、模型评价

接下来需要对模型进行预测并评价。决定系数是反映模型拟合优度的重要的统计量。下面是决定系数的简单介绍。

假设一数据集包括y1,...,yn共n个观察值,相对应的模型预测值分别为f1,...,fn。定义残差ei = yi − fi,然后可以推得:

  • 平均观察值:$\overline{y}=\frac{1}{n} \sum_{i=1}^{n} y_{i}$
  • 总平方和(总波动):$S S_{\mathrm{tot}}=\sum_{i}\left(y_{i}-\overline{y}\right)^{2}$
  • 回归平方和:$S S_{\mathrm{reg}}=\sum_{i}\left(f_{i}-\overline{y}\right)^{2}$
  • 残差平方和:$S S_{\mathrm{res}}=\sum_{i}\left(y_{i}-f_{i}\right)^{2}=\sum_{i} e_{i}^{2}$

由此,决定系数可定义为

$$
R^{2} \equiv 1-\frac{S S_{\mathrm{res}}}{S S_{\mathrm{tot}}}
$$

线性回归(右侧)的效果比起平均值(左侧)越好,决定系数的值就越接近于1。 蓝色正方形表示线性回归的残差的平方,红色正方形数据表示对于平均值的残差的平方

从上述定义可以看出,R2取值在0到1之间,且无单位,其数值大小反映了回归贡献的相对程度,即在因变量Y的总变异中回归关系所能解释的百分比。R2是最常用于评价回归模型优劣程度的指标,R2越大(接近于1),所拟合的回归方程越优。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#训练数据的预测值
y_train_pred = model.predict(X_train)
#绘制最佳拟合线:标签用的是训练数据的预测值y_train_pred
plt.plot(X_train, y_train_pred, color='black', linewidth=3, label="best line")

#测试数据散点图
plt.scatter(X_train, Y_train, color="blue", label="train data")
plt.scatter(X_test, Y_test, color='red', label="test data")

#添加图标标签
plt.legend(loc=2)
plt.xlabel("Hours")
plt.ylabel("Score")
#显示图像
plt.show()
1
2
score = model.score(X_test,Y_test)
print(score)
1
0.8161998892697531

三、sklearn实现多元线性回归

上面我们了解了线性回归相关的分析流程,接下来对多元线性回归进行分析。

1、数据导入

下面我们从http://www-bcf.usc.edu/~gareth/ISL/Advertising.csv 来下载数据集Advertising.csv,其数据描述了一个产品的销量与广告媒体的投入之间影响。

1
2
3
4
5
6
7
8
9
10
11
12
13
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression

adv_data = pd.read_csv("../Advertising.csv")

#清洗不需要的数据
new_adv_data = adv_data.iloc[:,1:]

print(new_adv_data.head())
1
2
3
4
5
6
      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
1
print('\nShape:',new_adv_data.shape)
1
Shape: (200, 4)
  • 标签值(sales):

    • Sales:对应产品的销量
  • 特征值(TV,Radio,Newspaper):

    • TV:对于一个给定市场中单一产品,用于电视上的广告费用(以千为单位)
    • Radio:在广播媒体上投资的广告费用
    • Newspaper:用于报纸媒体的广告费用

在这个案例中,通过不同的广告投入,预测产品销量。因为响应变量是一个连续的值,所以这个问题是一个回归问题。数据集一共有200个观测值,每一组观测对应一个市场的情况。

接下来我们对数据进行描述性统计,以及寻找缺失值(缺失值对模型的影响较大,如发现缺失值应替换或删除),且利用箱图来可视化查看数据集,在描述统计之后对数据进行相关性分析,以此来查找数据中特征值与标签值之间的关系。

2、数据描述

1
2
#数据描述
new_adv_data.describe()
TVradionewspapersales
count200.000000200.000000200.000000200.000000
mean147.04250023.26400030.55400014.022500
std85.85423614.84680921.7786215.217457
min0.7000000.0000000.3000001.600000
25%74.3750009.97500012.75000010.375000
50%149.75000022.90000025.75000012.900000
75%218.82500036.52500045.10000017.400000
max296.40000049.600000114.00000027.000000

3、缺失值检验

1
2
#缺失值检验
new_adv_data[new_adv_data.isnull()==True].count()
1
2
3
4
5
TV           0
radio 0
newspaper 0
sales 0
dtype: int64

4、箱形图

1
2
new_adv_data.boxplot()
plt.show()

5、相关系数

1
new_adv_data.corr()
TVradionewspapersales
TV1.0000000.0548090.0566480.782224
radio0.0548091.0000000.3541040.576223
newspaper0.0566480.3541041.0000000.228299
sales0.7822240.5762230.2282991.000000

从corr表中看出,TV特征和销量是有比较强的线性关系的,而Radio和Sales线性关系弱一些但是也是属于强相关的,Newspaper和Sales线性关系更弱。

接下来建立散点图来查看数据间的关系,这里使用seaborn的pairplot来绘画3种不同的因素对标签值的影响。

1
2
3
# 通过加入一个参数kind='reg',seaborn可以添加一条最佳拟合直线和95%的置信带。
sns.pairplot(new_adv_data, x_vars=['TV','radio','newspaper'], y_vars='sales', size=7, aspect=0.8,kind = 'reg')
plt.show()

如上图所示,可以了解到不同的因素对销量的预测线(置信度= 95 %),也可以大致看出不同特征对于标签值的影响。

6、建立模型

在了解了数据的各种情况后需要对数据集建立模型,建立模型的第一步就是建立训练集与测试集,同样的我们将会使用sklearn中的train_test_split函数来创建。

1
2
3
4
X_train,X_test,Y_train,Y_test = train_test_split(new_adv_data.iloc[:,:3],new_adv_data.sales,train_size = 0.8,test_size = 0.2)

print("原始数据特征:",new_adv_data.iloc[:,:3].shape,",训练数据特征:",X_train.shape,",测试数据特征:",X_test.shape)
print("原始数据标签:",new_adv_data.sales.shape,",训练数据标签:",Y_train.shape,",测试数据标签:",Y_test.shape)
1
2
原始数据特征: (200, 3) ,训练数据特征: (160, 3) ,测试数据特征: (40, 3)
原始数据标签: (200,) ,训练数据标签: (160,) ,测试数据标签: (40,)

建立初步的数据集模型之后将训练集中的特征值与标签值放入LinearRegression()模型中且使用fit函数进行训练,在模型训练完成之后会得到所对应的线性回归方程式的截距intercept_与回归系数coef_

1
2
3
4
5
6
model = LinearRegression()
model.fit(X_train,Y_train)

a = model.intercept_ #截距
b = model.coef_ #回归系数
print("最佳拟合线:截距",a,",回归系数:",b)
1
最佳拟合线:截距 2.7283555357741935 ,回归系数: [ 0.04665892  0.19082797 -0.00099623]

这里我们得到的多元线性回归模型的函数为:y = 2.79 + 0.04 * TV + 0.187 * Radio - 0.002 * Newspaper

对于给定了Radio和Newspaper的广告投入,如果在TV广告上每多投入1个单位,对应销量将增加0.04711个单位。就是加入其它两个媒体投入固定,在TV广告上每增加1000美元,销量将增加47.11。但是大家注意这里的newspaper的系数居然是负数,所以我们可以考虑不使用newspaper这个特征。

7、模型评估

接下来对数据集进行预测与模型测评。同样使用predict与score函数来获取所需要的预测值与得分。

1
2
score = model.score(X_test,Y_test) 
print(score)
1
0.9028750733729708

下面简单看下这个训练好的模型的预测值,然后和实际值进行可视化对比:

1
2
Y_pred = model.predict(X_test)
print(Y_pred)
1
2
3
4
5
6
7
[ 5.88658387 19.04199354 19.52640033 17.70822778 10.71716318  5.2030417
24.26646035 12.76256323 9.82255698 18.52225153 21.3759422 7.45340113
4.28180262 19.73124396 20.7434211 20.46565901 8.02140753 21.32827151
15.59515058 8.76729192 18.44116227 10.52302626 13.16705012 15.18864077
9.9417295 15.2734619 12.1230943 10.46509366 11.64698928 13.9265036
15.06236293 6.66457005 12.6079046 18.00988685 17.10350626 9.24127309
11.36598772 18.937225 6.43101507 6.89520345]
1
2
3
4
5
6
plt.plot(range(len(Y_pred)),Y_pred,'b',label="predict")
plt.plot(range(len(Y_test)),Y_test,'r',label="test")
#显示图像
#添加图标标签
plt.legend(loc=2)
plt.show()

四、参考

赞赏一杯咖啡
0%