目录
一 前言
本文从Logistic回归的原理开始讲起，补充了书上省略的数学推导。本文可能会略显枯燥，理论居多，Sklearn实战内容会放在下一篇文章。自己慢慢推导完公式，还是蛮开心的一件事。
二 Logistic回归与梯度上升算法
Logistic回归是众多分类算法中的一员。通常，Logistic回归用于二分类问题，例如预测明天是否会下雨。当然它也可以用于多分类问题，不过为了简单起见，本文暂先讨论二分类问题。首先，让我们来了解一下，什么是Logistic回归。
1 Logistic回归
假设现在有一些数据点，我们利用一条直线对这些点进行拟合(该线称为最佳拟合直线)，这个拟合过程就称作为回归，如下图所示：
Logistic回归一种二分类算法，它利用的是Sigmoid函数阈值在[0,1]这个特性。Logistic回归进行分类的主要思想是：根据现有数据对分类边界线建立回归公式，以此进行分类。其实，Logistic本质上是一个基于条件概率的判别模型(Discriminative Model)。
所以要想了解Logistic回归，我们必须先看一看Sigmoid函数 ，我们也可以称它为Logistic函数。它的公式如下：
整合成一个公式，就变成了如下公式：
下面这张图片，为我们展示了Sigmoid函数的样子。
z是一个矩阵，θ是参数列向量(要求解的)，x是样本列向量(给定的数据集)。θ^T表示θ的转置。g(z)函数实现了任意实数到[0,1]的映射，这样我们的数据集([x0,x1,…,xn])，不管是大于1或者小于0，都可以映射到[0,1]区间进行分类。hθ(x)给出了输出为1的概率。比如当hθ(x)=0.7，那么说明有70%的概率输出为1。输出为0的概率是输出为1的补集，也就是30%。
如果我们有合适的参数列向量θ([θ0,θ1,…θn]^T)，以及样本列向量x([x0,x1,…,xn])，那么我们对样本x分类就可以通过上述公式计算出一个概率，如果这个概率大于0.5，我们就可以说样本是正样本，否则样本是负样本。
举个例子，对于"垃圾邮件判别问题"，对于给定的邮件(样本)，我们定义非垃圾邮件为正类，垃圾邮件为负类。我们通过计算出的概率值即可判定邮件是否是垃圾邮件。
那么问题来了！如何得到合适的参数向量θ?
根据sigmoid函数的特性，我们可以做出如下的假设：
代价函数：
我们要介绍如何拟合逻辑回归模型的参数𝜃。具体来说，我要定义用来 拟合参数的优化目标或者叫代价函数，这便是监督学习问题中的逻辑回归模型的拟合问题。
对于线性回归模型，我们定义的代价函数是所有模型误差的平方和。理论上来说，我们 也可以对逻辑回归模型沿用这个定义，但是问题在于，当我们将 
带入到这样定义了的代价函数中时，我们得到的代价函数将是一个非凸函数(non-convexfunction)。
这意味着我们的代价函数有许多局部最小值，这将影响梯度下降算法寻找全局最小值。
线性回归的代价函数为 :
h𝜃(𝑥)与 𝐶𝑜𝑠𝑡(h𝜃(𝑥),𝑦)之间的关系如下图所示:
这样构建的𝐶𝑜𝑠𝑡(h𝜃(𝑥),𝑦)函数的特点是:当实际的 𝑦 = 1 且h𝜃(𝑥)也为 1 时误差为 0， 当 𝑦 = 1 但h𝜃(𝑥)不为 1 时误差随着h𝜃(𝑥)变小而变大;当实际的 𝑦 = 0 且h𝜃(𝑥)也为 0 时 代价为 0，当𝑦 = 0 但h𝜃(𝑥)不为 0 时误差随着 h𝜃(𝑥)的变大而变大。
将构建的 𝐶𝑜𝑠𝑡(h𝜃(𝑥),𝑦)简化如下
𝐶𝑜𝑠𝑡(h𝜃(𝑥),𝑦) = −𝑦 × 𝑙𝑜𝑔(h𝜃(𝑥)) − (1 − 𝑦) × 𝑙𝑜𝑔(1 − h𝜃(𝑥))
带入代价函数得到:
其中，m为样本的总数，y(i)表示第i个样本的类别，x(i)表示第i个样本，需要注意的是θ是多维向量，x(i)也是多维向量。
在得到这样一个代价函数以后，我们便可以用梯度下降算法来求得能使代价函数最小的参数了。
怎么求解使J(θ)最小的θ值呢？因为是求最小值，所以我们需要使用梯度下降算法。如果面对的问题是求解使J(θ)最大的θ值，那么我们就需要使用梯度上升算法。面对我们这个问题，如果使J(θ) := -J(θ)，那么问题就从求极小值转换成求极大值了，使用的算法就从梯度下降算法变成了梯度上升算法，它们的思想都是相同的，学会其一，就也会了另一个。本文使用梯度上升算法进行求解。
2 梯度上升算法
说了半天，梯度上升算法又是啥？J(θ)太复杂，我们先看个简单的求极大值的例子。一个看了就会想到高中生活的函数：
来吧，做高中题。这个函数的极值怎么求？显然这个函数开口向下，存在极大值，它的函数图像为：
求极值，先求函数的导数：
令导数为0，可求出x=2即取得函数f(x)的极大值。极大值等于f(2)=4
但是真实环境中的函数不会像上面这么简单，就算求出了函数的导数，也很难精确计算出函数的极值。此时我们就可以用迭代的方法来做。就像爬坡一样，一点一点逼近极值。这种寻找最佳拟合参数的方法，就是最优化算法。爬坡这个动作用数学公式表达即为：
其中，α为步长，也就是学习速率，控制更新的幅度。效果如下图所示：
比如从(0,0)开始，迭代路径就是1->2->3->4->…->n，直到求出的x为函数极大值的近似值，停止迭代。我们可以编写Python3代码，来实现这一过程：
# -*- coding:UTF-8 -*-
"""
函数说明:梯度上升算法测试函数
求函数f(x) = -x^2 + 4x的极大值
"""
def Gradient_Ascent_test():
def f_prime(x_old): #f(x)的导数
return -2 * x_old + 4
x_old = -1 #初始值，给一个小于x_new的值
x_new = 0 #梯度上升算法初始值，即从(0,0)开始
alpha = 0.01 #步长，也就是学习速率，控制更新的幅度
presision = 0.00000001 #精度，也就是更新阈值
while abs(x_new - x_old) > presision:
x_old = x_new
x_new = x_old + alpha * f_prime(x_old) #上面提到的公式
print(x_new) #打印最终求解的极值近似值
if __name__ == '__main__':
Gradient_Ascent_test()
代码运行结果如下：
结果很显然，已经非常接近我们的真实极值2了。这一过程，就是梯度上升算法。那么同理，J(θ)这个函数的极值，也可以这么求解。公式可以这么写：
由上小节可知J(θ)为：
sigmoid函数为：
那么，现在我只要求出J(θ)的偏导，就可以利用梯度上升算法，求解J(θ)的极大值了。
那么现在开始求解J(θ)对θ的偏导，求解如下(数学推导)：
其中：
再由：
可得：
接下来，就剩下第三部分：
综上所述：
因此，梯度上升迭代公式为：
知道了，梯度上升迭代公式，我们就可以自己编写代码，计算最佳拟合参数了。
三 Python3实战
数据集已经为大家准备好，下载地址：https://github.com/Jack-Cherish/Machine-Learning/blob/master/Logistic/testSet.txt
这就是一个简单的数据集，没什么实际意义。让我们先从这个简单的数据集开始学习。先看下数据集有哪些数据：
这个数据有两维特征，因此可以将数据在一个二维平面上展示出来。我们可以将第一列数据(X1)看作x轴上的值，第二列数据(X2)看作y轴上的值。而最后一列数据即为分类标签。根据标签的不同，对这些点进行分类。
那么，先让我们编写代码，看下数据集的分布情况
# -*- coding: UTF-8 -*-
import matplotlib.pyplot as plt
import numpy as np
"""
函数说明:加载数据
Parameters:
无
Returns:
dataMat - 数据列表
labelMat - 标签列表
Author:
yangshaojun
Modify:
2020-02-01
"""
# 加载数据 特征变为3维 吴恩达通常令第一个特征为1 即w0
def loadDataSet():
dataMat = [] # 此列表存储特征
labelMat = [] # 此列表存标签
fr = open('testSet.txt') # 打开文件
for line in fr.readlines(): # 逐行读取
lineArr = line.strip().split('\t') # 去回车，放入列表
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])]) # 每个样本特征作为列表存储
labelMat.append(int(lineArr[2])) # 添加到标签列表中
fr.close() # 关闭文件
return dataMat, labelMat # 返回
"""
函数说明:绘制数据集
Parameters:
无
Returns:
无
Author:
yangshaojun
Modify:
2020-02-01
"""
def plotDataSet():
dataMat, labelMat = loadDataSet()
dataArr = np.array(dataMat) # 转换成numpy的array数组
n = np.shape(dataArr)[0] # 数据样本的数目
xcode1 = []; ycode1 = [] #正样本
xcode2 = []; ycode2 = [] #负样本
for i in range(n): #根据数据集标签进行分类
if int(labelMat[i]) == 1:
xcode1.append(dataArr[i, 1]); ycode1.append(dataArr[i, 2]) #1为正样本
else:
xcode2.append(dataArr[i, 1]); ycode2.append(dataArr[i, 2]) #0为负样本
fig = plt.figure()
ax = fig.add_subplot(111) #添加subplot
ax.scatter(xcode1, ycode1, s=20, c='red', marker='s', alpha=0.5) #绘制正样本
ax.scatter(xcode2, ycode2, s=20, c='green', alpha=0.5) #绘制负样本
plt.title('DataSet') #绘制title
plt.xlabel('x') #绘制label
plt.ylabel('y')
plt.show() #显示
运行结果如下：
从上图可以看出数据的分布情况。假设Sigmoid函数的输入记为z，那么z=w0x0 + w1x1 + w2x2，即可将数据分割开。其中，x0为全是1的向量，x1为数据集的第一列数据，x2为数据集的第二列数据。另z=0，则0=w0 + w1x1 + w2x2。横坐标为x1，纵坐标为x2。这个方程未知的参数为w0，w1，w2，也就是我们需要求的回归系数(最优参数)。
2 训练算法
在编写代码之前，让我们回顾下梯度上升迭代公式：
将上述公式矢量化：
根据矢量化的公式，编写代码如下：
# -*- coding: UTF-8 -*-
import matplotlib.pyplot as plt
import numpy as np
"""
函数说明:加载数据
Parameters:
无
Returns:
dataMat - 数据列表
labelMat - 标签列表
Author:
yangshaojun
Modify:
2020-02-01
"""
# 加载数据 特征变为3维 吴恩达通常令第一个特征为1 即w0
def loadDataSet():
dataMat = [] # 此列表存储特征
labelMat = [] # 此列表存标签
fr = open('testSet.txt') # 打开文件
for line in fr.readlines(): # 逐行读取
lineArr = line.strip().split('\t') # 去回车，放入列表
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])]) # 每个样本特征作为列表存储
labelMat.append(int(lineArr[2])) # 添加到标签列表中
fr.close() # 关闭文件
return dataMat, labelMat # 返回
"""
函数说明:sigmoid函数
Parameters:
inX - 数据
Returns:
sigmoid函数
Author:
yangshaojun
"""
def sigmoid(inX):
return 1.0 / (1 + np.exp(-inX))
"""
函数说明:梯度上升算法
Parameters:
dataMatIn - 数据集
classLabels - 数据标签
Returns:
weights.getA() - 求得的权重数组(最优参数)
[[ 4.12414349]
[ 0.48007329]
[-0.6168482 ]]
Author:
yangshaojun
Modify:
2020-02-01
"""
def gradAscent(dataMatIn, classLabels):
dataMatrix = np.mat(dataMatIn) # 列表转换成numpy的mat 100行3列
labelMat = np.mat(classLabels).transpose() #转换成numpy的mat,并进行转置 100行1列
m, n = np.shape(dataMatrix) # 返回dataMatrix的大小。m为行数,n为列数。
alpha = 0.001 #移动步长,也就是学习速率,控制更新的幅度。
maxCycles = 500 # 最大迭代次数
weights = np.ones((n, 1)) # 创建一个3行1列的矩阵 矩阵元素为1
for k in range(maxCycles):
h = sigmoid(dataMatrix * weights) #梯度上升矢量化公式
error = (labelMat - h)
weights = weights + alpha * dataMatrix.transpose() * error
return weights.getA() # 将矩阵转换为数组，返回权重数组
if __name__ == '__main__':
dataArr, labelMat = loadDataSet()
weights = gradAscent(dataArr, labelMat)
print(weights)运行结果如图所示：
可以看出，我们已经求解出回归系数[w0,w1,w2]。
通过求解出的参数，我们就可以确定不同类别数据之间的分隔线，画出决策边界。
# -*- coding: UTF-8 -*-
import matplotlib.pyplot as plt
import numpy as np
"""
函数说明:加载数据
Parameters:
无
Returns:
dataMat - 数据列表
labelMat - 标签列表
Author:
yangshaojun
Modify:
2020-02-01
"""
# 加载数据 特征变为3维 吴恩达通常令第一个特征为1 即w0
def loadDataSet():
dataMat = [] # 此列表存储特征
labelMat = [] # 此列表存标签
fr = open('testSet.txt') # 打开文件
for line in fr.readlines(): # 逐行读取
lineArr = line.strip().split('\t') # 去回车，放入列表
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])]) # 每个样本特征作为列表存储
labelMat.append(int(lineArr[2])) # 添加到标签列表中
fr.close() # 关闭文件
return dataMat, labelMat # 返回
"""
函数说明:sigmoid函数
Parameters:
inX - 数据
Returns:
sigmoid函数
Author:
yangshaojun
"""
def sigmoid(inX):
return 1.0 / (1 + np.exp(-inX))
"""
函数说明:梯度上升算法
Parameters:
dataMatIn - 数据集
classLabels - 数据标签
Returns:
weights.getA() - 求得的权重数组(最优参数)
[[ 4.12414349]
[ 0.48007329]
[-0.6168482 ]]
Author:
yangshaojun
Modify:
2020-02-01
"""
def gradAscent(dataMatIn, classLabels):
dataMatrix = np.mat(dataMatIn) # 列表转换成numpy的mat 100行3列
labelMat = np.mat(classLabels).transpose() #转换成numpy的mat,并进行转置 100行1列
m, n = np.shape(dataMatrix) # 返回dataMatrix的大小。m为行数,n为列数。
alpha = 0.001 #移动步长,也就是学习速率,控制更新的幅度。
maxCycles = 500 # 最大迭代次数
weights = np.ones((n, 1)) # 创建一个3行1列的矩阵 矩阵元素为1
for k in range(maxCycles):
h = sigmoid(dataMatrix * weights) #梯度上升矢量化公式
error = (labelMat - h)
weights = weights + alpha * dataMatrix.transpose() * error
return weights.getA() # 将矩阵转换为数组，返回权重数组
"""
函数说明:绘制数据集
Parameters:
weights - 权重参数数组
Returns:
无
Author:
yangshaojun
Modify:
2020-02-01
"""
def plotBestFit(weights):
dataMat, labelMat = loadDataSet()
dataArr = np.array(dataMat) #转换成numpy的array数组
n = np.shape(dataArr)[0] # 返回样本数
xcode1 = []; ycode1 = [] #正样本
xcode2 = []; ycode2 = [] #负样本
for i in range(n): #根据数据集标签进行分类
if int(labelMat[i]) == 1:
xcode1.append(dataArr[i, 1]); ycode1.append(dataArr[i, 2]) #1为正样本
else:
xcode2.append(dataArr[i, 1]); ycode2.append(dataArr[i, 2]) #0为负样本
fig = plt.figure()
ax = fig.add_subplot(111) #添加subplot
ax.scatter(xcode1, ycode1, s=30, c='red', marker='s') #绘制正样本
ax.scatter(xcode2, ycode2, s=30, c='green') #绘制负样本
x = np.arange(-3.0, 3.0, 0.1)
y = (-weights[0] - weights[1] * x) / weights[2] # 这里的y的值即x2
ax.plot(x, y)
plt.xlabel('X1');
plt.ylabel('X2');
plt.show()
这个分类结果相当不错，从上图可以看出，只分错了几个点而已。但是，尽管例子简单切数据集很小，但是这个方法却需要大量的计算(300次乘法)。因此下篇文章将对改算法稍作改进，从而减少计算量，使其可以应用于大数据集上。
四 总结
Logistic回归的一般过程：
收集数据：采用任意方法收集数据。
准备数据：由于需要进行距离计算，因此要求数据类型为数值型。另外，结构化数据格式则最佳。
分析数据：采用任意方法对数据进行分析。
训练算法：大部分时间将用于训练，训练的目的是为了找到最佳的分类回归系数。
测试算法：一旦训练步骤完成，分类将会很快。
使用算法：首先，我们需要输入一些数据，并将其转换成对应的结构化数值；接着，基于训练好的回归系数，就可以对这些数值进行简单的回归计算，判定它们属于哪个类别；在这之后，我们就可以在输出的类别上做一些其他分析工作。
其他：
Logistic回归的目的是寻找一个非线性函数Sigmoid的最佳拟合参数，求解过程可以由最优化算法完成。
本文讲述了Logistic回归原理以及数学推导过程。
下篇文章将讲解Logistic回归的改进以及Sklearn实战内容。
939
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
