预测数值型数据：回归

2023-10-27

本文传送机

用线性回归找到最佳拟合直线

局部加权线性回归

用线性回归找到最佳拟合直线

线性回归

优点：结果易于理解，计算上不复杂

缺点：对非线性的数据拟合不好

适用数据类型：数值型和标称型数据

回归方程为： $\hat{Y}=Xw$ ，其中 $w$ 为回归系数

回归的一般方法

收集数据：采用任意方法收集数据

准备数据：回归需要数值型数据，标称型数据将被转换为二值型数据

分析数据：绘制出数据的可视化二维图将有助于对数据做出理解和分析，在采用压缩减法求得新回归系数之后，可以将新拟合线在图上作为对比

训练算法：找到回归系数

测试算法：使用 $R^2$ 或者预测值和数据的拟合度来分析模型的效果

使用算法：使用回归，可以在给定输入的时候预测出一个数值，这是对分类方法的提升，因为这样可以预测连续型数值而不仅仅是离散的类别标签。

现在的问题是，给定 $\bold{X}$ 和对应的 $y$ ，怎样才能找到 $w$ 呢？一个常用的方法就是找出使误差最小的 $w$ 。这里的误差是指预测值 $\hat{y}$ 和真实值之间的差值，使用丐武训哈的简单累加将使得正差值和负差值相互抵消，所以我们采用平方误差。

平方误差为：

$\varepsilon =\sum_{i=1}^m(y_i-x_i^Tw)^2$

写成矩阵运算的形式为：

$\varepsilon =(Y-\Bold{X}w)^T(Y-\Bold{X}w)=w^TX^TXw-2w^TX^TY+Y^TY$

为求 $\varepsilon$ 的极小值，可对 $\varepsilon$ 进行求导， $\varepsilon$ 的极值点一定出现在 $\frac{\partial \varepsilon }{\partial w}=0$ 的位置：

$\frac{\partial \varepsilon }{\partial w}=2X^TXw-2X^TY=0$

可得 $\hat{w}=(X^TX)^{-1}X^TY$

这里的 $\hat{w}$ 便是回归系数的最优估计向量。

代码：

import numpy as np
import matplotlib.pyplot as plt


def load_dataset(filename):
    feature_number = len(open(filename).readline().split('\t'))-1
    data_mat = []
    label_mat = []
    file = open(filename)
    for line in file.readlines():
        line_array = []
        curr_line = line.strip().split('\t')
        for i in range(feature_number):
            line_array.append(float(curr_line[i]))
        data_mat.append(line_array)
        label_mat.append(curr_line[-1])
    return data_mat, label_mat


def linear_regress(x_array, y_array):
    x_mat = np.mat(x_array, dtype=np.float64)
    y_mat = np.mat(y_array, dtype=np.float64).T
    xTx = x_mat.T * x_mat
    if np.linalg.det(xTx) == 0.0:
        print("This matrix is singular, cannot do inverse")
        return
    w_hat = xTx.I * (x_mat.T * y_mat)
    return w_hat


if __name__ == '__main__':
    x_array, y_array = load_dataset('ex0.txt')
    w_hat = linear_regress(x_array, y_array)
    x_mat = np.mat(x_array, dtype=np.float64)
    y_mat = np.mat(y_array, dtype=np.float64)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(x_mat[:, 1].flatten().A[0], y_mat.T[:, 0].flatten().A[0])
    x_copy = x_mat.copy()
    x_copy.sort(0)
    y_hat = x_copy * w_hat
    ax.plot(x_copy[:, 1], y_hat, c="red")
    plt.show()

结果：

那么，如何判断这些模型的好坏呢？

有种方法可以计算预测序列 $\hat{Y}$ 和真实值 $Y$ 之间的匹配程度，那就是计算这两个序列的相关系数（numpy中提供了corrcoef方法）：

print(np.corrcoef(y_hat.T, y_mat))

结果：

[[1.         0.13653777]
 [0.13653777 1.        ]]

结果的含义是： $\hat{Y}$ 和自己的匹配最完美，而和 $Y$ 的相关系数为0.1365

局部加权线性回归

线性回归的一个问题是有可能出现欠拟合现象，因为它求得是具有最小均方误差得无偏估计。显而易见，如果模型欠拟合将不能取得最好得预测结果。所以有些方法允许在估计中引入一些偏差，从而降低预测得均方误差。

