机器学习之Logistic回归与Python实现

logistic回归是一种广义的线性回归，通过构造回归函数，利用机器学习来实现分类或者预测。

一 Logistic回归概述

Logistic回归的主要思想是，根据现有的数据对分类边界建立回归公式，从而实现分类（一般两类）。“回归”的意思就是要找到最佳拟合参数，其中涉及的数学原理和步骤如下：
（1）需要一个合适的分类函数来实现分类【单位阶跃函数、Sigmoid函数】
（2）损失函数（Cost函数）来表示预测值（h(x)）与实际值(y)的偏差(h−y),要使得回归最佳拟合，那么偏差要尽可能小（偏差求和或取均值）。
（3）记J(ω)表示回归系数为ω时的偏差，那么求最佳回归参数ω就转换成了求J(ω)的最小值。【梯度下降法】
所以，接下来就围绕这几个步骤进行展开。

1.1 分类函数

假设要实现二分类，那么可以找一个函数，根据不同的特征变量，输出0和1，并且只输出0和1，这种函数在某个点直接从0跳跃到1，如：
 
但是这种函数处理起来，稍微有点麻烦，我们选择另外一个连续可导的函数，也就是Sigmoid函数,函数的公式如下：

这个函数的特点是，当x=0.5时，h(x)=0.5，而x越大，h(x)越接近1，x越小，h(x)越接近0。函数图如下：

这个函数很像阶跃函数，当x>0.5，就可以将数据分入1类；当x<0.5，就可以将数据分入0类。

确定了分类函数，接下来，我们将Sigmoid函数的输入记为z，那么

向量x是特征变量，是输入数据，向量w是回归系数是特征
之后的事情就是如何确定最佳回归系数ω(w0,w1,w2,...,wn)

1.2 Cost函数

现有

对于任意确定的x和w,有：

这个函数可以写成：

取似然函数：

求对数似然函数：

因此，就构造得到了函数J(w)来表示预测值与实际值的偏差，于是Cost函数也可以写成：

所以，我们可以用J(w)来表示预测值与实际值的偏差，也就是Cost函数，接下里的任务，就是如何让偏差最小，也就是J(w)最大

Question:为什么J(w)可以表示预测值与实际值的大小，为什么J(w)最大表示偏差最小。

我们回到J(w)的推导来源，来自
P(y=1|x,w)=hw(x)和P(y=0|x,w)=1−hw(x)，
那么显然有
当x>0,此时y=1，1/2<hw(x)<1，所以P(y=1|x,w)=hw(x)>1/2
当x<0,此时y=0，0<hw(x)<1/2，所以P(y=0|x,w)=1−hw(x)>1/2
所以，无论如何取值，P(y=0|x,w)都大于等于1/2,P(y=0|x,w)越大，越接近1，表示落入某分类概率越大，那么分类越准确，预测值与实际值差异就越小。
所以P(y=0|x,w)可以表示预测值与实际值的差异，且P(y=0|x,w)越大表示差异越小，所以其似然函数J(w)越大，预测越准确。

所以，接下来的任务，是如何求解J(w)最大时的w值，其方法是梯度上升法。

1.3 梯度上升法求J(w)最大值

梯度上升法的核心思想是：要找某个函数的最大值，就沿着这个函数梯度方向探寻，如果梯度记为∇，那么函数f(x,y)的梯度是：

梯度上升法中，梯度算子沿着函数增长最快的方向移动（移动方向），如果移动大小为α（步长），那么梯度上升法的迭代公式是：

问题转化成：

首先，我们对J(w)求偏导：

在第四至第五行的转换，用到的公式是：

将求得的偏导公式代入梯度上升法迭代公示：

可以看到，式子中所有函数和输入的值，都是已知的了。接下来，可以通过Python实现Logistic回归了。

二、Python算法实现

2.1 梯度上升法求最佳回归系数

首先，数据取自《机器学习实战》中的数据，部分数据如下：

-0.017612   14.053064   0
-1.395634   4.662541    1
-0.752157   6.538620    0
-1.322371   7.152853    0
0.423363    11.054677   0
0.406704    7.067335    1

