机器学习之回归模型

基本形式

线性模型(linear model)就是试图通过属性的线性组合来进行预测的函数，基本形式如下：

f(x)=wTx+b

许多非线性模型可在线性模型的基础上通过引入层结构或者高维映射（比如核方法）来解决。线性模型有很好的解释性。

线性回归

线性回归要求均方误差最小:

(w∗,b∗)=argmin∑i=1m(f(xi)−yi)2

均方误差有很好的几何意义，它对应了常用的欧式距离(Euclidean distance)。基于均方误差最小化来进行模型求解称为最小二乘法(least square method)，线性回归中，最小二乘发就是试图找到一条直线，使得所有样本到直线的欧式距离之和最小。

我们把上式写成矩阵的形式：

w∗=argmin(y−Xw)T(y−Xw)

这里我们把b融合到w中，X中最后再加一列1。为了求最小值，我们对w求导并令其为0：

2XT(Xw−y)=0

当XTX为满秩矩阵(full-rank matrix)时是可逆的。此时：

w=(XTX)−1XTy

令xi=(xi,1)，可以得到线性回归模型：

f(xi)=xTi(XTX)−1XTy

最小二乘法的代码如下：

def standRegres(xArr,yArr):    xMat = mat(xArr); yMat = mat(yArr).T    xTx = xMat.T*xMat    if linalg.det(xTx) == 0.0:        print "This matrix is singular, cannot do inverse"        return    ws = xTx.I * (xMat.T*yMat)    return ws

然而，现实情况是很多时候，XTX并不满秩，这样w的解可能有多个解，它们都能使得均方误差最小化。碰到数据特征比样本数还要多的情况，一种方案是引入正则化(regularization)项。另一种方案是可以“缩减系数”，比如用岭回归(ridge regression)、LASSO法（该方法效果好但是计算复杂），前向逐步回归（可以得到与LASSO差不多的效果，但更容易实现）、LAR、PCA回归等。

局部加权线性回归

我们在介绍岭回归之前先看一下局部加权线性回归。这首针对线性回归容易出现“欠拟合”现象提出的，因为它求的是具有最小均方误差的无偏估计。如果模型欠拟合将不能取得最好的效果。所以有些方法在估计中引入一些偏差，从而降低均方误差。其中一个方法就是局部加权线性回归（Locally Wegihted Linear Regression,LWLR），我们给待预测点附近的每个点赋予一定的权重。回归系数w形式如下：

w=(XTWX)−1XTWy

w就是用来给每个数据点赋权重的矩阵。
LWLWR使用核方法来赋权重，最常用的就是高斯核：

w(i,i)=exp(|xi−x|−2k2)

这样构建了一个只含对角元素的权重矩阵w，并且点x与x(i)越近，w(i,i)将会越大。k决定了对附近的点赋多大权重。
函数代码如下：

def lwlr(testPoint,xArr,yArr,k=1.0):    xMat = mat(xArr); yMat = mat(yArr).T    m = shape(xMat)[0]    weights = mat(eye((m)))    for j in range(m): #next 2 lines create weights matrix        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

文末给出的下载链接中regreesion.py文件中，有个线性加权回归例子。

岭回归

就是在XTX加上一个λI使得矩阵非奇异，这样XTX+λI就可以求逆了。于是：

w=(XTX+λI)−1XTy

岭回归最先用来处理特征数多于样本数的情况，现在也用于在估计中加入偏差，从而得到更好的估计。这里通过引入λ来限制所有w之和，通过引入该惩罚项，能够减少不重要的参数。这种方法叫做“缩减”(shrinkage)。我们通过预测误差最小化来得到λ，即通过选取不同的λ来重复训练测试过程，最终得到一个使预测误差最小的λ。
相关代码如下：

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 wsdef ridgeTest(xArr,yArr):    xMat = mat(xArr); yMat=mat(yArr).T    yMean = mean(yMat,0)    yMat = yMat - yMean     #to eliminate X0 take mean off of Y    #regularize X's    xMeans = mean(xMat,0)   #calc mean then subtract it off    xVar = var(xMat,0)      #calc variance of Xi then divide by it    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

ridgeRegres用于计算回归系数，ridgeTest用于在一组λ上测试结果。
同样在使用特征时，都要先进行归一化处理。

LASSO

LASSO是least absolute shrinkage and selection operator的简称。我们注意到增加如下约束，最小二乘法回归会得到与岭回归一样的公式：

∑k=1nw2k≤λ

上式限定了所有回归系数的平方和不大于λ。在用最小二乘法回归时，当两个或更多特征相关时，可能会得到一个很大的正系数和一个很大的负系数，正因为上面约束，岭回归避免了这个问题。
LASSO的约束为：

∑k=1n|wk|≤λ

虽然用绝对值取代了平方和，但结果却差别很大。在λ足够小的时候，一些系数会因此被迫缩减为0，这样结果就有很好的稀疏性，这个特性可以帮助我们更好地理解数据。两个看起来相差无几，但细微的变化却极大增加了计算复杂度。为了在此约束下求解，需要使用二次规划算法。

前向逐步回归

前向逐步回归可以得到与LASSO差不多的结果，但更加简单，是一种贪心算法，每一步都尽可能减少误差。一开始，所有权重设为1，然后每一步所做的决策是对某个权重增加或减少一个很小的值。
代码如下：

