机器学习实践——Logistic回归

利用Logistic回归进行分类的主要思想：根据现有数据对分类边界线建立回归公式，以此进行分类。

二值型输出分类器的数学原理：理想的函数是能够接受所有的输入，然后预测出类别。例如两个类的情况下，函数的输出是0或是1，这就需要Sigmoid函数。

                                          

因此，为了实现回归分类器，我们可以把每个特征都乘上一个回归系数，然后把所有的结果值相加，总和带入函数中，进而得到一个0-1之间的值，大于0.5就归为1，否则归为0。所以现在目的就是通过最优化方法求的最佳回归系数。

梯度上升法：

思想：想找到一个函数的最大值就应该在该函数的梯度上查找。

梯度上升的迭代公式：

W: = W+a*Grad f(x)

这个公式将一直被迭代，直至达到某个停止条件

 

训练算法：使用梯度上升找到最佳参数

#加载数据，将文本数据转换为列表形式
def loadDataSet():
    dataMat = [];labelMat = []
    fr = open('Ch05/testSet.txt')
    for line in fr.readlines():
        line = line.strip().split()
        #第一列数据全为设置为1.0
        dataMat.append([1.0,float(line[0]),float(line[1])])
        labelMat.append(int(line[2]))
    return dataMat,labelMat

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

#输入数据集以及类标签，输出每个特征的类标签
def gradAscent(dataMatIn,classLabels):
    #将数据集转换为矩阵形式，为了接下来的矩阵计算
    dataMatrix = mat(dataMatIn)
    #print dataMatrix
   # print shape(dataMatrix)
    #transpose函数将矩阵进行转置，这里是将行向量转为列向量
    labelMat = mat(classLabels).transpose()
    m,n =shape(dataMatrix)
    #这是步数
    alpha = 0.001
    #迭代的次数
    maxCycles = 500
    #将初始的最优的系数设为1，n对征应数据集的特数
    weights  = ones((n,1))
    #这里的迭代逐步求梯度是用到数学推导的公式
    #不理解的小伙伴可以去查阅资料
    for k in range(maxCycles):
        h = sigmoid(dataMatrix * weights)
        error = (labelMat - h)
        weights = weights+alpha * dataMatrix.transpose() * error
    return weights

我们通过matplotlib中的方法将数据可视化

分析数据:画出决策边界

#画出决策边界
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:
            #类标号为0的坐标，第一个特征时x坐标，第二个特征为有坐标
            xCord1.append(dataArr[i][1]);yCord1.append(dataArr[i][2])
        else:
            # 类标号为1的坐标
            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和y的值
    x = arange(-3.0,3.0,0.1)
    """这里最大的疑惑都在y上，x的取值是（-3,3）
    想想公式：w0*x0+w1*x1+w2*x2=Z
    因为开始定义x0 =1 x1就是x的值，x2就是y,z是两个类别的分界z=0
    所以带入推导得：y=(-w0-w1*x)/w2
    """
    y = (-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x,y)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

d,l = loadDataSet()
w = gradAscent(array(d),l)
plotBestFit(w.getA())

实验结果：

这种梯度上升算法在计算矩阵乘法的时候效率太低，如果数据少还可以，但现实中数据通常很大，所以进行该算法的改进

#随机梯度上升算法，就是基于增量式的学习方法，online learning
#将逐个样本点进行迭代最优回归系数
def stocGradAscent0(dataMatrix,classLabels):
    m,n = shape(dataMatrix)
    alpha = 0.01
    weights = 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

实验结果：

这种方法将每个样本点进行迭代，是一种在线学习的方法，但是实验结果有待改善。所有介绍另一种改进算法策略，将迭代与随机样本点更新结合在一起

 

