预测数值型数据：回归

项目代码

回归

回归一词是由达尔文的表兄弟Francis Galton发明的，Galton根据上一代的豌豆种子（双亲）的尺寸来预测下一代豌豆种子（孩子）的尺寸。Galton在大量对象上应用了回归分析，甚至包含人的身高。Galton注意到，如果双亲的高度比平均高度高，他们子女也倾向于比平均高度高，但尚不及双亲。孩子的高度向着平均高度回退(回归)。Galton在多项研究上都注意到这个现象，所以尽管这个英文单词跟数值预测没有任何关系，但这种研究方法仍被称作回归。
本节中的回归的目的是预测数值型的目标值。
假设方程 $y=w_1{x}_1+w_2{x}_2+...+w_n{x}_n$，我们可以称该方程为回归方程，其中$w$称作回归系数，求这些回归系数的过程就是回归。

线性回归

假设输入数据为矩阵$X$，其对应值为$y$，回归系数向量为$w$。对于给定数据$x_i$，其预测值可以表示为 ${\hat{y}}_i=x_i{}^{T}w$ 。
我们的目的是求出向量$w$，这里用的方法是求预测值 $\hat{y}$ 和真实值 $y$ 之间的差值，采用平均误差，并计算使误差总和最小的向量$w$。
平方误差可以表示为：${\sum }_{i=1}^m(y_i-{\hat{y}}_i{)}^2$ $=$ ${\sum }_{i=1}^m(y_i-x_i{}^{T}w{)}^2$
上式用矩阵形式可以表示为：$(y-Xw{)}^T(y-Xw)$，对$w$求导，得到$X^T(y-Xw)$，另其等于0，则$\hat{w}=(X^{T}X{)}^{-1}{X}^{T}y$，$\hat{w}$ 表示当前估计出的 $w$ 的最优解。
上式中存在 $(X^{T}X{)}^{-1}$ ，因此需要对矩阵求逆，因此该方程只在逆矩阵存在的情况下适用。
计算下图数据的回归系数：

加载数据及求解代码：

import numpy as np
import matplotlib.pyplot as plt
import numpy.linalg as nlinalg

def read_data(path):
    data = np.loadtxt(path)
    return data[:, :2], data[:, -1]

def stand_reg(x,y):
    x_mat = np.matrix(x)
    y_mat = np.matrix(y).T
    xTx = x_mat.T*x_mat
    if nlinalg.det(xTx) == 0:
        print('This matrix is singular, cannot do inverse')
        return
    ws = xTx.I*(x_mat.T*y_mat)
    return ws
 
 if __name__ == '__main__':

    data, la = read_data('./data/ex0.txt')
    ws = stand_reg(data,la)

    plt.scatter(data[:,1], la)
    y_hat = np.matrix(data)*ws
    print('corrcoef:', np.corrcoef(y_hat.T,np.matrix(la))) # 计算预测值与实际值之间的相关系数，确定结果好坏
    plt.plot(data[:,1], y_hat, color='r')
    plt.show()

结果如图：

局部加权线性回归

线性回归有可能出现欠拟合现象，因为它求的是具有最小均方误差的无偏估计。因此有些方法允许在估计中引入一些偏差，从而降低预测的均方误差。
局部加权线性回归（Locally Weighted linear Regression，LWLR）就是其中一个方法，其给待预测点附近的每个点赋予一定的权重。
该算法解出的回归系数$w$形式如下：$\hat{w}=(X^{T}WX{)}^{-1}{X}^{T}Wy$，其中 $W$ 是一个矩阵，用来给每个数据点赋予权重。
LWLR使用‘核’来对附近的点赋予更高的权重。核的类型可以自由选择，做常用的就是高斯核：$w(i,i)=exp(\frac{|x^{(i)}-x|}{-2k^2})$，其中$k$需要用户指定。
代码实现如下：

import numpy as np
import numpy.linalg as nlinalg


def lwlr_(test_data, x, y, k=1.0):
    x_mat = np.matrix(x)
    y_mat = np.matrix(y).T
    n = x.shape[0]
    weights = np.matrix(np.eye((n)))
    for j in range(n):
        diff_mat = test_data - x_mat[j, :]
        weights[j, j] = np.exp(diff_mat * diff_mat.T / (-2 * k ** 2))
    xTx = x_mat.T * (weights * x_mat)
    if nlinalg.det(xTx) == 0:
        print('This matrix is singular, cannot do inverse')
        return
    ws = xTx.I * (x_mat.T * (weights * y_mat))
    return test_data*ws


def test_lwlr(testArr, x, y, k=1.):
    m = testArr.shape[0]
    y_hat = np.zeros(m)
    for i in range(m):
        y_hat[i] = lwlr_(testArr[i], x, y, k)
    return y_hat
