机器学习笔记（九）：多元线性回归 | 凌云时刻

凌云时刻 · 技术

导读：这篇笔记主要介绍线性回归算法，对在第一篇笔记中介绍过的线性回归算法进行实现。kNN算法主要解决的是分类问题，并且它的结果不具备良好的解释性。线性回归算法主要解决回归问题，它的结果具有良好的可解释性，和kNN算法的介绍过程一样，线性回归算法也蕴含了机器学习中很多重要的思想，并且它是许多强大的非线性模型的基础。

作者 | 计缘

来源 | 凌云时刻（微信号：linuxpk）

多元线性回归

上面我们讲了简单线性回归，也就是只关注样本数据的一个特征，这一节我们来看看多元线性回归，既关注样本数据的多个特征。

 

像上图中展示的，简单线性回归我们只关注   一个特征，假如   是一个特征向量，那此时就变成了关注多个特征的多元线性回归问题，如下图所示：

此时相对与简单线性回归的公式   ，因为   成为了特征行向量，所以我们将aa也看作特征系数列向量，那么公式展开后可以写成这样：

   

‍‍我们按惯例，将多元线性回归公式中的特征系数称为   :

   

那么描述每行样本数据中每个特征和其预测值的公式为：

   

我们对上面的公式再进一步做一下处理，将截距   也乘以一个特征   ，不过该特征值恒等于1:

   

我们将上面公式用矩阵的方式表示为（用到矩阵相乘的知识，1行n列矩阵乘以n行1列矩阵为1行1列矩阵，既一个标量）：

   

我们可以直接将上面公式表示为：

   

我们知道简单线性回归就是使   尽可能小，那多元线性回归也是一样的，只不过多元线性回归中的   的公式不一样而已，将上面的公式代入：

   

将上面公式里的平方展开：

   

矩阵的相乘求和可以转换为如下形式（矩阵或向量的相同元素相乘再求和等于该向量或矩阵的转置乘以该向量或矩阵）：

   

所以多元线性回归就是求使上面公式尽可能小的   列向量。

上面的公式出现了矩阵乘积转置，看看下面这张手记就能明白如何转换矩阵乘积转置了：

 

如此，上面的公式可以转换为：

   

要求   Floss函数的最小值既对该函数中的每项   求偏导数：

   

这里用到了矩阵求导的知识，具体可见下面的手记：

上面两张手记简要介绍了矩阵求导的一种方法，根据求导的类型，然后确定是分子布局还是分母布局，最后查阅常用的矩阵求导表确定结果。

•      ：分子是标量，分母是列向量，是标量/向量类型，并且是分母布局，查表满足   条件，故可得结果为0。

•    分子是标量，分母是列向量，是标量/向量类型，并且是分母布局，查表满足   条件，故可得结果为   。

•    ：分子是矩阵，分母是列向量，是矩阵/向量类型，并且是分母布局，但是根据矩阵乘积转置的规则，可以将分子转换为   ，和第二项偏导表达式相同，故结果为   

  。

•    ：分子为标量，分布为列向量，是标量/向量类型，并且是分母布局，因为   还是矩阵，所以查表满足   条件，故可得结果为   。

 

所以   函数的最小值为：

    

   

便可求得   的值为：

   

 

根据逆矩阵的转换规则可得：

   

上面的公式就是多元线性回归的正规方程解 (Normal Equation)。

 使用正规方程解实现多元线性回归

在实现之前，我们先明确一下   ，它是一个列向量：

   

其中   是多元线性方程的截距（intercept），    θ1到   才是系数（coefficients）。在PyCharm中，我们创建一个类LinearRegression，然后构造函数如下：

 

class LinearRegression:

	def __init__(self):
		# 截距theta0
		self.intercept_ = None
		# 系数，theta1 ... thetaN
		self.coef_ = None
		# theta列向量
		self._theta = None

	def __repr__(self):
		return "LinearRegression()"

然后我们来看训练的过程，既fit方法，这里因为是用正规方程解实现的，所以我们将训练方法的名称取为fit_normal()：

 

def fit_normal(self, X_train, y_train):
	# 根据训练数据集X_train，y_train训练LinearRegression模型
	assert X_train.shape[0] == y_train.shape[0], \
		"特征数据矩阵的行数要等于样本结果数据的行数"

	# 计算X_b矩阵，既将X_train矩阵前面加一列，元素都为一
	X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
	# 实现正规方式解
	self._theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y_train)

	# 取到截距和系数
	self.intercept_ = self._theta[0]
	self.coef_ = self._theta[1:]

	return self

