对于一个样本点 $x^{(i)}$ ，我们根据直线方程都可以得到一个预测值 $\hat{y}^{(i)}=a x^{(i)}+b$ 。直线方程能够尽可能地拟合真实情况，也就是说真实值 $y^{(i)}$ 和预测值 $\hat{y}^{(i)}$ 之间的差距应尽量小。所有的样本的误差都小，说明直线方程拟合性好。而一条直线的具体表示又是由其斜率a和截距b来决定的。因此，这个问题就变成了当a与b取什么值时，样本误差最小。

由于描述样本误差的指标有很多，因此我们需要先挑选一个合适的。通常来说，为了防止出现正误差值和负误差值相抵的情况，使用绝对值来表示误差，即 $\left|y^{(i)}-\hat{y}^{(i)}\right|$ 。可是，通过一个绝对值函数来求最小值往往会存在一个问题，那就是不能求导。因为绝对值函数不是一个处处可导的函数。为了保证函数可导从而方便通过导数找到最小值，我们自然想到使用 $\left(y^{(i)}-\hat{y}^{(i)}\right)^{2}$ 来替换绝对值。因此，此时对于所有的n个样本而言，误差函数就变成了

$\sum_{i=1}^{n}\left(y^{(i)}-\hat{y}^{(i)}\right)^{2}=\sum_{i=1}^{n}\left(y^{(i)}-a x^{i}-b\right)^{2}$

此时的误差函数也被称之为损失函数或者目标函数，能令损失函数最小的参数a和b也就是最佳直线的斜率和截距。这时，可以用令偏导数为0的方式求出满足条件的a和b。具体推导过程见2.5。结果如下：

$a=\frac{\sum_{i=1}^{n}\left(x^{i}-\bar{x}\right)\left(y^{i}-\bar{y}\right)}{\sum_{i=1}^{n}\left(x^{(i)}-\bar{x}\right)^{2}}=\frac{\bar{xy}-\bar{x}\bar{y}}{\bar{x^2}-\bar{x}^2}$

$b=\bar{y}-a \bar{x}$

2.3 千奇百怪的单样本损失函数

首先需要注意，这部分讲的损失函数不同于2.2中的简单线性回归模型的损失函数。这里的损失函数是对于单个样本而言的。损失函数是衡量预测模型预测期望结果表现的指标。损失函数越小，模型的鲁棒性越好[1]。常见的单样本损失函数有以下几种：

2.3.1 0-1损失函数

分类问题常用的损失函数，分类正确结果为0，错误为1。

