线性回归模型

2019-07-13 本文已影响33人 To_QT

线性模型.png

1.表达形式：

一个示例有 $d$ 个属性： $\vec{x}=(x_1, x_2,...,x_d)$ ，其表达形式为：
$f(x)=w^Tx+b \tag{1.1}$
其中 $w=(w_1, w_2, ..., w_d)$

2.标准线性回归

令预测值与真实值的均方根误差最小化
$\begin{align*} (w^*,b^*)=\mathop{\arg\min}_{(w,b)} \sum_{i=1}^{m}(f(x_i)-y_i)^2 \tag{2.1} \end{align*}$

2.1 线性回归（二元一次函数，此时 $w$ 是一个数）

2.1.1均方误差最小化

$\begin{align*} (w^*,b^*) =\mathop{\arg\min}_{(w,b)} \sum_{i=1}^{m}(w·x_i + b-y_i)^2 \tag{2.2} \end{align*}$

2.1.2最小二乘参数估计

令： $E_{(w,b)}=\sum_{i=1}^{m}(w·x_i+b-y_i)^2$ ，对 $w,b$ 求导，令导数为0，解得 $w,b$ 分别为：
$w=\frac {\sum_{i=1}^{m}y_i(x_i-\bar{x})}{\sum _{i=1}^{m} {x^2_i}- \frac {1}{m}(\sum_{i=1}^{m}x_i)^2} \tag{2.3}$

$b=\frac {1}{m}\sum_{i=1}^{m}(y_i-wx_i) \tag{2.4}$

2.2 多元线性回归（此时 $w$ 是一个向量）

令设训练样本共有 $m$ 个，每个样本有 $d$ 个属性，把 $\boldsymbol{b}$ 与 $\boldsymbol{x}$ 写到一起，记为 $\boldsymbol{X}$ ，则 $\boldsymbol{X}$ 为：

根据向量求模长的公式，可写为

令 $E_{\boldsymbol{\hat{w}}}=(\boldsymbol{y}-\boldsymbol{X}·\boldsymbol{w^T})^T ·(\boldsymbol{y}-\boldsymbol{X}·\boldsymbol{w^T})$

考虑矩阵求导：

矩阵乘法

所以，对 $E_{\boldsymbol{\hat{w}}}$ 求导，有

令 $\frac {\partial{E_{\boldsymbol{\hat{w}}}}}{\partial{\hat{w}}}=0$ ，考虑下列两种情况：

当 $\boldsymbol{X}$ 为满秩矩阵时：
$\boldsymbol{\hat{w}^*}=\boldsymbol{(X^TX)^{-1}X^Ty}$ ，令 $\boldsymbol{\hat{x}}=[\boldsymbol{x}, 1]$ ，则学习到的模型为：
$\begin{align} f(\boldsymbol{\hat{x}_i})=\boldsymbol{\hat{x}_i}\boldsymbol{(X^TX)^{-1}}\boldsymbol{X^Ty} \tag{2.6} \end{align}$
当 $\boldsymbol{X}$ 为非满秩矩阵时（数据条数少于未知数个数）：
在这样的情况下，可能有多个解（ $\boldsymbol{\hat{w}}$ ）此时需要考虑正则化

3 局部加权线性回归

由于线性回归求的是具有最小均方误差的无偏估计，因此通常伴随着欠拟合的现象，一些方法会在估计方法中加入偏差，降低预测的均方误差。

3.1 基本思想