其中的一个方法是局部加权线性回归（Locally Weighted Linear Regression, LWLR），在该算法中，我们给待测点附近的每一个点赋予一定的权重，然后在这个子集上基于最小均方差来进行普通的回归。与kNN一样，这种算法每次预测均需要事先选取出对应的数据子集。该算法决出的回归系数 $\hat{w}$ 的形式如下：

$\hat{w}=(X^TWX)^{-1}X^TWY$

其中 $W$ 是一个矩阵，用来给每个数据点赋予权重。

LWLR（局部加权线性回归）使用“核”（与支持向量机中的核类似）来对附近点赋予更高的权重。核的类型可以自由选择，最常用的核就是高斯核，高斯核对应的权重如下：

$w(i,i)=\exp \left(\frac{\left | x^{(i)}-x \right | ^ 2} {-2k^2}\right)$

这样就构建了一个只含有对角元素的权重矩阵 $W$ ，并且点 $x$ 与 $x^{(i)}$ 越近， $w(i,i)$ 将会越大。其中 $k$ 是一个需要用户指定的参数，它决定了对附近的点赋予多大的权重，这也是LWLR唯一需要考虑的参数。

代码：

import numpy as np
import matplotlib.pyplot as plt


def load_dataset(filename):
    feature_number = len(open(filename).readline().split('\t'))-1
    data_mat = []
    label_mat = []
    file = open(filename)
    for line in file.readlines():
        line_array = []
        curr_line = line.strip().split('\t')
        for i in range(feature_number):
            line_array.append(float(curr_line[i]))
        data_mat.append(line_array)
        label_mat.append(curr_line[-1])
    return data_mat, label_mat


def lwlr(test_point, x_array, y_array, k=1.0):
    x_mat = np.mat(x_array, dtype=np.float64)
    y_mat = np.mat(y_array, dtype=np.float64).T
    m = np.shape(x_mat)[0]
    weights = np.mat(np.eye(m))
    for j in range(m):
        diff_mat = test_point - x_mat[j, :]
        weights[j, j] = np.exp(diff_mat * diff_mat.T / (-2.0 * k ** 2))
    xTx = x_mat.T * (weights * x_mat)
    if np.linalg.det(xTx) == 0.0:
        print("This matrix is singular, cannot do inverse")
        return
    w_hat = xTx.I * (x_mat.T * weights * y_mat)
    return test_point * w_hat


def lwlr_test(test_array, x_array, y_array, k=1.0):
    m = np.shape(test_array)[0]
    y_hat = np.zeros(m)
    for i in range(m):
        y_hat[i] = lwlr(test_array[i], x_array, y_array, k)
    return y_hat


if __name__ == '__main__':
    x_array, y_array = load_dataset('ex0.txt')
    x_array = np.array(x_array, dtype=np.float64)
    y_array = np.array(y_array, dtype=np.float64)
    y_hat = lwlr_test(x_array, x_array, y_array, 0.003)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(x_array[:, 1], y_array)
    sorted_index = x_array[:, 1].argsort(0)
    sorted_x_array = x_array[sorted_index]
    sorted_y_hat = y_hat[sorted_index]
    ax.plot(sorted_x_array[:, 1], sorted_y_hat, c="red")
    plt.show()

结果：

通过缩减系数来“理解”数据

如果特征比样本点还多（n>m），也就是说输入数据的矩阵 $X$ 不是满秩矩阵。非满秩矩阵在求逆时会出现问题。

为了解决这个问题，统计学家引入了岭回归（ridge regression）的概念，同时还有lasso法（这种方法效果好但计算复杂）；第二种缩减的方法为前向逐步法，可以得到和lasso差不多的效果，且更容易实现。

岭回归

简单来说，岭回归就是在矩阵 $X^TX$ 上加一个 $\lambda I$ 从而使得矩阵变得非奇异，这样便能对 $X^TX+\lambda I$ 求逆，其中 $\lambda$ 是一个用户自定义参数。在这种情况下，回归系数的计算公式为：

$\hat{w}=(X^TX+\lambda I)^{-1}X^TY$

岭回归最先用来处理特征数多余样本数的情况，现在也用于在估计中加入偏差，从而得到更好的估计。这里通过引入 $\lambda$ 来限制所有的 $w$ 之和，通过引入该惩罚项，能够减少不重要的参数，这个技术在统计学中称为缩减（shrinkage）

缩减方法可以去掉不重要的参数，因此能更好得理解数据。此外，与简单的线性回归相比，缩减法能够取得更好的预测效果。

