5.机器学习实战之Logistic回归

返回目录

上一篇：朴素贝叶斯

1. 简单理论介绍

1.1  Sigmoid函数

    在线性代数中都知道有一个单位阶跃函数，但是该函数在跳跃点上从0瞬间跳跃到1，这个瞬间跳跃过程有时很难处理，幸好有一个函数具有类似的性

质，而且更容易处理，这就是Sigmoid函数。Sigmoid函数的计算公式如下：

 

输入一个z，得到的结果总在0~1之间。

下面给出Sigmoid函数的图形：

为了实现logistic回归分类器，我们可以在每个特征上乘以一个回归系数，然后把所有的结果值相加，将这个总和结果代入sigmoid函数中，进而得到一个范围在0~1之间的数值。这样，z可以表示为  ， x 为特征，也可以表示成向量形式，即，任何大于0.5的数据被分入1类，小于0.5的数据被归入0类。所以，logistic回归也可以被看成是一种概率估计。

根据图形可以得到，如果z > 0, h(z) > 0.5，分到1类；如果z < 0, h(z) < 0.5，分到0类；如果z = 0，即=0，称为决策边界。

    确定了分类器的函数形式后，现在要求的就是最佳的回归系数是多少，也就是说向量是多少？

1.2 Logsic回归的优缺点

优点：计算代价不高，易于理解和实现

缺点：容易欠拟合，分类精度可能不高

2. 算法实现 

2.1 梯度上升法

梯度上升法基于的思想是要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。该公式将一直被迭代执行，直到达到某个停止条件为止，

比如迭代次数达到某个指定值或者算法达到某个可以允许的误差范围。

准备数据：

def loadDataSet():
    dataMat = []; labelMat = []
    path = r'D:\machine learning\machine learning in action\ch05\data\testSet.txt'
    fr = open(path)
    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

梯度上升法：

import numpy as np

def sigmoid(inX):
    return 1.0/(1 + np.exp(-inX))
    
def gradAscent(dataMatIn,classLabels):
    dataMatrix = np.mat(dataMatIn)# MxN,M表示训练样本个数，这里M=100，N表示特征个数，这里N=3
    labelMatrix = np.mat(classLabels).transpose()# 转置后变为Nx1
    m,n = np.shape(dataMatrix)# 
    alpha = 0.001# 步长
    maxCycles = 500# 最大迭代次数
    weights = np.ones((n,1))# Nx1矩阵
    
    for k in range(maxCycles):
        # sigmoid函数的参数为dataMatrix与weights相乘后的结果，注意不是元素相乘
        h = sigmoid(dataMatrix*weights)# 得到模型的预测值，Mx1
        error = (labelMatrix - h)# 使用真实值减去预测值得到误差，Mx1
        weights = weights + alpha * dataMatrix.transpose() * error
    
    return weights

注意：numpy包下的两个array相乘时，做的是元素相乘，两个matrix相乘时，做的是矩阵相乘。

测试：

dataArr,labelMat = loadDataSet()
weights = gradAscent(dataArr,labelMat)
print weights

# out
'''
[[ 4.12414349]
 [ 0.48007329]
 [-0.6168482 ]]
'''

2.2 画出决策边界

%matplotlib inline
def plotBestFit(weights):
    import matplotlib.pyplot as plt
    dataMat,labelMat = loadDataSet()
    dataArr = np.array(dataMat)
    m = np.shape(dataArr)[0]
    xcord1 =[];ycord1 = []
    xcord2 =[];ycord2 = []
    
    for i in range(m):
        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 = np.arange(-3.0,3.0,0.1)
    if isinstance(weights, np.matrixlib.defmatrix.matrix):
        y = np.array((-weights[0]-weights[1]*x)/weights[2])[0]
    elif isinstance(weights ,np.ndarray):
        y = (-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x,y)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

plotBestFit(weights)

2.3 随机梯度上升法

    梯度上升法在每次更新回归系数时都需要遍历整个数据集，如果数据集的特征特别多，那么该方法计算速度就会很慢。一种改进方法是一次只用一个样本点来更新回归系数，该方法成为随机梯度上升法。由于可以在新的样本到来时对分类器进行增量式更新，因而随机梯度上升法是一个在线学习算法。与“在线学习”相对应，一次处理所有数据被称为“批处理”。

随机梯度上升法的伪代码：

所有回归系数初始化为1
对数据集中每个样本
    计算该样本的梯度
    使用alpha x gradient更新回归系数值
返回回归系数值

具体实现：

def stocGradAscent0(dataMatrix, classLabels):  
    dataMatrix = np.array(dataMatrix)  
    m,n = np.shape(dataMatrix)  
    alpha = 0.1  
    weights = np.ones(n)  
    for i in range(m):  
        h = sigmoid(sum(dataMatrix[i] * weights))  
        error = classLabels[i] - h  
        weights = weights + alpha * error * dataMatrix[i]  
    return weights 

随机梯度上升法和梯度上升法的区别：

1) 后者的变量h和误差error都是向量，而前者都是数值

2) 前者没有矩阵的转换过程，所有变量的数据类型都是numpy数组

测试：

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

2.4 改进随机梯度上升法

    在上节使用随机梯度上升法得到了最佳拟合直线图，但是发现拟合效果并不如梯度上升法那么完美，为此可以进行一些改进。

# numInter表示迭代次数，默认为150次
def stocGradAscent1(dataMatrix, classLabels, numInter = 150):  
    dataMatrix = np.array(dataMatrix)  
    m,n = np.shape(dataMatrix)  
    weights = np.ones(n)  
    for j in range(numInter):  
        dataIndex = range(m)  
        for i in range(m):  
            alpha = 4 / (1.0+j+i) + 0.01# alpha值每次迭代时都进行调整，基本上是趋向于减小的  
            randIndex = int(random.uniform(0, len(dataIndex))) # 随机选取更新 
            h = sigmoid(sum(dataMatrix[randIndex] * weights))  
            error = classLabels[randIndex] - h  
            weights = weights + alpha * error * dataMatrix[randIndex]  
            del[dataIndex[randIndex]]  
    return weights  

测试：

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

该算法主要改进了一下三个方面：

1) alpha在每次迭代的时候都会调整，这会缓解数据波动或者高频波动。虽然alpha会随着迭代次数不断减小，但永远不会减小到0，这是因为存在一个

0.01常数项。这样做的原因是为了保证多次迭代之后新数据仍然具有一定的影响。在降低alpha的函数中，alpha每次减少1/ (j + i)，其中j是迭代次数，i

是样本点的下标。这样当j << max(i)时，alpha就不是严格下降的。

2) 通过随机选取样本来更新回归系数，这样可以减少周期性波动。

3) 增加了一个迭代次数作为第三个参数。默认为150次。

下一篇： 支持向量机（SVM）