机器学习算法之回归详解
导语
回归：从一组数据出发，确定某些变量之间的定量关系式；即建立数学模型并估计未知参数。
回归的目的是预测数值型的目标值，它的目标是接受连续数据，寻找最适合数据的方程，并能够对特定值进行预测。这个方程称为回归方程，而求回归方程显然就是求该方程的回归系数，求这些回归系数的过程就是回归。
我们先来看看回归跟分类的区别：(引用 - 走刀口的回答 - 知乎
https://www.zhihu.com/question/21329754/answer/17901883)
分类和回归的区别在于输出变量的类型。
定量输出称为回归，或者说是连续变量预测；
定性输出称为分类，或者说是离散变量预测。
举个例子：
预测明天的气温是多少度，这是一个回归任务；
预测明天是阴、晴还是雨，就是一个分类任务。
普通线性回归
我们做线性回归是为了得到一个最优回归系数向量w使得当我们给定一个x能够通过y=xw预测y的值。假定输入数据存放在矩阵 X 中，而回归系数存在在向量 w 中。那么对于给定的数据X1，预测结果将会通过Y1=XT1w给出。
那么怎样的w才是最优的呢？在标准线性回归中我们需要找到是误差最小的w, 即预测的y值与真实的y值之间的差值，为了避免简单累加造成的正负差值相互抵消，这里采用了平方误差:
平方误差
用矩阵可表示为：(y-Xw)T(y-Xw)，)因为要求函数的极小值，对w求导可得:
对w求导
使其等于0，便可求出W的最优解:
求出w
上述求w最优解的方法也叫做最小二乘法。
需要注意的是，上述公式中包含(XTX)-1，也就是需要对矩阵求逆，因此这个方程只有在逆矩阵存在的时候有用，所以我们写代码时必须需要事先确定矩阵是否可逆。
此外，我们知道线性回归的方程的一般形式为：y=wx+b；即存在一定的偏移量b，于是，我们可以将回归系数和特征向量均增加一个维度，将函数的偏移值b也算作回归系数的一部分，占据回归系数向量中的一个维度，比如w=(w1,w2,...,wN,b)。相应地，将每条数据特征向量的第一个维度设置为1.0，即所有特征向量的第一个维度值均为1.0,这样，最后得到的回归函数形式为y=w[0]+w[1]x1+...+w[N]xN.
Python代码实现
# 普通线性回归
def standRegress(xArr, yArr):
# 用mat函数转换为矩阵之后可以才进行一些线性代数的操作
xMat = mat(xArr); yMat = mat(yArr).T
# XtX
xTx = xMat.T * xMat
# 判断矩阵是否可逆
if linalg.det(xTx) == 0:
print("This matrix is singular, cannot do inverse")
return
# 根据公式计算w
w_hat = xTx.I * (xMat.T * yMat)
return w_hat
线性回归拟合
相关系数(Correlation Coefficient)计算：
几乎任意数据集都可以用上述方法建立模型，怎么判断这些模型的好坏呢？为了计算预测值序列和真实值序列的匹配程度，可以计算出这两个序列的相关系数。
计算公式：
相关系数
公式就是用X、Y的协方差除以X的标准差和Y的标准差。协方差是衡量两个变量变化趋势是否相似的一种方法，是同向变化(同时变大或变小)还是反向变化(一个变大一个变小), 同向或者反向的程度如何，计算公式如下:
协方差计算
在Python中，NumPy库提供了相关系数的计算方法：
corrcoef(yEstimate, yActual)来计算相系数，该方法返回一个矩阵，该矩阵包所有两两组合的相关系数。
从上面的结果图可以看出，线性回归的一个问题是可能会欠拟合，标准的线性回归是一种最小均方差的无偏差估计，在计算所有点的时候都是无偏差的计算误差，如果针对不同的点能够对误差进行调整便可以一定程度上避免标准线性回归带来的欠拟合现象。接下来介绍的局部加权线性回归就是采取这种方法对值进行预测。
局部加权线性回归
在该算法中，我们给带预测点附近的每个点赋予一定的权重，越靠近预测点的数据点分配的权重越高，于分类算法kNN一样，此算法每次预测均需事先选取出对应的数据子集，算法代价高。我们用θ表示回归系数，w表示权重, 那么平方误差的表达式就变成:
包含权重的平方误差
用矩阵可表示为：
矩阵表示
对θ求导，令其等于0求极值得：
回归系数求解
其中的W是一个矩阵，用来给每个数据点赋予权重。那我们怎么计算这个W呢？
LWLR使用“核”（与支持向量机中的核类似）来对附近的点赋予更高的权重。核的类型可以自由选择，最常用的是高斯核，高斯核对应的权重如下：
高斯核权重
通过公式可以看到如果xi距离x的距离越小，W(i)就会越大，其中参数k决定了权重的大小。k越大权重的差距就越小，k越小权重的差距就很大，仅有局部的点参与进回归系数的求取，其他距离较远的权重都趋近于零。如果k去进入无穷大，所有的权重都趋近于1，W也就近似等于单位矩阵，局部加权线性回归变成标准的无偏差线性回归，会造成欠拟合的现象；当k很小的时候，距离较远的样本点无法参与回归参数的求取，会造成过拟合的现象。
Python代码实现
# 局部加权线性回归
def lwlr(testPoint, xArr, yArr, k=1.0):
# 先把数组转为矩阵，方便进行线性代数操作
xMat = mat(xArr); yMat = mat(yArr).T
# shape():返回矩阵的维数，是个数组
numLine = shape(xMat)[0]
# eye(num): 生成一个对角矩阵(对角线全为1，其它为0)
weights = mat(eye(numLine))
#计算权重
for i in range(numLine):
# x(i) - x
diffMat = testPoint - xMat[i, :]
weights[i, i] = exp(diffMat * diffMat.T / (-2 * k ** 2))
xTx = xMat.T * (weights * xMat)
if linalg.det(xTx == 0.0): #矩阵不可逆
print("This matrix ix singular, cannot do inverse")
return
# 根据公式求出回归系数w
w = xTx.I * xMat.T * weights * yMat
return testPoint * w
def lwlrTest(testArr, xArr, yArr, k=1):
numLines = shape(testArr)[0]
# 初始化y值，全为0
yHat = zeros(numLines)
for i in range(numLines):
yHat[i] = lwlr(testArr[i], xArr, yArr, k)
return yHat
k=0.1 太大，欠拟合
k=0.01 很好的拟合效果
k=0.003 太小，过拟合
局部加权线性回归也存在一个问题，对于每一个要预测的点，都要重新依据整个数据集计算一个线性回归模型出来，增加了很大的计算量，使得算法代价极高。从上图可以看出当k=0.01时可以得到一个很好的预测。但是下图展示了当k=0.01时(假定我们预测的是x = 0.5这个点)，每个点所占的权重：
每个点的权重图
可以发现除了预测点附近，大多数数据点的权重都接近于零，如果避免这些计算将可以减少程序运行时间，从而缓解因计算量增加所带来的问题。
岭回归
如果数据的特征比样本点还多应该怎么办？或者说输入的数据矩阵不是满秩矩阵，是否还能用前面介绍的普通回归或者局部加权回归来做预测？我们在求回归系数的时候需要算出(XTX)-1，当矩阵不是满秩矩阵的话是不能进行求逆运算的，所以之前的回归方法就不能用了，为了解决这个问题，我们需要对最初的标准线性回归做一定的变化使原先无法求逆的矩阵变得非奇异，使得问题可以稳定求解。我们可以通过缩减系数的方式来处理这些问题,例如岭回归和LASSO.
简单的说，岭回归就是在矩阵(XTX)上加上一个λI从而使得矩阵非奇异，进而能对(XTX) + λI求逆。其中矩阵I是一个单位矩阵，此时回归系数的计算公式将变为：
岭回归系数计算公式
Python代码实现
# 岭回归
def ridgeRegress(xMat, yMat, lam=0.2):
xTx = xMat.T * xMat
# print(shape(xTx))
# print(shape(xMat))
temp = xTx + eye(shape(xTx)[0]) * lam
if linalg.det(temp) == 0.0:
print("This matrix is singular, cannot do reverse")
return
w = temp.I * (xMat.T * yMat)
return w
很明显，我们要找到使预测误差最小的λ，不同的λ可以得到不同的参数w，因此我们可以改变λ的值来得到岭回归系数的变化，通过岭迹图我们可以观察较佳的λ取值：
def ridgeTest(xArr, yArr):
xMat = mat(xArr)
yMat = mat(yArr).T
xMean = mean(xMat, 0)
yMean = mean(yMat, 0)
yMat = yMat - yMean
#数据标注化：所有特征减去各自的均值并除以方差
xMat = (xMat - xMean) / var(xMat, 0)
numTestPts = 30
wMat = zeros((numTestPts, shape(xMat)[1]))
for i in range(numTestPts):
ws = ridgeRegress(xMat, yMat, exp(i-10))
wMat[i, :] = ws.T
return wMat
上述代码可以生成一个30行的回归系数，将其绘制出来可得到岭迹图：
岭迹图
该图绘出了回归系数与log(λ)的关系，在最左边λ系数最小时，可以得到所有系数的原始值(与标准线性回归相同); 而在右边，系数全部缩减为0, 从不稳定趋于稳定；为了定量的找到最佳参数值，还需要进行交叉验证。要判断哪些变量对结果的预测最具影响力，可以观察他们的系数大小即可。
中心化和标准化
在回归问题和一些机器学习算法中通常要对原始数据进行中心化和标准化处理，也就是需要将数据的均值调整到0，标准差调整为1, 计算过程很简单就是将所有数据减去平均值后再除以标准差，上述代码中就是先对数据进行了标准化处理，之所以需要进行中心化其实就是个平移过程，将所有数据的中心平移到原点。而标准化则是使得所有数据的不同特征都有相同的尺度Scale, 这样在使用梯度下降法以及其他方法优化的时候不同特征参数的影响程度就会一致了。
交叉验证
为了定量的找到岭回归最佳参数值λ，还需要进行交叉验证。一般地，交叉验证就是通过把数据集分成若干等份，每一份轮流作测试集，其余数据作训练集进行轮流训练与测试。如果把数据集分成n份，就叫n-folds交叉验证，这样会得到n个错误率，然后取这n个的平均值作为最终的结果。
# 交叉验证
def crossValidation(xArr, yArr, numSplit):
# 误差矩阵，每行有30个lam得到的结果
errorMat = zeros((numSplit, 30))
data_len = len(yArr)
indexList = list(range(len(yArr)))
for i in range(numSplit):
trainX = []
trainY = []
testX = []
testY = []
# 打乱列表顺序，获得随机效果
random.shuffle(indexList)
# 划分出训练集和测试集
for j in range(data_len):
if j < 0.9 * data_len:
trainX.append(xArr[indexList[j]])
trainY.append(yArr[indexList[j]])
else:
testX.append(xArr[indexList[j]])
testY.append(yArr[indexList[j]])
# 求出此次的回归系数 30 * 8
wMat = ridgeTest(trainX, trainY)
# 求每一组回归系数的误差
for group in range(30):
testXMat = mat(testX)
trainXMat = mat(trainX)
trainXMean = mean(trainX, 0)
trainXVar = var(trainXMat, 0)
# 训练集做了标准化处理，测试集也要做相同处理
testXMat = (testXMat - trainXMean) / trainXVar
trainYMat = mat(trainY)
trainYMean = mean(trainYMat, 0)
yEst = testX * mat(wMat[group, :]).T + trainYMean
errorMat[i, group] = rssError(yEst.T.A, array(testY))
#计算每个lamda的平均误差
meanErrors = mean(errorMat, 0)
minErrors = float(min(meanErrors))
bestWeightIndex = nonzero(meanErrors == minErrors)
print(bestWeightIndex[0])
LASSO回归
LASSO(The Least Absolute Shrinkage and Selection Operator)是另一种缩减方法，将回归系数收缩在一定的区域内。LASSO的主要思想是构造一个一阶惩罚函数获得一个精炼的模型, 通过最终确定一些变量的系数为0进行特征筛选。
LASSO的惩罚项为:
LASSO惩罚项
而岭回归的惩罚项：
岭回归的惩罚项
我们可以看到，它俩唯一的不同点在于，这个约束条件使用绝对值代替了平方和。虽然约束形式只是稍作变化，但是相比岭回归可以直接通过矩阵运算得到回归系数相比，这细微的变化却极大的增加了计算复杂度。下面介绍更为简单的方法来得到与之差不多的结果。
前向逐步回归
前向逐步回归是一种贪心算法，即每一步都尽可能的减少误差。一开始，所有的权重都设为1，然后每一步所做的决策是对某个权重增加或减少一个很小的值。
伪代码如下：
数据标准化，使其分布满足0均值和单位方差
在每轮迭代过程中:
设置当前最小误差lowestError为正无穷
对每个特征:
增大或缩小:
改变一个系数得到一个新的w
计算新W下的误差
如果误差Error小于当前最小误差lowestError:
设置Wbest等于当前的W
将w设置为新的Wbest
python代码实现
def stageWise(xArr, yArr, eps=0.01, numIt=100):
xMat = mat(xArr)
yMat = mat(yArr).T
# 标准化处理
yMat -= mean(yMat, 0)
xMat = (xMat - mean(xMat, 0)) / var(xMat, 0)
ws = zeros((shape(xMat)[1], 1))
returnMat = zeros((numIt, shape(xMat)[1]))
# 每轮迭代
for i in range(numIt):
lowestError = inf
# 对每个特征
for j in range(shape(xMat)[1]):
# 增大缩小
for sign in [-1, 1]:
wsTest = ws.copy()
wsTest[j] += eps * sign
yEst = xMat * wsTest
rssE = rssError(yMat.A, yEst.A)
if rssE < lowestError:
lowestError = rssE
wsMax = wsTest
ws = wsMax.copy()
returnMat[i, :] = ws.T
return returnMat
我们取步长为0.01，迭代200次来看看结果：
eps = 0.01， numIt = 200
我们可以看到结果中w1和w6都是0，这说明他们对目标值的预测没有任何影响，也就是说这些特征可能是不需要的。另外，一段时间后系数就已经饱和并在特征值之间来回震荡，这是应为步长太大的缘故。我们可以看到第一个权重在0.04和0.05之间来回震荡。说明最优的权重应该是在0.04-0.05之间，这时我们应该把步长调小。但如果迭代次数过多，逐步线性回归算法就会与常规的最小二乘法效果类似。使用0.005的步长并经过1000次迭代后的结果图：
回归系数与迭代次数关系
逐步线性回归的好处在于它可以帮助人们理解现有的模型并作出改进。当构建了一个模型后，可以运行该算法找出重要的特征，这样就有可能及时停止对那些不重要特征的收集。最后，如果用于测试，该算法每100次迭代就可构建出一个模型，可以使用上文所述的交叉验证法来比较这些模型，最终选择使误差最小的模型。
