目录
本节涉及的相关代码和数据
预测数值型数据：回归
本章内容：
① 线性回归
② 局部加权线性回归
③ 岭回归和逐步线性回归
④ 预测鲍鱼年龄和玩具售价
1 利用线性回归找到最佳拟合直线
线性回归的相关特点主要是
优点：结果易于理解，计算上不复杂
缺点：对非线性数据不友好
适用数据类型：数值型和标称型
回归的一般方法：
① 收集数据：采用任意方法收集数据
② 准备数据：回归需要数值型数据，标称型数据将会被转为二值型数据
③ 分析数据：绘出数据的可视化二维图将有助于对数据做出理解和分析，在采用缩减法求得新回归系数之后，可以将新拟合曲线绘在图上作为对比
④ 训练算法：找到回归系数
⑤ 测试算法：使用R2或者预测值和数据的拟合度，来分析模型的结果
⑥ 使用算法：使用回归，可以在给定输入的时候预测出一个数值，这是对分类方法的提升，因为这样可以预测连续性数据而不仅仅是离散的数据标签
采用最小平方误差
-
# 标准回归函数和数据导入函数
-
from numpy
import *
-
# 加载数据
-
def
loadDataSet(
fileName):
-
# print(fileName)
-
numFeat=
len(
open(fileName).readline().split(
'\t'))-
1
-
dataMat=[]
-
labelMat=[]
-
fr=
open(fileName)
-
for line
in fr.readlines():
-
lineArr=[]
-
curLine=line.strip().split(
'\t')
-
for i
in
range(numFeat):
-
lineArr.append(
float(curLine[i]))
-
# print(lineArr)
-
# print(lineArr)
-
dataMat.append(lineArr)
-
labelMat.append(
float(curLine[-
1]))
-
# print(dataMat)
-
return dataMat,labelMat
-
-
# 线性回归函数
-
def
standRegres(
xArr,yArr):
-
# 数组转为矩阵
-
xMat=mat(xArr)
-
yMat=mat(yArr).T
-
xTx=xMat.T*xMat
-
# linalg.det()函数可以直接用来计算行列式
-
if linalg.det(xTx)==
0.0:
-
print(
"this matrix is singular,cannot do inverse")
-
return
-
# matrix.I获得与给定矩阵相同大小的乘法逆
-
ws=xTx.I*(xMat.T*yMat)
-
# 返回权重
-
return ws
调用上述函数：
-
xArr,yArr=loadDataSet(
'ex0.txt')
-
# xArr
-
# 回归系数
-
ws=standRegres(xArr,yArr)
-
ws
输出得到的结果为：
将得到的权重系数，画图表示为：
-
# 绘图
-
def
plotData(
xArr,yArr,ws):
-
import matplotlib.pyplot
as plt
-
xMat = mat(xArr)
-
yMat = mat(yArr)
-
figure = plt.figure()
-
ax = figure.add_subplot(
111)
-
# 取第二个特征绘图
-
# flatten()函数转化成一维数组
-
# matrix.A属性返回矩阵变成的数组，和getA()方法一样
-
# 绘制点
-
ax.scatter(xMat[:,
1].flatten().A[
0], yMat.T[:,
0].flatten().A[
0])
-
# 返回给定数据的数组形式的拷贝
-
xCopy = xMat.copy()
-
# 升序排序
-
xCopy.sort(
0)
-
print (ws.shape)
-
yHat = xCopy * ws
# yHat 表示拟合直线的纵坐标，用回归系数求出
-
ax.plot(xCopy[:,
1], yHat, c =
'green')
-
plt.show()
-
-
-
plotData(xArr,yArr,ws)
得到的输出图像为：
可以看到该拟合，虽然在数据计算上结果较好，但可能单纯的线性拟合并不能够得到很好的结果，因为原始数据很明显有一个波浪形的分布，而不单单是线性分布。
计算该拟合结果的相关系数为：
-
# 计算相关系数
-
def
calcCorrel(
xArr,yArr,ws):
-
xMat = mat(xArr)
-
yMat = mat(yArr)
-
yHat=xMat*ws
-
# 转置yHat向量，以保证都是行向量
-
# corrcoef函数计算皮尔逊相关系数
-
correl=corrcoef(yHat.T,yMat)
-
return correl
-
-
calcCorrel(xArr,yArr,ws)
得到的结果为：
2 局部加权线性回归
线性回归的一个问题就是有可能出现欠拟合的现象，因为他求的是具有最小均方误差的无偏估计。显而易见，如果模型欠拟合将不能取得最好的预测效果。所以有些方法允许在估计中引入一些偏差，从而降低预测的均方误差。
其中一个方法就是局部加权线性回归（LWLR）。我们给待预测点附近的每个点赋予一定的权重，然后在这个基础上基于最小均方差来进行普通的回归。与KNN一样，这种算法每次预测均需事先选取处对应的数据子集。LWLR使用‘核’来对附近的点赋予更高的权重。核的类型可以自由选择，常用的有高斯核。
-
def
lwlr(
testPoint,xArr,yArr,k=1.0):
-
xMat=mat(xArr)
-
yMat=mat(yArr).T
-
m=shape(xMat)[
0]
-
# eye()函数返回一个二维数组，对角线为1，其余地方为0
-
weights=mat(eye((m)))
-
for j
in
range(m):
-
diffMat=testPoint-xMat[j,:]
-
weights[j,j]=exp(diffMat*diffMat.T/(-
2.0*k**
2))
-
xTx=xMat.T*(weights*xMat)
-
if linalg.det(xTx)==
0.0:
-
print(
"this matrix is singular,cannot do inverse")
-
return
-
ws=xTx.I*(xMat.T*(weights*yMat))
-
-
return testPoint*ws
-
-
def
lwlrTest(
testArr,xArr,yArr,k=1.0):
-
m=shape(testArr)[
0]
-
yHat=zeros(m)
-
for i
in
range(m):
-
yHat[i]=lwlr(testArr[i],xArr,yArr,k)
-
return yHat
调用上述函数：
lwlr(xArr[0],xArr,yArr)得到的输出结果为：
画出该拟合的结果曲线为：
-
def
plotDatalwlr(
xArr,yArr,k=0.01):
-
yHat=lwlrTest(xArr,xArr,yArr,k)
-
xMat=mat(xArr)
-
# 对样本x升序排序并返回索引
-
srtInd=xMat[:,
1].argsort(
0)
-
xSort=xMat[srtInd][:,
0,:]
-
import matplotlib.pyplot
as plt
-
fig=plt.figure()
-
ax=fig.add_subplot(
111)
-
ax.plot(xSort[:,
1],yHat[srtInd])
-
ax.scatter(xMat[:,
1].flatten().A[
0],mat(yArr).T.flatten().A[
0],s=
2,c=
'red')
-
plt.show()
-
xArr,yArr=loadDataSet(
'ex0.txt')
-
plotDatalwlr(xArr,yArr)
得到的拟合图像为：
局部加权线性回归也存在一个问题，即增加了计算量，因为他对每个点做预测时都必须使用整个数据集。如果避免这些计算将可以减少程序运行时间，从而减缓因计算量增加带来的问题。
3 示例：预测鲍鱼的年龄
-
# 计算误差
-
def
rssError(
yArr,yHatArr):
-
# **表示幂运算，前面为底数，后面为指数
-
return ((yArr-yHatArr)**
2).
sum()
-
-
abX,abY=loadDataSet(
'abalone.txt')
-
-
plotDatalwlr(abX[
0:
99],abY[
0:
99],
0.1)
-
plotDatalwlr(abX[
0:
99],abY[
0:
99],
1)
-
plotDatalwlr(abX[
0:
99],abY[
0:
99],
10)
-
-
# 训练集
-
yHat01_=lwlrTest(abX[
0:
99],abX[
0:
99],abY[
0:
99],
0.1)
-
yHat1_=lwlrTest(abX[
0:
99],abX[
0:
99],abY[
0:
99],
1)
-
yHat10_=lwlrTest(abX[
0:
99],abX[
0:
99],abY[
0:
99],
10)
-
-
err01_=rssError(abY[
0:
99],yHat01_.T)
-
err1_=rssError(abY[
0:
99],yHat1_.T)
-
err10_=rssError(abY[
0:
99],yHat10_.T)
-
-
print(err01_,err1_,err10_)
-
-
# 测试集
-
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)
-
-
err01=rssError(abY[
100:
199],yHat01.T)
-
err1=rssError(abY[
100:
199],yHat1.T)
-
err10=rssError(abY[
100:
199],yHat10.T)
-
-
print(err01,err1,err10)
得到的输出结果为：
可以看到，使用较小的核可以获得较低的误差，但是使用最小的核可能会造成过拟合，对新数据不一定能够达到最好的效果
接下来再和简单的线性回归做个比较
-
ws=standRegres(abX[
0:
99],abY[
0:
99])
-
yHat=mat(abX[
100:
199])*ws
-
rssError(abY[
100:
199],yHat.T.A)
得到结果为：
简单线性回归达到了与局部加权线性回归类似的效果，这也表明：必须在未知数据上比较效果才能够选取到最佳模型
4 缩减系数来“理解”数据
如果数据的特征比样本点还多，也就是说输入数据的矩阵X不是满秩矩阵。非满秩矩阵在求逆时会出现问题，因此，引入了“岭回归”的概念，这就是第一种缩减方法，接着是lasso方法，效果很好但复杂，最后介绍了第二种缩减方法，称为前向逐步回归，可以得到与lasso差不多的效果，也更容易实现。
4.1 岭回归
岭回归就是在矩阵XTX上加入一个λI从而使矩阵非奇异，进而就能对矩阵求逆，其中矩阵I时一个单位矩阵，λ由用户自定义。
缩减方法可以去掉不重要的参数，因此能够更好的理解数据，此外，与简单的线性回归相比，缩减法能够取得更好的效果。
-
# 用于计算回归系数
-
# 实现给定系数下的岭回归
-
def
ridgeRegres(
xMat,yMat,lam=0.2):
-
xTx=xMat.T*xMat
-
denom=xTx+eye(shape(xMat)[
1])*lam
-
if linalg.det(denom)==
0.0:
-
print(
"this matrix is singlar,cannot do inverse")
-
return
-
ws=denom.I*(xMat.T*yMat)
-
return ws
-
-
# 在一组λ上测试结果
-
def
ridgeTest(
xArr,yArr):
-
xMat=mat(xArr)
-
yMat=mat(yArr).T
-
yMean=mean(yMat,
0)
-
yMat=yMat-yMean
-
xMeans=mean(xMat,
0)
-
# var()函数求样本方差的无偏估计值，如果参数是1，就是有偏估计
-
xVar=var(xMat,
0)
-
xMat=(xMat-xMeans)/xVar
-
numTestPts=
30
-
wMat=zeros((numTestPts,shape(xMat)[
1]))
-
for i
in
range(numTestPts):
-
ws=ridgeRegres(xMat,yMat,exp(i-
10))
-
wMat[i,:]=ws.T
-
return wMat
-
-
abX,abY=loadDataSet(
'abalone.txt')
-
ridgeWeights=ridgeTest(abX,abY)
-
-
# 绘制出回归系数与log（λ）的关系
-
import matplotlib.pyplot
as plt
-
fig=plt.figure()
-
ax=fig.add_subplot(
111)
-
ax.plot(ridgeWeights)
-
plt.show()
得到的结果为：
4.2 lasso
lasso对回归系数做出了限定，使用绝对值取代了平方和，在λ足够小的时候，一些系数会因此被迫缩减到0，这个特性将有助于理解数据
4.3 前向逐步回归
属于一种贪心的算法，每一步都尽可能的减小误差。一开始所有的权重都设置为1，然后每一步所作的决策是对某个权重增加或减少一个很小的值
伪代码：
数据标准化，使其分布满足0均值和单位方差
在每轮迭代中：
设置当前最小误差lowestError为正无穷
对每个特征：
增大或减小：
改变一个系数得到一个新的w
计算新w下的误差
如果误差Error小于当前最小误差lowestError：
设置Wbest等于当前w
将w设置为新的Wbest
-
# 实现矩阵归一化
-
def
regularize(
xMat):
-
inMat = xMat.copy()
-
# 得到平均数，压缩行，对每一列求平均值
-
inMeans = mean(inMat,
0)
-
-
inVar = var(inMat,
0)
-
inMat = (inMat - inMeans)/inVar
-
return inMat
-
-
# 前向逐步线性回归
-
def
stageWise(
xArr,yArr,eps=0.01,numIt=100):
-
xMat=mat(xArr)
-
yMat=mat(yArr).T
-
# mean()函数求平均值
-
yMean=mean(yMat,
0)
-
yMat=yMat-yMean
-
xMat=regularize(xMat)
-
m,n=shape(xMat)
-
returnMat=zeros((numIt,n))
-
ws=zeros((n,
1))
-
wsTest=ws.copy()
-
wsMax=ws.copy()
-
for i
in
range(numIt):
-
print(ws.T)
-
lowestError=inf
-
for j
in
range(n):
-
for sign
in [-
1,
1]:
-
wsTest=ws.copy()
-
wsTest[j]+=eps*sign
-
yTest=xMat*wsTest
-
rssE=rssError(yMat.A,yTest.A)
-
if rssE<lowestError:
-
lowestError=rssE
-
wsMax=wsTest
-
ws=wsMax.copy()
-
returnMat[i,:]=ws.T
-
return returnMat
调用函数：
-
xArr,yArr=loadDataSet(
'abalone.txt')
-
stageWise(xArr,yArr,
0.01,
200)
得到的输出结果为：
与最小二乘法的结果进行对比得到
-
# 最小二乘法的结果
-
xMat=mat(xArr)
-
yMat=mat(yArr).T
-
xMat=regularize(xMat)
-
yM=mean(yMat,
0)
-
yMat=yMat-yM
-
weights=standRegres(xMat,yMat.T)
-
weights.T
结果为：
5 权衡偏差与方差
任何时候，一旦发现模型和测量值直接存在差异，就说明出现了误差。当考虑模型中的“噪声”或者说误差时，必须考虑其来源。可能会因为对复杂的模型进行简化，这将导致在模型和测量值之间出现“噪声”或者“误差”，若无法理解数据的真实生成过程，也会导致差异的发生。
一般认为，误差一般由三个部分组成：偏差、测量误差和随机噪声
方差是可以度量的，如果从总体数据中取一个随机样本集并用线性模型拟合，将会得到一组回归系数，同理，再取出另一组随机样本集并拟合，将会得到另一组回归系数。这些系数之间的差异大小也就是模型反差大小的反映。
6 示例：预测乐高玩具套装的价格
①收集数据：用google Shopping的APi收集数据
②准备数据：从返回的json数据中抽取价格
③分析数据：可视化并观察数据
④训练算法：构建不同的模型，采用逐步线性回归和直线的线性回归模型
⑤测试算法：使用交叉验证来测试不同的模型，分析哪个效果更好
⑥使用算法：这次联系的目的就是生成数据模型
6.1 收集数据：使用google购物的API
-
from bs4
import BeautifulSoup
-
def
scrapePage(
retX, retY, inFile, yr, numPce, origPrc):
-
# 打开并读取HTML文件
-
with
open(inFile, encoding=
'utf-8')
as f:
-
html = f.read()
-
# beautiful函数将复杂HTML文档转换成一个复杂的树形结构
-
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.99
调用上面函数：
-
lgX=[]
-
lgY=[]
-
setDataCollect(lgX,lgY)
得到结果：
6.2 训练算法：建立模型
为上面收集到的数据建立模型。构建出来的模型可以对售价做出预测
-
# 创建一个全为1 的矩阵
-
lgX1=mat(ones((
63,
5)))
-
print(lgX1[
0])
-
# 将原数据矩阵lgX复制到新数据矩阵lgX1的第一列到第五列
-
lgX1[:,
1:
5]=mat(lgX)
-
# 确认数据复制的正确性
-
print(lgX[
0])
-
print(lgX1[
0])
-
ws=standRegres(lgX1,lgY)
-
ws
得到的输出结果为：
交叉测试验证岭回归
-
# 交叉测试验证岭回归
-
# 交叉验证的次数
-
def
crossValidation(
xArr,yArr,numVal=10):
-
m=
len(yArr)
-
indexList=
list(
range(m))
-
errorMat=zeros((numVal,
30))
-
for i
in
range(numVal):
-
trainX=[]
-
trainY=[]
-
testX=[]
-
testY=[]
-
random.shuffle(indexList)
-
# 将数据分为测试集和训练集
-
for j
in
range(m):
-
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)
-
for k
in
range(
30):
-
matTestX=mat(testX)
-
matTrainX=mat(trainX)
-
meanTrain=mean(matTrainX)
-
varTrain=var(matTrainX,
0)
-
matTeatX=(matTestX-meanTrain)/varTrain
-
yEst=matTestX*mat(wMat[k,:]).T+mean(trainY)
-
errorMat[i,k]=rssError(yEst.T.A,array(testY))
-
meanErrors=mean(errorMat,
0)
-
minMean=
float(
min(meanErrors))
-
bestWeights=wMat[nonzero(meanErrors==minMean)]
-
xMat=mat(xArr)
-
yMat=mat(yArr).T
-
meanX=mean(xMat,
0)
-
varX=var(xMat,
0)
-
unReg=bestWeights/varX
-
# 岭回归的最佳模型是
-
print(
"the best model from Ridge Regression is :\n",unReg)
-
# 常数项
-
print(
"with constant term:",-
1*
sum(multiply(meanX,unReg))+mean(yMat))
-
-
crossValidation(lgX,lgY,
10)
得到的输出结果为：
这些系数是经过不同程度的缩减得到的
7 本章结束
转载：https://blog.csdn.net/weixin_46182244/article/details/127908310
查看评论