通过给待预测点附近的每个点赋予一定的权重，然后在这个子集上进行普通的回归（设计代价函数时，待预测点附近的点拥有更高的权重，权重随着距离的增大而缩减——这也就是名字中“局部”和“加权”的由来。）
那么，原本均方误差最小化的表达式为（公式2.2）：
$\begin{align*} (w^*,b^*) =\mathop{\arg\min}_{(w,b)} \sum_{i=1}^{m}(\sum_{j=1}^{p}w_j·x_{ij} + b-y_i)^2 \tag{2.2} \end{align*}$
现在，我们把每一个点都赋予一定的权重（用 $\theta_i$ 表示第 $i$ 个点的权重），目标是：离真实值越近，赋予的权重越大。因此，修改公式为：
$\begin{align*} (w^*,b^*) =\mathop{\arg\min}_{(w,b)} \sum_{i=1}^{m} \theta_i·(\sum_{j=1}^{p}w_j·x_{ij} + b-y_i)^2 \tag{3.1} \end{align*}$
公式(3.1)对权重 $w$ 求导并令其为0：
$\begin{align*} \frac {\partial (w^*,b^*)}{\partial w^*} &= -2 \boldsymbol{X^T \theta (y-XW) }=0 \\ \rightarrow \boldsymbol{X^TWy} &= \boldsymbol{X^T \theta XW} \\ \rightarrow \boldsymbol{W} &= \boldsymbol{(X^T \theta X)^{-1} X^T \theta y} \tag{3.2} \end{align*}$

3.2 权重怎么取？

首先明确一点：局部加权线性回归是一个非参数（non-parametric）算法。之前学习的（不带权）线性回归算法是有参数（parametric）算法，因为它有固定的有限数量的，能够很好拟合数据的参数（权重）。一旦我们拟合出权重并存储了下来，也就不需要再保留训练数据样本来进行更进一步的预测了。相比而言，用局部加权线性回归做预测，我们需要保留整个的训练数据，每次预测得到不同的权重，即参数不是固定的。

术语 “非参数” 粗略意味着：我们需要保留用来代表假设 h的内容，随着训练集的规模变化是呈线性增长的。

3.3 机器学习-核函数（核模型）

局部加权线性回归使用“核”函数对附近的点赋予更高的权重，目前常用的类型就是高斯核，具体表达式为：
$\begin{align*} w(i,i) &= exp(\frac{|x^{(i)}-x|^2 } {-2k^2} ) \tag{3.3} \end{align*}$

使用高斯核函数具有以下特征：

构建了一个只含有对角元素的权重矩阵 $w$ ，当点 $x$ 与 $x^{(i)}$ 越接近，则 $w(i,i)$ 将会越大，随着样本点与待预测点距离的递增，权重将会以指数级衰减。
$k$ 能够控制衰减的速度。若 $k$ 越大，则权重的宽度越大，衰减的速度越慢，容易导致欠拟合；若 $k$ 越小，则权重的宽度越窄，衰减的速度越快，容易造成过拟合；<当 $k=1$ 时，可认为是最小二乘法拟合的结果>

通常情况下，通过交叉验证或者网格搜索等方法确定最佳的k值。

权重宽度与K值的关系

为什么高斯核是对角矩阵？
因为在矩阵中，只有对角线上的元素值表示的是 $x_1^2,x_2^2,...,x_n^2$ 的系数。

上述两个方法：标准线性回归和局部加权线性回归都有一个重要基础——特征要比样本少，即输入数据的矩阵必须要是满秩矩阵。那么，当矩阵为非满秩矩阵时，都哪些处理方法？

4 岭回归

在标准线性回归中，权重的值为：
$\boldsymbol {\hat{w}^*} = \boldsymbol{(X^TX)^{-1}X^Ty}$

而岭回归的方法是在矩阵 $X^TX$ 上加一个 $\lambda I$ ，即从对 $X^TX$ 的逆转化为对 $X^TX+\lambda I$ 求逆，公式可写为：

$\begin{align*} \boldsymbol{\hat{w}^*}=\boldsymbol{(X^TX+\lambda I)^{-1}X^Ty} \tag{4.1} \end{align*}$
对比公式(2.2)可知，岭回归其实是优化下面这个问题：
$\begin{align*} (w^*,b^*) =\mathop{\arg\min}_{(w,b)} \sum_{i=1}^{m}(\sum_{j=1}^{p}w_j·x_{ij} + b-y_i)^2 + \lambda \sum_{j=1}^{p}w_j^2 \tag{4.2} \end{align*}$

