最小二乘法多元线性回归（岭回归实现）

采用岭回归的方式实现，数据集采用sklearn库中的波士顿房价数据集。

多元线性回归岭回归实现

一、理论分析

  岭回归的相关推导已经有很多了，这里不再赘述，这里主要分析在波士顿房价数据集中的一些分析。
由岭回归的损失函数：

进一步推导可得知参数

  （其中λ是系数，Ip是单位矩阵，矩阵的阶数和X^TX的阶数相同。
  对于波士顿房价数据集，输入X是有13项特征值，所以XTX是(13,13)的矩阵，故Ip是(13,13)的单位矩阵；系数λ取exp（-2），我查阅了一些资料，感觉系数用e^x函数比较适合，经过我的测试，发现当x=-2时拟合效果比较好。

二、实现代码

  代码包含两部分，一是岭回归的训练模型、预测及评估模块，全部封装在类LinRegression中，第二部分是调用LinRegression类训练、预测、评估，同时调用sklearn库中的LinearRegression包、Ridge包，与我编写LinRegression类做对比。

第一部分（LinRegression类）

import numpy as np
from sklearn.metrics import r2_score


class LinRegression:
    # 初始化
    def __init__(self):
        self.coef_ = None  # 系数（theta0~1）
        self.interception_ = None  # 截距（theta0/b 数）
        self._theta = None  # 全部的theta（向量）

    # 根据训练集训练模型
    def train(self, X_train, y_train):
        X_b = np.hstack([np.ones([len(X_train), 1]), X_train])  # 在训练集前添加一列全一（便于截距b计算）
        self._theta = np.linalg.inv((X_b.T.dot(X_b)) + np.exp(-2) * np.eye(X_b.shape[1])).dot(X_b.T).dot(y_train)  # 利用公式计算系数
        self.interception_ = self._theta[0]
        self.coef_ = self._theta[1:]
        return self

    # 根据输入的测试集预测
    def predict(self, X_predict):
        X_b = np.hstack([np.ones((len(X_predict), 1)), X_predict])
        y_predict = X_b.dot(self._theta)
        return y_predict

    # 用r2_score评估模型
    def score(self, X_test, y_test):
        y_predict = self.predict(X_test)
        score = r2_score(y_test, y_predict)
        return score

第二部分（main，调用类并做对比）

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from LinearRegression import LinRegression
from matplotlib import rcParams
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.metrics import r2_score

boston = datasets.load_boston()  # 导入波士顿房价数据集
X = boston.data
y = boston.target
X = X[y < 50.0]  # 选取房价小于50的数据集
y = y[y < 50.0]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=666)


# 方法1：岭回归
reg = LinRegression()
reg.train(X_train, y_train)  # 使用训练集训练模型
y_predict = reg.predict(X_test)  # 将测试集输入模型并输出
score = reg.score(X_test, y_test)  # 使用r2_score方法评价训练效果

'''
# 方法2：调用LinearRegression包
reg = LinearRegression()
reg.fit(X_train, y_train)
y_predict = reg.predict(X_test)
score = r2_score(y_test, y_predict)


# 方法3：调用Ridge包
reg = Ridge()
reg.fit(X_train, y_train)
y_predict = reg.predict(X_test)
score = r2_score(y_test, y_predict)
'''
print("score:", score)

# 测试集输出与模型输出图形化显示
rcParams['font.sans-serif'] = 'SimHei'
fig = plt.figure(figsize=(10, 6))
plt.plot(range(y_test.shape[0]), y_test, color='blue', linewidth=1.5, linestyle='-')
plt.plot(range(y_test.shape[0]), y_predict, color='red', linewidth=1.5, linestyle='-.')
plt.legend(['真实值', '预测值'])
plt.show()

三、代码解读及结果分析

1.代码解读

1）LinRegression类
  代码与公式的对应：
  第一行是对输入X进行处理，在X前添加一列全一，这样是为了计算截距b；
  第二行是对公式的实现，使用numpy的函数即可，其中

是一个单位矩阵乘系数，代表λIp

2）main，调用类
  数据集选用
  导入波士顿房价数据集后对数据集进行筛选，只选取y<50的数据；因为如果将全部数据集导入，再训练模型和测试，会发现最后结果比较差，使用r2_score函数计算的结果只有0.6左右。然后我查阅资料后发现，y>=50的数据在使用岭回归的时候比较差，所以我舍弃了这部分数据，最后r2_score函数计算的结果有0.81左右

2.结果分析

  代码中我实现了三种方式的最小二乘法，第一种是自己编写底层代码的岭回归，第二、三种是分别调用LinearRegression、Ridge来实现模型训练与测试。使用同样的训练集和测试集，最后对比三种方式的结果。

1）岭回归
r2_score结果：0.8140904410855684

2）LinearRegression
r2_score结果：0.8129794056212811

3）Ridge
r2_score结果：0.8090190492614575

  通过上述的对比可以发现，三者相差并不大，甚至Ridge（sklean库的岭回归实现）还略微不如线性回归和我自己实现的岭回归，我觉得可能一是数据集比较多，特征比较多，所以三者训练相差不大，二是Ridge选取的系数和我的不同。