《机器学习实战》逻辑回归算法简单Python实现

Logistic回归

假设有一些数据点，对这些数据点用一条直线进行拟合，这个拟合的过程称为回归。

使用逻辑回归进行二分类的思想是：根据现有的数据对分类边界线建立回归公式依此进行分类！

原理公式——求最佳的特征参数θ

分类可知真实值只有两个取值0和1，那么需要将计算得到的函数值转换为在0—1之间的数值，于是用了Sigmoid函数：

g(z)=1/(1+e^-z)

将参数换成一个线性回归公式得到

h(x)=g(θ.T·X）=1/（1+e^-θ.T·X）

01分类二项分布得到数据为正集的概率表达式为

对这个概率求其最大似然函数，并对最大似然函数求其对数将连乘转为求和

目的希望最大似然函数越大越好，求使该最大似然函数取得最大值的参数θ。可以使用梯度下降算法，但梯度下降算法求得的是最小值，因此将其取负号获得代价函数使用梯度下降求使其最小的参数θ（逆向思维我们就获得了使似然函数最大的参数θ）！

梯度下降算法的代价函数找到了，通过梯度下降的原理知求θ迭代更新的公式需要代价函数的偏导数，下面进行求J（θ）的偏导

Python代码实现

运行环境在Anaconda下的JupyterNotebook,使用《机器学习实战》第五章中的数据集。因为Python版本等原因，代码在原书的实现代码上做出了一些改动。

先实现数据集的读取，得到特征变量和标记

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
def loadDataSet():          #数据的读取，从txt文档中读入          
    dataMat = []   
    labelMat = []
    fr = open('logistic_02.txt')
    for line in fr.readlines():
        lineArr = line.strip().split() #split将一个字符串进行切片，默认分隔符为空字符
        if(len(lineArr)==0):  #line从文档头开始一行行读取，必须加上判断条件在读到末尾的下一行时跳出，否则会出现列表越界！
            break
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])]) #特征值x0,x1,x2，其中x0为偏移量设为1
        labelMat.append(int(lineArr[2]))                            #真实值（标记，分为0和1）
    return dataMat, labelMat

训练算法核心是实现以梯度下降算法为基础的参数θ更新。θ的迭代更新是需要一个停止终点的，方法有三种：1.迭代的次数限制 2.按照损失函数的精度 3.按照梯度下降的精度。这里使用最简单的第一种方法。

def sigmoid(z):                         
    return 1.0/(1+np.exp(-z))                                           
#返回最佳回归系数theta
def gradAscent(dataMatIn, classLabels):
    #将得到的列表转换为numpy矩阵类型
    dataMatrix = np.mat(dataMatIn) 
    labelMat = np.mat(classLabels).transpose()
    # m->样本数 n->特征数
    m= dataMatrix.shape[0]
    n=dataMatrix.shape[1]
    alpha = 0.001   #步长，学习率
    maxCycles = 500 #限制迭代次数
    #初始化theta参数，每个维度均为1.0
    theta = np.ones((n,1))
    for i in range(maxCycles):
        h = sigmoid(dataMatrix * theta)  #得到h(x)
        error = (h-labelMat)             #得到h(x)-y
        theta = theta - alpha * dataMatrix.transpose() *error  
        #dataMatrix.transpose()*error是代价函数的偏导,返回得到的特征参数θ0，θ1，θ2
    return np.array(theta)

使用得到的参数画出决策边界，进行直观比较

def plotBestFit(theta):              
    dataMat, labelMat = loadDataSet()
    dataArr = np.array(dataMat)
    # n->样本数
    n = dataArr.shape[0]
    #xcord1,ycord1代表正例特征
    #xcord2,ycord2代表负例特征
    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 = 'r', marker = 'o')  #画出各自的散点图
    ax.scatter(xcord2, ycord2, s = 30, c = 'b',marker='x')
    #设定边界直线x和y的值
    x = range(-3,3)
    """
    决策边界直线方程式
    w0*x0+w1*x1+w2*x2=f(x)
    f(x)=0是两个类别的分界处
    所以： w0+w1*x+w2*y=0 => y = (-w0-w1*x)/w2    
    """
    y = (-theta[0] - theta[1] * x) / theta[2]
    ax.plot(x, y)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

看出该特征参数已经比较准确了，可以使用测试集矩阵与参数矩阵乘法运算得到预测的真实值，该数值的大小反应该数据为正集的概率大小。

训练算法的优化

上面的训练算法称为批量梯度下降算法，是一种批处理（Batch learning)，因为每一次更新参数θ时都会使用整个数据集，当数据集大时计算量会十分庞大。该算法在接收了大量数据集后开始运算并且不再继续接收其它数据学习改进，因此是一种线下学习（offline learning）的算法。

因此可以使用随机梯度下降算法，它是一种在线学习的算法，关于在线学习这一段英文解释的很好：

In online learning, you train the system incrementally by feeding it data instances sequentially, either individually or by small groups called mini-batches. Each learning step is fast and cheap, so the system can learn about new data on the fly, as it arrives.Online learning is great for systems that receive data as a continuous flow and need to adapt to change rapidly or autonomously. It is also a good option if you have limited computing resources.

使用改进的随机梯度下降算法

与批量梯度下降不同，随机梯度下降每次都使用一个随机的样本来更新回归系数