$L(y, f(x))=\left\{\begin{array}{ll} 1, & y \neq f(x) \\ 0, & y=f(x) \end{array}\right.$

2.3.2 平方损失函数

用来描述回归问题，为真实值与预测值差值的平方。误差值越大、惩罚力度越强，即对表现差的值更敏感。

$L(y, f(x))=(y-f(x))^{2}$

2.3.3 绝对损失函数

可以用来描述回归问题，为真实值与预测值差的绝对值。

$L(y, f(x))=|y-f(x)|$

2.3.4 对数损失函数

可以用来描述分类问题。该损失函数用到了极大似然估计的思想。P(Y|X)通俗的解释就是：在当前模型的基础上，对于样本X，其预测值为Y，也就是预测正确的概率。由于概率之间的同时满足需要使用乘法，为了将其转化为加法，我们将其取对数。最后由于是损失函数，所以预测正确的概率越高，其损失值应该是越小，因此再加个负号取个反。

$L(y, f(x))=-\log P(y | x)$

另外，需要说明的是，对于多样本的对数损失函数需要假设所有样本是相互独立的。

2.4 侧重不同的多样本风险函数

损失函数是针对于单个样本的，但是一个训练数据集中存在多个样本，多个样本会给出多个损失。此时应该如何进行选择呢？为了解决这个问题，于是引出了风险函数的概念[1]。

2.4.1 期望风险

期望风险是损失函数的期望，用来表达理论上模型f(x)关于联合分布P(x, y)的平均意义下的损失。又叫期望损失/风险函数。

$R_{exp}(f)=E_{p}(L(y, f(x)))=\int_{x\times y} L(y, f(x)) P(x, y)dxdy$

2.4.2 经验风险

模型f(x)关于训练数据集的平均损失，称为经验风险或经验损失。其含义为：模型关于训练集的平均损失。具体来说就是将每个样本的损失加起来，然后平均一下。

$R_{\mathrm{emp}}(f)=\frac{1}{n} \sum_{i=1}^{n} L\left(y_{i}, f\left(x_{i}\right)\right)$

可以观察到该公式与2.2中简单线性回归模型的损失函数之间的关系。事实上，经验风险最小的模型为最优模型。在训练集上最小经验风险最小，也就意味着预测值和真实值尽可能接近，模型的效果越好。因此，公式含义还可以解释为训练样本集中对数损失函数平均值的最小值。

$\min _{f \in F} \frac{1}{n} \sum_{i=1}^{n} L\left(y_{i}, f\left(x_{i}\right)\right)$

期望风险是模型关于联合分布的期望损失，经验风险是模型关于训练样本数据集的平均损失。根据大数定律，当样本容量趋于无穷时，经验风险趋于期望风险。因此很自然地想到用经验风险去估计期望风险。

2.4.3 结构风险

由于训练样本个数有限，用经验风险去估计期望风险容易出现过拟合问题。因此需要对经验风险其进行修正。而结构风险指的是在经验风险上加上一个正则化项或者惩罚项。这样便可以实现修正。

$\min _{f \in F} \frac{1}{n} \sum_{i=1}^{n} L\left(y_{i}, f\left(x_{i}\right)\right)+\lambda J(f)$

2.5 最小二乘法

最小二乘法源于法国数学家阿德里安的猜想：对于测量值来说，让总的误差的平方最小的就是真实值。这是基于，如果误差是随机的，应该围绕真值上下波动。

意思就是说，如果误差是随机的，那么能让总的误差最小的测量值就可以认为是真实值。

高斯证明过：如果误差的分布是服从正态分布的，那么最小二乘法得到的值就最有可能是真实的值。

那，什么是“二乘”呢？其实就是“平方”的意思。是不是有些恍然大悟，明白了最小二乘法用的是什么误差计算方式。其实就是2.2中提到的真实值与预测值差的平方。

现在，我们承接2.2中的问题：如何利用偏导数为0求出令损失函数最小的a和b？具体过程参考如下：

首先，求解b：

$\begin{aligned} \frac{\partial L(a, b)}{\partial b}&=\sum_{i=1}^{n} 2\left(y^{(i)}-a x^{(i)}-b\right)(-1) \\ &=-2\sum_{i=1}^{n}\left(y^{(i)}-a x^{(i)}-b\right) \\ &=-2(\sum_{i=1}^{n} y^{(i)}-a \sum_{i=1}^{n} x^{(i)}-n b) \\ &=0 \end{aligned}$

$n b=\sum_{i=1}^{n} y^{(i)}-a \sum_{i=1}^{n} x^{(i)}$

这里我们用样本均值来代替求和，即用 $n \bar{x}$ 和 $n \bar{y}$ 来替换 $\sum x^{(i)}$ 和 $\sum y^{(i)}$ 。可得，

$b=\bar{y}-a \bar{x}$

然后，求解a：

$\begin{aligned} \frac{\partial L(a, b)}{\partial a} &=\sum_{i=1}^{n} 2\left(y^{(i)}-a x^{(i)}-b\right)\left(-x^{(i)}\right) \\ &=-2(\sum_{i=1}^{n}\left(y^{(i)}-a x^{(i)}-b\right) x^{(i)})\\ &=0 \end{aligned}$

将 $b=\bar{y}-a \bar{x}$ 代入到上式可得，

$\sum_{i=1}^{n}\left(y^{(i)}-a x^{(i)}-\bar{y}+a \bar{x}\right) x^{(i)}=0$

打开括号，整理得，

$a=\frac{\sum_{i=1}^{n}\left(x^{(i)} y^{(i)}-x^{(i)} \bar{y}\right)}{\sum_{i=1}^{n}\left(\left(x^{(i)}\right)^{2}-\bar{x} x^{(i)}\right)}$

用同上的替换方式去除求和符号得，

$a=\frac{\bar{xy}-\bar{x}\bar{y}}{\bar{x^2}-\bar{x}^2}$

2.6 简单线性回归的实现及应用（以波士顿房价数据为例）

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split

boston = datasets.load_boston()
X = boston.data[:, 12]
y = boston.target

X = X[y < 50].reshape(-1, 1)
y = y[y < 50].reshape(-1, 1)
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state = 2020)