def stageWise(xArr,yArr,eps=0.01,numIt=100):    xMat = mat(xArr); yMat=mat(yArr).T    yMean = mean(yMat,0)    yMat = yMat - yMean     #can also regularize ys but will get smaller coef    xMat = regularize(xMat)    m,n=shape(xMat)    returnMat = zeros((numIt,n)) #testing code remove    ws = zeros((n,1)); wsTest = ws.copy(); wsMax = ws.copy()    for i in range(numIt):#could change this to while loop        #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

逐步线性回归主要优点在于它可以帮助人们理解现有模型并做出改进。当构建一个模型后，运行该算法找出重要的特征，这样就可以及时停止对那些不重要特征的收集。当应用缩减方法时，模型增加了偏差(bias)，也减少了方差。

对数几率回归（逻辑回归）

上面是回归任务，如果我们需要的是分类任务呢？我们只需要找一个单调可微函数将分类任务的真实标记与线性回归的预测值就可以了。

二分类最理想的是“单位阶跃函数”(unit-step function)。但是此函数不连续，于是我们希望找一个单调可微的替代函数：

y=11+e−z

logistic function是一种Sigmoid函数，我们将z代入线性回归模型：

y=11+e−wTx+b

变换一下：

lny1−y=wTx+b

如果y是正样本的可能性，那么1-y就是负样本的可能性，两者比值称为几率，左边的意义就是对数几率。
实际上上式是在用线性回归的模型预测结果去逼近真实标记的对数几率。虽然名字叫“回归”，却是一种分类方法。这种方法优点很多，它直接对分类可能性进行建模，无需事先假设数据分布。此外，对率函数是任意阶可导的凸函数，利用现有的数值优化算法可以直接得到最优解。
下面来看以下具体推导：
假设y=1是二分类中正样本，

lnp(y=1|x)p(y=0|x)=wTx+b

p(y=1|x)=ewTx+b1+ewTx+b

p(y=0|x)=11+ewTx+b

于是，我们采用极大似然来估计w和b:

l(w,b)=∑i=1mlnp(yi|xi;w,b)

即令每个样本属于其真实标记的概率越大越好。
我们可以再转换成最小值问题，用经典的数值优化算法如梯度下降法(gradient descent method)，牛顿法（Newton method）解决。

在此处我们简便起见，直接采用梯度上升算法。我们知道梯度算子总是指向函数值最快的方向。我们每次向增长最快的方向移动一个值，称为步长，记为a，所以梯度上升算法的公式为：

w=w+a∂f(w)∂w

该公式将一直迭代执行，直到达到某个停止条件（比如迭代次数）为止。
实现代码为：(当然之前数据先标准化)

def gradAscent(dataMatIn, classLabels):    dataMatrix = mat(dataMatIn)                 labelMat = mat(classLabels).transpose()     m,n = shape(dataMatrix)    alpha = 0.001    maxCycles = 500    weights = ones((n,1))    for k in range(maxCycles):                      h = sigmoid(dataMatrix*weights)        error = (labelMat - h)              #vector subtraction        weights = weights + alpha * dataMatrix.transpose()* error    return weights

代码中求导不好操作，这里我们用差分可以起到同样的效果。
梯度上升算法在每次更新回归系数时要遍历整个数据集，改进方法是一次仅用一个样本点来更新回归系数，称为“随机梯度上升算法”，这是一种在线算法。上述代码仍有改进之处，比如a在每次迭代时候都会调整，这会缓解数据波动或者高频波动。a会随着得到次数不断减小，但永远不会减到0，这样多次迭代之后新数据仍具有一定的影响。
改进后的随机梯度上升代码如下：

def stocGradAscent(dataMatrix, classLabels, numIter=150):    m,n = shape(dataMatrix)    weights = ones(n)   #initialize to all ones    for j in range(numIter):        dataIndex = range(m)        for i in range(m):            alpha = 4/(1.0+j+i)+0.0001    #apha decreases with iteration, does not             randIndex = int(random.uniform(0,len(dataIndex)))#go to 0 because of the constant            h = sigmoid(sum(dataMatrix[randIndex]*weights))            error = classLabels[randIndex] - h            weights = weights + alpha * error * dataMatrix[randIndex]            del(dataIndex[randIndex])    return weights

这里用个例子来使用Logistic回归来预测患有疝病的马存活问题。但是我们发现了样本数据缺失问题，通常的解决方法有：

• 使用可用特征均值填补
• 使用特殊值来填补缺失值，如-1
• 忽略有缺失值的样本
• 使用相似样本的均值填补
• 使用另外的机器学习算法预测缺失值
• …

这里我们全部用0来代替缺失值，这样在更新时不会影响系数。即特征对应0时，系数将保持0.5。我们在数据中还发现了标签缺失，这个在强化学习中可以通过更靠近某些标签就将其标记为那个标签，在本例中由于此类较少，我们直接丢弃。

线性判别分析

（Linear Discriminant Analysis,LDA）与LDA(Latent Dirichlet Allocation)算法不要搞混了哦，虽然简写都是一样的。
LDA的思想很简单：设法将样例投影到一条直线上，使得同类样例的投影点近可能接近，异类样例的投影尽可能原理，在对新样本进行分类时，将其投影到这条直线上，再根据投影点位置来确定分类。

多分类学习

二分类问题推广到多分类主要用拆解法，主要策略有：
- 一对一
- 一对多
- 多对多
这样那个容易碰到类别不平衡问题，就是某类的数据量特别大或者特别少。

最后就是本文数据集和代码的下载地址啦，请点击这里。