def GradAscent1(dataArr, classLabels, numIter = 150):
    dataMatrix=np.mat(dataArr)    #将得到的列表转换为矩阵
    #m为样本数n为特征数
    m=dataMatrix.shape[0]
    n=dataMatrix.shape[1]
    theta = np.ones((n,1))        #初始化特征参数矩阵
    # 随机梯度, 缺省则默认循环150次
    for j in range (numIter):
        dataIndex = range(m)
        for i in range(m):    #随机选取一个数据来更新theta，使用过后从在本次迭代中删除不再使用
            #alpha会随迭代次数不断减小，步长不断减小移动更精确
            alpha = 4 / (4.0+i+j) + 0.01
            # random.uniform(x, y) 方法将返回随机[x,y]范围内的一个实数。
            randIndex = int(rd.uniform(0,len(dataIndex)))
            #随机选取一个更新特征参数θ
            h = sigmoid(np.dot(dataMatrix[randIndex],theta))
            error = h-classLabels[randIndex]
            theta = theta - alpha  * dataMatrix[randIndex].T*error
            del (list(dataIndex)[randIndex])   #从列表中删除该值
    return theta

相比于批量梯度下降，随机梯度下降主要在两处进行了改进：

第一处是学习率在每次迭代的时候都会调整，会随着迭代次数不断变小，这样可以缓解数据波动或高频波动。

第二处是通过随机选取样本来更新回归系数，这样可以减少周期性的波动。

dataMat,labelMat=loadDataSet()
theta=GradAscent1(dataMat,labelMat)
print(theta)
theta=theta.tolist()
mytheta=np.array([[theta[0][0]],[theta[1][0]],[theta[2][0]]])
plotBestFit(mytheta)

可视化得到的数据结果：

看起来结果还不差，但是此次只在数据集上进行了150次的运算。

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

使用Logistic回归分类来预测患有疝病的马的存活问题，数据训练集中每个样本都有多个特征以及它的标记值，通过该训练集可以使用逻辑回归算法来训练得到特征参数从而进行预测，但是该数据集的问题是存在30%的值是缺失的。

在准备数据训练集中很重要的一项就是处理数据中可能存在的缺失值，在《机器学习实战》中给出了一些可选的做法：

• 使用可用特征的均值来填补缺失值；
• 使用特殊值来填补缺失值，如-1；
• 忽略有缺失值的样本；
• 使用相似样本的均值添补缺失值；
• 使用另外的机器学习算法预测缺失值。

对于该数据集的Logistic回归我们可以将缺失的数据全部处理为零：其一在梯度下降中更新公式theta=theta-alpha*datamatrix.T*error如果缺失值为零则alpha*datamatrix.T*error也为零对于回归系数不做出更新；其二是sigmoid(0)值为0.5，对于数据的分类预测不具有偏向性。

注意对于应该如何处理缺失的数据是没有标准答案的，其取决于实际应用中的需求。

import numpy as np
import random as rd
def sigmoid(z):
    return 1/(1.0+np.exp(-z))
def GradAscent(datamat,classlabels):
    datamatrix=np.mat(datamat)        #将得到的列表数据转为矩阵
    labelmat=np.mat(classlabels).transpose()
    m=datamatrix.shape[0]             #m->样本数 n->特征数
    n=datamatrix.shape[1]
    alpha=0.01         
    maxx=500
    theta=np.ones((n,1)) 
    for i in range(maxx):             #因为数据集本身很小，使用批量梯度下降即可
        h=sigmoid(datamatrix*theta)
        error=h-labelmat
        theta=theta-alpha*datamatrix.transpose()*error
    return np.array(theta)            #返回theta参数矩阵
#预测函数，把得到的特征参数与测试集数据运算，然后以0.5为界进行分类
def predict(theta,datamatrix):       
    y=sigmoid(sum(datamatrix*theta))
    y=y.tolist()
    flag=(y[0]>0.5)
    if flag:
        return 1
    else:
        return 0
def colictest():
    ftrain=open('预处理的疝气症病马训练集.txt')
    ftest=open('疝气症病马测试集.txt')
    trainset=[];trainlabels=[]
    for line in ftrain.readlines():    #从训练集逐行读取数据
        curline=line.strip().split()
        if(len(curline)==0):
            break
        linearr=[]
        for i in range(len(curline)-1):
            linearr.append(float(curline[i]))
        trainset.append(linearr)
        trainlabels.append(float(curline[-1]))
    theta=GradAscent(trainset,trainlabels)   #把得到的两个数据列表传入梯度下降计算theta
    errorcount = 0;numtestvec=0.0     #分别为预测错误个数和测试集的个数
    for line in ftest.readlines():    #逐个读取测试集进行测试
        curline=line.strip().split()
        if(len(curline)==0):
            break
        numtestvec+=1.0
        linearr=[]
        for i in range(len(curline)):
            linearr.append(float(curline[i]))
        if int(predict(theta,np.array(linearr)))!=int(curline[-1]):
            errorcount+=1
    errorrate=(float(errorcount)/numtestvec)*100
    print("测试集的错误率为：%f%%" %errorrate)

可以看出测试集的错误率在30%附近，这个结果已经比较好了，因为训练数据集中就有30%的数据缺失。当然，再通过调整学习率步长和迭代次数还可能找到错误率更小的特征参数。