机器学习之回归算法—Logistic回归

Logistic回归是机器学习中的一种分类模型,由于其算法的简单和高效,在实际应用中非常广泛。例如,在实际生活中,当判断类似是否为垃圾邮件、病人是否得了癌症等以二值型(即1或0)输出的分类时,我们可以使用基于对数几率回归的Logistic回归。

注意:logistic又翻译为”逻辑回归”,不过其本意应为“对数几率回归”,与”逻辑”一词并无关系

Logistic回归

什么是回归:假设现在有一些数据点,我们用一条直线对其进行拟合(该线称为最佳拟合直线),这个拟合过程就称为回归

Sigmoid函数

由于通常的线性回归会拟合出一条直线,如果我们想要做的是分类任务该怎么办?只需找一个单调可微函数将该分类任务的真实标记y与线性回归模型的预测值连接起来。例如,在两个类的情况下,函数输出0或1。这里,我们要将线性回归产生的预测实值转化为0/1值,通常我们会想到”单位阶跃函数”。

单位阶跃函数的原理是,若预测值z大于零就判为正例,小于零就判为反例。但该函数不连续,而且判断方式不具有普遍性,于是我们希望能找到在一定程度上近似单位阶跃函数的”替代函数”,Sigmoid函数(也称Logistic函数)正是这样的一种替代函数。

Sigmoid函数的计算公式如下:

该函数的变化曲线图如下:

可以看出,该曲线将z值转化为一个接近0或1的y值,并且在输出z=0时变化很陡,保证了阶跃函数的特征同时又具有连续性。

Sigmoid函数的输入记为z,一般是线性回归函数,即z=w0x0+w1x1+w2x2+…….+wnxn,在这里,我们采用向量的写法:

整合到Sigmoid函数中:

可以看出,为了实现Logistic分类器,我们可以在每个特征上乘以一个回归系数,然后把所有的结果值相加,将这个总和代入Sigmoid函数中。那么问题来了,如何寻找最佳回归系数θ呢?

最佳回归系数

我们将Sigmoid函数中的y视为类后验概率估计p(y=1|x),进行如下假设:

然后将它们合二为一:

合并出来的Cost,我们称之为代价函数(Cost Function)。当y等于1时,(1-y)对应的项为0,;当y等于0时,y对应的那项为0。为了简化问题,我们对整个表达式求对数(将指数函数对数化是处理数学问题常见的方法)

这个代价函数,是针对一个样本而言的,给定一个样本,我们就可以通过这个代价函数求出样本所属类别的概率。这时,我们可以通过最大似然法来估计,即令每个样本所属真实类别越大越好。所以,得到如下公式:

综上所述,满足J(θ)的最大的θ值即是我们需要求解的模型。

梯度上升算法

由于需要求出最大值,所以我们要用到最优化理论的一些知识。下面首先介绍梯度上升这一优化算法,梯度上升算法基于的思想是:要找到某函数的最大值,最好的方法打是沿着该函数的梯度方向探寻。如果梯度记为▽,则函数f(x,y)的梯度有下式表示:

这个梯度表示要沿x方向移动∂f(x,y)/∂x,沿y方向移动∂f(x,y)/∂y。梯度算子总是指向函数值增长最快的方向,但在真实环境下并不容易找出最快的方向,此时我们可以用迭代的方法来做,就像爬坡一样,一点一点逼近极值。梯度上升迭代算法如下:

我们现在要求J()的极值,所以建立如下迭代算法:

那么,现在我只要求出J(θ)的偏导,就可以利用梯度上升算法,求解J(θ)的极大值了。

所以,梯度上升迭代公式为:


基于《机器学习实战》这本书,这里打算将实现Logistic回归分为以下几个步骤,分别为:收集并分析数据、梯度上升算法、绘制数据集及决策边界、分析迭代算法以及改进梯度上升算法几个步骤完成。

准备数据

实验数据有两个特征,可以在二维坐标用x,y表示,数据集一共三列,前两列表示数据特征,第三列表示标签,下面用代码将其提取。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
# 加载数据
def loadDataSet():
dataMat =[] # 存放数据
labelMat = [] # 存放标签
fr = open('testSet.txt','r')
for line in fr.readlines():
lineArr = line.strip().split()
dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])
labelMat.append(int(lineArr[2]))
return dataMat,labelMat
#绘制数据集
def plotDataSet():
dataMat,labelMat = loadDataSet()
dataArr = np.array(dataMat) # 构建成numpy的array数组
n = dataArr.shape[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=20, c='red', marker='s',alpha=.5)
ax.scatter(xcord2, ycord2, s=20, c= 'green', alpha=.5)
plt.title('DataSet')
plt.xlabel('x'); plt.ylabel('y')
plt.show()

运行结果如下显示:

由于该数据集有两个特征,那么Sigmoid函数的输入z=w0x0+w1x1+w2x2,其中,x1为第一列数据,x2为第二列数据,为了方便计算,x0为全是1的列向量(这里可以理解为方便向量计算的一种方式),而w0、w1、w2就是我们要求的回归系数(最优参数)。

训练算法

首先,我们将梯度上升迭代公式矢量化:

