Logistic回归实战

本文借鉴于《机器学习实战》，本文对其进行了代码更新与注释。
以二分类为例初步介绍Logistc回归

基于Logistic回归和Sigmoid函数的分类

1、函数公式

2、函数曲线图

3、曲线图分析

Sigmoid（x）函数曲线，其变动范围于（0，1）之间，随着x的增大，对应的Sigmoid值逼近于1，随着x的减小，Sigmoid值逼近于0。在x=0处，函数Sigmoid(x)近似于跳跃函数。
4、Sigmoid函数与Logistic回归的关系

因此为了实现Logistic回归分类器，可以在每一个特征上都乘以回归系数，所有的结果值相加，将这个和带入Sigmoid函数，进而得到一个位于（0，1）范围内的值，任何大于0.5的数据被分入1类，否则分入0类。即

Sigmoid(W0X0+W1X1+…+Wn*Xn)？>0.5:1，

Wi表示特诊Xi上的回归系数

Sigmoid(W0X0+W1X1+…+Wn*Xn)=0.5为两个类别的分界处，即

W0X0+W1X1+…+Wn*Xn=0为决策边界。

Sigmoid函数的公式可得：

其中向量X是分类器的输入数据，向量w为我们要找到的最佳参数（回归系数向量）

####基于最优化方法的最佳回归系数确定
1、梯度上升法
基本思想：要找到某函数的最大值，最好的方法就是沿着该函数的梯度方向寻找。梯度公式如下：

现有梯度上升算法的迭代公式，其中alpha代表移动步长

迭代条件终止条件：迭代次数达到某个指定值、或算法达到某个允许的误差范围之内。
####算法实现
现有数据testSet.txt，其中含有100个样本点，每个点包含两个数据值特征：X1和X2，以及类标签（0、1）

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

（1）批处理梯度上升算法伪代码
每个回归系数初始化为1
重复R次：
        计算整个数据集的梯度
        使用alpha*gradient更新回归系数的向量
返回回归系数

import numpy as np
from numpy import *
from math import *
def loadDataSet(filepath):
    dataMat=[]; labelMat=[]
    fr=open(filepath)
    for line in fr.readlines():
        lineArr=line.strip().split()
        #添加特征X0，并将值设为1.0,便于计算
        dataMat.append([1.0,np.float(lineArr[0]),np.float(lineArr[1])])
        labelMat.append(int(lineArr[2]))
    return dataMat,labelMat
#Sigmoid函数，inX为矩阵
def sigmoid(inX):
    m,n=np.shape(inX)
    returnX=np.ones((m,n))
    for i in range(m):
        for j in range(n):
            returnX[i][j]=1.0/(returnX[i][j]+exp(-inX[i][j]))
    return returnX

#Logistic回归梯度上升优化算法，遍历500次数据集的批量上升优化算法
def gradAscent(dataMatIn, classLabels):
    #(100,3)矩阵
    dataMatrix = mat(dataMatIn)
    #(1,100)矩阵
    labelMat = mat(classLabels).transpose()
    m,n = shape(dataMatrix)
    alpha = 0.001  #alpha为向目标移动的步长
    maxCycles = 500  #迭代次数
    weights = ones((n,1)) #初始化回归系数为1
    for k in range(maxCycles):
        h = sigmoid(dataMatrix*weights) #整个数据集的梯度，大于0.5的数据划分到1类，小于0.5的数据划分到0类
        #计算真实类别与预测类别的差值，以该差值的方向调整回归系数
        error = (labelMat - h)
        weights = weights + alpha * dataMatrix.transpose()* error
    return weights

（2）在线学习随机梯度上升算法伪代码
每个回归系数初始化为1
对数据集中每个样本：
        计算该样本的梯度
        使用alpha*gradient更新回归系数值
返回回归系数值

#只遍历一次数据集的在线学习梯度上升算法
def sigmoid0(inX):
    return 1.0/(1+exp(-inX))
def stocGradAscent0(dataMatrix,classLabels):
    #dataMatrix为ndarray类型（100，3）
    dataMatrix=array(dataMatrix)
    m,n=shape(dataMatrix)
    alpha=0.01 #alpha为向目标移动的步长
    weights=ones(n)
    for i in range(m):
        h=sigmoid0(sum(dataMatrix[i]*weights))  #h梯度为数值
        error=float(classLabels[i]-h) #数值计算
        weights=weights+alpha*error*dataMatrix[i]
    return weights