fig, axs = plt.subplots(2, 1)

# 手动计算
a = (np.mean(X_train * y_train) - np.mean(X_train) * np.mean(y_train)) / \
    (np.mean(X_train ** 2) - np.mean(X_train) ** 2)
b = np.mean(y_train) - a * np.mean(X_train)
print(a, b)

axs[0].scatter(X_test, y_test, c = 'b')
axs[0].plot(X_test, a * X_test + b, c = 'r')

# 调用sklearn模块
lr = LinearRegression()
lr.fit(X_train, y_train)
predictions = lr.predict(X_test)
print(lr.coef_[0], lr.intercept_)

axs[1].scatter(X_test, y_test, c = 'b')
axs[1].plot(X_test, predictions, c = 'r')

plt.show()

3. 多元线性回归

3.1 定义及损失函数

相对于简单线性回归，多元线性回归则是不再只考虑单一样本特征与样本标签的关系。此时，我们可以把含有m个特征的样本n个样本构成的总体X写成如下的矩阵形式：

$X=\left(\begin{array}{ccccc} 1 & x_{1}^{(1)} & x_{2}^{(1)} & \ldots & x_{m}^{(1)} \\ 1 & x_{1}^{(2)} & x_{2}^{(2)} & \ldots & x_{m}^{(2)} \\ \ldots & & & & \ldots \\ 1 & x_{1}^{(n)} & x_{2}^{(n)} & \ldots & x_{m}^{(n)} \end{array}\right)$

其中，每一行的数据都是同一个样本的各个特征。另外，还可以发现第一列全部为1，这样做的目的是为了统一格式，将常数项系数纳入到X中可以得到更加简练的结果[3]。

而此时某样本的预测值表达式如下

$\hat{y}_{(i)}=\theta_{0}+\theta_{1} x_{1}^{(i)}+\theta_{2} x_{2}^{(i)}+\ldots+\theta_{n} x_{n}^{(i)}$

如果将样本总体考虑在内的话，可以写成如下的矩阵形式

$\hat{y}=X\cdot\theta$

此时，我们仍然使用2.2中的损失函数的概念，则在多元线性回归中，写成矩阵形式的损失函数表达式为

$\left(y-X \cdot \theta\right)^{T}\left(y - X \cdot \theta\right)$

同理，这里一样可以先打开括号然后再利用导数求最小值。只不过这里涉及到矩阵和向量的求导。相关的求导规则可以参考《the Matrix Cookbook》[5]。

这里略过过程，最后的结果是

$\theta=\left(X^{T} X\right)^{-1} X^{T} y$

关于该最优参数解的形式的进一步讨论见3.3。

3.2 多元线性回归的实现及应用（以波士顿房价数据为例）

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split

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

X = X[y < 50]
y = y[y < 50]

# 手动计算
# 在数据第一列上添加1
X = np.concatenate((np.array([1] * X.shape[0]).reshape(-1, 1), X), axis = 1)
# 计算参数theta
theta = np.dot(np.linalg.inv(np.dot(X.T, X)), X.T).dot(y)
coef_ = theta[1:]
intercept_ = theta[0]
print(coef_, intercept_)

# 调用sklearn模块
lr = LinearRegression()
lr.fit(X, y)
predictions = lr.predict(X)
print(lr.coef_, lr.intercept_)