机器学习--线性回归-蒲公英云

标准线性回归：

局部加权线性回归：

线性回归的一个问题是有可能出现欠拟合现象，因为它求的是具有小均方误差的无偏估计。显而易见，如果模型欠拟合将不能取得好的预测效果。所以有些方法允许在估计中引入一些偏差，从而降低预测的均方误差。

其中的一个方法是局部加权线性回归（Locally Weighted Linear Regression，LWLR）。在该方法中，我们给待预测点附近的每个点赋予一定的权重。与kNN一样，这种算法每次预测均需要事先选取出对应的数据子集。该算法解除回归系数W的形式如下：

其中W是一个矩阵，这个公式跟我们上面推导的公式的区别就在于W，它用来给每个店赋予权重。

LWLR使用”核”（与支持向量机中的核类似）来对附近的点赋予更高的权重。核的类型可以自由选择，最常用的核就是高斯核，高斯核对应的权重如下：

岭回归：

简单说来，岭回归就是在普通线性回归的基础上引入单位矩阵。回归系数的计算公式变形如下：

公式中，矩阵I是一个mxm的单位矩阵，加上一个λI从而使得矩阵非奇异，进而能对矩阵求逆。

岭回归最先用来处理特征数多于样本数的情况，现在也用于在估计中加入偏差，从而得到更好的估计。这里通过引入λ来限制了所有w之和，通过引入该惩罚项，能够减少不重要的参数，这个技术在统计学中也可以叫做缩减（shrinkage）。

缩减方法可以去掉不重要的参数，因此能更好地裂解数据。此外，与简单的线性回归相比，缩减法能够取得更好的预测效果。

逐步线性回归：

前向逐步线性回归算法属于一种贪心算法，即每一步都尽可能减少误差。我们计算回归系数，不再是通过公式计算，而是通过每次微调各个回归系数，然后计算预测误差。那个使误差最小的一组回归系数，就是我们需要的最佳回归系数。

#

小结：

回归系数ws是一个1*n(n是特征数)的matrix

为了更好计算，一般在特征集前加一列值为1的常数项

缩减方法（逐步线性回归或岭回归），就是将一些系数缩减成很小的值或者直接缩减为0。这样做，就增大了模型的偏差（减少了一些特征的权重），通过把一些特征的回归系数缩减到0，同时也就减少了模型的复杂度。消除了多余的特征之后，模型更容易理解，同时也降低了预测误差。但是当缩减过于严厉的时候，就会出现过拟合的现象，即用训练集预测结果很好，用测试集预测就糟糕很多。

岭回归 lasso 向前逐步回归都是缩减法，要标准化数据：
对于训练特征集(trainX)采用训练特征集减去均值并除以方差，(matTrainX-meanTrainX)/varTrainX
对于训练标签集(trainY)采用训练标签集减去均值，(trainY-meanTrainY)，
上两步直接把trainX trainY传入ridgeTest即可，ridgeTest会采用上述方法进行标准化，得到回归系数wMat；
对于测试特征集(testX)采用测试特征集减去训练特征集均值并除以训练特征集方差，matTestX=(matTestX-meanTrainX)/varTrainX
对于测试标签集，通过测试特征集乘以对应的回归系数mat(wMat[k,:])并加上训练标签集的均值得到 matTestX*mat(wMat[k,:]).T+mean(trainY)

一般线性回归都会在特征集前加一列为1的常数项，
必须这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。
如果需要处理的问题是动态变化的，那么可以适当加大上述常数项，来确保新的值获得更大的回归系数。
也可以在运算公式里加，这个例子就是在运算公式里加。

numpy用到的内容(个别的)

nonzero() nonzero函数是numpy中用于得到数组array中非零元素的位置（数组索引）的函数。

参考：https://blog.csdn.net/lilong117194/article/details/78283358

multiply() 矩阵对应元素相乘

参考：https://blog.csdn.net/u012609509/article/details/70230204/ 注意这篇文章是以np.array 做的例子

numpy array 对象—索引，以数组为索引