训练的实现过程很简单，首先求出X_b，也就是将X_train矩阵前面填加元素为1的一列，这里我们使用到了Numpy的np.hstack方法，也就是在水平方向组合矩阵，另外使用了np.ones这个快捷创建元素为1的矩阵的方法。然后实现了上文中推导出的正规方程解，其中np.linalg.inv是对矩阵求逆的方法。

然后来看预测的过程：

 

# 给定待预测数据集X_predict，返回表示X_predict的结果向量
	def predict(self, X_predict):
		assert self.intercept_ is not None and self.coef_ is not None, \
		"截距和系数都不为空，表示已经经过了fit方法"
		assert X_predict.shape[1] == len(self.coef_), \
		"要预测的特征数据集列数要与theta的系数数量相等"

		X_b = np.hstack([np.ones((len(X_predict), 1)), X_predict])
		return X_b.dot(self._theta)

首先有两个断言增加健壮性，然后同样是计算出X_b，最后乘以训练出的   就得到了预测的结果。

最后同样增加一个评分的方法，使用我们之前实现的   方法来计算评分：

 

# 根据测试数据集X_test和y_test确定当前模型的准确度
	def score(self, X_test, y_test):
		y_predict = self.predict(X_test)
		return r2_score(y_test, y_predict)

 在Jupyter Notebook中使用LinearRegression

我们使用Scikit Learn中提供的波士顿房价数据来验证我们实现的基于正规方程解的多元线性回归。

 

import numpy as np
from sklearn import datasets

boston = datasets.load_boston()
X = boston.data
y = boston.target

X = X[y < 50.0]
y = y[y < 50.0]

# 拆分训练数据和测试数据
from myML.modelSelection import train_test_split
X_train, y_train, X_test, y_test = train_test_split(X, y, seed=123)

# 调用我们实现的类
from myML.LinearRegression import LinearRegression
reg = LinearRegression()
reg.fit_normal(X_train, y_train)

reg.coef_
# 结果
array([ -1.02162853e-01,   2.51989824e-02,  -4.45146218e-02,
		-1.71466010e-01,  -1.17374943e+01,   4.01098742e+00,
		-2.93108282e-02,  -1.12226201e+00,   2.34868501e-01,
		-1.30958162e-02,  -8.32044126e-01,   8.27684963e-03,
		-3.18223340e-01])

reg.intercept_
# 结果
29.513591626754184

reg.score(X_test, y_test)
# 结果
0.80030886154058956

我们再来看看Scikit Learn中提供的线性回归的使用方式：

 

from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(X_train, y_train)

lin_reg.coef_
# 结果
array([ -1.02162853e-01,   2.51989824e-02,  -4.45146218e-02,
		-1.71466010e-01,  -1.17374943e+01,   4.01098742e+00,
		-2.93108282e-02,  -1.12226201e+00,   2.34868501e-01,
		-1.30958162e-02,  -8.32044126e-01,   8.27684963e-03,
		-3.18223340e-01])

lin_reg.intercept_
# 结果
29.513591626752053

lin_reg.score(X_test, y_test)
# 结果
0.80030886154069325

我们可以看到使用Scikit Learn中的线性回归和我们实现的线性回归结果是一样的，但这里要注意的是Scikit Learn中的线性回归实现方式并不是正规方程解，只是因为样本数据量比较少，所以得到了相同的结果。下一篇笔记会介绍另外一种实现多元线性回归的方式。

线性回归的可解释性

我们将全量的波士顿房价进行预测，求出系数来看一下：

 

from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(X, y)
lin_reg.coef_
# 结果
array([ -1.05574295e-01,   3.52748549e-02,  -4.35179251e-02,
		 4.55405227e-01,  -1.24268073e+01,   3.75411229e+00,
		-2.36116881e-02,  -1.21088069e+00,   2.50740082e-01,
		-1.37702943e-02,  -8.38888137e-01,   7.93577159e-03,
		-3.50952134e-01])

我们看到系数有正有负，表示了正相关和负相关，既正系数越大的特征对结果的影响越大，负系数越小对结果影响越大，我们来看看影响波士顿房价的特征都是什么，先将系数按索引从小到大排序：

 

np.argsort(lin_reg.coef_)
# 结果
array([ 4,  7, 10, 12,  0,  2,  6,  9, 11,  1,  8,  3,  5])

# 按索引顺序排列出对应的房价特征
boston.feature_names[np.argsort(lin_reg.coef_)]
# 结果
array(['NOX', 'DIS', 'PTRATIO', 'LSTAT', 'CRIM', 'INDUS', 'AGE', 'TAX',
	   'B', 'ZN', 'RAD', 'CHAS', 'RM'],
	  dtype='