机器学习实战【9】（岭回归和Lasso）

本博客记录《机器学习实战》（MachineLearningInAction）的学习过程，包括算法介绍和python实现。

岭回归

岭回归（Ridge Regression），是一种线性回归方法，在最小二乘法的基础上加入一个正则项，以解决样本数少于特征数以及特征存在线性相关（多重共线性）时产生的问题。

最小二乘法

多元线性模型y=Xβ+ϵ中，ϵ为残差项，残差项越小，模型也就越拟合数据。最小二乘法就是求解使得残差平方和最小的参数的方法。残差平方和RSS为：

Q(β)=(y−Xβ)T(y−Xβ)

把这个式子求导计算一波就可以推出β的最小二乘估计：

β^=(XTX)−1XTy

如果X存在线性相关关系，或者p小于n，这时逆矩阵就无法求解。

岭回归

在最小二乘估计的基础上，岭回归增加了一项，称为岭回归估计：

β^=(XTX+λI)−1XTy

此时求解的式子变成了：

β^=argmin{(∑i=1nyi−β0−∑j=1pXijβj)2+λ∑j=1pβ2j}

其中新增的最后一项称为惩罚函数。进一步，这个问题又可以转化成：

β^=argmin(∑i=1nyi−β0−∑j=1pXijβj)2s.t.∑j=1pβ2j<t

相当于在最小二乘法的基础上添加了对参数β的约束。岭参数λ接近0时，岭回归变为最小二乘法；岭参数λ变大时，参数β趋近于0。

LASSO

LASSO算法与岭回归类似，但是其惩罚函数是一阶的，即：

β^=argmin{(∑i=1nyi−β0−∑j=1pXijβj)2+λ∑j=1p|βj|}

相比岭回归，LASSO更易求解，并且更容易筛选变量。

python实现：

程序中计算了从小到大30个λ参数的结果，通过交叉验证的方法，取使得最终的平方误差最小的一组参数作为结果。需要注意的是数据集在处理之前要先进行标准化以去除公式中的常数项。

# 数据标准化
def regularize(xMat):  # regularize by columns
inMat = xMat.copy()
inMeans = np.mean(inMat, 0)  # calc mean then subtract it off
inVar = np.var(inMat, 0)  # calc variance of Xi then divide by it
inMat = (inMat - inMeans) / inVar
return inMat
# 岭回归
def ridgeRegres(xMat, yMat, lam=0.2):
xTx = xMat.T * xMat
denom = xTx + np.eye(np.shape(xMat)[1]) * lam
if np.linalg.det(denom) == 0.0:
print("This matrix is singular, cannot do inverse")
return
ws = denom.I * (xMat.T * yMat)
return ws
# 岭回归测试，返回一系列（30个）lambda下计算出的回归参数，从e-10到e20
def ridgeTest(xArr, yArr):
xMat = np.mat(xArr)
yMat = np.mat(yArr).T

# 归一化数据
yMean = np.mean(yMat, 0)
yMat = yMat - yMean  # to eliminate X0 take mean off of Y
xMat = regularize(xMat)

numTestPts = 30
wMat = np.zeros((numTestPts, np.shape(xMat)[1]))
for i in range(numTestPts):
ws = ridgeRegres(xMat, yMat, np.exp(i - 10))
wMat[i, :] = ws.T
return wMat
# 交叉验证
def crossValidation(xArr, yArr, numVal=10):
m = len(yArr)
indexList = range(m)
errorMat = np.zeros((numVal, 30))  # create error mat 30columns numVal rows
for i in range(numVal):
trainX = []
trainY = []
testX = []
testY = []
# 随机取90%的数据训练，剩下10%的数据测试
np.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]])

# 计算30个岭回归lambda参数下的平方误差
wMat = ridgeTest(trainX, trainY)  # get 30 weight vectors from ridge
for k in range(30):  # loop over all of the ridge estimates
matTestX = np.mat(testX)
matTrainX = np.mat(trainX)
meanTrain = np.mean(matTrainX, 0)
varTrain = np.var(matTrainX, 0)
matTestX = (matTestX - meanTrain) / varTrain  # regularize test with training params
yEst = matTestX * np.mat(wMat[k, :]).T + np.mean(trainY)  # test ridge results and store
errorMat[i, k] = rssError(yEst.T.A, np.array(testY))
print(errorMat[i, k])

# 取多次迭代的误差平均，使用平均误差最小的岭回归参数计算出的结果
meanErrors = np.mean(errorMat, 0)  # calc avg performance of the different ridge weight vectors
minMean = float(min(meanErrors))
bestWeights = wMat[np.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 = np.mat(xArr)
yMat = np.mat(yArr).T
meanX = np.mean(xMat, 0)
varX = np.var(xMat, 0)
unReg = bestWeights / varX
print("the best model from Ridge Regression is:\n", unReg)
print("with constant term: ", -1 * sum(np.multiply(meanX, unReg)) + np.mean(yMat))