a=x[np.array([0,1]),np.array([3,2])] #获取(0,3)元素和(1,2)元素，两个数组构成坐标对

参考：https://www.jianshu.com/p/acfb224c150c?from=timeline

代码部分

from numpy import *
from bs4 import BeautifulSoup
import random 
'''
ws是一个1*n(n是特征数)的matrix
为了更好计算，一般在特征集前加一列值为1的常数项
岭回归 lasso 向前逐步回归都是缩减法，要标准化数据：
对于训练特征集(trainX)采用训练特征集减去均值并除以方差，(matTrainX-meanTrainX)/varTrainX
对于训练标签集(trainY)采用训练标签集减去均值，(trainY-meanTrainY)，
上两步直接把trainX trainY传入ridgeTest即可，ridgeTest会采用上述方法进行标准化，得到回归系数wMat；
对于测试特征集(testX)采用测试特征集减去训练特征集均值并除以训练特征集方差，matTestX=(matTestX-meanTrainX)/varTrainX
对于测试标签集，通过测试特征集乘以对应的回归系数mat(wMat[k,:])并加上训练标签集的均值得到 matTestX*mat(wMat[k,:]).T+mean(trainY)
'''
#格式化特征数据和标签
def loadDataSet(filename):
    #定义测试数据和标签数据
    dataMat=[]; labelMat=[]
    #定义特征值的长度，测试文件中前两列数据是特征值
    numFeat=len(open(filename).readline().split('\t'))-1
    #定义文件
    fr=open(filename)
    '''
    遍历文件内容，
    定义一个lineArr数组，以便保存每行数据，
    把文件的每一行放在curLine数组中，
    截取curLine中的前两个数据，添加到lineArr中，
    把lineArr添加到dataMat中，把curLine数组中最后一个数据添加到labelMat中
    '''
    for line in fr.readlines():
        lineArr=[]
        curLine=line.strip().split('\t')
        for i in range(numFeat):
            lineArr.append(float(curLine[i]))
        dataMat.append(lineArr)
        labelMat.append(float(curLine[-1]))
    return dataMat,labelMat        
'''
标准线性回归
'''    
def standRegression(xArr,yArr):
    #把x，y数组转换成矩阵，由于y是一维数组，所以需要把y转置，即mat.T
    xMat=mat(xArr);yMat=mat(yArr).T
    xTx=xMat.T*xMat
    #判断xTx矩阵是否可逆
    if linalg.det(xTx)==0.0:
        print('This matrix is singular, cannot do inverse')
        return
    ws=xTx.I*(xMat.T*yMat)
    return ws    
'''
局部加权线性回归 
k是计算权重的一个常数，k值为1时，系数与标准的线性回归一样，k值太小时会造成过拟合
'''
def lwlr(testpoint,xArr,yArr,k=1.0):
    #把特征数据和标签数据转成矩阵，标签数据yMat是行矩阵需要转置
    xMat=mat(xArr); yMat=mat(yArr).T
    #获得文件中数据行数
    m=shape(xArr)[0]
    #创建只含对角元素的权重矩阵weights
    weights=mat(eye((m)))
    '''
    通过LWLR公式计算权重矩阵weights,
    测试点与每个特征数据通过计算，得到这个点对应的权重
    '''
    for j in range(m):
        diffMat=testpoint-xMat[j,:]
        weights[j,j]=exp(diffMat*diffMat.T/(-2.0*k**2))
    #还是利用平方误差计算出局部加权后的回归系数ws
    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
#lwlr测试
def lwlrTest(testArr, xArr,yArr,k=1.0):
    #获得测试数据的行数m
    m=shape(testArr)[0]
    #创建存放预测值的一维数组
    yHat=zeros(m)
    #遍历测试数据通过lwlr方法计算出与其对应的预测值
    for i in range(m):
        yHat[i]=lwlr(testArr[i],xArr,yArr,k)
    return yHat            
#通过平方差计算预测值与真实值的误差
def rssError(yArr,yHatArr):
    error=((yArr-yHatArr)**2).sum()    
    return error