if __name__ == '__main__':
    data, la = read_data('./data/ex0.txt')
    # ws = stand_reg(data,la)
    #
    # plt.scatter(data[:,1], la)
    # y_hat = np.matrix(data)*ws
    # print('corrcoef:', np.corrcoef(y_hat.T,np.matrix(la))) # 计算预测值与实际值之间的相关系数，确定结果好坏
    # plt.plot(data[:,1], y_hat, color='r')
    # plt.show()

    # lwlr
    y_hat = test_lwlr(data, data, la)
    plt.scatter(data[:, 1], la)
    x_mat = np.matrix(data)
    Ind = x_mat[:,1].argsort(0)
    data_sort = x_mat[Ind][:,0,:]
    plt.plot(data_sort[:, 1], y_hat[Ind], color='r')
    plt.show()

实验结果如图，k值依次为0.003，0.01，1：

岭回归

数据集中样本特征比样本点还多的，即输入矩阵$X$不是满秩矩阵，非满秩矩阵求逆时会出现问题。
岭回归（ridge regression）就是在矩阵 $X^{T}X$ 上加一个$\lambda I$ 从而使矩阵非奇异，进而能够对$X^{T}X$+ $\lambda I$ 求逆。其中矩阵 $I$是一个$m\ast m$ 的单位矩阵。$\lambda$为用户自定义参数。
回归系数的形式变为：$\hat{w}=(X^{T}X+\lambda I{)}^{-1}{X}^{T}y$
岭回归最先用来处理特征数多于样本数的情况，现在也用于在估计中加入偏差，从而得到更好的估计，引入$\lambda$来限制所有$w$ 之和，通过引入该惩罚项，能够减少不重要的参数，这在统计学中叫做缩减（shrinkage）。
代码实现：

import numpy as np
import numpy.linalg as nlinalg


def ridge_regress(x_mat, y_mat, lam=1.0):

    xTx = x_mat.T * x_mat
    denom = xTx + np.eye(x_mat.shape[1])*lam
    if nlinalg.det(denom) == 0:
        print('This matrix is singular, cannot do inverse')
        return
    ws = denom.I * (x_mat.T*y_mat)
    return ws


def test_ridge_regress(x, y):
    x_mat = np.matrix(x)
    y_mat = np.matrix(y).T
    # 数据标准化
    y_mean = np.mean(y_mat, 0)
    y_mat = y_mat-y_mean
    x_means = np.mean(y_mat,0)
    x_vars = np.var(y_mat,0)
    x_mat = (x_mat-x_means)/x_vars

    nums = 30
    w_mat = np.zeros((nums, x_mat.shape[1]))
    for i in range(nums):
        ws = ridge_regress(x_mat, y_mat, np.exp(i-10))
        w_mat[i,:] = ws.T
    return w_mat
if __name__ == '__main__':
    # data, la = read_data('./data/ex0.txt')
    # ws = stand_reg(data,la)
    #
    # plt.scatter(data[:,1], la)
    # y_hat = np.matrix(data)*ws
    # print('corrcoef:', np.corrcoef(y_hat.T,np.matrix(la))) # 计算预测值与实际值之间的相关系数，确定结果好坏
    # plt.plot(data[:,1], y_hat, color='r')
    # plt.show()

    # lwlr
    # y_hat = test_lwlr(data, data, la)
    # plt.scatter(data[:, 1], la)
    # x_mat = np.matrix(data)
    # Ind = x_mat[:,1].argsort(0)
    # data_sort = x_mat[Ind][:,0,:]
    # plt.plot(data_sort[:, 1], y_hat[Ind], color='r')
    # plt.show()

    # ridge
    data, la = read_data('./data/abalone.txt')
    ridge_weights = test_ridge_regress(data,la)
    plt.plot(ridge_weights)
    plt.show()

当$\lambda$ 不同时的实验结果图：

$\lambda$ 非常小时，系数与普通回归一样。而$\lambda$ 非常大时，所有回归系数缩减为0。

lasso

岭回归的约束条件如下：
${\sum }_{k=1}^n{w}_k^2<=\lambda$，等于添加了L2正则化。
而lasso的约束条件为：
${\sum }_{k=1}^n|w_k|<=\lambda$，等于添加了L1正则化.
lasso增加了计算复杂度，为了解出这个回归系数，需要用到二次规划算法。

前向逐步回归

（上图来源机器学习实战）
逐步线性回归的好处在于可以帮助人们理解现有模型的并作出改进。当构建一个模型后，可以运行该算法找出重要特征，这样就有可能及时停止对那些不重要特征的收集。