logistic回归算法原理及python实现

1 logistic回归与sigmoid函数

考虑如下线性函数：

$$\begin{equation} y=\pmb{w}^\rm{T}\pmb{x}+b \tag{1} \end{equation}$$
输出 $y$为连续的实值，如何让输出成为二值来完成二分类任务？即$y \in \{0,1\}$,最理想的是单位阶跃函数即：
$$\begin{equation} y= \left\{ \begin{array}{lll} 0,\qquad z<0& \\ 0.5,\quad z=0 & \\ 1,\qquad z>0&\\ \end{array} \tag{2} \right. \end{equation}$$
但是，单位阶跃函数不连续，不利于求解权值，构建模型。于是sigmoid函数（对数几率函数，logistic function）出现了，他单调可微，并且形似阶跃函数，其公式描述如下所示:
$$\begin{equation} y=\frac{1}{1+e^{-(\pmb{w}^\rm{T}\pmb{x}+b )}}\tag{3} \end{equation}$$
sigmoid曲线如下图所示：

令$y$表示当输入为$x$时，输出为正例的概率（可能性），即$y=P(Y=1|X=x)$,$1-y$表示当输入为$x$时，输出为反例的概率（可能性），即$1-y=P(Y=0|X=x)$，两者的比值$\frac{y}{1-y}$称为几率（odds），对其取对数即达到对数几率，所以logistic回归又称为对数几率回归。因此根据（2）可得对数几率回归（logistic回归）公式如下所示：

$$\begin{equation} \log \frac{y}{1-y}=\pmb{w}^\rm{T}\pmb{x}+b \tag{4} \end{equation}$$

2 二项逻辑斯蒂回归模型与参数估计

由式（3）可得二项逻辑斯蒂回归模型如下所示：

$$\begin{equation} \log \frac{P(Y=1|X=x)}{1-P(Y=1|X=x)}=\pmb{w}^\rm{T}\pmb{x} \tag{5} \end{equation}$$
其中， $\pmb{w}=(w_1,w_2,...,w_m,b)^\rm{T}$， $\pmb{x}=(x_1,x_2,...,x_m,1)$
学习模型的关键是对权值 $\pmb{w}$的学习，已知的是训练样本即输入及其对应的标签，利用已知输入样本来如何学习权值？该学习过程可以转化为带约束的最优化问题，或者以极大似然函数为目标函数（策略）并使用梯度上升或者牛顿法等最优化算法。下边以后者为学习方法，来求解逻辑斯蒂回归模型。
极大似然函数的假设是：训练样本出现的概率最大。换句话所就是，有些事情具有多种可能，而其中一种可能值出现，其他可能值未出现，在这个过程中，出现的可能值具有较大概率，所以才会出现。
一种学习方法的假设很重要，合理、科学的假设代表了学习方法的正确方向，在该假设条件下，得出的模型往往能够达到预期效果。
设训练样本 $\{\pmb X,\pmb{y}\}$, $\pmb{X}=\{\pmb{x_j}\},\pmb{x_j} \in R^n,\pmb{y}\in R^n,y_i\in \{0,1\},i=1,2,...,n,j=1,2,...,m$则逻辑斯蒂回归输出 $\hat{y}=\frac{1}{1+e^{-(\pmb{w}^\rm{T}\pmb{x})}}\in (0,1)$为区间在0和1的连续实值（表示概率）。则 样本的似然函数为：
$$\begin{equation} l(\pmb{w})=\prod\limits_{i=1}^{n}\hat{y}_i^{y_i}(1-\hat{y}_i)^{(1-y_i)} \tag{6} \end{equation}$$
对数似然函数为：
$$\begin{equation} L(\pmb{w})=\sum\limits_{i=1}^{n}(y_i\log\hat{y}_i+(1-y_i)log(1-\hat{y_i}))\\ \quad =\sum\limits_{i=1}^{n}(y_i\log\frac{\hat{y}_i}{(1-\hat{y_i})}+log(1-\hat{y_i}))\\ \qquad =\sum\limits_{i=1}^{n}(y_i\pmb{w}^\rm{T}\pmb{x}_i-log(1+e^{(\pmb{w}^\rm{T}\pmb{x}_i)}))\tag{7} \end{equation}$$
则逻辑斯蒂回归模型学习可转化为如下最优化问题：
$$\begin{equation} \max\limits_{\pmb{w}}L(\pmb{w}) \qquad \tag{8} \end{equation}$$
采用梯度上升算法来求解函数的最大值（梯度下降求解函数的最小值）：
对式（7）对权值求偏导得如下公式：
$$\begin{equation} \nabla_{\pmb{w}}=\left[ \begin{matrix} \nabla_{w_1}\\ \nabla_{w_2}\\ \vdots\\ \nabla_{w_m}\\ \end{matrix} \right] =\frac{\partial L(\pmb{w})}{\partial \pmb{w}} =\sum\limits_{i=1}^{n}(y_i\pmb{x}_i-\frac{1}{1+e^{(\pmb{w}^\rm{T}\pmb{x}_i)}}e^{(\pmb{w}^\rm{T}\pmb{x}_i)}\pmb{x}_i)\\ =\sum\limits_{i=1}^{n}(y_i-\frac{1}{1+e^{-(\pmb{w}^\rm{T}\pmb{x}_i)}})\pmb{x}_i=\sum\limits_{i=1}^{n}(y_i-\hat{y}_i)\pmb{x}_i\\ =\pmb{X}^\rm{T}(\pmb{y}-\hat{\pmb{y}}) \tag{9} \end{equation}$$
在此需注意到： $\pmb{y}-\hat{\pmb{y}}$为误差向量，梯度上升算法的迭代公式如下所示：
$$\begin{equation} \pmb{w}:= \pmb{w}+\alpha \nabla_{\pmb{w}}\tag{10} \end{equation}$$
其中， $\alpha$为步长因子，需人为给定， $\pmb{w}$的初始值一般设置为 $[-0.01,0.01]$之间。梯度下降算法为： $\pmb{w}:= \pmb{w}-\alpha \nabla_{\pmb{w}}$。

3 利用python根据梯度上升算法和随机梯度上升算法实现二项逻辑斯蒂回归

梯度上升算法在每次更新回归系数时都要遍历整个数据集，当数据集有亿级样本时，计算复杂度将会很大，为解决此问题，提出了随机梯度上升算法，该算法每次只用一个样本来更新回归系数，而且可以在新的样本来临时对分类器进行增量式的更新，因此随机梯度上升算法是一个在线学习算法。与在线学习相对应的是一次处理所有数据集的称为“批处理”。下图为梯度上升算法、随机梯度上升算法以及改进的随机梯度上升算法的流程图。

说明：
1）学习步长一般设置为alpha=0.01，过大会影响分类精度
2）改进的随机梯度算中学习步长更新公式alpha = 4/(1+i+j) + 0.01中的常数项可根据具体情况进行修改。
3）改进的随机梯度上升算法具有随机性，每次得到的分类模型精度不确定

3.1 python实现程序

1. 测试数据生成程序

# -*- coding: utf-8 -*-
"""
Created on Wed Jun 07 14:33:08 2017
生成test数据集
@author: liujiping
"""

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(1)

y1_x1_mu,y1_x1_sigma = 0,2 #均值与标准差
y1_x2_mu,y1_x2_sigma = 3,3.5
y1_x1 = np.random.normal(y1_x1_mu,y1_x1_sigma,50)
y1_x2 = np.random.normal(y1_x2_mu,y1_x2_sigma,50)
y1 = np.ones((50,1))

y0_x1_mu,y0_x1_sigma = 0,1 #均值与标准差
y0_x2_mu,y0_x2_sigma = 10,3
y0_x1 = np.random.normal(y1_x1_mu,y1_x1_sigma,50)
y0_x2 = np.random.normal(y0_x2_mu,y0_x2_sigma,50)
y0 = np.zeros((50,1))

x1 = list(y1_x1)
x1.extend(list(y0_x1))

x2 = list(y1_x2)
x2.extend(list(y0_x2))

y = list(y1)
y.extend(list(y0))
#------保存数据 ------
with open("logisticData.txt","w") as f:#使用with不需要f.close()
    for i in range(100):
        write_str = '%f %f %d\n'%(x1[i],x2[i],y[i])
        f.write(write_str)
#------画出散点图  ------       
line1,= plt.plot(x1[:50],x2[:50],'ro',label = 'class1')
line2, = plt.plot(x1[50:],x2[50:],'b^',label ='class0')
plt.xlabel('x1')
plt.ylabel('x2')
plt.legend(handles=[line1,line2],loc = 2)
plt.show()        

1. logistic回归实现程序

import numpy as np
import matplotlib.pyplot as plt
def loadDataSet():
    '''
    输入：无
    功能：读取txt中数据并将输入和标签分别存储
    输出：输入数据dataMat,标签数据labelMat
    '''
    dataMat = [];labelMat = []
    with open('logisticData.txt') as f:
        for line in f.readlines():
            lineList = line.strip().split()
            dataMat.append([1.0,float(lineList[0]),float(lineList[1])])
            labelMat.append(int(lineList[2]))
    return dataMat,labelMat
def sigmoid(z):
    '''
    输入：加权的输入数据w*x
    功能：执行sigmoid变换
    输出：sigmoid变换值，值域（0,1）
    '''
    return 1.0/(1+np.exp(-z))

def gradientAscent(dataMat,labelMat):
    '''
    输入：输入特征矩阵，标签列表
    功能：批处理梯度上升算法，更新权重
    输出：权重向量
    '''
    dataMatrix = np.mat(dataMat)
    labelMatrix = np.mat(labelMat).transpose()
    n,m = np.shape(dataMatrix)
    alpha = 0.01#梯度算法的步长，可以控制收敛的速度以及最后模型的精度
    maxCycles = 500#批处理，权值跟新的最大次数
    weights = np.ones((m,1))*0.01 #初始化权值，权值个数等于特征个数（包括常数项1）
    for k in range(maxCycles):
        predictLabel = sigmoid(dataMatrix*weights)
        error = (labelMatrix - predictLabel)
        #计算梯度
        gradient = dataMatrix.transpose() * error
        #更新权重
        weights = weights +alpha * gradient
    return weights,error
def plotDataSet(dataMat):
    '''
    输入：从txt文本文件里读取的输入数据
    功能：画出前两个特征的二维图
    输出：散点图
    '''
    x1 = np.mat(dataMat)[:,1]
    x2 = np.mat(dataMat)[:,2]
    line1,= plt.plot(x1[:50],x2[:50],'ro',label = 'class1')
    line2, = plt.plot(x1[50:],x2[50:],'b^',label ='class0')
    plt.xlabel('x1')
    plt.ylabel('x2')
    plt.legend(handles=[line1,line2],loc = 2)
    #plt.show() 

def plotBestFit(dataMat,weights):
    '''
    输入：输入数据,权值矩阵
    功能：画出前两个特征的二维图及分类曲线
    输出：散点图
    '''    
    plt.figure()
    plotDataSet(dataMat)
    x = np.mat(np.arange(-4.0,4.0,0.1))
    y = (-weights[0]-weights[1] * x)/weights[2]
    plt.plot(x.transpose(),y.transpose())
    plt.show()

def logisticTest(weights,testData):
    '''
    输入：权值，测试数据
    功能：利用训练的数据计算测试数据对应的标签
    输出：测试数据的标签
    '''
    n,m = np.shape(np.mat(testData))
    testLabel = np.zeros((n,1))
    for i in range(n):
        testLabel[i] = weights[0]+weights[1]*testData[i][0]+weights[2]*testData[i][1]
        if testLabel[i] >= 0.5: testLabel[i] = 1
        else:  testLabel[i] = 0
    return testLabel

def stoGradAscent(dataMat,labelMat):
    '''
    输入：输入特征矩阵，标签列表
    功能：随机梯度上升算法，更新权重，只用了一遍数据集
    输出：权重向量
    '''
    dataMatrix = np.mat(dataMat)
    n,m = np.shape(dataMatrix)
    alpha = 0.01
    weights =np.mat(np.ones((m,1)))
    for i in range(n):
        predictlabel = sigmoid(dataMatrix[i] * weights)
        error = labelMat[i] - predictlabel
        #计算梯度
        gradient = np.multiply(dataMatrix[i],error)
        #更新权重
        weights = weights + alpha * gradient.transpose()        
    return weights

def improvedStoGradAscent(dataMat,labelMat,numInter = 150):   
    '''
    输入：输入特征矩阵，标签列表,最大迭代次数（决定了所有训练样本被使用多少次）
    功能：改进的随机梯度上升算法，更新权重，随机选取训练样本中的数据
    输出：权重向量
    '''
    dataMatrix = np.mat(dataMat)
    n,m = np.shape(dataMatrix)
    weights =np.mat(np.ones((m,1)))    
    for j in range(numInter):
        dataIndex = range(n)
        for i in range(n):
            #修改学习步长，缓解数据波动,由于常数项的存在，alpha不是严格下降的
            #alpha =  0.01
            alpha = 2/(1.0+j+i) + 0.0001
            #获得随机样本索引
            randIndex = int(np.random.uniform(0,len(dataIndex)))
            predictlabel = sigmoid(dataMatrix[randIndex] * weights)
            error = labelMat[randIndex] - predictlabel
            gradient = np.multiply(dataMatrix[randIndex],error)
            weights = weights + alpha * gradient.transpose()

            del dataIndex[randIndex]
    return weights

if __name__ == '__main__':
    dataMat,labelMat = loadDataSet()
    #plt.figure(1)
    #plotDataSet(dataMat)  
    #------批处理梯度上升算法------
    #weights,error = gradientAscent(dataMat,labelMat)
    #-----随机梯度上升算法---
    #weights = stoGradAscent(dataMat,labelMat)
    #-----改进的随机梯度上升算法---
    weights = improvedStoGradAscent(dataMat,labelMat,numInter = 150)    
    plotBestFit(dataMat,weights)
    testLabel = logisticTest(weights,[[0,0],[0,10]])

1. 程序运行结果
   1）批处理梯度上升算法结果图及测试数据分类结果
   
   

测试分类结果

testLabel
Out[2]: 
array([[ 1.],
       [0.]])

2）随机梯度上升算法结果图及测试数据分类结果

测试数据分类结果

testLabel
Out[4]: 
array([[ 1.],
       [0.]])

3）改进的随机梯度上升算法结果图及测试数据分类结果

测试数据分类结果

testLabel
Out[9]: 
array([[ 1.],
       [0.]])