'''
岭回归 
lam是引入的惩罚项，
lam越小，惩罚项越小，计算出的回归系数与标准的线性回归系数越接近
lam越大，惩罚项越大，计算出的回归系数越趋于0
矩阵I是一个n*n(特征数*特征数)的单位矩阵
'''
def ridgeRegress(xMat,yMat,lam=0.2):
    xTx=xMat.T*xMat
    #通过numpy.eye创建一个单位矩阵
    denom=xTx+eye(shape(xMat)[1])*lam
    if linalg.det(xTx)==0.0:
        print('This matrix is singular, cannot do inverse')
        return
    ws=denom.I*(xMat.T*yMat)
    return ws        
'''
ridgeTest是对岭回归做测试，分别用30个lam值，看其对回归系数的影响
'''
def ridgeTest(xArr,yArr):
    #把特征数据和标签数据转成矩阵，标签数据yMat是行矩阵需要转置
    xMat=mat(xArr); yMat=mat(yArr).T
    #抵消xMat标准化数据后，对yMat的影响
    yMean=mean(yMat,0)
    yMat=yMat-yMean
    #对xMat标准化，可以使用下面的regularize方法
    xMeans=mean(xMat,0)
    xVar=var(xMat,0)
    xMat=(xMat-xMeans)/xVar
    #定义要使用30个lam值做测试
    numTestPts=30
    #创建存放不同lam值，计算出来的回归系数ws的矩阵
    wMat=zeros((numTestPts,shape(xMat)[1]))
    #传入30个lam值，计算出每个lam对应的回归系数w，并把其填入到wMat中
    for i in range(numTestPts):
        ws=ridgeRegress(xMat,yMat,exp(i-10))
        wMat[i,:]=ws.T
    return wMat    
'''
标准化数据
所有特征减去各自的均值并除以方差
mean(mat,0)压缩行，求每列的均值
var(mat,0)压缩行，求每列的方差
'''
def regularize(xMat):
    InMat=xMat.copy()
    xMeans=mean(InMat,0)
    xVar=var(InMat,0)
    InMat=(InMat-xMeans)/xVar
    return InMat
'''
使用岭回归算法
'''
# def ridgeAbalone(testPoint, xArr, yArr,lam=0.2):
#     xMat=mat(xArr); yMat=mat(yArr).T; testMat=mat(testPoint)
#     #标准化xArr,yArr,testPoint
#     yMean=mean(yMat,0)
#     yMat=yMat-yMean
#     xMat=regularize(xMat)
#     testMat=regularize(testPoint)
#     m=shape(testMat)[0]
#     yHat=zeros(m)
#     #使用岭回归计算出回归系数ws
#     ws=ridgeRegress(xMat,yMat,lam)
#     #计算预测值yHat，因为之前是用标准化后的数据，所以要加上均值
#     yHat=testMat*ws+mean(yArr)
#     return yHat
'''
向前逐步回归
数据标准化，使其分布满足0均值和单位方差
每次迭代过程中：
    设置当前最小误差lowestError为正无穷
    对每个特称：
        增大或减小：
            改变一个系数得到一个新的w
            计算新的wTest下的误差
            如果当前误差Error小于当前最小误差lowestError:那么将wBest设置为wTest
        将w设置为新的wBest
'''    
def stageWise(xArr,yArr,eps=0.01,numIt=100):
    #同样需要对数据做标准化
    xMat=mat(xArr);yMat=mat(yArr).T
    yMean=mean(yArr,0)
    yMat=yMat-yMean
    xMat=regularize(xMat)
    #获取特征数据的行，列 m,n
    m,n=shape(xMat)
    #定义要返回的回归系数矩阵，里面存放着每一次迭代的系数
    returnMat=zeros((numIt,n))
    #定义回归系数ws 用于测试(增大或减小)的wsTest 每次迭代得到的wsMax
    ws=zeros((n,1)); wsTest=ws.copy(); wsMax=ws.copy()
    '''
    通过增大或减小ws，得到预测值yTest,并计算出误差值
    对比当前误差rssE与当前最小误差lowestError
    如果当前误差rssE小于当前最小误差lowestError，则将wMax设置为wTest
    '''
    #进行numIt次迭代
    for i in range(numIt):
        print(ws.T)
        #设当前最小误差是正无穷大
        lowestError=inf
        #遍历每个特征计算该次迭代的回归系数
        for j in range(n):
            #对wsTest进行加减比较
            for sign in [-1,1]:
                #拷贝当前ws值赋给wsTest
                wsTest=ws.copy()
                wsTest[j]+=eps*sign
                #通过wsTest计算出对应的预测值yTest
                yTest=xMat*wsTest
                #计算yTest的误差rssE
                rssE=rssError(yMat.A,yTest.A)
                #如果当前误差rssE小于当前最小误差lowestError，则将wMax设置为wTest
                if rssE<lowestError:
                    lowestError=rssE
                    wsMax=wsTest
        #把拷贝wsMax的值赋给ws            
        ws=wsMax.copy()
        #把该次迭代得到的ws.T(因为得到的ws是1*n的矩阵，所以要转置)添加到回归系数矩阵returnMat
        returnMat[i,:]=ws.T
    return returnMat                
