机器学习实战 -- 预测数值型数据: 回归
标签： 机器学习实战 -- 预测数值型数据: 回归 代码人生博客 51CTO博客
2023-07-20 18:24:17 94浏览
学习内容
线性回归
局部加权线性回归
岭回归和逐步线性回归
预测鲍鱼年龄和玩具售价
回归的概念在中学阶段已经接触过了,这里简单介绍 .
形如 y = ax + b 这样一个函数方程 , 我们已经知道 大量的 (x,y) , 让我找到一组合适的参数 a,b ,使得我们预测的函数和 实际的函数误差很小。简单来说 a ,b 叫做 回归系数 , 求 a , b 的过程就是回归。
那么应该怎么从一大堆数据里求出回归方程呢? 假定输入数据放在 矩阵 X 中 , 回归系数放在 向量W 中 .对于给定的数据 X1 , 预测结果通过 Y1 = X1^T w 给出。我们考虑找出使 预测的 yhat 与 实际的 y 之间的差距最小,也即误差最小。我们用平方误差计算。
通过变换,并对 w 求导 ,得
标准回归函数和数据导入函数
# regression.py
import numpy as np
# 加载数据 生成数据和 标签
def loadDataSet(filename) :
dataMat =[]
labelMat = []
with open(filename,'r') as fr :
for line in fr.readlines() :
curLine = line.strip().split('\t')
lineArr = [float(item) for item in curLine]
dataMat.append(lineArr[:-1])
labelMat.append(lineArr[-1])
return dataMat,labelMat
# 这个函数来求最佳 w
def standRegres(xArr,yArr) :
xMat = np.matrix(xArr)
yMat = np.matrix(yArr).T
xTx = np.dot(xMat.T,xMat)
# 如果该矩阵的行列式为0 , 不存在逆
if np.linalg.det(xTx) == 0.0 :
print("This matrix is singular , cannot do inverse ")
return None
ws = xTx.I *(xMat.T *yMat)
return ws>>> xArr, yArr = regression.loadDataSet('ex0.txt')
>>> xArr[0:2]
[[1.0, 0.067732], [1.0, 0.42781]]>>> ws = regression.standRegres(xArr,yArr)
>>> ws
matrix([[3.00774324],
[1.69532264]])>>> import numpy as np
>>> xMat = np.matrix(xArr)
>>> yMat = np.matrix(xArr)
>>> yHat = np.dot(xMat,ws)
>>> yHat
matrix([[3.12257084],
[3.73301922],
[4.69582855],
[4.25946098],
[4.67099547],
... ..........fig = plt.figure()
ax = fig.add_subplot(111)
# 画出实际的散点图
ax.scatter(xMat[:, 1].A, yMat.A)
ax.set_xlabel("$x$")
ax.set_ylabel("$y$")
# 画出回归方程
xCopy = xMat.copy()
xCopy.sort(0)
yHat = np.dot(xCopy, ws)
ax.plot(xCopy[:, 1], yHat)
plt.show()
ex0.txt的数据集和它的最佳拟合直线
计算预测值和真实值的相关系数
In [1]: import regression
In [2]: import numpy as np
In [3]: xArr,yArr = regression.loadDataSet('ex0.txt')
In [4]: xMat = np.matrix(xArr)
In [5]: yMat = np.matrix(yArr)
In [6]: ws = regression.standRegres(xArr,yArr)
In [7]: yHat = np.dot(xMat,ws)
In [8]: print(yHat.shape)
(200, 1)
In [9]: print(yMat.shape)
(1, 200)
In [10]: print(np.corrcoef(yHat.T,yMat))
[[1. 0.98647356]
[0.98647356 1. ]]结果是 副对角线是 : 0.98 , 可以看出相关性很大 ,
局部加权线性回归
为了减缓欠拟合现象,提供模型的泛化能力. 采用局部加权线性回归 。
在最小均方误差（lowest mean-squared error）意义下，线性回归实现对统计量的无偏估计，但线性回归倾向于欠拟合（underfit）。通过增加偏置项，线性回归能够进一步减小均方误差。
*局部加权线性回归（Locally weighted linear regression，LWLR）*通过对感兴趣数据点附近的数据赋予权值，减少均方误差：
其中，W \mathbf{W}W为数据点的权值矩阵。LWLR通过类似“核”技巧的方式，使邻近数据点（nearby points）相比其它数据点或的更大的权重。常用核为高斯（Gaussian）核函数：
该方程定义了一个对角（diagonal）权值矩阵W，数据点x(i) 与x 越近，其权值w(i,i) 越大。预定义的参数k 用于指定x的影响范围（权值衰减速度）。
# regression.py
def lwlr(testPoint,xArr,yArr,k=1.0) :
xMat = np.matrix(xArr)
yMat = np.matrix(yArr).T
m = xMat.shape[0]
# 创建对角矩阵
weights = np.matrix(np.eye(m))
for j in range(m) :
diffMat = testPoint - xMat[j,:]
weights[j,j] = np.exp(diffMat*diffMat.T/(-2.0*k**2))
xTx = xMat.T * weights *xMat
if np.linalg.det(xTx) == 0.0 :
print("This matrix is singular , cannot do inverse ")
return None
ws = xTx.I * xMat.T * weights *yMat
return ws
def lwlrTest(testArr,xArr,yArr,k=1.0) :
m = np.shape(testArr)[0]
yHat = np.zeros(m)
for i in range(m):
yHat[i] = np.dot(testArr[i],lwlr(testArr[i],xArr,yArr,k))
return yHat单点进行估计 :
>>> yArr[0]
3.176513
>>> regression.lwlr(xArr[0],xArr,yArr,1.0)
matrix([[3.00714496],
[1.69638809]])
>>> regression.lwlr(xArr[0],xArr,yArr,0.001)
matrix([[2.85149878],
[5.1712412 ]])为了得到数据集所有点的估计, 调用 lwlrTest()函数
>>> yHat = regression.lwlrTest(xArr,xArr,yArr,0.003)
>>> yHat[:10]
array([3.20200665, 3.75940186, 4.53670134, 4.25050564, 4.56094936,
3.93721635, 3.53392289, 3.15405352, 3.12604366, 3.14881027])matplotlib 绘制图像 :
# test.py
yHat =regression.lwlrTest(xArr,xArr,yArr,1)
xMat = np.matrix(xArr)
#
srtInd = xMat[:,1].argsort(0)
# print("srtInd : " , srtInd)
xSort = xMat[srtInd][:,0,:]
# print("xSort : ",xSort)
# 第一个图
fig = plt.figure(figsize=(5,8))
ax = fig.add_subplot(311)
ax.plot(xSort[:,1] ,yHat[srtInd] )
ax.scatter(xMat[:,1].A ,yArr ,s=2 ,c ='red')
# print("xMAt " , xMat[:,1].A)
# 第二个图
yHat =regression.lwlrTest(xArr,xArr,yArr,0.01)
ax = fig.add_subplot(312)
ax.plot(xSort[:,1],yHat[srtInd])
ax.scatter(xMat[:,1].A,yArr,s =2 ,c ="red")
yHat =regression.lwlrTest(xArr,xArr,yArr,0.003)
ax = fig.add_subplot(313)
ax.plot(xSort[:,1],yHat[srtInd])
ax.scatter(xMat[:,1].A,yArr,s=2,c="red")
plt.show()
示例 : 预测鲍鱼的年龄
import numpy as np
import regression
abX, abY = regression.loadDataSet("abalone.txt")
yHat01 = regression.lwlrTest(abX[0: 99], abX[0: 99], abY[0: 99], 0.1)
yHat1 = regression.lwlrTest(abX[0: 99], abX[0: 99], abY[0: 99], 1)
yHat10 =regression.lwlrTest(abX[0: 99], abX[0: 99], abY[0: 99], 10)
print(regression.rssError(abY[0 : 99], yHat01.T))
print(regression.rssError(abY[0 : 99], yHat1.T))
print(regression.rssError(abY[0 : 99], yHat10.T))
56.84000194514532
429.89056187017104
549.1181708826828预测新数据 :
abX, abY = loadDataSet("abalone.txt")
abY = np.array(abY)
yHat01 = lwlrTest(abX[100 : 199], abX[0 : 99], abY[0 : 99], 0.1)
yHat1 = lwlrTest(abX[100 : 199], abX[0 : 99], abY[0 : 99], 1)
yHat10 = lwlrTest(abX[100 : 199], abX[0 : 99], abY[0 : 99], 10)
print(rssError(abY[100 : 199], yHat01.T))
print(rssError(abY[100 : 199], yHat1.T))
print(rssError(abY[100 : 199], yHat10.T))40659.27596541478
573.5261441898057
517.5711905382693在新数据(测试集)这时 k = 0.1 时 误差却很大, 说明 k = 0.1 时过拟合现象严重。
再来和简单的线性回归作比较 :
ws = regression.standRegres(abX[0 : 99], abY[0 : 99])
yHat = np.dot(np.matrix(abX[100 : 199]), ws)
print(regression.rssError(abY[100 : 199], yHat.A.T))
# 误差 518.6363153245542缩减系数来"理解" 数据
如果数据的特征比样本点还多增么办? 这时就不能用线性回归和之前的方法来做测试了.这时因为
不在可逆。
矩阵X 不是满秩矩阵. 这时引入了岭回归 。
岭回归
简单来说,岭回归就是在矩阵
上加一个
从而使得矩阵非奇异,进而对
求逆。 I 是 m × m 的单位矩阵 。在这种情况下,回归系数的计算公式将变成 :
岭回归最先用来处理特征数多于样本数的情况,现在也用于在估计中加入偏差,从而得到更好的估计。这里引入
来限制了所有 w 之和 ,通过引入该该惩罚项,防止过拟合,提供模型的泛化能力。
# regression.py
# 岭回归
def ridgeRegres(xMat, yMat, lam=0.2):
m, n = xMat.shape
xTx = np.dot(xMat.T, xMat)
denom = xTx + lam * np.eye(n)
if np.linalg.det(denom) == 0:
print("This matrix is singular, can not do inverse")
return None
ws = denom.I * np.dot(xMat.T, yMat)
return ws
def ridgeTest(xArr, yArr):
xMat = np.matrix(xArr)
yMat = np.matrix(yArr).T
# 标准化
m, n = xMat.shape
yMat = yMat - yMat.mean()
xMat = (xMat - xMat.mean(axis=0)) / xMat.var(axis=0)
numTestPts = 30
wMat = np.zeros((numTestPts, n))
for i in range(numTestPts):
ws = ridgeRegres(xMat, yMat, np.exp(i - 10))
wMat[i, :] = ws.T
return wMatabX, abY = regression.loadDataSet("abalone.txt")
ridgeWeights = regression.ridgeTest(abX, abY)
log_lam = [item - 10 for item in range(30)]
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(log_lam, ridgeWeights)
ax.set_xlabel("log($\\lambda$)")
plt.show()
该图绘出了回归系数与 log(
) 的关系。在最左边,即
最小时,可以得到所有系数的原始值(与线性回归一致);而在右边,系数全部缩减为0; 在中间的某值将取得最好的预测效果。为了定量地找到最佳参数值,还需要进行交叉验证。另外,要判断哪些变量对结果预测最具影响力,观察上图他们对应的系数大小就可以。
lasso
在增加如下约束时,普通的最小二乘法回归会得到与岭回归一样的公式 :
上式限定了所有回归系数的平方和要小于等于
。
与岭回归类似,另外一个缩减方法 lasso 也对回归系数做了限定 ,约束如下 :
对于这种缩减方法,它的一个特点就是, 在
足够小的时候,一些回归系数会被迫缩减为0 ,也就是抑制了一些特征的表达,因为这些特征对于我们预测并不会起到决定性的作用,我们把它丢弃,这样计算复杂度也就降低了。(具体可参照 https://baijiahao.baidu.com/s?id=1621054167310242353&wfr=spider&for=pc 查看L1正则和 L2 正则的优劣)
前向逐步回归
前向逐步回归算法可以得到与 lasso 差不多的效果,但是更加简单,他属于一种贪心算法,即每一步都尽可能减少误差。
伪代码:
1. 数据标准化 ,使其分布满足0均值 和 1 方差
2. FOR迭代:
3. SET lowestError = inf
4. FOR 遍历特征:
5. FOR 增加、减小特征权值:
6. 修改该特征的权值，得到新的权值向量W
7. 计算新特征向量的误差Error
8. IF Error < lowestERROR:
9. SET Wbest = W
10. 更新：SET W = Wbest#前向逐步线性回归
def stageWise(xArr, yArr, eps=0.01, numIt=100):
xMat = np.matrix(xArr)
yMat = np.matrix(yArr).T
m, n = xMat.shape
# 标准化
xMat = (xMat - xMat.mean(axis=0)) / xMat.var(axis=0)
yMat -= yMat.mean()
returnMat = np.zeros((numIt, n))
ws = np.zeros((n, 1))
wsMax = ws.copy()
# 对于每一轮
for i in range(numIt):
# print(ws.T)
# 设置当前最小误差为 inf
lowestError = np.inf
# 对于每一个特征
for j in range(n):
# 增大或者缩小
for sign in [-1, 1]:
# 改变一个系数
wsTest = ws.copy()
wsTest[j] += eps * sign
yTest = np.dot(xMat, wsTest)
# 计算当前新W的误差
rssE = rssError(yMat.A, yTest.A)
#
if rssE < lowestError:
# 更新最小误差
lowestError = rssE
wsMax = wsTest
ws = wsMax.copy()
returnMat[i, :] = ws.T
return returnMat>>> import regression
>>> xArr,yArr = regression.loadDataSet('abalone.txt')
>>> regression.stageWise(xArr,yArr,0.01,200)
array([[ 0. , 0. , 0. , ..., 0. , 0. , 0. ],
[ 0. , 0. , 0. , ..., 0. , 0. , 0. ],
[ 0. , 0. , 0. , ..., 0. , 0. , 0. ],
...,
[ 0.05, 0. , 0.09, ..., -0.64, 0. , 0.36],
[ 0.04, 0. , 0.09, ..., -0.64, 0. , 0.36],
[ 0.05, 0. , 0.09, ..., -0.64, 0. , 0.36]])值得注意的是有些参数是0 (w1 , w6 ) ,这表明他们对目标值的影响较小可以忽略不计,也就是说这些特征可能并不需要。
下面试着用更小的步长和更多的步数 :
>>> regression.stageWise(xArr,yArr,0.001,5000)
array([[ 0. , 0. , 0. , ..., 0. , 0. , 0. ],
[ 0. , 0. , 0. , ..., 0. , 0. , 0. ],
[ 0. , 0. , 0. , ..., 0. , 0. , 0. ],
...,
[ 0.043, -0.011, 0.12 , ..., -0.963, -0.105, 0.187],
[ 0.044, -0.011, 0.12 , ..., -0.963, -0.105, 0.187],
[ 0.043, -0.011, 0.12 , ..., -0.963, -0.105, 0.187]])接下来 把这些结果和最小二乘法进行比较
>>> xMat = np.matrix(xArr)
>>> yMat = np.matrix(yArr).T
>>> xMat = (xMat - xMat.mean(axis=0))/xMat.var(axis=0)
>>> yMat -=yMat.mean(axis=0 )
>>> weights = regression.standRegres(xMat,yMat.T)
>>> weights.T
matrix([[ 0.0430442 , -0.02274163, 0.13214087, 0.02075182, 2.22403814,
-0.99895312, -0.11725427, 0.16622915]])wMat = stageWise(xArr, yArr, 0.005, 1000)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(list(range(1000)), wMat)
plt.show()
使用0.005的eps 值,并经过1000次迭代
逐步线性回归的主要优点在于他可以帮助人们理解现有的模型并作出改进。当构建看一个模型后,可以运行该算法找出重要的特征,这样就可能及时停止对那些不重要特征的收集。
当应用缩减方法时,模型也就增加了偏差,同时却减少了模型的方差.
方差指的是模型之间的差异,而偏差指的是模型预测值和数据之间的差异。
实例 : 预测乐高玩具套装的价格
(1) 收集数据 从 html 网页上爬取
(2) 准备数据
(3) 分析数据 观察数据
(4) 训练数据 构建不同的模型 ,采用逐步线性回归和 线性回归
(5) 测试数据 使用交叉验证来测试不同的模型,分析哪个模型好
(6)使用数据
乐高数据集
链接：https://pan.baidu.com/s/1dLn8gqc61OjtEIxtaI1lAw
提取码：ypp8
收集数据
由于书籍中采用的 爬虫爬取,但现在网站进不去,我们采用的是解析 html 文件 ,从中获取我们想要的数据。
# lego.py
# -*-coding:utf-8 -*-
import numpy as np
from bs4 import BeautifulSoup
import random
# 购物信息的获取函数
def scrapePage(retX, retY, inFile, yr, numPce, origPrc):
'''
函数说明:从页面读取数据，生成retX和retY列表
:param retX: 数据X
:param retY: 数据Y
:param inFile: HTML文件
:param yr: 年份
:param numPce: 乐高部件数目
:param origPrc: 原价
:return:
'''
# 打开并读取HTML文件
with open(inFile, encoding='utf-8') as f:
html = f.read()
soup = BeautifulSoup(html)
i = 1
# 根据HTML页面结构进行解析
currentRow = soup.find_all('table', r = "%d" % i)
while(len(currentRow) != 0):
currentRow = soup.find_all('table', r = "%d" % i)
title = currentRow[0].find_all('a')[1].text
lwrTitle = title.lower()
# 查找是否有全新标签
if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1):
newFlag = 1.0
else:
newFlag = 0.0
# 查找是否已经标志出售，我们只收集已出售的数据
soldUnicde = currentRow[0].find_all('td')[3].find_all('span')
if len(soldUnicde) == 0:
print("商品 #%d 没有出售" % i)
else:
# 解析页面获取当前价格
soldPrice = currentRow[0].find_all('td')[4]
priceStr = soldPrice.text
priceStr = priceStr.replace('$','')
priceStr = priceStr.replace(',','')
if len(soldPrice) > 1:
priceStr = priceStr.replace('Free shipping', '')
sellingPrice = float(priceStr)
# 去掉不完整的套装价格
if sellingPrice > origPrc * 0.5:
print("%d\t%d\t%d\t%f\t%f" % (yr, numPce, newFlag, origPrc, sellingPrice))
retX.append([yr, numPce, newFlag, origPrc])
retY.append(sellingPrice)
i += 1
currentRow = soup.find_all('table', r = "%d" % i)
def setDataCollect(retX, retY):
"""
函数说明:依次读取六种乐高套装的数据，并生成数据矩阵
"""
scrapePage(retX, retY, './lego/lego8288.html', 2006, 800, 49.99) #2006年的乐高8288,部件数目800,原价49.99
scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99) #2002年的乐高10030,部件数目3096,原价269.99
scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99) #2007年的乐高10179,部件数目5195,原价499.99
scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99) #2007年的乐高10181,部件数目3428,原价199.99
scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99) #2008年的乐高10189,部件数目5922,原价299.99
scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99) #2009年的乐高10196,部件数目3263,原价249.99import lego
import regression
lgX = [] ; lgY = []
lego.setDataCollect(lgX,lgY)
print(lgX)
2006 800 0 49.990000 85.000000
2006 800 0 49.990000 102.500000
2006 800 0 49.990000 77.000000
商品 #4 没有出售
2006 800 0 49.990000 162.500000
2002 3096 0 269.990000 699.990000
2002 3096 0 269.990000 602.000000
2002 3096 0 269.990000 515.000000
2002 3096 0 269.990000 510.000000
2002 3096 0 269.990000 375.000000
2002 3096 1 269.990000 1050.000000
2002 3096 0 269.990000 740.000000
2002 3096 1 269.990000 759.000000
2002 3096 0 269.990000 730.000000
2002 3096 1 269.990000 750.000000
....... . ...........训练数据
# lego.py
lgX = [] ; lgY = []
setDataCollect(lgX,lgY)
# print("lgx : \n " , lgX)
# print("lgy : \n" ,lgY)
# print(np.shape(lgX))
lgX1 = np.matrix(np.ones((63,5)))
lgX1[:,1:5] = np.matrix(lgX)
print(lgX1[0])
ws = regression.standRegres(lgX1,lgY)
print(ws)
# 检查一下模型是否有效
print()
print(lgX1[0]*ws , " " , "实际: " , lgY[0])
print(lgX1[-1]*ws," " , "实际: " , lgY[-1])
print(lgX1[43]*ws," " , "实际: " , lgY[43])输出 :
[[1.000e+00 2.006e+03 8.000e+02 0.000e+00 4.999e+01]]
[[ 5.53199701e+04]
[-2.75928219e+01]
[-2.68392234e-02]
[-1.12208481e+01]
[ 2.57604055e+00]]
[[76.07418864]] 实际: 85.0
[[431.17797682]] 实际: 331.51
[[516.20733116]] 实际: 530.0可以看到模型还算有效.但是模型本身并不能令人满意。他对于数据拟合得很好 .
下面使用缩减法,岭回归进行一次实验 。
# regression
def crossValidation(xArr, yArr, numVal=10):
m = len(yArr) # 统计样本个数
indexList = list(range(m)) # 生成索引值列表
errorMat = np.zeros((numVal, 30)) # create error mat 30columns numVal rows
for i in range(numVal): # 交叉验证numVal次
trainX = [];
trainY = [] # 训练集
testX = [];
testY = [] # 测试集
random.shuffle(indexList) # 打乱次序
for j in range(m): # 划分数据集:90%训练集，10%测试集
if j < m * 0.9:
trainX.append(xArr[indexList[j]])
trainY.append(yArr[indexList[j]])
else:
testX.append(xArr[indexList[j]])
testY.append(yArr[indexList[j]])
wMat = ridgeTest(trainX, trainY) # 获得30个不同lambda下的岭回归系数
for k in range(30): # 遍历所有的岭回归系数
matTestX = np.mat(testX);
matTrainX = np.mat(trainX) # 测试集
meanTrain = np.mean(matTrainX, 0) # 测试集均值
varTrain = np.var(matTrainX, 0) # 测试集方差
matTestX = (matTestX - meanTrain) / varTrain # 测试集标准化
yEst = matTestX * np.mat(wMat[k, :]).T + np.mean(trainY) # 根据ws预测y值
errorMat[i, k] = rssError(yEst.T.A, np.array(testY)) # 统计误差
meanErrors = np.mean(errorMat, 0) # 计算每次交叉验证的平均误差
minMean = float(min(meanErrors)) # 找到最小误差
bestWeights = wMat[np.nonzero(meanErrors == minMean)] # 找到最佳回归系数
xMat = np.mat(xArr);
yMat = np.mat(yArr).T
meanX = np.mean(xMat, 0);
varX = np.var(xMat, 0)
unReg = bestWeights / varX # 数据经过标准化，因此需要还原
print('%f%+f*年份%+f*部件数量%+f*是否为全新%+f*原价' % (
(-1 * np.sum(np.multiply(meanX, unReg)) + np.mean(yMat)), unReg[0, 0], unReg[0, 1], unReg[0, 2], unReg[0, 3]))regression.crossValidation(lgX,lgY,10)
59856.168464-29.822292*年份+0.000961*部件数量+1.851745*是否为全新+2.011681*原价总结
与分类一样,回归也是预测目标值的过程,只不过不同之处在于, 前者预测连续型变量,而后者预测离散型变量。
在回归中求特征对应的最佳回归系数是最小化均方误差。给定 输入矩阵 X ,如果
可逆,回归法直接用就可以;否则
考虑岭回归 或者 lasso 。
好博客就要一起分享哦！分享海报
此处可发布评论
评论（0）展开评论
展开评论
您可能感兴趣的博客
