统计学习方法|逻辑斯蒂原理剖析及实现
阅读数：12,208前言
《统计学习方法》一书在前几天正式看完，由于这本书在一定程度上对于初学者是有一些难度的，趁着热乎劲把自己走过的弯路都写出来，后人也能走得更顺畅一点。
以下是我的github地址，其中有《统计学习方法》书中所有涉及到的算法的实现，也是在写博客的同时编写的。在编写宗旨上是希望所有人看完书以后，按照程序的思路再配合书中的公式，能知道怎么讲所学的都应用上。（有点自恋地讲）程序中的注释可能比大部分的博客内容都多。希望大家能够多多捧场，如果可以的话，为我的github打颗星，也能让更多的人看到。
github：
相关博文：
正文
逻辑斯蒂回归的直观理解
我们还记不记得感知机一章中给出的这张图？还记得感知机是怎么工作的吗？
感知机生成一个超平面（在二维中表现为一条直线），这个超平面可以将图中的两类点区分开。像下面这张图，图中的所有直线都可以作为感知机的超平面（因为它们都能把训练集中的点完美地分开）。
但是感知机存在一个很重要的问题，我们只用sign(w*x+b)输出的+1和-1来判断点的类别是不是太简单了？是不是有点生硬的感觉。
这么简单的判别方式真的会很有效吗？
当然了，虽然我们已经程序测试过正确率很高，但总是让人有点担心是否在很多情况下都能很好地工作。事实上我们从小到大一直会听到一些升学考试差一分两分的例子，那么差一分和高一分的学生真的就是天壤之别吗？感知器也是如此：在超平面左侧距离0.001的点和右边0.001距离的点输出就是+1和-1这天壤之别的差距真的合适吗？
此外我们也知道机器学习中通常会对目标函数进行微分进而梯度下降，但我们看上面这张图。很明显x=0是跳跃间断点，因此sign是一条不光滑的函数，没有办法进行微分。emmm……咦？那我记得感知器用了梯度下降，它是怎么去梯度下降的？
我们回忆一下感知器的梯度下降方式，确实用了梯度下降，但是发现没有，它把sign的壳子去掉了，对sign内部进行梯度下降，相对于其他能直接微分的算法来说，感知器的这种方式确实有点不太好。
因此既然我们要对感知器进行优化，那么上文的这两个主要的缺陷咱们得想想法子看看能不能弥补。我们首先试试解决第一个问题：
1.怎么解决极小距离差别带来的+1和-1的天壤之别？
在逻辑斯蒂回归中大致思想与感知器相同，但在输出+1与-1之间存在一些差别。在朴素贝叶斯中我们提到过P(Y|X)，它表示在给定X条件下，Y发生的概率。逻辑斯蒂也使用了同样的概念，它使用p(Y=-1|X)和p(Y=1|X)来表示该样本分别为-1或1的概率（实际上逻辑斯蒂并非强制要求标签必须为1或-1，可以用任意标签来表示）。这样当再出现样本X1距离为-0.001时，可能P(Y=1|X1)=0.49，P(Y=0|X1)=0.51,那么我们觉得X1为0的概率更大一点，但我们同时也清楚程序可能并不太确定X1一定为0。
使用概率作为输出结果使得样本在距离很小的差别下不再强制地输出+1和-1这两种天壤之别的结果，而是通过概率的方式告诉你结果可能是多少，同时也告诉你预测的不确信程度。这样子看起来让人比较安心一点不是吗？
2.怎么让最终的预测式子可微呢？
虽然无法微分并不会阻碍感知器的正常工作（事实上只是避开了sign），但对于很多场合都需要微分的机器学习来说，能找到一个可以微分的最终式子是很重要的。为了解决第一个问题我们提出了一种概率输出模型，那么感知器的sign也需要被随之替换为一种能输出概率大小的函数。具体函数在下文中详细讲解，其中值得高兴的是，我们找到的概率输出式子是平滑的，可微的，所以第二个问题也就迎刃而解了。
逻辑斯蒂回归的数学角度（配合《统计学习方法》食用更佳）
与感知器一样，我们先不管内部怎么工作，最好能够得到一个函数f(x)，我们将样本x放进去以后，它告诉我们属于1的概率是多少。比如说f(y=1|x)=0.2，f(y=0|x)=0.8，我们就知道该样本的标签大概率是0。总结一下也就是分别比较两种类别的概率大小，概率大的那一方则作为预测类别进行输出。f函数的定义如下图所示。
下图就是f(wx+b)的图像，我们假设一下点X在超平面右边，那么距离应当为正，如果距离无穷大时，exp(w*x + b)无穷大，P(Y=1|x)=1，也就是概率为1，表极其确信。如果wx+b是一个很接近0的数，exp(w*x + b)接近1，P(Y=1|x)=0.5，表示两边都有可能，不太确定。我们对于每一个样本点分别计算属于两类的概率，概率较大的一方作为预测的输出。
既然模型有了，那么就需要计算参数以此来得到划分的超平面。在计算参数上一般使用最大似然函数，关于最大似然函数的介绍可以在网上进行相关浏览。由于我们使用的sigmod函数很光滑，是可微函数，因此可以直接对似然函数求导，进行梯度下降找到参数的最优点。
关于似然函数关于W的求导式可以看下图：
结语：逻辑斯蒂在思想上与感知器极为相近，很多细节问题在感知器中已有说明因此在本文中并未详细展开。如果对部分地方有疑问可查看感知器一章节，此外也可以评论或邮箱的形式联系我。
贴代码：（建议去本文最上方的github链接下载，有书中所有算法的实现以及详细注释）
# coding=utf-8
# Author:Dodo
# Date:2018-11-27
# Email:lvtengchao@pku.edu.cn
# Blog:www.pkudodo.com
'''
数据集：Mnist
训练集数量：60000
测试集数量：10000
------------------------------
运行结果：
正确率：98.91%
运行时长：59s
'''
import time
import numpy as np
def loadData(fileName):
'''
加载Mnist数据集
:param fileName:要加载的数据集路径
:return: list形式的数据集及标记
'''
# 存放数据及标记的list
dataList = []; labelList = []
# 打开文件
fr = open(fileName, 'r')
# 将文件按行读取
for line in fr.readlines():
# 对每一行数据按切割福','进行切割，返回字段列表
curLine = line.strip().split(',')
# Mnsit有0-9是个标记，由于是二分类任务，所以将标记0的作为1，其余为0
# 验证过<5为1 >5为0时正确率在90%左右，猜测是因为数多了以后，可能不同数的特征较乱，不能有效地计算出一个合理的超平面
# 查看了一下之前感知机的结果，以5为分界时正确率81，重新修改为0和其余数时正确率98.91%
# 看来如果样本标签比较杂的话，对于是否能有效地划分超平面确实存在很大影响
if int(curLine[0]) == 0:
labelList.append(1)
else:
labelList.append(0)
#存放标记
#[int(num) for num in curLine[1:]] -> 遍历每一行中除了以第一哥元素（标记）外将所有元素转换成int类型
#[int(num)/255 for num in curLine[1:]] -> 将所有数据除255归一化(非必须步骤，可以不归一化)
dataList.append([int(num)/255 for num in curLine[1:]])
# dataList.append([int(num) for num in curLine[1:]])
#返回data和label
return dataList, labelList
def predict(w, x):
'''
预测标签
:param w:训练过程中学到的w
:param x: 要预测的样本
:return: 预测结果
'''
#dot为两个向量的点积操作，计算得到w * x
wx = np.dot(w, x)
#计算标签为1的概率
#该公式参考“6.1.2 二项逻辑斯蒂回归模型”中的式6.5
P1 = np.exp(wx) / (1 + np.exp(wx))
#如果为1的概率大于0.5，返回1
if P1 >= 0.5:
return 1
#否则返回0
return 0
def logisticRegression(trainDataList, trainLabelList, iter = 200):
'''
逻辑斯蒂回归训练过程
:param trainDataList:训练集
:param trainLabelList: 标签集
:param iter: 迭代次数
:return: 习得的w
'''
#按照书本“6.1.2 二项逻辑斯蒂回归模型”中式6.5的规则，将w与b合在一起，
#此时x也需要添加一维，数值为1
#循环遍历每一个样本，并在其最后添加一个1
for i in range(len(trainDataList)):
trainDataList[i].append(1)
#将数据集由列表转换为数组形式，主要是后期涉及到向量的运算，统一转换成数组形式比较方便
trainDataList = np.array(trainDataList)
#初始化w，维数为样本x维数+1，+1的那一位是b，初始为0
w = np.zeros(trainDataList.shape[1])
#设置步长
h = 0.001
#迭代iter次进行随机梯度下降
for i in range(iter):
#每次迭代冲遍历一次所有样本，进行随机梯度下降
for j in range(trainDataList.shape[0]):
#随机梯度上升部分
#在“6.1.3 模型参数估计”一章中给出了似然函数，我们需要极大化似然函数
#但是似然函数由于有求和项，并不能直接对w求导得出最优w，所以针对似然函数求和
#部分中每一项进行单独地求导w，得到针对该样本的梯度，并进行梯度上升（因为是
#要求似然函数的极大值，所以是梯度上升，如果是极小值就梯度下降。梯度上升是
#加号，下降是减号）
#求和式中每一项单独对w求导结果为：xi * yi - (exp(w * xi) * xi) / (1 + exp(w * xi))
#如果对于该求导式有疑问可查看我的博客 www.pkudodo.com
#计算w * xi，因为后式中要计算两次该值，为了节约时间这里提前算出
#其实也可直接算出exp(wx)，为了读者能看得方便一点就这么写了，包括yi和xi都提前列出了
wx = np.dot(w, trainDataList[j])
yi = trainLabelList[j]
xi = trainDataList[j]
#梯度上升
w += h * (xi * yi - (np.exp(wx) * xi) / ( 1 + np.exp(wx)))
#返回学到的w
return w
def test(testDataList, testLabelList, w):
'''
验证
:param testDataList:测试集
:param testLabelList: 测试集标签
:param w: 训练过程中学到的w
:return: 正确率
'''
#与训练过程一致，先将所有的样本添加一维，值为1，理由请查看训练函数
for i in range(len(testDataList)):
testDataList[i].append(1)
#错误值计数
errorCnt = 0
#对于测试集中每一个测试样本进行验证
for i in range(len(testDataList)):
#如果标记与预测不一致，错误值加1
if testLabelList[i] != predict(w, testDataList[i]):
errorCnt += 1
#返回准确率
return 1 - errorCnt / len(testDataList)
if __name__ == '__main__':
start = time.time()
# 获取训练集及标签
print('start read transSet')
trainData, trainLabel = loadData('../Mnist/mnist_train.csv')
# 获取测试集及标签
print('start read testSet')
testData, testLabel = loadData('../Mnist/mnist_test.csv')
# 开始训练，学习w
print('start to train')
w = logisticRegression(trainData, trainLabel)
#验证正确率
print('start to test')
accuracy = test(testData, testLabel, w)
# 打印准确率
print('the accuracy is:', accuracy)
# 打印时间
print('time span:', time.time() - start)
13条评论