'''
函数说明:从页面读取数据，生成retX和retY列表
    Parameters:
       retX - 数据X
       retY - 数据Y
       inFile - HTML文件
       yr - 年份
       numPce - 乐高部件数目
       origPrc - 原价
   Returns:
       无
'''
def scrapePage(retX, retY, inFile, yr, numPce, origPrc):    
    # 打开并读取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):
    #2006年的乐高8288,部件数目800,原价49.99
    scrapePage(retX, retY, './lego/lego8288.html', 2006, 800, 49.99)
    #2002年的乐高10030,部件数目3096,原价269.99
    scrapePage(retX, retY, './lego/lego10030.html', 2002, 3096, 269.99)
    #2007年的乐高10179,部件数目5195,原价499.99
    scrapePage(retX, retY, './lego/lego10179.html', 2007, 5195, 499.99)
    #2007年的乐高10181,部件数目3428,原价199.99
    scrapePage(retX, retY, './lego/lego10181.html', 2007, 3428, 199.99)
    #2008年的乐高10189,部件数目5922,原价299.99
    scrapePage(retX, retY, './lego/lego10189.html', 2008, 5922, 299.99)
    #2009年的乐高10196,部件数目3263,原价249.99
    scrapePage(retX, retY, './lego/lego10196.html', 2009, 3263, 249.99)        
