Logistic虽然不是十大经典算法之一,但却是数据挖掘中常用的有力算法,所以这里也专门进行了学习,以下内容皆为亲自实践后的感悟和总结(Logistic原理、代码实现和优化、真实样例数据、sklearn实现)。为了记录的比较清楚,所以内容可能有点多,但都比较浅显,下面进入正文。(运算数据来自机器学习第5章)
Logistic原理:
大体的思路:为了实现Logistic回归分类器,我们可以在每个特征上都乘以一个回归系数,然后把所有的结果值相加,将这个总和带入Sigmoid函数中,进而得到一个范围在0~1之间的数值。任何大于0.5的数据被分入1类,小于0.5即归为0类。所以Logistic回归也可以被看成是一种概率估计。
问题来了,如何求这个回归系数呢?下面围绕的关键点也就是求这个最佳回归系数。
1:逻辑回归并非回归算法,而是分类算法,”回归“一词来源于最佳拟合。(个人理解——拟合:调整分类边界)
用一条直线对一些数据点进行拟合(该线称为最佳拟合直线),而拟合过程称作回归。
2:Logistic分类的思想:根据现有数据对分类边界线建立回归公式。
由于Logistic属于分类算法,所以我们需要一个能够进行分类的函数,在单位阶跃函数和sigmoid函数之间,我们选择了更容易理解和处理的sigmoid函数。当sigmoid值>0.5为一类,<0.5为一类,0.5为分界线。
(在z=0时sigmoid=0.5,当横坐标够大时,在0处该函数很像阶跃函数)
所以现在思路明确了,只需要用测试数据求得sigmoid的值就能根据该值对测试数据分类了。如何求得sigmoid函数的值呢?(求唯一变量z)下面我在纸上写了逻辑回归求sigmoid最主要的思路(包括用梯度算法求最佳回归系数)
如果对Logistic一无所知,那么上面的步骤你可能只是有了大概的思路,而不知道如何去求对应值。下面具体对上面的步骤解释一下:
- 观察sigmoid公式得知只需求Z这个未知值
- Z=W0X0+W1X1+W2X2+...+WnXn (Wi:最佳回归系数 Xi:特征值) :表示将两个对应数值相乘后求和得到Z
- 每个特征对应的最佳回归系数我们使用梯度上升的优化方法来求(梯度即函数单调方向,也是函数变化最快的方向)
- 梯度只是方向,而α是步长,两者相乘来求移动距离
- 最后使用梯度算法迭代公式,来不断更新W
- 其中的梯度也就是损失函数在此处的导数,下面代码可能更好解释理解
(梯度具体说明文章:https://blog.csdn.net/qq_36523839/article/details/83153859)
代码实现和优化:
(说明:我用pandas数据结构也能实现,代码量也少,但相对与numpy理解复杂一点,所以这里还是用作者的numpy来实现)
样本示例(三列,前两列为点的x、y值,最后一列表示分类):
代码一(最原始没有优化过的梯度更新算法):
import numpy as np
import matplotlib.pyplot as plt
# 读取文本数据,返回数据集和目标值
def loadDataSet():
dataMat = []
labelMat = []
fr = open('testSet.txt')
for line in fr.readlines():
lineArr = line.strip().split()
# 该数据集中,添加了一列并初始化为1,便于后续的计算,但是在其他数据集中,一般没有必要添加
dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])
labelMat.append(int(lineArr[2]))
return dataMat,labelMat
# 运算的核心函数Sigmoid
def sigmoid(inX):
return 1.0 / (1+np.exp(-inX))
# 核心函数alpha*gradient更新回归系数
def gradAscent(dataMatIn, classLabels):
# 为了便于计算,mat将两个列表数据转换为numpy的数组
dataMatrix = np.mat(dataMatIn)
labelMat = np.mat(classLabels).transpose() # transpose矩阵转置
m,n = np.shape(dataMatrix)
alpha = 0.001 # 设置步长
maxCycles = 500 # 设置循环次数
weights = np.ones((n,1)) # 初始化每个特征的回归系数为1
for k in range(maxCycles):
# 得到每一行的sigmoid值 (两个相乘得到z值)
h = sigmoid(dataMatrix*weights) # 矩阵相乘 sigmoid(sum(每个特征*每个系数)) 行*列,由于是矩阵相乘,所以在相乘时便求和了
# 用一直更新的回归系数求预测值,然后与真实值之间求误差,误差越来越小,则系数更新变化也越来越小,最后趋于稳定
error = (labelMat - h) # 每行分类与对应sigmoid相减 误差值越来越小了
# 数据集转置*误差 每条样本*该样本的误差值
weights = weights + alpha * dataMatrix.transpose() * error
return weights
dataArr,labelMat = loadDataSet()
weights = gradAscent(dataArr,labelMat)
print(weights)
重点说明更新原理:梯度上升法本质是通过梯度联系自变量与因变量,使得自变量的变化能定向的改变因变量。这里我们讨论error的大小与回归系数weights的关系,自变量事实上是weights,因变量是error,我们将自变量weights结合error的大小进行变化(error的作用与梯度类似)可以使得因变量error绝对值最小化;例如,若error>0,则weights变大,由于error与weights成反比,所以error变小,即趋近于0,绝对值变小;若error<0,则weights变小,导致下次迭代后error变大,即error绝对值变小。∴该迭代结果是使error的绝对值最小化 。
weights=weights+alpha∗error∗dataMat
添加一个可视化的函数:
def plotBestFit(weights):
dataMat,labelMat=loadDataSet()
dataArr = np.array(dataMat)
n = np.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 = np.arange(-3.0, 3.0, 0.1)
y = (-weights[0]-weights[1]*x)/weights[2] #由sigmoid=0这个分界点 使用唯一非常量Z 来构造函数
ax.plot(x, y)
plt.xlabel('X1')
plt.ylabel('X2')
plt.show()
plotBestFit(weights.getA())
该函数为画最佳拟合直线的函数,调用该函数只需传入我们获得了的回归系数weights,函数内读取文本,然后画出所有点
重点说一下这条线如何画:当sigmoid=0.5,此时为0、1类的分界点,而等于0.5表示Z=0,所以根据Z值的公式Z=sum(各回归系数*各特征) ——>0=W0X0+W1X1+W2X2 由于本数据中X0=1 所以推导出y = (-weights[0]-weights[1]*x)/weights[2]
改进:梯度上升算法在每次更新回归系数时都需要遍历整个数据集,该方法在处理100个左右的数据集时尚可,但是如果由数十亿样本和大量特征,那么该方法的计算复杂度就太高了。一种改进方法是一次仅用一个样本点来更新回归系数,该方法成为随机梯度上升算法。
代码二(用stoGradAscent0函数替代gradAscent函数):
def stocGradAscent0(dataMatrix, classLabels):
dataMatrix = np.array(dataMatrix) #将列表转换格式
m,n = np.shape(dataMatrix)
alpha = 0.01
weights = np.ones(n)
for i in range(m):
h = sigmoid(sum(dataMatrix[i]*weights))
error = classLabels[i] - h
weights = weights + alpha * error * dataMatrix[i]
return weights
随机梯度上升算法与梯度上升算法很相似,区别在于:第一,后者的h和error都是向量,前者全是数值;第二,前者没有矩阵转换的过程,所有变量的数据类型都是Numpy数组。所以可以根据随机梯度上升的运算来理解系数的更新过程。
看看绘图的效果(plotBestFit(weights)):
效果一般,较差,但是前者是在整个数据集上迭代了500次才有了好的结果,所以这里我们再优化算法,使它在整个数据集上迭代200次。
代码三(用stoGradAscent0函数替代gradAscent函数):
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
dataMatrix = np.array(dataMatrix) #将列表转换格式
m,n = np.shape(dataMatrix)
weights = np.ones(n)
for j in range(numIter):
dataIndex = list(range(m))
for i in range(m):
alpha = 4/(1.0+j+i)+0.0001 #变化的alpha
randIndex = int(np.random.uniform(0,len(dataIndex))) #随机样本下标
h = sigmoid(sum(dataMatrix[randIndex]*weights))
error = classLabels[randIndex] - h
weights = weights + alpha * error * dataMatrix[randIndex]
del(dataIndex[randIndex]) #从整个数据中去除已经使用过的样本
return weights
第二次优化有两点:第一,alpha每次迭代都会调整,这会缓解数据波动或高频波动,此外alpha可能减小,但有一个常数项,所以该值会一直存在影响力;第二,通过随机取样来更新回归系数,减少周期性的波动。
我绘制了优化过后的拟合直线 以及 三个特征的变化(趋于稳定的波动):
1
2
根据第二次的优化,随机梯度上升算法已经能以更快的时间更好的划分数据了。
真实样例数据:
以上为Logistic的原理说明和样例代码的实现,下面使用优化过的随机梯度上升算法处理一个具体的数据集——《病马的死亡率》
数据样例:(100个样本,每个样本20个特征,缺失值以0填充这样不影响回归系数)
import numpy as np
# 读取文本数据,返回数据集和目标值
def loadDataSet():
dataMat = []
labelMat = []
fr = open('testSet.txt')
for line in fr.readlines():
lineArr = line.strip().split()
# 该数据集中,添加了一列并初始化为1,便于后续的计算,但是在其他数据集中,一般没有必要添加
dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])
labelMat.append(int(lineArr[2]))
return dataMat,labelMat
# 运算的核心函数Sigmoid
def sigmoid(inX):
return 1.0 / (1+np.exp(-inX))
# 核心函数alpha*gradient更新回归系数
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
m,n = np.shape(dataMatrix)
weights = np.ones(n)
for j in range(numIter):
dataIndex = list(range(m))
for i in range(m):
alpha = 4/(1.0+j+i)+0.0001
randIndex = int(np.random.uniform(0,len(dataIndex)))
h = sigmoid(sum(dataMatrix[randIndex]*weights))
error = classLabels[randIndex] - h
weights = weights + alpha * error * dataMatrix[randIndex]
del(dataIndex[randIndex])
return weights
# 根据sigmoid来判断分类
def classifyVector(inX, weights):
prob = sigmoid(sum(inX*weights))
if prob > 0.5: return 1.0
else: return 0.0
# 封装函数 获取文本并调用函数
def colicTest():
frTrain = open('horseColicTraining.txt')
frTest = open('horseColicTest.txt')
trainingSet = []
trainingLabels = []
for line in frTrain.readlines():
currLine = line.strip().split('\t')
lineArr =[]
for i in range(21):
lineArr.append(float(currLine[i]))
trainingSet.append(lineArr)
trainingLabels.append(float(currLine[21]))
trainWeights = stocGradAscent1(np.array(trainingSet), trainingLabels, 1000) # 获得训练集的回归系数,测试集使用
errorCount = 0
numTestVec = 0.0
for line in frTest.readlines():
numTestVec += 1.0
currLine = line.strip().split('\t')
lineArr =[]
for i in range(21):
lineArr.append(float(currLine[i]))
if int(classifyVector(np.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.0
for k in range(numTests):
errorSum += colicTest()
print("after %d iterations the average error rate is: %f" % (numTests, errorSum/float(numTests)))
multiTest()
由于有30%的数据缺失,所以错误率稍高一点。
sklearn实现:
由于篇幅太多,这里仅仅使用sklearn简单实现,如果需要深入了解逻辑回归sklearn的实现方式可以参考官方文档或查询相对专业一点的文献。
import pandas as pd
from sklearn.linear_model import LogisticRegression
train = pd.read_table('horseColicTraining.txt',header=None) # 读取文件 DataFrame
test = pd.read_table('horseColicTest.txt',header=None)
train_label = train[21] # 两个数据集作同样的操作,获得数据集和标签
train.drop([21],axis=1,inplace=True)
test_label = test[21]
test.drop([21],axis=1,inplace=True)
lr_model = LogisticRegression() #调用模型,但是并未经过任何调参操作,使用默认值
lr_model.fit(train,train_label) #训练模型
print(lr_model.score(test,test_label)) #获取测试集的评分
未经调参,直接使用sklearn的逻辑回归,准确率还可以,所以在正确使用sklearn并调节好参数时,sklearn是一个很强大的工具。
借用作者的总结:
参考资料:https://blog.csdn.net/lvsolo/article/details/50966012 回归系数更新的分析来自于此
参考书籍:《机器学习实战》