先定义函数来获取数去，然后定义分类函数Sigmoid函数，最后利用梯度上升法求解回归系数w。
建立一个logRegres.py文件，输入如下代码：

from numpy import *
#构造函数来获取数据
def loadDataSet():
    dataMat=[];labelMat=[]
    fr=open('machinelearninginaction/Ch05/testSet.txt')
    for line in fr.readlines():
        lineArr=line.strip().split()
        dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])#特征数据集，添加1是构造常数项x0
        labelMat.append(int(lineArr[-1]))#分类数据集
    return dataMat,labelMat

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

def gradAscent(dataMatIn,classLabels):
    dataMatrix=mat(dataMatIn) #(m,n)
    labelMat=mat(classLabels).transpose() #转置后(m,1)
    m,n=shape(dataMatrix)
    weights=ones((n,1)) #初始化回归系数，(n,1)
    alpha=0.001 #定义步长
    maxCycles=500 #定义最大循环次数
    for i in range(maxCycles):
        h=sigmoid(dataMatrix * weights) #sigmoid 函数
        error=labelMat - h #即y-h，（m,1）
        weights=weights + alpha * dataMatrix.transpose() * error #梯度上升法
    return weights

在python命令符中输入代码对函数进行测试：
In [8]: import logRegres
   ...:
In [9]: dataArr,labelMat=logRegres.loadDataSet()
   ...:
In [10]: logRegres.gradAscent(dataArr,labelMat)
    ...:
Out[10]:
matrix([[ 4.12414349],
        [ 0.48007329],
        [-0.6168482 ]])

于是得到了回归系数。接下来根据回归系数画出决策边界wTx=0
定义作图函数：