'''
必须这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。
如果需要处理的问题是动态变化的，那么可以适当加大上述常数项，来确保新的值获得更大的回归系数。
'''
'''
交叉验证测试岭回归
Parameters:
        xArr - x数据集
        yArr - y数据集
        numVal - 交叉验证次数
'''
def crossValidation(xArr,yArr,numVal=10):
    #统计样本个数
    m=len(yArr)
    #生成索引值列表
    indexList=list(range(m))
    #创建存放误差的矩阵，行数是要交叉验证的次数numVal，列数是岭回归所用lam的次数，岭回归测试方法中(ridgeTest)用了30个lam值，所以这里也是30
    errorMat=zeros((numVal,30))
    '''
    进行numVal次交叉验证，
    通过random.shuffle对索引列表indexList中的元素随机排序，
    样本集的90%作为训练集，10%作为测试集
    把随机筛选出来的测试集传入岭回归测试函数ridgeTest中，得到存放30组回归系数array wMat
    在每次交叉验证中，
    遍历30组回归系数，分别得到30个测试标签，
    通过rssError函数计算误差，并把他们添加到该次交叉验证对应的误差矩阵errorMat里
    '''
    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)
            meanTrainX=mean(matTrainX,0)
            varTrainX=var(matTrainX,0)
            matTestX=(matTestX-meanTrainX)/varTrainX
            yEst=matTestX*mat(wMat[k,:]).T+mean(trainY)
            errorMat[i,k]=rssError(array(testY),yEst.T.A)
    #计算误差平均值        
    meanErrors=mean(errorMat,0)
    #选取误差平均值中最小的误差
    minMean=float(min(meanErrors))
    '''
    通过nonzero返回最小误差minMean在平均误差meanErrors中的索引index，
    这个索引index是误差矩阵errorMat的列索引，比如index=1，也就是errorMat第二列，
    该列对应的lam值是一样的，所以最佳回归系数就是这个lam值在wMat里对应的那行，
    即wMat[index]，也就是 wMat[nonzero(meanErrors == minMean)]
    '''
    bestWeights = wMat[nonzero(meanErrors == minMean)]
    xMat = mat(xArr); yMat = mat(yArr).T
    meanX = mean(xMat,0); varX = var(xMat,0)
    #数据经过标准化，因此需要还原，unReg就是还原后的最佳回归系数
    unReg = bestWeights/varX
    print('the best model from Ridge Regression is:\n',unReg)
    '''
    -1*sum(multiply(meanX,unReg))+mean(yMat)) 相当于在特征集前加一列值为1的常数项,
    一般线性回归都会在特征集前加一列为1的常数项，
    必须这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。
    如果需要处理的问题是动态变化的，那么可以适当加大上述常数项，来确保新的值获得更大的回归系数。
    也可以在运算公式里加，这个例子就是在运算公式里加。
    multiply()是计算矩阵对应元素相乘，所以这两个矩阵的维度是一样的
    '''
    print('with constant term:', -1*sum(multiply(meanX,unReg))+mean(yMat))
    print('%f%+f*年份%+f*部件数量%+f*是否为全新%+f*原价' % ((-1 *sum(multiply(meanX,unReg)) + mean(yMat)),\
            unReg[0,0], unReg[0,1], unReg[0,2], unReg[0,3]))

分析数据代码—用wlr(局部加权线性回归)

import matplotlib.pyplot as plt
from numpy import *
import regression
'''
调用regression中的方法分别得到k值为1.0，0.03，0.01的预测值
'''
xArr,yArr=regression.loadDataSet('ex0.txt')
yHat1=regression.lwlrTest(xArr,xArr,yArr,1.0)
yHat2=regression.lwlrTest(xArr,xArr,yArr,0.03)
yHat3=regression.lwlrTest(xArr,xArr,yArr,0.01)
#把x，y数组转换成矩阵
xMat=mat(xArr); yMat=mat(yArr)
#通过argsort方法，以xMat第二行正序排列，得到其index值srtInd
srtInd=xMat[:,1].argsort(0)
#因为xMat[srtInd]是三维数组，通过截取获得排列后的xSort，这里也可以用 xSort=xMat.sort(0)
xSort=xMat[srtInd][:,0,:]
'''
用plt创建3行1列的图板，
分别显示k值为1.0，0.03，0.01时的图例
'''
fig=plt.figure()
ax1=fig.add_subplot(311)
ax1.scatter(xMat[:,1].flatten().A[0],yMat.flatten().A[0],s=2,c='red')
ax1.plot(xSort[:,1].flatten().A[0],yHat1[srtInd].T.flatten())
ax2=fig.add_subplot(312)
ax2.scatter(xMat[:,1].flatten().A[0],yMat.flatten().A[0],s=2,c='red')
ax2.plot(xSort[:,1].flatten().A[0],yHat2[srtInd].T.flatten())
ax3=fig.add_subplot(313)
ax3.scatter(xMat[:,1].flatten().A[0],yMat.flatten().A[0],s=2,c='red')
ax3.plot(xSort[:,1].flatten().A[0],yHat3[srtInd].T.flatten())
plt.show()

使用岭回归，分析lam值对回归系数的影响代码

import matplotlib.pyplot as plt
from numpy import *
import regression
abX,abY = regression.loadDataSet("abalone.txt")
ridgeWeights = regression.ridgeTest(abX,abY)
fig = plt.figure()
ax = fig.add_subplot(111)
axisX = [(exp(i - 10)) for i in range(30)]
#log(λ)作为横坐标
ax.plot(axisX, ridgeWeights)
#ax.plot(ridgeWeights)
plt.show()