根据矢量化的公式,编写如下代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
def sigmoid(inX):
return 1.0/(1+np.exp(-inX))
# 梯度上升算法
def gradAscent(dataMatIn, classLabels):
dataMatrix = np.mat(dataMatIn)
labelMat = np.mat(classLabels).transpose() # 转换成numpy的mat并进行转置
m, n = np.shape(dataMatrix) # 返回dataMatrix矩阵的大小,m是行数,n是列数
alpha = 0.01 # 移动步长,也就是学习速率
maxCycles = 500 # 最大迭代次数
weights = np.ones((n,1))
weights_array = np.array([])
for k in range(maxCycles): # 梯度上升矢量化公式
h = sigmoid(dataMatrix*weights)
error = (labelMat - h)
weights = weights + alpha * dataMatrix.transpose() * error
weights_array = np.append(weights_array,weights)
weights_array = weights_array.reshape(maxCycles, n)
return weights.getA(),weights_array # 将矩阵转化为数组,并返回权重数组
# 绘制基于回归系数的数据集及决策边界
def plotBestFit(weights):
dataMat, labelMat = loadDataSet()
dataArr = np.array(dataMat)
n = dataArr.shape[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=20, c='red', marker='s', alpha=.5) # 绘制正样本
ax.scatter(xcord2, ycord2, s=20, c='green',alpha=.5) # 绘制负样本
x = np.arange(-3.0, 3.0, 0.1)
y = (-weights[0] - weights[1] * x) / weights[2]
ax.plot(x, y)
plt.title('BestFit')
plt.xlabel('X1')
plt.xlabel('y2')
plt.show()

运行结果如下:

[[ 12.78439352]
[ 1.11820764]
[ -1.74082686]]

可以看出,梯度上升迭代算法已经计算出最佳回归系数[w0,w1,w2],然后将拟合直线在图中显示出来,分类效果还算不错。

改进的随机梯度上升算法

从刚才的算法中我们可以看出,每次计算最佳回归系数,都要遍历整个数据集,迭代500次就要遍历500次,当处理这种小型的数据集还可以,但是如果有成千上万的特征和样本,那么该方法的计算复杂度就太高了。一种改进方法是使用改进的随机梯度上升算法。由于可以在新样本到来时对分类器进行增量更新,因而随机梯度上升算法是一个在线学习算法,下面给出代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# 改进的随机梯度上升算法
def stocGradAscentl(dataMatrix, classLabels, numIter=150):
m, n = np.shape(dataMatrix)
weights = np.ones(n)
weights_array = np.array([])
for j in range(numIter):
dataIndex = list(range(m))
for i in range(m):
alpha = 4/(1.0 + j + i) + 0.01 # 降低alpha的大小,每次减小1/(j+1)
randIndex = int(random.uniform(0,len(dataIndex))) # 随机选取样本
h = sigmoid(sum(dataMatrix[randIndex] * weights)) # 选择随机选取的一个样本,记作h
error = classLabels[randIndex] - h # 计算误差
weights = weights + alpha * error * dataMatrix[randIndex] # 更新回归系数
weights_array = np.append(weights_array, weights, axis=0) # 添加回归系数到数组中
del(dataIndex[randIndex]) # 删除已经使用的样本
weights_array = weights_array.reshape(numIter*m, n) # 改变维度
return weights, weights_array

从该程序可以看出,改进的随机梯度上升算法相比梯度上升算法有一下几个特点:

  • 变量h和误差都是向量,前者都是数值
  • 没有矩阵的转化过程,所有变量的数据类型都是numpy数组
  • alpha在每次迭代时都会调整,但永远不会减小到零,会缓解数据波动
  • 随机选取样本来更新回归系数,每次随机从列表中选取一个值,然后从列表中删除这个值(再进行下一次迭代),这种方法可以减少周期性波动。

综上所述,使用改进的梯度上升算法可以减少数据集遍历的次数,使回归系数收敛更快,降低算法复杂度。

从疝气病症状预测病马的死亡率

处理数据的缺失值

在实际情况中,数据常常并不是完整的。但数据相当昂贵,扔掉或重新获取都是不可取的。所以必须彩英一些方法来解决问题:

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

对于本算法,numpy数组不允许数据包含缺失值。这里选择实数0来替换所有缺失值,恰好能适用与Logistic回归,由于Sigmoid(0)=0.5,即它对任何结果的预测不具有任何倾向性,因此不会对误差向造成任何影响。

用Logistic进行分类

使用Logistic回归方法进行分类并不需要做太多工作,所需做的只是吧数据集上的每个特征向量乘以最优化方法得来的回归系数,然后再将该成绩结果求和,最后输入到Sigmoid函数中即可。如果对应的Sigmoid值大于0.5就预测类别标签为1,否则为0。下面看代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
"""
说明:分类函数
基本思想:把数据集的每个特征向量诚意最优化方法得来的最优回归系数,
然后在将该乘积结果求和,最后输入到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(len(currLine)-1):
lineArr.append(float(currLine[i]))
trainingSet.append(lineArr) # 相当于将数据的内容按每行一个list(lineArr)的方式存放到一个总的list(trainSet)中
trainingLabels.append(float(currLine[-1]))
trainWeights = stocGradAscentl(np.array(trainingSet), trainingLabels, 500)
errorCount = 0; numTestVec = 0.0
for line in frTest.readlines(): # 测试数据
numTestVec += 1.0
currLine = line.strip().split('\t')
lineArr = []
for i in range(len(currLine) - 1):
lineArr.append(float(currLine[i]))
if int(classifyVector(np.array(lineArr),trainWeights)) != int(currLine[-1]): # 见classifyVector函数说明
errorCount += 1
errorRate = (float(errorCount)/numTestVec) * 100
print("测试集错误率为: %.2f%%" % errorRate)
if __name__ == '__main__':
colicTest()

运行结果如下:

可以看出,这种随机梯度算法错误率还是蛮高的。事实上,这个结果并不算太差,因为还有30%的数据缺失。如果调整最优化算法的迭代次数,错误率还可以下降到20%左右。在Sklearn中,可以通过调整最优化算法来降低错误率,以后打算用一篇文章单独讲一下用Sklearn实现各种机器学习算法。

Powered by Hexo and Hexo-theme-hiker

Copyright © 2016 - 2019 Dreaouth All Rights Reserved.

访客数 : | 访问量 :