#改进的随机梯度上升算法，将之前两种方法结合，
# 即采用迭代，并且采用随机的样本点更新回归系数
def stoGradAscent1(dataMatrix,classlabels,numIter = 150):
    m,n = shape(dataMatrix)
    weights = ones(n)
    #迭代的次数
    for j in range(numIter):
        dataIndex = range(m)
        for i in range(m):
            #每次都调整步数
            alpha = 4/(1.0+j+i)+0.1
            #随机的生成一个数，以此找到样本点
            randInx = int(random.uniform(0,len(dataIndex)))
            h = sigmoid(sum(dataMatrix[randInx]*weights))
            error = classlabels[randInx] - h
            weights = weights + alpha * error * dataMatrix[randInx]
            #书中这种删除方法还是会让数据重复出现，并且有些样本点并不能参与回归系数的更新
            del(dataIndex[randInx])
    return weights

实验结果：

这次的改动主要有两点
第一：是alpha值的在每次迭代都会调整，虽然alpha值会随着迭代一直减少，但永远不会减少到0.这样做的原因是为了保证在多次迭代之后新数据仍然有一定的影响力。

第二：这里不再是逐个样本点更新回归系数而是采用随机样本点进行更新回归系数，这种方法随机从列表中选出一个值，这里有必要说一下，书中给出的随机选取样本点的代码是错误的，不但有些样本点会重复出现，甚至有些样本点并不能参与回归系数的更新中。在这按照这种随机法给出我的实验验证：

# encoding: utf-8
from numpy import  *
dataIndex =range(10)
dataSet = ['1','2','3','4','5','6','7','8','9','10']
for i in range(10):
  randInx = int(random.uniform(0,len(dataIndex)))
  #打印数据集，看是否所有的样本点都参与更新回归
  print dataSet[randInx]
  del dataIndex[randInx]

实验结果:

从实验结果中可以看出，1重复了3次，3也重复了，并且6,8,10并未出现，由此可以看出这种随机方法的错误之处。

 

示例：从疝气病症预测病马的死亡率

首先要解决的就是处理数据中的缺失值，这里选择将缺失值置为0。因为这恰好适用于Logistic回归，回归系数的更新公式：

weights = weights + alpha * error * dataMatrix[randInx]

这里将 dataMatrix[randInx] 的值为0，那么weights将不做更新，并不影响结果

def classifyVecor(inX,weights):
    prob = sigmoid(sum(inX*weights))
    if prob > 0.5:
       return 1.0
    else:
        return 0.0

def colicTest():
    frTrain = open('Ch05/horseColicTraining.txt')
    frTest = open('Ch05/horseColicTest.txt')
    trainingSet = [];trainingLabel = []
    #获得训练集
    for line in frTrain.readlines():
        Currentline = line.strip().split()
        lineArr = []
        #一共有20个特征，1个类标签
        for i in range(21):
            lineArr.append(float(Currentline[i]))
        trainingSet.append(lineArr)
        trainingLabel.append(float(Currentline[21]))
        #通过梯度上升求得回归系数
    trainWeights = stoGradAscent1(array(trainingSet),trainingLabel,500)
    errorCount = 0;numTestVec = 0.0
    #将检验集转换为数组
    for line in frTest.readlines():
        numTestVec += 1.0
        curline = line.strip().split('\t')
        lineArr = []
        for i in range(21):
            lineArr.append(float(curline[i]))
            #计算失败的次数
        if int(classifyVecor(array(lineArr),trainWeights))!=int(curline[21]):
            errorCount += 1.0
    errorRate = float(errorCount)/float(numTestVec)
    print "the error rate of this test is:%f"%errorRate
    return errorRate
#多次试验求平均错误率
def multiTest():
    numTests = 10;errorSum = 0.0
    for k in range(numTests):
        errorSum += colicTest()
    print "经过了%d次迭代后，平均错误率是：%f"%(numTests,errorSum/float(numTests))

multiTest()

实验结果：

 

总结：

Logistic回归的目的是寻找一个非线性的函数Sigmoid的最佳拟合参数，求解的过程可以由最优化算法来解决，在最优化算法中，有采用梯度上升算法，进而为了节省计算机资源从而采用在线学习的随机梯度上升算法， 可以随着数据的改变而改变参数值，而不需要读取整个数据集来进行批处理运算