logistic回归算法原理及python实现

1 logistic回归与sigmoid函数

考虑如下线性函数：

y=wwTxx+b(1)

输出y为连续的实值，如何让输出成为二值来完成二分类任务？即y∈{0,1},最理想的是单位阶跃函数即：

y=⎧⎩⎨⎪⎪0,z<00.5,z=01,z>0(2)

但是，单位阶跃函数不连续，不利于求解权值，构建模型。于是sigmoid函数（对数几率函数，logistic function）出现了，他单调可微，并且形似阶跃函数，其公式描述如下所示:

y=11+e−(wwTxx+b)(3)

sigmoid曲线如下图所示：



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

logy1−y=wwTxx+b(4)

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

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

logP(Y=1|X=x)1−P(Y=1|X=x)=wwTxx(5)

其中，ww=(w1,w2,...,wm,b)T，xx=(x1,x2,...,xm,1)

学习模型的关键是对权值ww的学习，已知的是训练样本即输入及其对应的标签，利用已知输入样本来如何学习权值？该学习过程可以转化为带约束的最优化问题，或者以极大似然函数为目标函数（策略）并使用梯度上升或者牛顿法等最优化算法。下边以后者为学习方法，来求解逻辑斯蒂回归模型。

极大似然函数的假设是：训练样本出现的概率最大。换句话所就是，有些事情具有多种可能，而其中一种可能值出现，其他可能值未出现，在这个过程中，出现的可能值具有较大概率，所以才会出现。

一种学习方法的假设很重要，合理、科学的假设代表了学习方法的正确方向，在该假设条件下，得出的模型往往能够达到预期效果。

设训练样本{XX,yy},XX={xjxj},xjxj∈Rn,yy∈Rn,yi∈{0,1},i=1,2,...,n,j=1,2,...,m则逻辑斯蒂回归输出y^=11+e−(wwTxx)∈(0,1)为区间在0和1的连续实值（表示概率）。则样本的似然函数为：

l(ww)=∏i=1ny^yii(1−y^i)(1−yi)(6)

对数似然函数为：

L(ww)=∑i=1n(yilogy^i+(1−yi)log(1−yi^))=∑i=1n(yilogy^i(1−yi^)+log(1−yi^))=∑i=1n(yiwwTxxi−log(1+e(wwTxxi)))(7)

则逻辑斯蒂回归模型学习可转化为如下最优化问题：

maxwwL(ww)(8)

采用梯度上升算法来求解函数的最大值（梯度下降求解函数的最小值）：

对式（7）对权值求偏导得如下公式：

∇ww=⎡⎣⎢⎢⎢⎢⎢∇w1∇w2⋮∇wm⎤⎦⎥⎥⎥⎥⎥=∂L(ww)∂ww=∑i=1n(yixxi−11+e(wwTxxi)e(wwTxxi)xxi)=∑i=1n(yi−11+e−(wwTxxi))xxi=∑i=1n(yi−y^i)xxi=XXT(yy−yy^)(9)

在此需注意到：yy−yy^为误差向量，梯度上升算法的迭代公式如下所示：

ww:=ww+α∇ww(10)

其中，α为步长因子，需人为给定，ww的初始值一般设置为[−0.01,0.01]之间。梯度下降算法为：ww:=ww−α∇ww。

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

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



说明：

1）学习步长一般设置为alpha=0.01，过大会影响分类精度

2）改进的随机梯度算中学习步长更新公式alpha = 4/(1+i+j) + 0.01中的常数项可根据具体情况进行修改。

3）改进的随机梯度上升算法具有随机性，每次得到的分类模型精度不确定

3.1 python实现程序

测试数据生成程序

# -*- 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()

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）批处理梯度上升算法结果图及测试数据分类结果



测试分类结果

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

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



测试数据分类结果

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

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



测试数据分类结果

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