概述
在这一部分内容中,我们将首次接触最优化算法。其实我们在日常生活中也会经常发现最优化问题,比如如何在最短时间内从A点到达B点?如何投入最少的时间获得最大的效益等等。可见,求解最优化问题将会是我们遇见的很多问题都变得简洁,下面我们将介绍几个最优化算法,并利用他们训练拿出一个非线性函数用于分类。
假如二维平面内有一些数据点,现在我们需要用一条直线来对这些点进行拟合,这种不断的求取最佳拟合曲线的过程就是一个最优化回归过程。Logistic回归进行分类的主要思想就是:根据现有的数据对分类的边界线建立回归公式来进行分类,既要找到最佳拟合参数集,训练这个分类器就是不断地寻找最佳拟合参数,使用最优化算法。
Logistic回归的一般过程:
- 收集数据:采用任意方法收集数据。
- 准备数据:因为需要进行的是距离计算,因此要求数据类型为数值型。结构化数据最佳。
- 分析数据:采用任意方法对数据进行分析。
- 训练算法:大部分时间将用于训练,训练的目的是为了找到最佳的分类回归系数。
- 测试算法:一旦训练步骤完成,分类将会很快
- 使用算法:首先:我们需要输入一些数据,并将其转换成对应的结构化数值,接着基于训练好的回归系数就可以对这些数值进行简单的回归计算,判定他们属于那个类别,在这之后,我们就可以在输出的类别上做一些其他分析工作。
我们用到的最优化算法包括基本梯度上升和改进的随机梯度上升算法,这些最优化算法将用于分类器的训练。最后的实例将会预测一匹病马是否能被治愈。
基于logistic回归和Sigmoid函数的分类
我们理想中的回归函数应该是能接受所有输入,然后输出对应的类别,假如只有两个类别,那么就输出0和1.Sigmoid函数就有着类似的性质,它的计算公式如下所示:
下图给出了该函数的曲线图。当x为0时,函数值为0.5,随着x的值逐渐增大,对应的函数值无限接近于1,同样被,当x逐渐减小,函数值将逼近于0。我们可以想象一下,当横坐标尺度足够大,那么Sigmoid函数看起来很像是一个阶跃函数。
当我们在实现Logistic回归分类器时,我们可以在每个特征上都乘以一个回归系数,然后将所有的结果值相加,最后将将这个总和代入Sigmoid函数,进而得到一个在0-1范围内的结果,当结果大于0.5,在数据被分为1类,小于0.5被分为0类。所以,Logistic回归其实也是一种概率估计。
现在的问题就是如何确定最佳回归系数的大小。
基于最优化方法的最佳回归系数确定
Sigmoid函数的输入记为Z,由下面的公式得出:
采用向量的写法形式,上述公式可以写成,向量x是分类器的输入数据,向量w就是最佳回归系数,为了寻找这个最佳参数,需要用到最优化理论的一些知识。
梯度上升法
梯度上升法的思想就是,如果想要寻找某个函数的最大值,最好的方法就是沿着该函数的梯度方向探寻。如果梯度记为∇,那么而函数f(x,y)的梯度由下面这个式子表示:
这个梯度意味着要沿着x的方向移动,沿着y的方向移动。其中,函数必须要在待计算的点上有定义并且可微分。
从途中我们可以看到,梯度算子总是指向函数增长最快的方向。这里所说的是移动方向,并未提到移动量的大小。这个量值称为步长,记作α,梯度上升算法的迭代公式如下:
这个公式将被一直迭代执行,直到达到某个停止条件为止,比如迭代次数达到某个指定值或者算法达到某个可以允许的误差范围。
梯度上升算法是用来求取函数的最大值,而梯度下降算法则是用来求取函数的最小值。
接下来我们来看一个Logistic分类器的应用例子,首先进行数据集的可视化,添加如下代码:
import matplotlib.pyplot as plt
import numpy as np
def showDataSet():
fr = open("testSet.txt")
filelines = fr.readlines()
dataMat = np.zeros((len(filelines),2))
ClassLabelColors = []
index = 0
for each in filelines:
line = each.strip()
ListFromline = line.split('\t')
dataMat[index,:] = ListFromline[0:2]
if ListFromline[-1] == '1':
ClassLabelColors.append('red')
else:
ClassLabelColors.append('black')
index += 1
plt.scatter(x=dataMat[:,0],y=dataMat[:,1],color=ClassLabelColors,s=15,alpha=0.5)
if __name__ == '__main__':
showDataSet()
执行程序可以看到数据集的可视化结果:
训练集中包含了100个样本点,每个点的横纵坐标就是他的特征,在此数据集上,我们将使用梯度上升法找到最佳回归参数,也就是拟合出Logistic回归模型的最佳参数。
梯度上升法的伪代码如下:
- 每个回归系数初始化为1
- 重复R次
- 计算整个数据集的梯度
- 使用alpha * gradient更新回归系数的向量
- 返回回归系数
下面的代码是梯度上升算法的具体实现,创建一个logRegres.py的文件,输入以下代码:
def loadDataSet():
'''
函数说明:读取数据集
Returns: dataMat - 包含了数据特征值的矩阵
dataLabel - 数据集的标签
'''
dataMat = [];dataLabel = []
fr = open("testSet.txt")
fileLines = fr.readlines()
for line in fileLines:
lineArr = line.strip().split('\t')
dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])
dataLabel.append(int(lineArr[2]))
return dataMat,dataLabel
def sigmoid(inx):
'''
函数说明:sigmoid函数
Parameters: inx - 经历一次前向传播还没有激活的值
Returns: 激活后的值
'''
return 1.0/(1+exp(-inx))
def gradAscent(dataMatIn,classLabels):
#将输入数据集转换为矩阵
dataMatrix = mat(dataMatIn)
LabelsMat = mat(classLabels).transpose()
m,n = shape(dataMatrix)
weights = ones((n,1)) #初始化系数矩阵,因为数据集为m行n列,所以权重矩阵为n行
maxIter = 500 #迭代次数
alpha = 0.001 #步长,也可以称作学习率
for i in range(maxIter):
z = sigmoid(dataMatrix * weights)
error = LabelsMat - z #计算误差
weights = weights + alpha * dataMatrix.transpose() * error #更新权重矩阵
return weights
if __name__=='__main__':
dataArr,labelMat = loadDataSet()
print(gradAscent(dataArr,labelMat))
执行上述程序,看看实际效果:
分析数据:画出决策边界
现在我们已经得到了回归系数,他确定了不同数据之间的分隔线,接下来的问题就是画出该分割线,添加以下代码:
# -*- coding: utf-8 -*-
"""
Created on Mon Nov 19 21:43:41 2018
@author: jacksong1996
"""
from numpy import *
import matplotlib.pyplot as plt
def loadDataSet():
'''
函数说明:读取数据集
Returns: dataMat - 包含了数据特征值的矩阵
dataLabel - 数据集的标签
'''
dataMat = [];dataLabel = []
fr = open("testSet.txt")
fileLines = fr.readlines()
for line in fileLines:
lineArr = line.strip().split('\t')
dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])
dataLabel.append(int(lineArr[2]))
return dataMat,dataLabel
def sigmoid(inx):
'''
函数说明:sigmoid函数
Parameters: inx - 经历一次前向传播还没有激活的值
Returns: 激活后的值
'''
return 1.0/(1+exp(-inx))
def gradAscent(dataMatIn,classLabels):
#将输入数据集转换为矩阵
dataMatrix = mat(dataMatIn)
LabelsMat = mat(classLabels).transpose()
m,n = shape(dataMatrix)
weights = ones((n,1)) #初始化系数矩阵,因为数据集为m行n列,所以权重矩阵为n行
maxIter = 500 #迭代次数
alpha = 0.001 #步长,也可以称作学习率
for i in range(maxIter):
z = sigmoid(dataMatrix * weights)
error = LabelsMat - z #计算误差
weights = weights + alpha * dataMatrix.transpose() * error #更新权重矩阵
return weights
def plotBestFit(weights):
dataMat,dataLabel = loadDataSet()
dataArray = array(dataMat)
xcord1 = [];ycord1 = []
xcord2 = [];ycord2 = []
n = shape(dataArray)[0]
for i in range(n):
if int(dataLabel[i]) == 0:
xcord1.append(dataArray[i,1])
ycord1.append(dataArray[i,2])
else:
xcord2.append(dataArray[i,1])
ycord2.append(dataArray[i,2])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xcord1,ycord1,color = 'red',marker = 's',s = 30)
ax.scatter(xcord2,ycord2,color = 'blue',s = 30)
x = arange(-3.0,3.0,0.1)
y = (-weights[0,0] - weights[1,0]*x)/weights[2,0] #0=w0X0+w1X1+w2X2
ax.plot(x,y)
plt.xlabel('X1');plt.ylabel('X2')
plt.show()
if __name__=='__main__':
dataArr,labelMat = loadDataSet()
weights = gradAscent(dataArr,labelMat)
plotBestFit(weights)
执行程序可以看到输出结果如下图所示,可以看到这个分类效果非常不错,仅从图上看只错分了几个点而已,尽管这个例子看起来很简单,但是却经过了大量的计算,下面我们将会对这个算法进行改进。
训练算法:随机梯度上升
因为梯度上升算法在每次更新回归系数时都需要遍历整个数据集,因此当数据集样本很大,特征成千上万时,那么该方法的计算成本就会大大加大。改进的方法时一次仅用一个样本点来更新回归系数,该方法称为随机梯度上升算法,由于可以在新样本到来时对分类器进行增量式更新,因而随机梯度上升算法是一个在线学习算法,而一次处理所有数据则被称为“批处理”。
随机梯度上升算法的伪代码如下所示:
- 所有回归系数初始化为1
- 对数据集中的每个样本
- 计算该样本的梯度
- 使用alpha x gradient更新回归系数值
- 返回回归系数值
具体代码如下所示:
from numpy import *
import matplotlib.pyplot as plt
import random
def loadDataSet():
'''
函数说明:读取数据集
Returns: dataMat - 包含了数据特征值的矩阵
dataLabel - 数据集的标签
'''
dataMat = [];dataLabel = []
fr = open("testSet.txt")
fileLines = fr.readlines()
for line in fileLines:
lineArr = line.strip().split('\t')
dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])
dataLabel.append(int(lineArr[2]))
return dataMat,dataLabel
def sigmoid(inx):
'''
函数说明:sigmoid函数
Parameters: inx - 经历一次前向传播还没有激活的值
Returns: 激活后的值
'''
return 1.0/(1+exp(-inx))
def gradAscent(dataMatIn,classLabels):
#将输入数据集转换为矩阵
dataMatrix = mat(dataMatIn)
LabelsMat = mat(classLabels).transpose()
m,n = shape(dataMatrix)
weights = ones((n,1)) #初始化系数矩阵,因为数据集为m行n列,所以权重矩阵为n行
maxIter = 500 #迭代次数
alpha = 0.001 #步长,也可以称作学习率
for i in range(maxIter):
z = sigmoid(dataMatrix * weights)
error = LabelsMat - z #计算误差
weights = weights + alpha * dataMatrix.transpose() * error #更新权重矩阵
return weights
def stocGradAscent(dataMatrix,classLabels,numIter = 150):
m,n = shape(dataMatrix)
weights = ones(n)
for i in range(numIter):
dataIndex = list(range(m))
for j in range(m):
alpha = 4/(1.0+i+j)+0.01
randindex = int(random.uniform(0,len(dataIndex)))
h = sigmoid(sum(dataMatrix[randindex]*weights))
error = classLabels[randindex] - h
weights = weights + [alpha*error*k for k in dataMatrix[randindex]]
del(dataIndex[randindex])
return weights
def plotBestFit(weights):
dataMat,dataLabel = loadDataSet()
dataArray = array(dataMat)
xcord1 = [];ycord1 = []
xcord2 = [];ycord2 = []
n = shape(dataArray)[0]
for i in range(n):
if int(dataLabel[i]) == 0:
xcord1.append(dataArray[i,1])
ycord1.append(dataArray[i,2])
else:
xcord2.append(dataArray[i,1])
ycord2.append(dataArray[i,2])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xcord1,ycord1,color = 'red',marker = 's',s = 30)
ax.scatter(xcord2,ycord2,color = 'blue',s = 30)
x = arange(-3.0,3.0,0.1)
#y = (-weights[0,0] - weights[1,0]*x)/weights[2,0] #0=w0X0+w1X1+w2X2
y = (-weights[0] - weights[1]*x)/weights[2]
ax.plot(x,y)
plt.xlabel('X1');plt.ylabel('X2')
plt.show()
if __name__=='__main__':
dataArr,labelMat = loadDataSet()
weights = stocGradAscent(dataArr,labelMat)
plotBestFit(weights)
执行程序,打印输出带分割线的图片:
可以看到这张图类似于上一张图。但是他的计算成本大大降低了,该算法改进的地方在于:1 学习率衰减,alpha会随着迭代次数不断减小,但不会减小到0,另外alpha也不是严格下降的,这种避免参数的严格下降也经常见于模拟退火算法等其他优化算法中。2 随机选取样本来更新回归系数,这种做法将减少周期性波动.下图显示了每次迭代时各个回归系数的变坏情况。
实例:从疝气病症预测病马的死亡率
疝气病时描述马胃肠痛的术语,本例的数据集中包含了医院检测马疝气病的一些指标,有些指标比较主观,有的指标难以测量,比如马的疼痛级别。除了部分指标较为主观和难以测量,还需要处理数据集中的缺失值情况,再利用Logistic回归和随机梯度上升算法预测病马的生死。
在处理数据的缺失值时,我们有以下可选的做法:
- 使用可用特征的均值来填补缺失值
- 使用特殊值来填补缺失值,如-1
- 忽略有缺失值的样本
- 使用相似样本的均值填补缺失值
- 使用另外的机器学习算法预测缺失值
在这里,我们选择用0来代替缺失值,这样做的直觉是因为0在回归系数更新时不会影响到系数的值。回归系数的更新公式如下:
又因为sigmoid(0) = 0.5,即他对结果的预测也不具有任何倾向性;此外,该数据集中的特征取值一般不为0,因此在某种意义上说他也满足“特殊”这个要求。
如果类别标签已经丢失,那么我们的简单做法是将该条数据丢弃,这是因为类别标签与特征不同,很难确定采用某个合适的值来替换
测试算法:用logistic回归进行分类:
添加以下代码:
def classifyVector(inx,weights):
'''
函数说明:分类函数
Parameter:inx - 输入的特征向量
weights - 回归系数
return: 类别标签
'''
prob = sigmoid(sum(inx*weights))
if prob > 0.5:
return 1.0
else:
return 0.0
def colicTest():
'''
函数说明:打开测试集和训练集,训练数据预测错误率
Parameter:无
Retrurn:errorRate - 错误概率
'''
frTrain = open("horseColicTraining.txt")
frTest = open("horseColicTest.txt")
trainList = [];trainLabels = []
for line in frTrain.readlines():
lineArr = []
currLine = line.strip().split('\t')
for i in range(21):
lineArr.append(float(currLine[i]))
trainList.append(lineArr)
trainLabels.append(float(currLine[21])) #处理训练集数据
weight = stocGradAscent(array(trainList),trainLabels,500) #随机梯度上升得到回归系数
errorCount = 0;numTest = 0.0
for line in frTest.readlines():
lineArr = []
numTest += 1.0
currLine = line.strip().split('\t')
for i in range(21):
lineArr.append(float(currLine[i])) #处理测试集数据
if int(classifyVector ( array(lineArr), weight)) != int( currLine[21] ):
errorCount += 1
errorRate = float(errorCount) / numTest #计算错误率
print("the test errorrate is %f" % errorRate)
return errorRate
def multiTest():
'''
函数说明:计算迭代numtests次后的错误率
'''
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)))
执行程序,打印预测错误率如下所示:
因为有百分之30的数据缺失,所以得到的平均错误率为36%左右还是可以接受的,如果调整迭代次数以及学习率,平均错误率还会下降,但是可能会造成过拟合。
小结
逻辑回归是非常常见的模型,训练速度很快,虽然使用起来没有支持向量机(SVM)那么占主流,但是解决普通的分类问题是足够了,训练速度也比起SVM要快不少。如果你要理解机器学习分类算法,那么第一个应该学习的分类算法个人觉得应该是逻辑回归。理解了逻辑回归,其他的分类算法再学习起来应该没有那么难了。
PS:
本篇文章所涉及到的数据集以及完整代码都存放在github中
欢迎各位大佬fock or star
https://github.com/AdventureSJ/ML-Notes