机器学习系列-Logistic Regression(1)

机器学习

这是记录自学的过程，目前的理论基础就是：大学高等数学+线性代数+概率论。编程基础：C/C++,python

在观看机器学习实战这本书，慢慢介入。相信有读过以上三门课的人完全可以开始自学机器学习了，当然我上面这三门课学的一般，所以你只知道有这么一个公式或名词，不懂可以百度之深究之。在写这篇文章的时候作者机器学习还没学完，故文章中的错误还请不吝指出。再次声明，系列文章只是分享学习过程，学习点滴，不能保证文章的技术含量。后续再技术的不断完善中，我会重新再捋一遍这个学习过程，纠正其中错误。

目前的学习方法如下：

理论：斯坦福Andrew NG machine learning，系列课程。

代码实战：引用书籍 机器学习实战（python）。

真枪实战： 研究一个案例，建模选择合适的机器学习算法，分析数据。

演变：基于C/C++ 在linux下实现部分机器学习过程。将在学习中期阶段去实现。

一、Logistic Regression和Sigmoid函数

              关于Logistic Regression的定义，网上颇多，这里不再赘述。按道理来说机器学习的步骤是分析一个案例，然后选择一个函数模型h(x),引用一段andrew的原话：

We could approach the classification problem ignoring the fact that y is

discrete-valued, and use our old linear regression algorithm to try to predict

y given x. However, it is easy to construct examples where this method

performs very poorly. Intuitively, it also doesn’t make sense for hθ(x) to take

values larger than 1 or smaller than 0 when we know that y ∈ {0, 1}.

              linear regression的例子，其实我们高中就有接触过，所谓的一次函数类似y=ax+b，就是线性回归的函数简单单参数模型，这里建模的目的是为了pridict y，显然y的值不可确定范围。但是Logistic Regression，它的预测结果范围是可知的，比如本篇案例中Y的值就有明确定义 y ∈ {0, 1}.那么就要重新寻找函数模型h(x).

              有些技术总是站在巨人的肩膀上，大学课程中显然已经有接触到符合上述定义的函数模型：

              至于有些人会问，初学者怎么会知道选择这样一个模型，那我只能说可能你对高等数学的内容有点忘了，选择这样一个函数模型原因在于，g(z)的取值范围在0~1之间，可以规定<0.5属于0，>0.5属于1.所以选择什么函数模型，是基于实际问题的分析，然后在数学模型中筛选出一种合适的函数。此函数学名 Sigmoid：是一种阶跃函数，当横坐标轴刻度足够大时，它的图像看上去就像一个阶跃函数：

二、最佳回归系数的理论推导

              如上我们选择Sigmoid函数。g(z)中的z可以理解为数据输入，数据可能不止一类，我们可以定义一下公式：

z = w0x0+w1x1+w2x2+w3x3+.......+wnxn.

上述，x0,x1,x2......为输入的数据，w0,w1,w2......为对应的参数，一组向量，也就是我们要通过训练大量数据得出的回归系数。为了便于理解线面选择一个有两个数据输入类型(X1,X2)的例子进一步分析。这里一般都假设x0=1.

              好了现在是确定了分类器的函数，接下的目的是如何预测Y呢？引入一个数据分析便于理解：

              上述数据集摘自机器学习实战中的数据片段，可以忽略他的实际意义，前面两列为输入数据X1，X2,后面一列为分类情况，可以看出有两个类别，基于这样的数据类型，我们就可以用上述Sigmoid函数模型来预测输入其他数据后的分类结果。

              接下来要做的事就是如何得出最佳回归系数，先了解一下理论知识。我们最终确定回归系数也是通过让h(x)不断接近y(数据集中实际的值)，然后选择最佳回归系数。这里引入一个cost function的概念：J(w) = h(x) - y. 我们训练数据的目的就是让J(w)达到最小值，到这里你可能注意到h(x)的模型中含有e，单纯的加减法只会让操作更加复杂，另外更为权威的解释就是此时J(w)是一个非凸函数，这意味着J(w)有许多局部最值，这样我们很难找到他的最值（通过梯度上升）。针对这种情况，数学中有另外一种方法来最小化它。

我们先对g(z)求导：

              g(z)的导函数看起来有点门道，当g(z)取0或者1时，g(z)有极值?也就是当z -> 正无穷：g(z) 接近1；z -> 负无穷：g(z) 接近0。(g(z) = h.(x))

