最小二乘法多元线性回归(岭回归实现)
采用岭回归的方式实现,数据集采用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选取的系数和我的不同。