[Python知识库]线性回归和Python代码实现

线性回归

针对特征$X=[x_1,x_2,...,x_m]$，真值为$y$。
如果使用线性回归模型描述真值与特征之间的关系，可以表达为
$$\hat{y}=w_1?x_1+w_2?x_2+...+w_m?x_m+b$$
其中，$\hat{y}$为预测值，$W=[w_1,w_2,...,w_m]$为权重系数，$b$为截距。

为了找到最佳的$W$和$b$，可以建立最优化模型，最小化MSE
$$min(y?\hat{y}{)}^2$$
MSE越小，说明预测值和真值的差异越小，预测模型越准确。

当数据集中包含多组数据时，即$y=[y_1,y_2,...,y_n]$，上式变为
$$min\sum _{i=1}^n(y_i?{\hat{y}}_i{)}^2$$

自编代码

本节考虑最简单的情况，$m=1$，即只有一个特征。
此时$\hat{y}$的表达式可以简化为
$$\hat{y}=w?x+b$$
最小化MSE的表达式变为
$$min\sum _{i=1}^n(y_i?(w?x_i+b){)}^2$$
需要注意的是，上式中$w$和$b$是未知量，$x_i,y_i$为已知量。
将上式按照$w$和$b$进行归并，得到
$$min(\sum _{i=1}^n{x}_i^2?w^2+n?b^2+\sum _{i=1}^{n}2x_i?wb?\sum _{i=1}^{n}2x_i{y}_i?w?\sum _{i=1}^{n}2y_i?b+\sum _{i=1}^n{y}_i^2)$$

从式中可以看出，该问题为一个二次规划问题，因此可以调用cvxopt包进行求解。
cvxopt的使用方法可以参考python求解二次规划问题。

以下代码实现了上述问题的优化求解。
由于${\sum }_{i=1}^n{y}_i^2$是常量，不影响$w$和$b$的最优化，所以代码中，没有包含此项。

from cvxopt import matrix, solvers
from sklearn.model_selection import train_test_split


def lr_by_self(X, y):
    # 数据集拆分
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)

    # 构造二次规划的P
    f1 = [2 * (X_train ** 2).sum(), 2 * X_train.sum()]
    f2 = [2 * X_train.sum(), 2 * len(X_train)]
    P = matrix([f1, f2])

    # 构造二次规划的q
    q = [-2 * (X_train.T * y_train).sum(), -2 * y_train.sum()]
    q = matrix([q])

    # 调用cvxopt求解二次规划
    result = solvers.qp(P, q)

    print('lr_by_self: w = {:.2f}, b = {:.2f}'.format(result['x'][0], result['x'][1]))

sklearn代码

在sklearn中，可以直接调用sklearn.linear_model中的LinearRegression模块，实现线性回归的目标。
需要注意的是，虽然在优化时，目标函数为MSE，但是其自带的.score()函数默认指标是r2。

from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score


def lr_by_sklearn(X, y):
    # 数据集拆分
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)
    lr = LinearRegression().fit(X_train, y_train)

    print('lr_by_sklearn: w = {:.2f}, b = {:.2f}'.format(lr.coef_[0], lr.intercept_))

    # 默认性能指标是r2
    print('training set score: {:.2f}'.format(lr.score(X_train, y_train)))
    print('test set score: {:.2f}'.format(lr.score(X_test, y_test)))

    # r2_score中，必须是真值在前，预测值在后
    print('training set r2 score: {:.2f}'.format(r2_score(y_train, lr.predict(X_train))))
    print('test set r2 score: {:.2f}'.format(r2_score(y_test, lr.predict(X_test))))

代码验证

使用wave数据集进行测试。

if __name__ == '__main__':
    X_arr, y_arr = wave()
    lr_by_self(X_arr, y_arr)
    lr_by_sklearn(X_arr, y_arr)

程序输出结果如下所示。
显然，自编代码和sklearn代码得到的$w$和$b$是一致的。
同时，sklearn代码中的.score()函数也被确认是r2。

lr_by_self: w = 0.39, b = -0.03
lr_by_sklearn: w = 0.39, b = -0.03
training set score: 0.67
test set score: 0.66
training set r2 score: 0.67
test set r2 score: 0.66