机器学习 —— 多元线性回归

Posted 2020-10-04 魔法少女小Q

什么是多元线性回归？

　　在线性回归分析中，如果有两个或两个以上的自变量，就称为多元线性回归（multivariable linear regression）。如果我们要预测一套房子的价格，影响房子价格的因素可能包括：面积、卧室数量、层数以及房龄，我们用x₁、x₂、x₃、x₄来代表这4个特征（features）。在这里，我们假设房子的价格和这4个特征值都是线性相关的，并且用h_θ来表示预测房子的价格，那么有： \[{h_\theta }(x) = {\theta _0} + {\theta _1}{x_1} + {\theta _2}{x_2} + {\theta _3}{x_3} + {\theta _4}{x_4}\]

其中θ为线性回归中的参数。我们把x、θ分别表示为列向量的形式：

\[\theta  = {\left( {{\theta _0},{\theta _1},{\theta _2},{\theta _3},{\theta _4}} \right)^T}\]

\[x = {\left( {1,{x_1},{x_2},{x_3},{x_4}} \right)^T}\]

那么我们可以发现 \[{h_\theta }(x) = {\theta ^T} \bullet x\]

 

如何实现多元线性回归？

对于拥有n个特征值、m个数据样本的数据，我们可以用一个m*(n+1)矩阵的形式表示出来，其中矩阵的每一行（即x⁽ⁱ⁾）为一个数据样本，每一列代表一个特征值，x_j⁽ⁱ⁾代表第i个数据样本中第j个特征值（x₀⁽ⁱ⁾=0）

\[X = \left( {\begin{array}{*{20}{c}}
1&{x_{^1}^{(1)}}&{x_2^{(1)}}& \cdots &{x_n^{(1)}} \\
1&{x_{^1}^{(2)}}&{x_2^{(2)}}& \cdots &{x_n^{(2)}} \\
1&{x_{^1}^{(3)}}&{x_2^{(3)}}& \cdots &{x_n^{(3)}} \\
\vdots & \vdots & \vdots & \cdots & \vdots \\
1&{x_{^1}^{(m)}}&{x_2^{(m)}}& \cdots &{x_n^{(m)}}
\end{array}} \right) = \left( {\begin{array}{*{20}{c}}
{{x^{(1)}}^T} \\
{{x^{(2)}}^T} \\
{{x^{(3)}}^T} \\
\vdots \\
{{x^{(m)}}^T}
\end{array}} \right)\]

样本的实际值（对于房价预测来说就是房价）用y来表示：

\[y = {\left( {{y_1},{y_2}, \ldots ,{y_m}} \right)^T}\]

线性回归参数θ为：

\[\theta  = {\left( {{\theta _0},{\theta _1}, \ldots ,{\theta _n}} \right)^T}\]

接下来我们定义代价函数J_θ（cost function）为：

\[{J_\theta }(X) = \frac{1}{{2m}}\sum\limits_{i = 1}^m {{{({h_\theta }({x^{(i)}}) - {y^{(i)}})}^2}}  = \frac{1}{{2m}}\sum\limits_{i = 1}^m {{{({\theta ^T} \bullet {x^{(i)}} - {y^{(i)}})}^2}}  = \frac{1}{{2m}}\sum\limits_{i = 1}^m {{{({x^{(i)T}} \bullet \theta  - {y^{(i)}})}^2}} \]

我们还可以把J_θ利用矩阵表达，获得更为简洁的形式

\[{J_\theta }(X) = \frac{1}{{2m}}{(X \bullet \theta  - y)^T} \bullet (X \bullet \theta  - y)\]

现在，我们只需要让代价函数J_θ最小，就能得到最优的θ参数。那么，要怎样才能使J_θ最小呢？有两个办法，一个是梯度下降法（gradient descent），一个是标准方程法（norm equation）。

