「Machine Learning」线性回归
1、初识线性回归
所谓线性回归（Linear Regression），其最本质的特点就是可以用来根据已有的数据探究一个（或者多个）自变量与因变量之间的线性关系，从而对未知自变量所对应因变量进行预测。以单个自变量为例：
上图中，黄色数据点为真实样本数据，坐标轴横轴表示自变量x，纵轴表示因变量y。观察发现黄色数据点之间具有很强的线性相关关系，那么我们就可以根据已知的数据点，拟合出来一条近似的线性直线（蓝色直线）。有了这条直线，我们就可以输入任意的自变量x来预测到未知的y。
举个栗子：我们想要探究修🐶运动距离与消耗馒头的关系。经过对比实验发现：
| 运动距离（km） | 消耗馒头（个） |
|---|---|
| 1 | 3 |
| 2 | 5 |
| 3 | 6.8 |
修🐶的运动距离与消耗馒头的数量之间好像存在一定的线性关系，那么我们根据这个已有的数据，如果能够得到线性方程$y=ax+b$中的$a$和$b$，可以大概拟合出来一个方程：
$$ bread = 2 × Distance_{dogs} + 1 $$
有了这个方程，我们就可以对修🐶每次运动结束后，需要补充的馒头数量有了一个大概的估算。某天修🐶跑了10公里，虽然我们以前并没有统计过修🐶跑10公里可以消耗多少馒头的数据，但是我们可以根据拟合的线性公式推测，修🐶大概需要消耗21个馒头。
上面所讲述的就一个线性回归问题。
简单来说，就是通过一组有线性关系的自变量和因变量，找出其拟合方程的的$a$和$b$，然后即可通过线性方程$y=ax+b$对未知自变量对应的因变量进行预测。
如果你也是接触生化专业的学生，那么一定对标准曲线非常熟悉，制作标准曲线的过程就是一个典型的线性回归问题。想要了解的同学可以参考我的Blog：Python数据分析——以我硕士毕业论文为例，里面有较为详细的讲解。
线性回归算法的优点：
- 解决回归问题；
- 思想简单，实现容易；
- 许多强大的非线性模型的基础；
- 结果具有很好的可解释性；
- 蕴含Machine Learning中的很多重要思想。
至此，如果你理解了线性回归的大概含义，那么可能会提出一个问题。得到线性方程$y=ax+b$中，如何保证所求得的$a$和$b$可以使得所预测的未知自变量所对应因变量与其真实值的误差最小呢？也就是说，我们该如何量化“足够靠近”？
这个时候，最小二乘法就要登场了。
最小二乘法的推导
损失函数
如上图所示，假设我们找到了最佳拟合方程：
$$ y=ax+b $$
则对应每一个样本点$x^{(i)}$的预测值可以写为：
$$ \hat y^{(i)} = a x^{(i)} + b $$
设样本点$x^{(i)}$的真实值为$y{(i)}$。
我们所希望的就是$\hat y^{(i)}$与$y{(i)}$之间的差值尽可能的小，两者之间的差值可以写为：
$$ y^{(i)} - \hat y^{(i)} $$
为了排除正负相消的情况，加上一个绝对值：
$$ |y^{(i)} - \hat y^{(i)}| $$
但是对于绝对值函数，可能会出现不可导的点（为了求出可以满足两者之间插值最小的$a$和$b$，我们需要借助于求导）。因此，我们可以将一个点预测值与真实值的误差写为：
$$ (y^{(i)} - \hat y^{(i)})^2 $$
那么一组数据中m个点，所有预测值与真实值的误差之和就是：
$$ \sum_{i=1}^{m}(y^{(i)} - \hat y^{(i)})^2 $$
带入$\hat y^{(i)} = a x^{(i)} + b$，可得：
$$ \sum_{i=1}^{m}(y^{(i)} - a x^{(i)} - b)^2 $$
在Machine Learning中，我们一般称上式为损失函数（loss function）或者效用函数（utility function），我们的目标就是要使得所求得的$a$和$b$可以使上式的结果尽可能的小。
通过分析问题，确定问题的损失函数或者效用函数；通过最优化损失函数或者效用函数，获得机器学习的模型。这近乎是所有参数学习算法的套路。
a和b的寻找
对于
$$ \sum_{i=1}^{m}(y^{(i)} - a x^{(i)} - b)^2 $$
而言，其中的$y$与$x$都是具体点所对应的值，因此上式中在求解的$a$与$b$的过程中看作为一个常数。
令：
$$ J(a,b)=\sum_{i=1}^{m}(y^{(i)} - a x^{(i)} - b)^2 $$
我们会发现$J(a,b)$是一个二元（$a$与$b$）二次（$a$的幂为2）函数，求二元函数$J(a,b)$的最小值，那我们只需要对其中的$a$与$b$分别求偏导即可：
$$ \frac{\partial J(a,b)}{\partial a}=0, \frac{\partial J(a,b)}{\partial b}=0 $$
其中对$b$求偏导为：
$$ \frac{\partial J(a,b)}{\partial b}= \sum_{i=1}^{m}2(y^{(i)} - a x^{(i)} - b)(-1) = 0 $$
$$ \\ \sum_{i=1}^{m}(y^{(i)} - a x^{(i)} - b) = 0 \\ \sum_{i=1}^{m}y^{(i)} - a\sum_{i=1}^{m}x^{(i)} - \sum_{i=1}^{m}b = 0 \\ \sum_{i=1}^{m}y^{(i)} - a\sum_{i=1}^{m}x^{(i)} - mb = 0 \\ \frac{\sum_{i=1}^{m}y^{(i)} - a\sum_{i=1}^{m}x^{(i)}}{m} = b \\ $$
而在上式中，如果令$y$与$x$的均值分别为的$\bar{y}$与$\bar{x}$，则：
$$ \bar{y} = \frac{\sum_{i=1}^{m}y^{(i)}}{m}, \bar{x} = \frac{\sum_{i=1}^{m}x^{(i)}}{m} $$
所以带入，即可求得：
$$ b=\bar{y}-a \bar{x} $$
然后再来求$a$：
$$ \frac{\partial J(a,b)}{\partial a} = \sum_{i=1}^{m}2(y^{(i)} - a x^{(i)} - b)(-x^{(i)}) = 0 \\ \sum_{i=1}^{m}(y^{(i)} - a x^{(i)} - b)x^{(i)} = 0 \\ \sum_{i=1}^{m}(y^{(i)} - a x^{(i)} - \bar{y} + a \bar{x})x^{(i)} = 0 \\ \sum_{i=1}^{m}(y^{(i)}x^{(i)} - a (x^{(i)})^2 - \bar{y}x^{(i)} + a \bar{x}x^{(i)}) = 0 \\ \sum_{i=1}^{m}(y^{(i)}x^{(i)} - \bar{y}x^{(i)}) - a\sum_{i=1}^{m}((x^{(i)})^2 - \bar{x}x^{(i)}) = 0 \\ \frac{\sum_{i=1}^{m}(y^{(i)}x^{(i)} - \bar{y}x^{(i)})}{\sum_{i=1}^{m}((x^{(i)})^2 - \bar{x}x^{(i)})} = a \\ $$
又因为：
$$ \sum_{i=1}^{m} x^{(i)} \bar{y}=\bar{y} \sum_{i=1}^{m} x^{(i)}=\bar{y} \cdot m \bar{x}=\bar{x}\sum_{i=1}^{m} y^{(i)} = \sum_{i=1}^{m}\bar{x} \cdot \bar{y} $$
所以带入即可得：
$$ \frac{\sum_{i=1}^{m}(y^{(i)}x^{(i)} - \bar{y}x^{(i)} -\bar xy^{(i)} + \bar x\bar y)} {\sum_{i=1}^{m}((x^{(i)})^2 - \bar{x}x^{(i)} -\bar xx^{(i)} + \bar x^2)} = a \\ \frac{\sum_{i=1}^{m}(x^{(i)} - \bar{x})(y^{(i)} - \bar{y})} {\sum_{i=1}^{m}(x^{(i)} - \bar{x})^2} = a \\ $$
至此，我们已经可以根据一组数据$y$与$x$得输入，来通过求和和求均值等操作来求出其拟合方程$y=ax+b$。
纯Python实现简单的线性回归
创建一组数组：
x = np.array([1., 2, 3, 4, 5])
y = np.array([1., 3, 2, 3, 5])绘出图形观察下数据：
plt.scatter(x, y)
plt.axis([0, 6, 0, 6])
plt.show()分别计算$a$与$b$的值：
# Calculate a.
d, m = 0, 0
for x_i, y_i in zip(x, y):
d += (x_i - x.mean()) * (y_i - y.mean())
m += (x_i - x.mean()) ** 2
a = d / m
# Calculate b.
b = y.mean() - a * x.mean()
a, b所以说$\hat y$就等于：
y_hat = a * x + b
y_hat查看拟合的直线：
plt.scatter(x, y)
plt.plot(x, y_hat, color='r')
plt.axis([0, 6, 0, 6])
plt.show()封装类对象
我们可以将上面所实现的简单线性回归模型封装为一个类：
class SimpleLinearRegression:
def __init__(self):
"""
Initialization parameters
"""
self.a_ = None
self.b_ = None
def fit(self, x_train, y_train):
"""
Fit the data and calculate self.a_ and self.b_.
:param x_train: data of x.
:param y_train: data of y.
:return: Object
"""
# Data type conversion.
if type(x_train) != np.ndarray: x_train = np.array(x_train)
if type(y_train) != np.ndarray: y_train = np.array(y_train)
# Check the length of data.
assert x_train.ndim == 1 == y_train.ndim, 'The dimension of x and y must be one!!!'
assert len(x_train) == len(y_train), 'The length of x and y must be equal!!!'
# Calculate a and b.
d, m = 0, 0
for x, y in zip(x_train, y_train):
d += (x - x_train.mean()) * (y - y_train.mean())
m += (x - x_train.mean()) ** 2
self.a_ = d / m
self.b_ = y_train.mean() - a * x_train.mean()
return self
def predict(self, x_test):
"""
Using input x_test to calculate y_test, y_test = self.a_ * x_test + self.b_
:param x_test: Input x_test data.
"""
assert self.a_ is not None, 'Please run object.fit() before predict!!!'
y_test = np.array([self.predict_(x) for x in x_test])
return y_test
def predict_(self, x):
return self.a_ * x + self.b_
def __repr__(self):
return 'SimpleLinearRegression'
demo_ob = SimpleLinearRegression()
# demo_ob.fit(x, y)
demo_ob.fit_vector(x, y)
print(demo_ob.a_, demo_ob.b_)
demo_ob.predict_(x)在上面的计算方法中，我们使用for循环计算d与m的值，这是一种十分低效的手段。我们可以借助Numpy的向量化计算来优化：
def fit_vector(self, x_train, y_train):
"""
Fit the data and calculate self.a_ and self.b_.
:param x_train: data of x.
:param y_train: data of y.
:return: Object
"""
# Data type conversion.
if type(x_train) != np.ndarray: x_train = np.array(x_train)
if type(y_train) != np.ndarray: y_train = np.array(y_train)
# Check the length of data.
assert x_train.ndim == 1 == y_train.ndim, 'The dimension of x and y must be one!!!'
assert len(x_train) == len(y_train), 'The length of x and y must be equal!!!'
# Calculate a and b.
x_mean, y_mean = np.mean(x_train), np.mean(y_train)
d = (x_train - x_mean).dot(y_train - y_mean)
m = (x_train - x_mean).dot(x_train - x_mean)
self.a_ = d / m
self.b_ = y_mean - self.a_ * x_mean
return self在上面的类中添加一个fit_vector()方法，与fit()不同的是，fit_vector()方法采用向量化运算求d与m的值。
查看两者的效率对比：
NUM = 100000
rand_x = np.random.random(size=NUM)
rand_y = rand_x * 2.0 + 3.0 + np.random.normal(size=NUM)查看使用for循环fit()方法的效率：
%%time
demo_ob = SimpleLinearRegression()
demo_ob.fit(rand_x, rand_y)
print(demo_ob.a_, demo_ob.b_)
# demo_ob.predict_(rand_x)1.9940105008360216 3.6060045680996287
CPU times: total: 7.64 s
Wall time: 22.4 s
查看使用向量化运算fit_vector()方法的效率：
%%time
demo_ob = SimpleLinearRegression()
demo_ob.fit_vector(rand_x, rand_y)
print(demo_ob.a_, demo_ob.b_)
# demo_ob.predict_(rand_x)1.994010500836032 3.6060045680996287
CPU times: total: 0 ns
Wall time: 1.84 ms
线性回归算法的评估
由上面的讲解过程可得，线性回归算法的误差产生于：
$$ \sum_{i=1}^{m}(y^{(i)}_{test} - \hat y^{(i)}_{test})^2 $$
由此产生了几种评价误差的指标。分别是：
1、均方误差（Mean Squared Error, MSE）
$$ \frac{1}{m}\sum_{i=1}^{m}(y^{(i)}_{test} - \hat y^{(i)}_{test})^2 $$
但是由于MSE的量纲为$y$的平方，并不与$y$统一。RMSE则解决了这个问题。
2、均方根误差（Root Mean Squared Error, RMSE）
$$ \sqrt{MSE_{test}} = \sqrt{\frac{1}{m}\sum_{i=1}^{m}(y^{(i)}_{test} - \hat y^{(i)}_{test})^2} $$
3、平均绝对误差（Mean Absolute Error, MAE）
$$ \frac{1}{m}\sum_{i=1}^{m}|y^{(i)}_{test} - \hat y^{(i)}_{test}| $$
在scikit-learn中可以直接调用类对象使用计算RMSE与MAE：
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_squared_error
mean_absolute_error(y_test, y_predict), mean_squared_error(y_test, y_predict)这个时候如果你学习过KNN或者其他算法，会想到能否也使用0\~1来表述线性回归算法的准确度呢？
当然是可以的：
$$ R^2 = 1-\frac{\sum_{i=1}^{m}(y^{(i)} - \hat y^{(i)})^2}{\sum_{i=1}^{m}(y^{(i)} - \bar y^{(i)})^2}= 1-\frac{\sum_{i=1}^{m}(y^{(i)} - \hat y^{(i)})^2/m}{\sum_{i=1}^{m}(y^{(i)} - \bar y^{(i)})^2/m} = 1-\frac{MSE(\hat y, y)}{Var(y)} $$
- $R^2<=1$
- $R^2$越大越好，当我们的预测模型不犯任何错误是，$R^2$得到最大值1；
- 当我们的模型等于基准模型时，$R^2$为0；
- 如果$R^2<0$，说明我们学习到的模型还不如基准模型，此时，很有可能我们的数据不存在任何线性关系。
在scikit-learn中可以直接调用类对象使用计算：
from sklearn.metrics import r2_score
r2_score(y_test, y_predict)实际上，$R^2<=1$也是scikit-learn中LinearRegression所封装的方法。
多元线性回归
算法讲解
上面我们讲述了单个自变量的线性回归，然后实际应用中，会有多个自变量引起应变量变化的情况。例如，房间的数量、距离市中心的距离、房子的年龄等多种因素共同决定着房子的售价。令自变量：
$$ X^{(i)} = (X_1^{(i)}, X_2^{(i)}, X_3^{(i)}, \dots , X_n^{(i)}) $$
所以说此时对应的$\hat y$就变成了：
$$ \hat y = \theta_0 + \theta_1 X_1 + \theta_2 X_2 + \dots + \theta_n X_n $$
这个时候的目标就变成了，找到一组$\theta_0, \theta_1, \theta_2, \dots , \theta_n$，使得
$$ \sum_{i=1}^{m}(y^{(i)} - \hat y^{(i)})^2 $$
足够小。
令$X_0=1$，则
$$ \begin{equation*} \begin{split} \hat y &= \theta_0 + \theta_1 X_1 + \theta_2 X_2 + \dots + \theta_n X_n \\ &= \theta_0 X_0 + \theta_1 X_1 + \theta_2 X_2 + \dots + \theta_n X_n \\ &= X^{(i)} \cdot \theta \end{split} \end{equation*} $$
所以说：
$$ X= \begin{bmatrix} 1& X_1^{(1)}& X_2^{(1)}& \dots& X_n^{(1)}& \\ 1& X_1^{(2)}& X_2^{(2)}& \dots& X_n^{(2)}& \\ \vdots& \vdots& \vdots& \ddots & \vdots& \\ 1& X_1^{(m)}& X_2^{(m)}& \dots& X_n^{(m)}& \end{bmatrix}, \theta = \begin{bmatrix} \theta_0\\ \theta_1\\ \theta_2\\ \dots\\ \theta_n \end{bmatrix}, \hat y = X \cdot \theta $$
优化目标就变成了使
$$ (y - X \cdot \theta)^T(y - X \cdot \theta) $$
尽量小，这个时候
$$ \theta = (X^TX)^{-1}X^Ty $$
这个$\theta$被称为多元线性回归的正规方程解（Normal Equation）。其中$\theta_0$是截距（Intercept），$\theta_1, \theta_2, \dots , \theta_n$是系数（Coefficients）。
具体的求解过程，参考「深度学习」PyTorch笔记-02-线性回归模型。
上述$\theta$的求解过程，使用线性代数的知识更容易理解，以一元回归方程为例：
$$ \begin{equation*} \begin{split} y &= X \begin{bmatrix} b\\ a \end{bmatrix} \\ X^T y &= X^TX \begin{bmatrix} b\\ a \end{bmatrix} \\ \begin{bmatrix} b\\ a \end{bmatrix} &= (X^TX)^{-1}X^Ty \end{split} \end{equation*} $$
类的封装
封装一个多元线性回归求解的类：
class MulLinearReg:
"""
The object can support the multiple linear regression
"""
def __init__(self):
self._theta = None
self.coef_ = None
self.interception_ = None
def fit_norm(self, X_train, y_train):
assert X_train.shape[0] == y_train.shape[0], \
'The size of X and y must be equal!!!'
X_b = np.hstack([np.ones([X_train.shape[0], 1]), X_train])
self._theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y_train)
self.interception_ = self._theta[0]
self.coef_ = self._theta[1:]
return self
def predict(self, X_predict):
assert self.interception_ is not None and self.coef_ is not None, \
'You must fit this object before predict!!!'
assert X_predict.shape[1] == len(self.coef_), \
'The feature of X_predict must be equal to X_train!!!'
X_b = np.hstack([np.ones([X_predict.shape[0], 1]), X_predict])
return X_b.dot(self._theta)
def __repr__(self):
return 'MulLinearReg()'在sk-learn中使用
from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(X_train, y_train)
y_predict = lin_reg.predict(X_test)
r2_score(y_test, y_predict)
lin_reg.score(X_test, y_test)