一、回归介绍
监督学习指的是有目标变量或预测目标的机器学习方法。回归与分类的不同，就在于其目标变量是连续数值型，最直接的办法是依据输入写出一个目标值的计算公式。假如你想要预测姐姐男友汽车的功率大小，可能会这么计算：
HorsePower = 0.0015 * annualSalary - 0.99 * hoursListeningToPublicRadio
这就是所谓的回归方程（regression equation），其中的0.0015和-0.99称作回归系数（regression weights），求这些回归系数的过程就是回归。一旦有了这些回归系数，再给定输入，做预测就非常容易了。具体的做法是用回归系数乘以输入值，再将结果全部加在一起，就得到了预测值。
说到回归，一般都是指线性回归（linear regression），所以本章里的回归和线性回归代表同一个意思。线性回归意味着可以将输入项分别乘以一些常量，再将结果加起来得到输出。需要说明的是，存在另一种称为非线性回归的回归模型，该模型不认同上面的做法，比如认为输出可能是输入的乘积。这样，上面的功率计算公式也可以写做：HorsePower = 0.0015 * annualSalary / hoursListeningToPublicRadio
二、用线性回归找到最佳拟合直线
全文代码在最后第八点里。
应该怎么从一大堆数据里求出回归方程呢？假定输入数据存放在矩阵X中，结果存放在向量y中：
而回归系数存放在向量w中：
那么对于给定的数据x1，即矩阵X的第一列数据，预测结果u1将会通过如下公式给出：
现在的问题是，手里有数据矩阵X和对应的标签向量y，怎么才能找到w呢？一个常用的方法就是找出使误差最小的w。这里的误差是指预测u值和真实y值之间的差值，使用该误差的简单累加将使得正差值和负差值相互抵消，所以我们采用平方误差。
平方误差和可以写做：
用矩阵表示还可以写做：
找到w，使平方误差和最小。因为我们认为平方误差和越小，说明线性回归拟合效果越好。
现在，我们用矩阵表示的平方误差和对w进行求导：
令上述公式等于0，得到：
w上方的小标记表示，这是当前可以估计出的w的最优解。从现有数据上估计出的w可能并不是数据中的真实w值，所以这里使用了一个"帽"符号来表示它仅是w的一个最佳估计。
值得注意的是，上述公式中包含逆矩阵，也就是说，这个方程只在逆矩阵存在的时候使用，也即是这个矩阵是一个方阵，并且其行列式不为0。
上述的最佳w求解是统计学中的常见问题，除了矩阵方法外还有很多其他方法可以解决。通过调用NumPy库里的矩阵方法，我们可以仅使用几行代码就完成所需功能。该方法也称作OLS， 意思是“普通小二乘法”（ordinary least squares）。
现有一数据集：
第一列都为1.0，即x0。第二列为x1，即x轴数据。第三列为x2，即y轴数据。首先绘制下数据，看下数据分布。
如何判断拟合曲线的拟合效果的如何呢？我们可以使用corrcoef方法计算这两个序列的相关系数来比较预测值和真实值的相关性。
可以看到，对角线上的数据是1.0，因为yMat和自己的匹配是完美的，而YHat和yMat的相关系数为0.98。
三、局部加权线性回归
线性回归的一个问题是有可能出现欠拟合现象，因为它求的是具有最小均方误差的无偏估计。显而易见，如果模型欠拟合将不能取得最好的预测效果。所以有些方法允许在估计中引入一些偏差，从而降低预测的均方误差。
中的一个方法是局部加权线性回归（Locally Weighted Linear Regression，LWLR）。在该方法中，我们给待预测点附近的每个点赋予一定的权重。与kNN一样，这种算法每次预测均需要事先选取出对应的数据子集。该算法解除回归系数W的形式如下：
其中 w 是一个矩阵，用来给每个数据点赋予权重。LWLR使用“核”（与支持向量机中的核类似）来对附近的点赋予更高的权重。核的类型可以自由选择，最常用的核就是高斯核，高斯核对应的权重如下：
这样就构建了一个只含对角元素的权重矩阵 w ，并且点 x 与 x(i) 越近， w(i,i) 将会越大。上述公式包含一个需要用户指定的参数 k ，它决定了对附近的点赋予多大的权重，这也是使用LWLR时唯一需要考虑的参数。
当k = 1.0时权重很大，如同将所有的数据视为等权重，得出的最佳拟合直线与标准的回归一致。使用k = 0.01得到了非常好的效果，抓住了数据的潜在模式。下图使用k = 0.003纳入了太多的噪声点，拟合的直线与数据点过于贴近。所以，图中的最下图是过拟合的一个例子，而最上图则是欠拟合的一个例子。
四、示例：预测鲍鱼的年龄
abalone.txt文件中记录了鲍鱼的年龄，这个数据来自UCI数据集合的数据。鲍鱼年龄可以从鲍鱼壳的层数推算得到。
最后一列代表的是鲍鱼的真实年龄，前面几列的数据是一些鲍鱼的特征，例如鲍鱼壳的层数等。
可以看到，当k=0.1时，训练集误差小，但是应用于新的数据集之后，误差反而变大了。这就是经常说道的过拟合现象。我们训练的模型，我们要保证测试集准确率高，这样训练出的模型才可以应用于新的数据，也就是要加强模型的普适性。可以看到，当k=1时，局部加权线性回归和简单的线性回归得到的效果差不多。这也表明一点，必须在未知数据上比较效果才能选取到最佳模型。那么最佳的核大小是10吗？或许是，但如果想得到更好的效果，应该用10个不同的样本集做10次测试来比较结果。
四、缩减系数来“理解”数据
如果数据的特征比样本点还多，就不可以使用线性回归和之前的方法来做预测了，这是因为在计算 (X^T X)^ -1的时候会出错。如果特征比样本点还多（n > m），也就是说输入数据的矩阵X不是满秩矩阵。非满秩矩阵在求逆时会出现问题。
为了解决这个问题，统计学家引入了岭回归（ridge regression）、lasso法和前向逐步回归。
4.1岭回归
岭回归即我们所说的L2正则线性回归，在一般的线性回归最小化均方误差的基础上增加了一个参数w的L2范数的罚项，从而最小化罚项残差平方和：
简单说来，岭回归就是在矩阵X^TX上加一个λI 从而使得矩阵非奇异，进而能对X^TX + λI求逆。其中矩阵 I 是一个m×m的单位矩阵，对角线上元素全为1，其他元素全为0。而λ是一个用户定义的数值，后面会做介绍。在这种情况下，回归系数的计算公式将变成：
岭回归最先用来处理特征数多于样本数的情况，现在也用于在估计中加入偏差，从而得到更好的估计。这里通过引入λ来限制了所有w之和，通过引入该惩罚项，能够减少不重要的参数，这个技术在统计学中也可以叫做缩减（shrinkage）。
上图绘制了回归系数与log(λ)的关系。在最左边，即λ最小时，可以得到所有系数的原始值（与线性回归一致）；而在右边，系数全部缩减成0；在中间部分的某个位置，将会得到最好的预测结果。想要得到最佳的λ参数，可以使用交叉验证的方式获得。
4.2前向逐步线性回归
前向逐步线性回归算法属于一种贪心算法，即每一步都尽可能减少误差。我们计算回归系数，不再是通过公式计算，而是通过每次微调各个回归系数，然后计算预测误差。那个使误差最小的一组回归系数，就是我们需要的最佳回归系数。
上图是反应了迭代次数与回归系数的关系曲线。可以看到，有些系数从始至终都是约为0的，这说明它们不对目标造成任何影响，也就是说这些特征很可能是不需要的。逐步线性回归算法的优点在于它可以帮助人们理解有的模型并做出改进。当构建了一个模型后，可以运行该算法找出重要的特征，这样就有可能及时停止对那些不重要特征的收集。
五、示例、预测乐高玩具套件的价格
乐高（LEGO）公司生产拼装类玩具，由很多大小不同的塑料插块组成。一般来说，这些插块都是成套出售，它们可以拼装成很多不同的东西，如船、城堡、一些著名建筑等。乐高公司每个套装包含的部件数目从10件到5000件不等。
一种乐高套件基本上在几年后就会停产，但乐高的收藏者之间仍会在停产后彼此交易。本次实例，就是使用回归方法对收藏者之间的交易价格进行预测。
5.1、获取数据
全部的数据集在各个网页当中，需要利用爬虫只是解析抓取出来。
抓取数据结果如下：
5.2、建立模型
抓取下来数据集后，接下来就是训练模型。首先需要添加全为0的特征X0列。因为线性回归的第一列特征要求都是1.0。然后使用最简单的普通线性回归进行建模。训练出来的线性回归模型如截图最下面所示：
我们使用岭回归，通过交叉验证，找到使误差最小的λ对应的回归系数。
这里随机选取样本，因为其随机性，所以每次运行的结果可能略有不同。不过整体如上图所示，可以看出，它与常规的最小二乘法，即普通的线性回归没有太大差异。我们本期望找到一个更易于理解的模型，显然没有达到预期效果。
现在，我们看一下在缩减过程中回归系数是如何变化的。
看运行结果的第一行，可以看到最大的是第4项，第二大的是第2项。
因此，如果只选择一个特征来做预测的话，我们应该选择第4个特征，也就是原始加个。如果可以选择2个特征的话，应该选择第4个和第2个特征。
这种分析方法使得我们可以挖掘大量数据的内在规律。在仅有4个特征时，该方法的效果也许并不明显；但如果有100个以上的特征，该方法就会变得十分有效：它可以指出哪个特征是关键的，而哪些特征是不重要的。
六、应用scikit-learn构建线性回归模型
sklearn.linear_model提供了很多线性模型，包括岭回归、贝叶斯回归、Lasso等。现在我们用sklearn实现岭回归。
我正则化项系数设为0.5，其余参数使用默认。可以看到，获得的结果与上面的结果类似。
七、小结：
与分类一样，回归也是预测目标值的过程。回归与分类的不同点在于，前者预测连续类型变量，而后者预测离散类型变量。
在回归方程里，求得特征对应的最佳回归系数的方法是最小化误差的平方和。给定输入矩阵X，如果X^TX的逆存在并可以求得的话，回归法都可以直接使用。数据集上计算出的回归方程并不一定意味着它是最佳的，可以使用预测值y_hat和原始值y的相关性来度量回归方程的好坏。
当数据的样本数比特征数还少时候，矩阵 X^TX的逆不能直接计算。即便当样本数比特征数多时， X^TX的逆仍有可能无法直接计算，这是因为特征有可能高度相关。这时可以考虑使用岭回归，因为当 X^TX的逆不能计算时，它仍保证能求得回归参数。
岭回归是缩减法的一种，相当于对回归系数的大小施加了限制。另一种很好的缩减法是lasso。Lasso难以求解，但可以使用计算简便的逐步线性回归方法来求得近似结果。
缩减法还可以看做是对一个模型增加偏差的同时减少方差。
八、贴上全部代码：
# -*- conding: utf-8 -*-
import matplotlib.pyplot as plt
from numpy import *
from matplotlib.font_manager import FontProperties
from bs4 import BeautifulSoup
def load_data_set(file_name):
"""
Function:
加载数据集
Parameters:
file_name - 文件名
Returns:
data_mat - 数据列表
label_mat - 标签列表
Modify:
2018-10-22
"""
num_feat = len(open(file_name).readline().split('\t')) - 1
data_mat = []
label_mat = []
fr = open(file_name)
for line in fr.readlines():
line_arr = []
cur_line = line.strip().split('\t')
for i in range(num_feat):
line_arr.append(float(cur_line[i]))
data_mat.append(line_arr)
label_mat.append(float(cur_line[-1]))
return data_mat, label_mat
def plot_data_set():
"""
Function:
绘制数据集
Parameters:
无
Returns:
Modify:
2018-10-22
"""
data_mat, label_mat = load_data_set('../Machine/machinelearninginaction/Ch08/ex0.txt')
n = len(data_mat)
xcord = []
ycord = []
for i in range(n):
xcord.append(data_mat[i][1])
ycord.append(label_mat[i])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xcord, ycord, s=20, c='blue', alpha=0.5)
plt.title('data_set')
plt.xlabel('X')
plt.show()
def stand_regres(x_arr, y_arr):
"""
Function:
计算回归系数w
Parameters:
x_arr - 数据列表
y_arr - 标签列表
Returns:
ws - 回归系数
Modify:
2018-10-22
"""
x_mat = mat(x_arr)
y_mat = mat(y_arr).T
xTx = x_mat.T * x_mat
if linalg.det(xTx) == 0.0:
print('This matrix is singular, cannot do inverse')
return
ws = xTx.I * (x_mat.T * y_mat)
return ws
def plot_regression():
"""
Function:
加载数据集
Parameters:
无
Returns:
dataMat - 数据列表
labelMat - 标签列表
Modify:
2018-10-22
"""
data_mat, label_mat = load_data_set('../Machine/machinelearninginaction/Ch08/ex0.txt')
ws = stand_regres(data_mat, label_mat)
x_mat = mat(data_mat)
y_mat = mat(label_mat)
x_coyp = x_mat.copy()
x_coyp.sort(0)
y_hat = x_coyp * ws
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x_coyp[:, 1], y_hat, c='red')
ax.scatter(x_mat[:, 1].flatten().A[0], y_mat.flatten().A[0], s=20, c='blue', alpha=0.5)
plt.title('data_set')
plt.xlabel('X')
plt.show()
def lwlr(test_point, x_arr, y_arr, k=1.0):
"""
Function:
使用局部加权线性回归计算回归系数w
Parameters:
test_point - 测试样本点
x_arr - x数据集
y_arr - y数据集
Returns:
ws - 回归系数
Modify:
2018-10-24
"""
x_mat = mat(x_arr)
y_mat = mat(y_arr).T
m = shape(x_mat)[0]
# 生成对角矩阵
weights = mat(eye((m)))
for j in range(m):
diff_mat = test_point - x_mat[j, :]
weights[j, j] = exp(diff_mat * diff_mat.T / (-2.0 * k ** 2))
x_t_x = x_mat.T * (weights * x_mat)
if linalg.det(x_t_x) == 0.0:
print('矩阵为奇异矩阵,不能求逆')
return
# 计算回归系数
ws = x_t_x.I * (x_mat.T * (weights * y_mat))
return test_point * ws
def lwlr_test(test_arr, x_arr, y_arr, k=1.0):
"""
Function:
局部加权线性回归测试
Parameters:
test_arr - 测试数据集
x_arr - x数据集
y_arr - y数据集
Returns:
dataMat - 数据列表
labelMat - 标签列表
Modify:
2018-10-24
"""
m = shape(test_arr)[0]
y_hat = zeros(m)
for i in range(m):
y_hat[i] = lwlr(test_arr[i], x_arr, y_arr, k)
return y_hat
def plot_lwlr_regression():
"""
Function:
绘制多条局部加权回归曲线
Parameters:
无
Returns:
无
Modify:
2018-10-24
"""
font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
x_arr, y_arr = load_data_set('../Machine/machinelearninginaction/Ch08/ex0.txt')
y_hat_1 = lwlr_test(x_arr, x_arr, y_arr, 1.0)
y_hat_2 = lwlr_test(x_arr, x_arr, y_arr, 0.01)
y_hat_3 = lwlr_test(x_arr, x_arr, y_arr, 0.003)
x_mat = mat(x_arr)
y_mat = mat(y_arr)
# 排序，返回索引值
srt_ind = x_mat[:, 1].argsort(0)
x_sort = x_mat[srt_ind][:, 0, :]
fig, axs = plt.subplots(nrows=3, ncols=1, sharex=False, sharey=False, figsize=(10, 8))
# 绘制回归曲线
axs[0].plot(x_sort[:, 1], y_hat_1[srt_ind], c='red')
axs[1].plot(x_sort[:, 1], y_hat_2[srt_ind], c='red')
axs[2].plot(x_sort[:, 1], y_hat_3[srt_ind], c='red')
# 绘制原数据点
axs[0].scatter(x_mat[:, 1].flatten().A[0], y_mat.flatten().A[0], s=20, c='blue', alpha=0.5)
axs[1].scatter(x_mat[:, 1].flatten().A[0], y_mat.flatten().A[0], s=20, c='blue', alpha=0.5)
axs[2].scatter(x_mat[:, 1].flatten().A[0], y_mat.flatten().A[0], s=20, c='blue', alpha=0.5)
# 设置标题,x轴label,y轴label
axs0_title_text = axs[0].set_title(u'局部加权回归曲线,k=1.0', FontProperties=font)
axs1_title_text = axs[1].set_title(u'局部加权回归曲线,k=0.01', FontProperties=font)
axs2_title_text = axs[2].set_title(u'局部加权回归曲线,k=0.003', FontProperties=font)
# plt.setp():设置图标实例的属性
plt.setp(axs0_title_text, size=8, weight='bold', color='red')
plt.setp(axs1_title_text, size=8, weight='bold', color='red')
plt.setp(axs2_title_text, size=8, weight='bold', color='red')
plt.xlabel('X')
plt.show()
def rss_error(y_arr, y_hat_arr):
"""
Function:
误差大小评价函数
Parameters:
y_arr - 真实数据
y_hat_arr - 预测数据
Returns:
误差大小
Modify:
2018-10-24
"""
return ((y_arr - y_hat_arr) ** 2).sum()
def ridge_regress(x_mat, y_mat, lam=0.2):
"""
Function:
岭回归
Parameters:
x_mat - x数据集
y_mat - y数据集
lam - 缩减系数
Returns:
ws - 回归系数
Modify:
2018-11-07
"""
x_t_x = x_mat.T * x_mat
denom = x_t_x + eye(shape(x_mat)[1]) * lam
if linalg.det(denom) == 0.0:
print('矩阵为奇异矩阵,不能转置')
return
ws = denom.I * (x_mat.T * y_mat)
return ws
def ridge_test(x_arr, y_arr):
"""
Function:
岭回归测试
Parameters:
x_mat - x数据集
y_mat - y数据集
lam - 缩减系数
Returns:
w_mat - 回归系数矩阵
Modify:
2018-11-07
"""
x_mat = mat(x_arr)
y_mat = mat(y_arr).T
# 行与行操作，求均值
y_mean = mean(y_mat, axis=0)
# 数据减去均值
y_mat = y_mat - y_mean
# 行与行操作，求均值
x_means = mean(x_mat, axis=0)
# 行与行操作，求方差
x_var = var(x_mat, axis=0)
# 数据减去均值除以方差实现标准化
x_mat = (x_mat - x_means) / x_var
# 30个不同的lambda测试
num_test_pts = 30
# 初始回归系数矩阵
w_mat = zeros((num_test_pts, shape(x_mat)[1]))
for i in range(num_test_pts):
# lambda以e的指数变化，最初是一个非常小的数
ws = ridge_regress(x_mat, y_mat, exp(i - 10))
w_mat[i, :] = ws.T
return w_mat
def plot_ridge_regress_mat():
"""
Function:
绘制岭回归系数矩阵
Parameters:
无
Returns:
无
Modify:
2018-10-24
"""
font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
ab_x, ab_y = load_data_set('../Machine/machinelearninginaction/Ch08/abalone.txt')
redge_weights = ridge_test(ab_x, ab_y)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(redge_weights)
ax_title_text = ax.set_title(u'log(lambada)与回归系数的关系', FontProperties=font)
ax_xlabel_text = ax.set_xlabel(u'log(lambada)', FontProperties=font)
ax_ylabel_text = ax.set_ylabel(u'回归系数', FontProperties=font)
plt.setp(ax_title_text, size=20, weight='bold', color='red')
plt.setp(ax_xlabel_text, size=10, weight='bold', color='black')
plt.setp(ax_ylabel_text, size=10, weight='bold', color='black')
plt.show()
def regularize(x_mat, y_mat):
"""
Function:
数据标准化
Parameters:
x_mat - x数据集
y_mat - y数据集
Returns:
in_x_mat - 标准化后的x数据集
in_y_mat - 标准化后的y数据集
Modify:
2018-11-11
"""
in_x_mat = x_mat.copy()
in_y_mat = y_mat.copy()
in_y_mean = mean(y_mat, 0)
in_y_mat = y_mat - in_y_mean
in_means = mean(in_x_mat, 0)
in_var = var(in_x_mat, 0)
in_x_mat = (in_x_mat - in_means) / in_var
return in_x_mat, in_y_mat
def stage_wise(x_arr, y_arr, eps=0.01, num_it=100):
"""
Function:
前向逐步线性回归
Parameters:
x_arr - x输入数据
y_arr - y预测数据
eps - 每次迭代需要调整的步长
num_it - 迭代次数
Returns:
return_mat - 迭代次数
Modify:
2018-10-24
"""
x_mat = mat(x_arr)
y_mat = mat(y_arr).T
x_mat, y_mat = regularize(x_mat, y_mat)
m, n = shape(x_mat)
# 初始化num_it次迭代的回归系数矩阵
return_mat = zeros((num_it, n))
# 初始化回归系数矩阵
ws = zeros((n, 1))
ws_test = ws.copy()
ws_max = ws.copy()
# 迭代num_it次
for i in range(num_it):
# 打印当前回归系数矩阵
print(ws.T)
lowest_error = float('inf')
# 遍历每个特征的回归系数
for j in range(n):
for sign in [-1, 1]:
ws_test = ws.copy()
# 微调回归系数
ws_test[j] += eps * sign
# 计算预测值
y_test = x_mat * ws_test
# 计算平方误差
rss_e = rss_error(y_mat.A, y_test.A)
# 如果误差更小，则更新当前的最佳回归系数
if rss_e < lowest_error:
lowest_error = rss_e
ws_max = ws_test
ws = ws_max.copy()
# 记录numIt次迭代的回归系数矩阵
return_mat[i, :] = ws.T
return return_mat
def plot_stage_wise():
"""
Function:
绘制前向逐步线性回归
Parameters:
无
Returns:
无
Modify:
2018-10-24
"""
font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
ab_x, ab_y = load_data_set('../Machine/machinelearninginaction/Ch08/abalone.txt')
stage_weights = stage_wise(ab_x, ab_y, 0.005, 1000)
print(stage_weights)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(stage_weights)
ax_title_text = ax.set_title(u'前向逐步回归:迭代次数与回归系数的关系', FontProperties=font)
ax_xlabel_text = ax.set_xlabel(u'迭代次数', FontProperties=font)
ax_ylabel_text = ax.set_ylabel(u'回归系数', FontProperties=font)
plt.setp(ax_title_text, size=15, weight='bold', color='red')
plt.setp(ax_xlabel_text, size=10, weight='bold', color='black')
plt.setp(ax_ylabel_text, size=10, weight='bold', color='black')
plt.show()
def scrape_page(ret_x, ret_y, in_file, yr, num_pce, orig_prc):
"""
Function:
从html页面读取数据，生成ret_x和ret_y列表
Parameters:
ret_x - 数据x
ret_y - 数据y
in_file - HTML文件
yr - 年份
num_pce - 乐高部件数目
orig_prc - 原价
Returns:
无
Modify:
2018-10-24
"""
# 打开并读取HTML文件
with open(in_file, encoding='utf-8') as f:
html = f.read()
soup = BeautifulSoup(html, 'lxml')
i = 1
# 根据HTML页面结构进行解析
current_row = soup.find_all('table', r='%d' % i)
while (len(current_row) != 0):
current_row = soup.find_all('table', r='%d' % i)
title = current_row[0].find_all('a')[1].text
lwr_title = title.lower()
# 查找是否有全新标签
if (lwr_title.find('new') > -1) or (lwr_title.find('nisb') > -1):
new_flag = 1.0
else:
new_flag = 0.0
# 查找是否已经标志出售，我们只收集已出售的数据
sold_unicde = current_row[0].find_all('td')[3].find_all('span')
if len(sold_unicde) == 0:
print('商品 #%d 没有出售' % i)
else:
sold_price = current_row[0].find_all('td')[4]
price_str = sold_price.text
price_str = price_str.replace('$', '')
price_str = price_str.replace(',', '')
if len(sold_price) > 1:
price_str = price_str.replace('Free shipping', '')
selling_price = float(price_str)
if selling_price > orig_prc * 0.5:
print("%d\t%d\t%d\t%f\t%f" % (yr, num_pce, new_flag, orig_prc, selling_price))
ret_x.append([yr, num_pce, new_flag, orig_prc])
ret_y.append(selling_price)
i += 1
current_row = soup.find_all('table', r='%d' % i)
def set_data_collect(ret_x, ret_y):
"""
Function:
依次读取六种乐高套装的数据，并生成数据矩阵
Parameters:
ret_x - 数据x
ret_y - 数据y
Returns:
无
Modify:
2018-10-24
"""
file_path = './Lego/'
scrape_page(ret_x, ret_y, file_path + 'lego8288.html', 2006, 800, 49.99)
scrape_page(ret_x, ret_y, file_path + 'lego10030.html', 2002, 3096, 269.99)
scrape_page(ret_x, ret_y, file_path + 'lego10179.html', 2007, 5195, 499.99)
scrape_page(ret_x, ret_y, file_path + 'lego10181.html', 2007, 3428, 99.99)
scrape_page(ret_x, ret_y, file_path + 'lego10189.html', 2008, 5922, 299.99)
scrape_page(ret_x, ret_y, file_path + 'lego10196.html', 2009, 3263, 249.99)
def use_stand_regres():
"""
Function:
使用简单的线性回归进行建模
Parameters:
ret_x - 数据x
ret_y - 数据y
Returns:
无
Modify:
2018-10-24
"""
lg_x = []
lg_y = []
set_data_collect(lg_x, lg_y)
# lg_y = mat(lg_y).T
data_num, feature_num = shape(lg_x)
lg_x1 = mat(ones((data_num, feature_num + 1)))
lg_x1[:, 1:5] = mat(lg_x)
ws = stand_regres(lg_x1, lg_y)
print('%f%+f*年份%+f*部件数量%+f*是否为全新%+f*原价' % (ws[0], ws[1], ws[2], ws[3], ws[4]))
def cross_validation(x_arr, y_arr, num_val=10):
"""
Function:
交叉验证岭回归
Parameters:
x_arr - x数据集
y_arr - y数据集
num_val - 交叉验证次数
Returns:
w_mat - 回归系数矩阵
Modify:
2018-10-24
"""
# 统计样本个数
m = len(y_arr)
# 生成索引值列表
index_list = list(range(m))
# create error mat 30columns num_val rows
error_mat = zeros((num_val, 30))
# 交叉验证num_val次
for i in range(num_val):
# 训练集
train_x = []
train_y = []
# 测试集
test_x = []
test_y = []
# 打乱次序
random.shuffle(index_list)
# 划分数据集:90%训练集，10%测试集
for j in range(m):
if j < m * 0.9:
train_x.append(x_arr[index_list[j]])
train_y.append(y_arr[index_list[j]])
else:
test_x.append(x_arr[index_list[j]])
test_y.append(y_arr[index_list[j]])
# 获得30个不同lambda下的岭回归系数
w_mat = ridge_test(train_x, train_y)
# 遍历所有的岭回归系数
for k in range(30):
# 测试集
mat_test_x = mat(test_x)
mat_train_x = mat(train_x)
# 测试集均值
mean_train = mean(mat_train_x, 0)
# 测试集方差
var_train = var(mat_train_x, 0)
# 测试集标准化
mat_test_x = (mat_test_x - mean_train) / var_train
# 根据ws预测y值
y_est = mat_test_x * mat(w_mat[k, :]).T + mean(train_y)
# 统计误差
error_mat[i, k] = rss_error(y_est.T.A, array(test_y))
# 计算每次交叉验证的平均误差
mean_errors = mean(error_mat, 0)
# 找到最小误差
min_mean = float(min(mean_errors))
# 找到最佳回归系数
best_weights = w_mat[nonzero(mean_errors == min_mean)]
x_mat = mat(x_arr)
y_mat = mat(y_arr).T
mean_x = mean(x_mat, 0)
var_x = var(x_mat, 0)
# 数据经过标准化，因此需要还原
un_reg = best_weights / var_x
print('%f%+f*年份%+f*部件数量%+f*是否为全新%+f*原价' % (
(-1 * sum(multiply(mean_x, un_reg)) + mean(y_mat)), un_reg[0, 0], un_reg[0, 1], un_reg[0, 2], un_reg[0, 3]))
if __name__ == '__main__':
# plot_data_set()
# plot_regression()
# data_mat, label_mat = load_data_set('../Machine/machinelearninginaction/Ch08/ex0.txt')
# ws = stand_regres(data_mat, label_mat)
# x_mat = mat(data_mat)
# y_mat = mat(label_mat)
# y_hat = x_mat * ws
# print(corrcoef(y_hat.T, y_mat))
# plot_lwlr_regression()
# ab_x, ab_y = load_data_set('../Machine/machinelearninginaction/Ch08/abalone.txt')
# print('训练集与测试集相同:局部加权线性回归,核k的大小对预测的影响:')
# y_hat_01 = lwlr_test(ab_x[0:99], ab_x[0:99], ab_y[0:99], 0.1)
# y_hat_1 = lwlr_test(ab_x[0:99], ab_x[0:99], ab_y[0:99], 1)
# y_hat_10 = lwlr_test(ab_x[0:99], ab_x[0:99], ab_y[0:99], 10)
# print('k=0.1时,误差大小为:', rss_error(ab_y[0:99], y_hat_01.T))
# print('k=1 时,误差大小为:', rss_error(ab_y[0:99], y_hat_1.T))
# print('k=10 时,误差大小为:', rss_error(ab_y[0:99], y_hat_10.T))
# print('')
# print('训练集与测试集不同:局部加权线性回归,核k的大小是越小越好吗？更换数据集,测试结果如下:')
# y_hat_1 = lwlr_test(ab_x[100:199], ab_x[0:99], ab_y[0:99], 0.1)
# y_hat_2 = lwlr_test(ab_x[100:199], ab_x[0:99], ab_y[0:99], 1)
# y_hat_3 = lwlr_test(ab_x[100:199], ab_x[0:99], ab_y[0:99], 10)
# print('k=0.1时,误差大小为:', rss_error(ab_y[100:199], y_hat_1.T))
# print('k=1 时,误差大小为:', rss_error(ab_y[100:199], y_hat_2.T))
# print('k=10 时,误差大小为:', rss_error(ab_y[100:199], y_hat_3.T))
# print('')
# print('训练集与测试集不同:简单的线性归回与k=1时的局部加权线性回归对比:')
# print('k=1时,误差大小为:', rss_error(ab_y[100:199], y_hat_2.T))
# ws = stand_regres(ab_x[0:99], ab_y[0:99])
# y_hat = mat(ab_x[100:199]) * ws
# print('简单的线性回归误差大小:', rss_error(ab_y[100:199], y_hat.T.A))
# plot_ridge_regress_mat()
# plot_stage_wise()
# lg_x = []
# lg_y = []
# set_data_collect(lg_x, lg_y)
# use_stand_regres()
# lg_x = []
# lg_y = []
# set_data_collect(lg_x, lg_y)
# cross_validation(lg_x, lg_y)
# lg_x = []
# lg_y = []
# set_data_collect(lg_x, lg_y)
# print(ridge_test(lg_x, lg_y))
from sklearn import linear_model
reg = linear_model.Ridge(alpha=0.5)
lg_x = []
lg_y = []
set_data_collect(lg_x, lg_y)
reg.fit(lg_x, lg_y)
print('%f%+f*年份%+f*部件数量%+f*是否为全新%+f*原价' % (reg.intercept_, reg.coef_[0], reg.coef_[1], reg.coef_[2], reg.coef_[3]))
网友评论