def ridge_regress(x_mat, y_mat, lam=0.2):
    xTx = x_mat.T * x_mat
    denom = xTx + np.eye(np.shape(x_mat)[1]) * lam
    if np.linalg.det(denom) == 0.0:
        print("This matrix is singular, cannot do inverse")
        return
    w_hat = denom.T * (x_mat.T * y_mat)
    return w_hat

lasso

不难证明，在增加如下约束时，普通的最小二乘法回归会得到与岭回归的一样的公式：

$\sum_{k=1}^n w_k^2 \leqslant \lambda$

上式限定了所有回归系数的平方和不能大于 $\lambda$ 。使用普通的最小二乘法回归在当两个或更多的特征相关时，可能会得出一个很大的正系数和很大的负系数。正是因为上述限制条件的存在，使用岭回归可以避免这个问题。

与岭回归类似，另一个缩减方法lasso也对回归系数做出了限定，对应的约束条件如下：

$\sum_{k=1}^n\left | w_k \right |\leqslant \lambda$

唯一的不同点在于，这个约束条件使用绝对值取代了平方和。虽然约束形式只是稍作变化，结果却大项径庭：在 $\lambda$ 足够小的时候，一些系数会因此被迫缩减到0，这个特性可以帮助我们更好地理解数据，这两个约束条件在公式上看起来相差无几，但细微得变化却极大地增加了计算复杂度（为了在这个新得约束条件下解出回归系数，需要使用二次规划算法）。下面将介绍一个更为简单的方法来得到结果，该方法叫做前向逐步回归。

前向逐步回归

前向逐步回归算法可以得到与lasso差不多的效果，但更加简单。它属于一种贪心算法，即每一步都尽可能减小误差。一开始，所有权重都设为1，然后每一步所做的决策是对某个权重增加或减少一个很小的值。

前向逐步回归伪代码

数据标准化，使其分布满足0均值和单位方差

在每轮迭代过程中：

设置当前最小误差lowest_error为正无穷

对每个特征：

增大或缩小：

改变一个系数得到一个新的 $w$

计算新 $w$ 下的误差error

如果误差error小于当前最小误差lowest_error：设置 $\hat{w}$ 等于当前的 $w$

将 $w$ 设置为新的 $\hat{w}$

代码：

import numpy as np


def load_dataset(filename):
    feature_number = len(open(filename).readline().split('\t'))-1
    data_array = []
    label_array = []
    file = open(filename)
    for line in file.readlines():
        line_array = []
        curr_line = line.strip().split('\t')
        for i in range(feature_number):
            line_array.append(float(curr_line[i]))
        data_array.append(line_array)
        label_array.append(curr_line[-1])
    return data_array, label_array

def regularize(x_mat):
    in_mat = x_mat.copy()
    in_means = np.mean(in_mat, 0)
    in_var = np.var(in_mat, 0)
    in_mat = (in_mat - in_means) / in_var
    return in_mat

def calc_rss_error(y_array, y_hat_array):
    return ((y_array - y_hat_array) ** 2).sum()

def stepwise(x_array, y_array, eps=0.01, iter_num=100):
    x_mat = np.mat(x_array, dtype=np.float64)
    y_mat = np.mat(y_array, dtype=np.float64).T
    # 将特征值按照均值为0，方差为1进行标准化处理
    y_mean = np.mean(y_mat, 0)
    y_mat = y_mat - y_mean
    x_mat = regularize(x_mat)
    m, n = np.shape(x_mat)
    return_mat = np.zeros((iter_num, n))
    w_hat = np.zeros((n, 1))
    w_hat_test = w_hat.copy()
    w_hat_max = w_hat.copy()
    for i in range(iter_num):
        print("[iteration  %d]" % i, w_hat.T)
        lowest_error = np.inf
        for j in range(n):
            for sign in [-1, 1]:
                w_hat_test = w_hat.copy()
                w_hat_test[j] += eps * sign
                y_test = x_mat * w_hat_test
                rss_error = calc_rss_error(y_mat.A, y_test.A)
                if rss_error < lowest_error:
                    lowest_error = rss_error
                    w_hat_max = w_hat_test
        w_hat = w_hat_max.copy()
        return_mat[i, :] = w_hat.T
    return return_mat

if __name__ == '__main__':
    x_array, y_array = load_dataset('abalone.txt')
    w_hat = stepwise(x_array, y_array, 0.01, 200)
    print("w_hat: ", w_hat)

CONTACT INFORMATION

E-mail: birdguan@seu.edu.cn

QQ: 46611253

本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)

Machine Learning