机器学习实战——岭回归、缩减法

一、基本知识

1、岭回归：从公式看，加入正则化项（2范数）。

回归系数的计算公式为：

问题引入：若给定数据集X，如果XTX的逆存在，可以使用常规的线性回归方法。但是，

（1）数据样本数比特征数少的情况，矩阵的逆不能直接计算；

（2）即使样本数多于特征数，若特征高度相关，XTX的逆依然无法计算。

此时，可以考虑岭回归。

另，岭回归是有偏估计回归方法，引入lamda来限制所有系数之和，通过引入该惩罚项（从需要最小化的平方误差函数可以看出），能够减少不重要的参数；它是缩减法的一种，增大模型的偏差，从而达到更好的估计效果，还可以用于挖掘大量数据的内容在规律，根据系数的大熊啊可以得到最重要的一个或几个特征。

书籍：《机器学习实战》中文版
IDE：PyCharm Edu 4.02
环境：Adaconda3  python3.6

from numpy import *
# 自适应数据加载函数
# 不必指定特征数目,
def loadDataSet(fileName):     #general function to parse tab -delimited floats
    numFeat = len(open(fileName).readline().split('\t'))-1   #get number of fields
    dataMat = [];labelMat = []
    with open(fileName) as fr:
        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       # 返回列表
# 岭回归
# 参数xMat和yMat是矩阵格式
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 singular, cannot do inverse")
        return
    ws = denom.I * (xMat.T * yMat)
    return ws
def ridgeTest(xArr,yArr):
    xMat = mat(xArr)
    yMat = mat(yArr).T
    # 数据标准化 可以使用regularize()函数
    yMean = mean(yMat,0)
    yMat = yMat-yMean
    xMeans = mean(xMat,0)
    xVar = var(xMat,0)
    xMat = (xMat-xMeans)/xVar
    numTestPts =30   #设置岭参数lam的取值范围
    wMat = zeros((numTestPts,shape(xMat)[1]))
    for i in range(numTestPts):
        ws = ridgeRegres(xMat,yMat,exp(i-10))
        wMat[i,:] = ws.T       #每一行是一个lam对应的回归系数
    return wMat
# 岭回归交叉验证
def rssError(yArr,yHatArr): #yArr and yHatArr both need to be arrays
    return ((yArr-yHatArr)**2).sum()
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):   #create training set based on first 90% of values in indexList
            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) #wMat每一行是一组回归系数；numTestPts=30每组系数的个数
        for k in range(30): #loop over all of the ridge estimates
            matTestX = mat(testX);matTrainX = mat(trainX)
            meanTrain = mean(matTrainX,0)
            varTrain = var(matTrainX,0)
            matTestX = (matTestX-meanTrain)/varTrain   #regularize test with training params
            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)]
    #can unregularize to get model
    #when we regularized we wrote Xreg = (x-meanX)/var(x)
    #we can now write in terms of x not Xreg:  x*w/var(x) - meanX/var(x) +meanY
    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))
# 鲍鱼数据集验证岭回归
abX,abY = loadDataSet('abalone.txt')
print(crossValidation(abX,abY))

运行结果：
the best model from Ridge Regression is:
 [[  0.05915723  -3.07364169  14.96227846  11.73052929   8.92101407
  -20.08279615  -9.68807653   9.33245255]]
with constant term:  3.14629495468
None

2、loss回归、前向逐步回归

loss回归（缩减法之一）：加入正则化项——绝对值形式。但计算复杂，因此使用计算简单的前向逐步回归，可以达到与loss回归相同的效果。

前向逐步回归属于一种贪心算法，每一步都尽可能减少误差。

# 前向逐步回归
def stageWise(xArr,yArr,eps=0.01,numIt=100):
    xMat = mat(xArr);yMat = mat(yArr).T
    # 数据标准化
    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();wsBest=ws.copy()
    for i in range(numIt):
        lowestError = inf
        for j in range(n):
            for sign in [-1,1]:
                wsTest = ws.copy()   #备份ws，以判断到底是增加还是减小系数
                wsTest[j] += eps*sign
                yTest = xMat*wsTest
                rssE = ((yMat.A-yTest.A)**2).sum()
                if rssE < lowestError:
                    lowestError = rssE
                    wsBest = wsTest
        ws = wsBest.copy()
        returnMat[i,:] = ws.T
    return returnMat
weightsStageWise = stageWise(abX,abY)
print(weightsStageWise)

说明：贪心算法在所有特征上运行两次for循环，分别计算增加或减少该特征对误差的影响。

即当一次外循环（i）结束后才改变系数ws，而内循环（j）的每个值记录在wsBest中但每次运行时仍是未改变的ws，在每一次外循环中（i）：

先判断第一个特征应该增加还是减少系数值，然后判断第二个特征，

此时第一个特征值并未改变，因此代码中每次循环前需要复制一份 wsTest = ws.copy()

关于引用、copy()的使用，参见：http://blog.csdn.net/ckzhb/article/details/78659449