def plotBestFit(weights):
    import matplotlib.pyplot as plt
    dataMat,labelMat=loadDataSet()
    n=shape(dataMat)[0]
    xcord1=[];ycord1=[]
    xcord2=[];ycord2=[]
    for i in range(n):
        if labelMat[i]==1:
            xcord1.append(dataMat[i][1])
            ycord1.append(dataMat[i][2])
        else:
            xcord2.append(dataMat[i][1])
            ycord2.append(dataMat[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,3,0.1)
    y=(-weights[0,0]-weights[1,0]*x)/weights[2,0] #matix
    ax.plot(x,y)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

在Python的shell中对函数进行测试：

In [11]: weights=logRegres.gradAscent(dataArr,labelMat)

In [12]: logRegres.plotBestFit(weights)
    ...:

2.2 算法改进

(1) 随机梯度上升
上述算法，要进行maxCycles次循环，每次循环中矩阵会有m*n次乘法计算，所以时间复杂度（开销）是maxCycles*m*n，当数据量较大时，时间复杂度就会很大。因此，可以是用随机梯度上升法来进行算法改进。

随机梯度上升法的思想是，每次只使用一个数据样本点来更新回归系数。这样就大大减小计算开销。
代码如下：

def stocGradAscent(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] #array 和list矩阵乘法不一样
    return weights

注意：gradAscent函数和这个stocGradAscent函数中的h和weights的计算形式不一样，因为
前者是的矩阵的计算，类型是numpy的matrix，按照矩阵的运算规则进行计算。
后者是数值计算，其类型是list，按照数值运算规则计算。

对随机梯度上升算法进行测试：

In [37]: dataMat,labelMat=logRegres.loadDataSet()
    ...:
In [38]: weights=logRegres.stocGradAscent(array(dataMat),labelMat)
    ...:
In [39]: logRegres.plotBestFit(mat(weights).transpose())
    ...:

输出的样本数据点和决策边界是：

（2）改进的随机梯度上升法

def stocGradAscent1(dataMatrix,classLabels,numIter=150):
    m,n=shape(dataMatrix)
    weights=ones(n)
    for j in range(numIter):
        dataIndex=list(range(m))
        for i in range(m):
            alpha=4/(1+i+j)+0.01#保证多次迭代后新数据仍然具有一定影响力
            randIndex=int(random.uniform(0,len(dataIndex)))#减少周期波动
            h=sigmoid(sum(dataMatrix[randIndex] * weights))
            error=classLabels[randIndex]-h
            weights=weights + alpha*dataMatrix[randIndex]*error
            del(dataIndex[randIndex])
    return weights 

在Python命令符中测试函数并画出分类边界：

In [188]: weights=logRegres.stocGradAscent1(array(dataMat),labelMat)
     ...:
In [189]: logRegres.plotBestFit(mat(weights).transpose())
     ...:

（3）三种方式回归系数波动情况
普通的梯度上升法：
 
随机梯度上升：
 
改进的随机梯度上升

评价算法优劣势看它是或否收敛，是否达到稳定值，收敛越快，算法越优。

三 实例

3.1 通过logistic回归和氙气病症预测马的死亡率

数据取自《机器学习实战》一书中的氙气病症与马死亡的数据，部分数据如下：

2.000000    1.000000    38.500000   66.000000   28.000000   3.000000    3.000000    0.000000    2.000000    5.000000    4.000000    4.000000    0.000000    0.000000    0.000000    3.000000    5.000000    45.000000   8.400000    0.000000    0.000000    0.000000
1.000000    1.000000    39.200000   88.000000   20.000000   0.000000    0.000000    4.000000    1.000000    3.000000    4.000000    2.000000    0.000000    0.000000    0.000000    4.000000    2.000000    50.000000   85.000000   2.000000    2.000000    0.000000
2.000000    1.000000    38.300000   40.000000   24.000000   1.000000    1.000000    3.000000    1.000000    3.000000    3.000000    1.000000    0.000000    0.000000    0.000000    1.000000    1.000000    33.000000   6.700000    0.000000    0.000000    1.000000

通过21个特征数据，来对结果进行分类和预测。

#定义分类函数，prob>0.5，则分入1，否则分类0
def classifyVector(inX,trainWeights):
    prob=sigmoid(sum(inX*trainWeights))
    if prob>0.5:return 1
    else : return 0

def colicTest():
    frTrain = open('machinelearninginaction/Ch05/horseColicTraining.txt')#训练数据
    frTest = open('machinelearninginaction/Ch05/horseColicTest.txt')#测试数据
    trainSet=[];trainLabels=[]
    for line in frTrain.readlines():
        currLine=line.strip().split('\t')
        lineArr=[]
        for i in range(21):
            lineArr.append(float(currLine[i]))
        trainSet.append(lineArr)
        trainLabels.append(float(currLine[21]))
    trainWeights=stocGradAscent1(array(trainSet),trainLabels,500)#改进的随机梯度上升法
    errorCount=0;numTestVec=0
    for line in frTest.readlines():
        numTestVec+=1
        currLine=line.strip().split('\t')
        lineArr=[]
        for i in range(21):
            lineArr.append(float(currLine[i]))
        if classifyVector(array(lineArr),trainWeights)!=int(currLine[21]):
            errorCount+=1
    errorRate=(float(errorCount)/numTestVec)
    print('the error rate of this test is :%f'%errorRate)
    return errorRate
def multiTest():#进行多次测试
    numTests=10;errorSum=0
    for k in range(numTests):
        errorSum+=colicTest()
    print('after %d iterations the average error rate is:%f'%(numTests,errorSum/float(numTests)))

在控制台命令符中输入命令来对函数进行测试：

In [3]: logRegres.multiTest()
G:\Workspaces\MachineLearning\logRegres.py:19: RuntimeWarning: overflow encountered in exp
  return 1/(1+exp(-inX))
the error rate of this test is :0.313433
the error rate of this test is :0.268657
the error rate of this test is :0.358209
the error rate of this test is :0.447761
the error rate of this test is :0.298507
the error rate of this test is :0.373134
the error rate of this test is :0.358209
the error rate of this test is :0.417910
the error rate of this test is :0.432836
the error rate of this test is :0.417910
after 10 iterations the average error rate is:0.368657

分类的错误率是36.9%。