J_θ在样本数据X确定时，是关于θ的一个函数，而θ是一个（n+1）维的列向量，也就是说J_θ其实是一个（n+1）元函数。学过微积分的人应该知道，要求一个多元函数的极值，可以令多元函数对每个未知元的偏导数为0，从而求出未知元，再代入函数中，便可求出函数的极值。同样地，我们现在要求J_θ的极值，那么就要先使J_θ关于θ的偏导数，或者说关于θ₀、θ₁、θ₂、...、θn的偏导数都为0。

J_θ对θ_j的偏导数为：

\[\frac{{\delta J}}{{\delta {\theta _j}}} = \frac{1}{m}\sum\limits_{i = 1}^m {({x^{(i)T}} \bullet \theta - {y^{(i)}})x_j^{(i)}} = \frac{1}{m}\left( {\begin{array}{*{20}{c}}
{x_j^{(1)}}&{x_j^{(2)}}& \cdots &{x_j^{(m)}}
\end{array}} \right) \bullet \left( {X\theta - y} \right)\]

则J_θ对θ的偏导数为：

\[\frac{{\delta J}}{{\delta \theta }} = \left( {\begin{array}{*{20}{c}}
{\frac{{\delta J}}{{\delta {\theta _0}}}} \\
{\frac{{\delta J}}{{\delta {\theta _1}}}} \\
\vdots \\
{\frac{{\delta J}}{{\delta {\theta _n}}}}
\end{array}} \right) = \frac{1}{m}\left( {\begin{array}{*{20}{c}}
1&1& \cdots &1 \\
{x_1^{(1)}}&{x_1^{(2)}}& \cdots &{x_1^{(m)}} \\
\vdots & \vdots &{}& \vdots \\
{x_n^{(1)}}&{x_n^{(2)}}&{}&{x_n^{(m)}}
\end{array}} \right)(X\theta - y) = \frac{1}{m}{X^T}(X\theta - y)\]

 

标准方程法（Normal Equation）

既然我们知道了J_θ偏导数的表达公式，那么我们只要让它等于0便可求出我们想要的θ，即令：

\[\begin{gathered}
\frac{{\delta J}}{{\delta \theta }} = \frac{1}{m}{X^T}(X\theta - y) = 0 \hfill \\
X\theta - y = 0 \hfill \\
X\theta = y \hfill \\
{X^T}X\theta = {X^T}y \hfill \\
\theta = {({X^T}X)^{ - 1}}{X^T}y \hfill \\
\end{gathered} \]

注意：在第四步中，两边左乘一个X^T是为了在左边获得一个方阵（X^TX），因为只有方阵才有逆矩阵。

在Python中，只需要两步便可实现这个算法：

1 from numpy.linalg import pinv
2 theta = pinv(X.T.dot(X)).dot(X.T).dot(y)

 

梯度下降法（Gradient Descent）

标准方程法虽然实现起来很简单，但是在面对十万甚至是百万级别的数据样本量时，求解一个如此大的矩阵的逆需要花费大量的时间，这时候就需要考虑梯度下降法了。梯度下降法的思想是一点点、逐步靠近极小值，通过多次迭代，不断更新θ的值，让代价函数J_θ收敛于极小值。迭代的公式为：

\[\theta  = \theta  - \alpha \frac{{\delta J}}{{\delta \theta }}\]

我们已经知道了如何求J_θ关于θ的偏导数，只要一直迭代下去，直到J_θ等于0或者足够的小（比如小于1e^-5），最后我们便得到了想要的θ值。

以下我用Python实现了用梯度下降法求θ：

 1 import numpy as np
 2 
 3 def partial_derivative(X, theta, y):
 4     derivative = X.T.dot(X.dot(theta) - y) / X.shape[0]
 5     return derivative
 6 
 7 
 8 def gradient_descent(X, y, alpha=0.1):
 9     theta = np.ones(shape=y.shape, dtype=float)
10     partial_derivative_of_J = partial_derivative(X, theta, y)
11     while any(abs(partial_derivative_of_J) > 1e-5):
12         theta = theta - alpha * partial_derivative_of_J
13         partial_derivative_of_J = partial_derivative(X, theta, y)
14     return theta