回归算法之线性回归

算法,机器学习 2017-12-21

起步

线性回归是利用数理统计和归回分析,来确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法。与之前的分类问题( Classification )不一样的是,分类问题的结果是离散型的;而回归问题中的结果是数值型的。

描述数据特征

数理统计中,常用的描述数据特征的有:

均值( mean ):

又称平均数或平均值,是计算样本中算术平均数:

\bar{x} = \frac{\sum_{i=1}^{n} x_i}{n}

中位数( median )

将数据中的各个数值按照大小顺序排列,居于中间位置的变量。

众数( mode )

数据中出现次数最多的数。

方差( variance )

这是一种描述离散程度的衡量方式:

s^2 = \frac{\sum_{i=1}^{n} (x_i - \bar{x})^2}{n}

标准差 (standard deviation) 将方差开方就能得到标准差。

简单线性回归

简单线性回归( Simple Linear Regression )包含一个自变量 x 和一个因变量 y ,也就是一元线性回归。被用来描述因变量(y)和自变量(X)以及偏差(error)之间关系的方程叫做回归模型,这个模型是:

y = \beta_0 + \beta_1x + \varepsilon

其中偏差 ε 满足正态分布的。因此它的期望值是 0 。ε的方差(variance)对于所有的自变量 x 是一样的。

因此式子也可以写成:

E(y) = \beta_0 + \beta_1x

其中,β0 是回归线的截距,β1 是回归线的斜率,E(y) 是在一个给定 x 值下 y 的期望值(均值)。

假设我们的估值函数(注意,是估计函数):

\hat{y} = b_0 + b_1x

β 和 b 有什么关系呢? β0 是表示回归方程中真实的值,而 b0 是对 β0 的估值,通过训练能让 b0 趋近于 β0 。从而提高估值函数准确度。

如何练处适合简单线性回归模型的最佳回归线?

Image.png

寻找的回归方程,是与真实值的离散程度最小的:

\min \sum_{i=0}^{n}(y_i - \hat{y}_i)^2

这就类似于高中学的最小二乘法了,利用那时候学到的公式套上去就能得到该方程了:

\begin{align*}
b_1 &= \frac{\sum(x_i - \bar{x})(y_i - \bar{y})}{\sum(x_i - \bar{x})^2} \\
b_0 &= \bar{y} - b_1\bar{x}
\end{align*}

最小二乘法的推导过程可以查看维基,写得非常清楚:https://zh.wikipedia.org/wiki/最小二乘法

对应的python代码:

class SimpleLinearRegression(object):
    """
    简单线性回归方程,即一元线性回归,只有一个自变量,估值函数为: y = b0 + b1 * x
    """
    def __init__(self):
        self.b0 = 0
        self.b1 = 0
    def fit(self, x: list, y: list):
        n = len(x)
        x_mean = sum(x) / n
        y_mean = sum(y) / n
        dinominator = 0
        numerator = 0
        for xi, yi in zip(x, y):
            numerator += (xi - x_mean) * (yi - y_mean)
            dinominator += (xi - x_mean) ** 2
        self.b1 = numerator / dinominator
        self.b0 = y_mean - self.b1 * x_mean

    def pridict(self, x):
        return self.b0 + self.b1 * x

多元线性回归

当自变量有多个时,回归模型就变成了:

E(y) = \beta_0 + \beta_1x_1 + ... + \beta_nx_n

估值函数的自变量也变多:

\hat{y} = b_0 + b_1x_1 + ... + b_nx_n

令矩阵:

X=\begin{bmatrix}
1 \\
x_i \\
... \\
x_n
\end{bmatrix}
B=\begin{bmatrix}
b_0 \\
b_1 \\
... \\
b_n
\end{bmatrix}

函数可以写成:

h(x) = B^TX

有训练数据:

D = {(x^{(1)},y^{(1)}), (x^{(2)},y^{(2)}), ... ,(x^{(m)},y^{(m)})}

损失函数( cost function ),寻找的目标函数应该让损失函数尽可能小,这个损失函数为:

\begin{align*}
J(B) &= \frac1{2m} \sum_{i=1}^m (h(x^{(i)}) - y^{(i)})^2 \\
&= \frac1{2m}(XB - y)^T(XB-Y)
\end{align*}

目的就是求解出一个使得代价函数最小的 B。

梯度下降求解估值函数

梯度下降算法是一种求局部最优解的方法,这部分算法可以查看https://www.cnblogs.com/pinard/p/5970503.html 了解下,这里不展开。B 由多个元素组成,对于损失函数可以求偏导如下:

\frac{\partial}{\partial B_j} J(B) = \frac1m\sum_{i=1}^m (h(x^{(i)}) - y^{(i)})x_j^{(i)}

更新B:

B_j : B_j - \alpha \frac{\partial}{\partial B_j} J(B)

其中 α 表示学习率,决定了在梯度下降迭代的过程中,每一步沿梯度负方向前进的长度。

最小二乘法求解估值函数

多元多项式的情况一样可以利用最小二乘法来进行求解,将 m 条数据都带入估值函数可以得到线性方程组:

\begin{align*}
&b_0 + b_1x_{11} + b_2x_{12} + ... + b_nx_{1n} = y_1 \\
&b_0 + b_1x_{21} + b_2x_{22} + ... + b_nx_{2n} = y_2 \\
&... \\
&b_0 + b_1x_{m1} + b_2x_{m2} + ... + b_nx_{mn} = y_m 
\end{align*}

向量化表示:

X=\begin{bmatrix}
1 & x_{11} & x_{12} & ... & x_{1n} \\
1 & x_{21} & x_{22} & ... & x_{2n} \\
... \\
1 & x_{m1} & x_{m2} & ... & x_{mn} \\
\end{bmatrix}

B=\begin{bmatrix}
b_0 \\
b_1 \\
... \\
b_n \\
\end{bmatrix}

Y=\begin{bmatrix}
y_0 \\
y_1 \\
... \\
y_n \\
\end{bmatrix}

因为要使:

\sum_{i=1}^m (y_i - \hat{y})^2

的值最小,则系数矩阵中,各个参数的偏导的值都会0,因此可以得到 n 个等式:

\left\{
\begin{array}{lr}
\frac{\partial}{\partial b_0} J = -2 \sum_{i=1}^m (y_i - b_0 - b_1x_{i1} - ... - b_nx_{in}) = 0 \\
\frac{\partial}{\partial b_1} J = -2 \sum_{i=1}^m (y_i - b_0 - b_1x_{i1} - ... - b_nx_{in})x_{i1} = 0 \\
... \\
\frac{\partial}{\partial b_n} J = -2 \sum_{i=1}^m (y_i - b_0 - b_1x_{i1} - ... - b_nx_{in})x_{in} = 0
\end{array}
\right.

得到正规方程组:

\begin{align*}
mb_0 + \sum_{i=1}^mx_{i1}b_1 + ... + \sum_{i=1}^mx_{in}b_n &= \sum_{i=1}^my_i \\
\sum_{i=0}^mx_{i1}b_0 + \sum_{i=1}^mx_{i1}^2b_1 + ... + \sum_{i=1}^mx_{in}b_n &= \sum_{i=1}^mx_{i1}y_i \\
... \\
\sum_{i=0}^mx_{in}b_0 + \sum_{i=1}^mx_{i1}b_1 + ... + \sum_{i=1}^mx_{in}^2b_n &= \sum_{i=1}^mx_{in}y_i \\
\end{align*}

写成矩阵表示为:

X^TXB = X^TY

解得:

B = (X^TX)^{-1}X^TY

有的解法会把矩阵 X 的转置写成 X' ,我比较喜欢用上标 T 表示。

对应的python代码:

import numpy as np
class MyLinearRegression(object):
    def __init__(self):
        self.b = []

    def fit(self, x: list, y: list):
        # 为每条训练数据前都添加 1
        tmpx = [[1] for _ in range(len(x))]
        for i, v in enumerate(x):
            tmpx[i] += v

        x_mat = np.mat(tmpx)
        y_mat = np.mat(y).T
        xT = x_mat.T
        self.b = (xT * x_mat).I * xT * y_mat

    def predict(self, x):
        return np.mat([1] + x) * self.b

测试

将我们书写的代码与 sklearn 的进行对比看看:

# coding: utf-8
from sklearn import linear_model

x = [
    [100.0, 4.0],
    [50.0, 3.0],
    [100.0, 4.0],
    [100.0, 2.0],
    [50.0, 2.0],
    [80.0, 2.0],
    [75.0, 3.0],
    [65.0, 4.0],
    [90.0, 3.0],
    [90.0, 2.0]
]

y = [9.3, 4.8, 8.9, 6.5, 4.2, 6.2, 7.4, 6.0, 7.6, 6.1]

test_row = [50, 3]
linear = MyLinearRegression()
linear.fit(x, y)
print(linear.predict(test_row)) # [[ 4.95830457]]

sk = linear_model.LinearRegression()
sk.fit(x, y)
print(sk.predict([test_row])) # [ 4.95830457]

预测结果一致,开心。

总结

在线性回归方程的求解中,因为梯度下降求得是局部最优解,且系数需要随机初始化,因此一般解法是采用最小二乘法。最小二乘法中,如果两个训练数据一致,那么训练出来的模型也会一模一样。


本文由 hongweipeng 创作,采用 知识共享署名 3.0,可自由转载、引用,但需署名作者且注明文章出处。

赏个馒头吧