岭回归的意义:
加上了L2范数（目标函数的惩罚函数），作用是确保权重值不会很大，起到收缩的作用。对于岭回归来说，随着 $\lambda$ 的增大，模型的方差会减少（ $X^TX+\lambda I$ 增大，( $X^TX+\lambda I)^{-1}$ 减小， $w$ 减小，因此偏差增大，因此 $\lambda$ 能够平衡偏差与方差的关系。），·

4.1 岭回归的几何意义

$\begin{align*} \left \{ \begin{matrix} (w^*,b^*) =\mathop{\arg\min}_{(w,b)} \sum_{i=1}^{m}(\sum_{j=1}^{p}w_j·x_{ij} + b-y_i)^2 \\ 附加约束条件： \lambda \sum_{j=1}^{p}w_j^2 \leq t & \end{matrix} \right. \tag{4.3} \end{align*}$

4.1.1 岭回归添加回归系数平方和的原因

为了解决多重共线性的麻烦。作者在《岭回归和LASSO回归的区别》中提到了一个很通俗的例子，一个家庭的收入和支出具有很强的共线性，但是在岭回归系数平方和的约束下，通过调整权重值，即使收入有很大的正数，支出为很小的负数，最后预测的结果也不会有较大的偏差，这就是岭回归系数平方和约束的作用。

4.1.2 岭回归系数的性质

岭回归系数是OLS估计的线性变换
岭回归系数是有偏的
当lambda>0时，岭回归系数具有压缩性

4.1.3 岭参数的选择

我们知道岭回归系数会随着 $\lambda$ 的变化而变化，为保证选择出最佳的岭回归系数，该如何确定这个 $\lambda$ 值呢？一般我们会选择定性的可视化方法和定量的统计方法。对这种方法作如下说明：
1）绘制不同 $\lambda$ 值与对应的 $\beta$ 值之间的折线图，寻找那个使岭回归系数趋于稳定的 $\lambda$ 值；同时与OLS相比，得到的回归系数更符合实际意义；
2）方差膨胀因子法，通过选择最佳的 $\lambda$ 值，使得所有方差膨胀因子不超过10；
3）虽然 $\lambda$ 的增大，会导致残差平方和的增加，需要选择一个 $\lambda$ 值，使得残差平方和趋于稳定（即增加幅度细微）。

5 lasso 回归

上文提到，岭回归增加的约束是 $\sum_{k=1}^{p}\beta_k^2 \leq t$ ，那么在lasso回归中，添加的约束条件为： $\sum_{k=1}^{p}|\beta_k| \leq t$

$\begin{align*} (w^*,b^*) =\mathop{\arg\min}_{(w,b)} \sum_{i=1}^{m}(\sum_{j=1}^{p}w_j·x_{ij} + b-y_i)^2 + \lambda \sum_{j=1}^{p} |w_j| \tag{5.1} \end{align*}$

若采用这种约束，当 $\lambda$ 足够小的时候，一些系数会因此衰减为0。

但是，在新的约束条件下，若要解出回归系数，则需要使用二次规划法算法

6 前向逐步回归

前向逐步回归算法利用了贪心算法的思想，具体算法伪代码如下：

前向逐步回归-machine learning in action

7 代码实现

完整代码在这里

class Regression(object):
    def __init__(self):
        self.ws = None

    def standRgressFit(self, x_data, y_data):
        '''
        标准线性回归
        :param data: 待拟合的数据
        :param x_name: 输入的特征名称
        :param y_name: 拟合对象的特征名称
        :return: 系数
        '''
        x_matrix = np.mat(x_data.values)
        y_matrix = np.mat(y_data.values)

        xTx = x_matrix.T * x_matrix
        # 如果是非满秩矩阵
        if np.linalg.det(xTx) == 0:
            print("This matrix is singular, cannot do inverse")
            return
        else:
            self.ws = xTx.I * (x_matrix.T * y_matrix)

    def lwlrFit(self, test_point, x_data, y_data, k=1):
        '''
        局部加权线性回归
        :param test_point: 待预测点的特征
        :param x_data: 训练集的特征数据
        :param y_data: 训练集的结果数据
        :param k: 权重的宽度
        :return: 加权线性回归的权重矩阵
        '''
        xMat = np.mat(x_data)
        yMat = np.mat(y_data).T
        m = np.shape(xMat)[0]
        # 创建一个m*m的单位矩阵
        weights = np.mat(np.eye((m)))
        # 对角矩阵赋予高斯核函数
        for j in range(m):  # next 2 lines create weights matrix
            diffMat = test_point - xMat[j, :]  #
            weights[j, j] = np.exp(diffMat * diffMat.T / (-2.0 * k ** 2))
        xTx = xMat.T * (weights * xMat)
        if np.linalg.det(xTx) == 0.0:
            print("This matrix is singular, cannot do inverse")
            return
        self.ws = xTx.I * (xMat.T * (weights * yMat))

    def ridgeFit(self, x_data, y_data, lam=0.2):
        '''
        岭回归拟合
        :param x_data: 训练数据的特征值
        :param y_data: 训练数据的y值
        :param lam: lambda的取值
        :return: 线性回归的参数
        '''
        x_matrix = np.matrix(x_data)
        y_matrix = np.matrix(y_data)

        denom = x_matrix.T * x_matrix + lam * np.eye(x_matrix.shape[1])
        if np.linalg.det(denom) == 0:
            print("This matrix is singular, cannot do inverse")
            return
        else:
            self.ws = denom.I * (x_matrix.T * y_matrix)

    def stageWiseFit(self, x_data, y_data, epsilon=0.01, iterations=200):
        '''
        计算前向逐步回归的参数
        :param x_data: 训练数据集
        :param y_data: 训练数据的输出值
        :param epsilon: 步长
        :param iterations: 迭代次数
        :return: 返回每一次迭代过程中的权重大小
        '''
        x_matrix = np.matrix(x_data)
        y_matrix = np.matrix(y_data)

        weight_num = x_data.shape[1]
        # 初始化权重
        self.ws = np.zeros((weight_num, 1))
        bset_w = self.ws.copy()
        # 初始化权重的记录矩阵
        # （用returnMat矩阵记录iterations次迭代中的权重
        returnMat = np.zeros((iterations, weight_num))  # testing code remove

        for i in range(iterations):
            lowest_error = np.inf
            # 修改第j个特征的权重
            for j in range(weight_num):
                # sign 决定方向，往哪边走
                for sign in [-1, 1]:
                    current_w = self.ws.copy()
                    current_w[j] += epsilon * sign
                    y_hat = x_matrix * current_w

                    # 修改第j个特征后，误差是否有所降低？
                    current_error = rssError(y_matrix.A, y_hat.A)

                    if lowest_error > current_error:
                        lowest_error = current_error
                        bset_w = current_w
            self.ws = bset_w
            returnMat[i, :] = self.ws.T
        return returnMat

    def predict(self, x_data):
        '''
        预测数据
        :param x_data: 预测的特征值
        :return: y模拟值
        '''
        x_matrix = np.mat(x_data.values)
        return x_matrix * self.ws

def rssError(yArr, yHatArr): #yArr and yHatArr both need to be arrays
    '''
    计算mse
    :param yArr: y的真实值
    :param yHatArr: y的模拟值
    :return: mse
    '''
    return ((yArr-yHatArr)**2).sum()

def corrCoef(y_hat, y):
    '''
    计算相关系数
    :param y_hat: 模拟值
    :param y: 真实值
    :return: 相关系数矩阵
    '''
    return np.corrcoef(y_hat.T, y.T)

def regularize(data):
    '''
    标准化数据
    :param data: 待标准化的dataframe
    :return: 标准化后的dataframe
    '''
    data = (data - data.mean()) / data.var()
    return data

参考文献

岭回归和LASSO回归的区别
统计学习方法
西瓜书
机器学习实战

线性回归模型

1.表达形式：