不妨假设：

              假设我们数据集都是互相独立的，那么就可以用极大值似然估计法来求出上面的最大值,其作用为：在已知试验结果（即是样本）的情况下，用来估计满足这些样本分布的参数，把可能性最大的那个参数作为真实的参数估计。下面是似然函数的推导：

          求似然函数的极大值的一个方法就是求导，当然还有其他方法andrew ng后面还介绍fisher scoreing，但是这个算法偏难，我准备放到后面去研究实现。上述中的θ就是我们要确定的回归系数，似然函数的极大值受θ
影响，所以在寻找极大值的过程中我们采用了一个算法：梯度上升。关于梯度的概念：寻找某个函数的最大值，最好的方法是沿着这个函数的梯度的方向探寻。梯度的简单理解就是对某个点的求导。

我们定义以下式子来更新我们回归系数：

θ:=θ+α∇θℓ(θ)

α学习速率，可以理解为步长，更新回归系数的速度，由我们自己初始化定义，往往选值偏小慢慢收敛。

对似然函数的求导得到梯度:

代入上式回归系数的更新式子中：

接下来会用代码的形式用数据集来训练出一组参数出来。

三、梯度上升算法的代码实现

以上的推导最终得到上面的一个式子，梯度上升算法的核心。可以看到影响回归系数的主要还是数据本身和预测值与实际值之间的误差。下面是梯度上升的伪代码:

每个回归系数初始化为1；

重复R次：

计算整个数据集的梯度；

使用alpha * gradient 来更新回归系数的值（有可能是向量）；

返回回归系数。

以下讲解机器学习实战中提供的源代码：

1.load数据from file。

def loadDataSet():    dataMat = []; labelMat = []    fr = open('testSet.txt')    for line in fr.readlines():        lineArr = line.strip().split()        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])        labelMat.append(int(lineArr[2]))    return dataMat,labelMat

2.sigmoid 函数

def sigmoid(inX):    return 1.0/(1+exp(-inX))

3.梯度上升算法:

def gradAscent(dataMatIn, classLabels):    dataMatrix = mat(dataMatIn)             #转化成矩阵    labelMat = mat(classLabels).transpose() #convert to NumPy matrix    m,n = shape(dataMatrix)#获取矩阵属性，m行n列    alpha = 0.001#init alpha    maxCycles = 500#define the recycle times    weights = ones((n,1))#initialize the 回归系数，为n行一列单位向量    for k in range(maxCycles):              #heavy on matrix operations        h = sigmoid(dataMatrix*weights)     #用初始值回归系数乘以数据，可以看到#这边是整个数据集乘以回归系数，矩阵相乘： m*n 乘以 n*1，得到m行一列的单位向量,#也就是得到预测值h(x).        error = (labelMat - h)              #误差值，相当于(y - h(x))。        weights = weights + alpha * dataMatrix.transpose()* error #更新回归系数，代码中#对应上述式子：θj := θj + α (y(i) − hθ(x(i))) * xj( i).其中error为(y(i) − hθ(x(i)))#dataMatrix为输入数据集x。transpose为矩阵转置，dataMatrix转置之后为n*m 与 m*1的error可以相乘    return weights

上述代码保存到logRegres.py中，当然中文注释好像python无法编码，需要删除。在python实践验证:

>>>import logRegres

>>>dataArr,labelMat = logRegres.loadDataSet()

>>>logRegres.gradAscent(dataArr,labelMat )

将会得到回归系数向量。

matrix：

([[ 4.12414349], ->X0
        [ 0.48007329], ->X1
        [-0.6168482 ]]) ->X2

不过我们只关心X1，X2。因为X0对应的数据是我们自己初始化的，为1.

下面代码是将整个训练情况和数据分布用图像画出来：

def plotBestFit(weights):    import matplotlib.pyplot as plt    dataMat,labelMat=loadDataSet()    dataArr = array(dataMat)    n = shape(dataArr)[0]     xcord1 = []; ycord1 = []    xcord2 = []; ycord2 = []    for i in range(n):        if int(labelMat[i])== 1:            xcord1.append(dataArr[i,1]); ycord1.append(dataArr[i,2])        else:            xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2])    fig = plt.figure()    ax = fig.add_subplot(111)    ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')    ax.scatter(xcord2, ycord2, s=30, c='green')    x = arange(-3.0, 3.0, 0.1)    y = (-weights[0]-weights[1]*x)/weights[2]    ax.plot(x, y)    plt.xlabel('X1'); plt.ylabel('X2');    plt.show()