（3）改进的在线学习梯度上升优化算法
使用样本随机选择和alpha动态减少机制的随机梯度上升优化算法

#改进的在线梯度上升优化算法，数据集的迭代次数增加
#使用样本随机选择和alpha动态减少机制的随机梯度上升优化算法
def stocGradAscent1(dataMatrix,classLabels,numIter=150):  #numIter为迭代次数
    dataMatrix=array(dataMatrix)
    m,n=shape(dataMatrix)
    weights=ones(n)
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            #alpha每次迭代时都需要调整
            alpha=4/(1.0+j+i)+0.01
            #随机选取样本
            randIndex=int(random.uniform(0,len(dataIndex)))
            #计算该样本的梯度
            h=sigmoid0(sum(dataMatrix[randIndex]*weights))
            error=classLabels[randIndex]-h
            weights=weights+alpha*error*dataMatrix[randIndex]
            del(dataIndex[randIndex])
    return weights

画出数据集和Logistic回归最佳拟合直线

#画出数据集和Logistic回归最佳拟合直线函数
def plotBestFit(weights,dataArr,labelMat):
    import matplotlib.pyplot as plt
    dataMat,labelMat=dataArr,labelMat
    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(-4.0, 4.0, 0.1)

    # 最佳拟合直线函数,x轴代表特征X1（回归系数weights[0]）,y轴代表特征X2（回归系数weights[2]），特征X0=1
    # 梯度函数取值为0.5时，拟合直线是两个类别的分界线，因此W0*X0+W1*X1+...+Wn*Xn=0
    y = (-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x, y)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

dataArr,labelMat=loadDataSet('data/practice5/testSet.txt')
Weights=stocGradAscent1(dataArr,labelMat,numIter=500)
plotBestFit(Weights,dataArr,labelMat)

输出结果为：

用Logistic回归进行分类

现有实例：使用Logistic回归预测患有疝气病的马是否存活
数据集：
测试数据：horseColicTest.txt
训练数据horseColicTraining.txt：299个样本和20个特征，若对应的Sigmoid函数值大于0.5就预测类别标签为1，否则为0。

#Logistic回归分类函数
def classifyVector(inX,weights):
    prob=sigmoid0(sum(inX*weights))
    if prob>0.5 :
        return 1
    else:
        return 0
def coliTest():
    frTrain=open('data/practice5/horseColicTraining.txt')
    frTest=open('data/practice5/horseColicTest.txt')
    trainingSet=[]
    trainingLabels=[]
    for line in frTrain.readlines():
        currLine=line.strip().split('\t')
        lineArr=[]
        for i in range(len(currLine)-1):
            lineArr.append(currLine[i])
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[-1]))
    trainWeights=stocGradAscent1(trainingSet,trainingLabels)
    #存储最佳回归系数
    import pickle
    file=open('data/practice5/trainWeights.txt',mode='wb')
    pickle.dump(trainWeights,file)
    file.close()
    errorCount=0
    numTestVec=0.0
    for line in frTest.readlines():
        numTestVec+=1.0
        currLine=line.strip().split('\t')
        lineArr=[]
        for i in range(len(currLine)-1):
            lineArr.append(currLine[i])
        if int(classifyVector(lineArr,trainWeights))!=int(currLine[-1]):
            errorCount+=1
    errorate=float(errorCount)/numTestVec
    return errorate
 #取十次测试的错误率的平均值
def multiTest():
    numTests=10
    errorSum=0.0
    for k in range(numTests):
        errorSum+=coliTest()
    return errorSum/float(numTests)
if __name__=='__main__':
    inX=[]
    print('Logistic回归分类器的错误率为：'.format(multiTest()))
    for i in range(21):
        value=input('请输入特征值{}'.format(i+1))
        inX.append(float(value))
    import pickle
    file=open('data/practice5/trainWeights.txt',mode='rb')
    trainWeights=pickle.load(file)
    file.close()
    if int(classifyVector(inX,trainWeights))==1:
        print('该马仍存活')
    else:
        print('该马未存活')