接着输入：

>>> weights = logRegres.gradAscent(dataArr,labelMat)
>>> logRegres.plotBestFit(weights.getA())
可以得到如图：

按照书中的描述，这个分类方法效果相当不错，图上可以看出只分错了四个点。但是尽管例子简单且数据集很小，但这个方法却需要大量的计算：一次迭代需要300次乘法,这种算法还被叫做Batchgradient ascent ，大批量梯度算法。下面介绍书中另外一个改进算法，随机梯度上升。

四、优化算法：随机梯度上升

          前面也说过在每一次迭代中有300次乘法，就是因为遍历了整个数据集（数据集为300*3），代价颇高。改进算法为：每次迭代只用一个样本点来更新回归系数，该方法称为随机梯度上升算法，英文：Stochastic gradient ascent。

          随机梯度上升算法伪代码如下：

每个回归系数初始化为1；

对数据中每个样本：

计算整个数据集的梯度；

使用alpha * gradient 来更新回归系数的值（有可能是向量）；

返回回归系数。

以下是具体代码：

def stocGradAscent0(dataMatrix, classLabels):    m,n = shape(dataMatrix)    alpha = 0.01    weights = ones(n)   #初始化回归系数为{1,1,1} 为numpy数组，非矩阵    for i in range(m):        h = sigmoid(sum(dataMatrix[i]*weights))#dataMatrix[i]类似{x,y,z}为数据集中的一个数据，#与weights相乘之后再经过sum，结果为数值 x*1+y*1+z*1。        error = classLabels[i] - h#也是数值操作        weights = weights + alpha * error * dataMatrix[i] #每次用一个数据来更新回归系数    return weights

输入以下命令：

>>>dataArr,labelMat=logRegres.loadDataSet()
>>> weights = logRegres.stocGradAscent0(array(dataArr),labelMat)
>>> logRegres.plotBestFit(weights)

得到优化后的算法结果:

          单单从结果上来看，这个算法最后的效果不如人意，没有第一个算的准确度高，有差不多1/3的数据预测错误。关于如何判断一个优化算法优劣的可靠方法，书中有详细介绍这里不再赘述。

          我们直接来看又一次改进后的算法：

        1.之前我们的alpha是保持不变的，这边每次迭代的时候都会适当调整alpha，调整的思想：alpha还是有一个最基本的值比如0.01，在此基础上加上一个单调递减的式子（4/(1+i+j)）,j是迭代次数，为什么选择4了，我是这么理解的：

          i是数据下标，例子中例子中为0-300，蛮算为300，j是迭代次数可能值不确定但肯定不会太小，这样就是1+i+j > 300 ,这样的话3/(1+i+j)  < 0.01,把其中的3改为4，那这个式子就比0.01大了，最终可以理解为alpha加上一个以alpha为最大值单调递减的式子，这样子的行为为了保证常数alpha的值对数据还是有一定的影响。不然的话你想要是以0.5,0.08这样的最值单调减，可能就会使0.01失去意义。以上只是个人见解，最终还是要以实际数据为主，选择合适的式子、参数来确定。

        这个式子可能不会严格递减当j<<max(i)。使alpha随着迭代次数不断减小，但永远不会减到0，最小就是初始值0.01。

     2.通过随机获取样本来更新回归系数，更新之后再删除此样本，然后在进行下一次迭代这种方法也是遍历了所有的数据集，只是遍历的顺序是随机的，可以减少周期性波动，你想有时候可能一连串的数据波动很大，随机获取就可以减免这种情况，当然也会有随机到不好的时候。

         详细见以下代码:

def stocGradAscent1(dataMatrix, classLabels, numIter=150):    m,n = shape(dataMatrix)    weights = ones(n)   #initialize to all ones    for j in range(numIter): #initialize value is 150        dataIndex = range(m) #300        for i in range(m):            alpha = 4/(1.0+j+i)+0.001    #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

代码理解具体可见上述两点分析。

最后看一下这个算法实验结果：

>>>dataArr,labelMat=logRegres.loadDataSet()
>>> weights = logRegres.stocGradAscent1(array(dataArr),labelMat)
>>> logRegres.plotBestFit(weights)

可见这个比较接近第一次梯度上升算法的结果，但是使用的迭代次数和乘法次数都更少。

有兴趣的还可以通过调整alpha和迭代次数来看看调整后的结果。