1. 逻辑回归的模型函数
前面我们讲了线性回归模型，将线性模型用于回归问题中。这篇我们讲一下线性模型用于分类任务。在二分类问题中，对于线性回归所产生的预测值
我们需将这个预测z转化为0/1值，最理想的是“单位阶跃函数”，即若预测z大于零就判为1，若预测z小于零则判为反例，预测值为临界值则可以任意判别。但是由于单位阶跃函数不连续，我们希望可以找到一个连续的，同时在一定程度上近似“单位阶跃函数“的函数，那么逻辑斯特函数(logistic function)则是一个很不错的替代函数：
如下图所示，这是一种"Sigmoid"函数，它z值转化为一个接近0或1的y值。
我们将
代入上式，最终得到逻辑回归的模型函数，如下
作者个人理解：“单位阶跃函数”是将预测值直接转换为分类，而“逻辑斯特函数”是将预测值转换为0-1之间的一个数y，这个值正好可以当成某一类别的概率值，从而实现分类。
2. 逻辑回归的目标函数
逻辑回归的模型函数已经知道，求逻辑回归的模型函数，其实就是确定该模型函数theta这个参数的值。在讲如何求这个参数之前，需要求得逻辑回归的目标函数，而在求解目标函数的过程中，需先讲解一下以下知识点：
一个事件的几率（odds）：指该事件发生与不发生的概率比值，若事件发生概率p，那么事件发生的几率就是
那么该事件的对数几率（log odds,亦称logit）就是：
将式子
转化为如下
也就是说，输出y的对数几率是由输入x的线性函数表示的模型，这就是逻辑回归模型。当
的值越接近正无穷，y的值是越接近1的，那么我们将y视为样本x为正例，即y=1的概率。那么上述公式可转化为如下
从而得
进一步可转化为
在统计学中，常常使用极大似然估计法来求解，即找到一组参数，使得在这组参数下，我们的数据的似然度（概率）最大，似然函数如下：
取对数似然函数得
我们要最大化上式，也就等价于最小化
上式即可当做逻辑回归的目标函数
在逻辑回归模型中，我们最大化似然函数和最小化对数似然损失函数实际上是等价的。
3. 求解目标函数
上一块知识讲解中已经得到了目标函数，这一块知识讲解中，主要是讲解如何得到目标函数的最小值。自然而来的还是想到了梯度下降算法。我们求上述目标函数的梯度得
其中
那么对于
的一次迭代如下:
从而推出对于theta的一次迭代如下（并进行了向量化）
其中
上面将求和公式转化为了向量的形式，具体推导的过程已省略，可自行推导，也可以参考文末的梯度下降算法的向量化推导的参考链接。
4. 代码的实现
上面讲了逻辑回归的理论，下面我们通过代码来实现逻辑回归，例子主要是机器学习实战书上的例子。首先来看一下部分数据集（数据集的地址见文末给的链接地址），如下
-0.017612 14.053064 0
-1.395634 4.662541 1
-0.752157 6.538620 0
-1.322371 7.152853 0
...
-2.168791 0.143632 1
1.388610 9.341997 0
0.317029 14.739025 0我们需要从文件中读取数据，并返回读取之后的列表结果，实现的代码如下
'''
从文件中读取数据
'''
def loadDataSet(filePath):
dataList = []
labelList = []
f = open(filePath)
for line in f.readlines():
lineArr = line.strip().split()
dataList.append([1.0, float(lineArr[0]), float(lineArr[1])])
labelList.append(int(lineArr[2]))
return dataList, labelList之后我们采用梯度下降算法进行求解，在迭代过程中使用向量化，实现的代码如下：
'''
sigmoid函数
'''
def sigmoid(inX):
return 1.0 / (1 + np.exp(-inX))
'''
梯度下降函数
'''
def gradAscent(dataList, labelList):
dataMat = np.mat(dataList) # x数据转化为矩阵
labelMat = np.mat(labelList).T # y数据转换为矩阵并转置为列向量
m, n = np.shape(dataMat) # 返回矩阵的大小
alpha = 0.001 # 步长
maxCycles = 500 # 迭代次数
weights = np.ones((n, 1)) # 权重列向量
# 梯度下降算法
for k in range(maxCycles):
h = sigmoid(dataMat * weights)
error = h - labelMat
weights = weights - alpha * dataMat.T * error
return weights通过函数gradAscent可以得到一组回归系数，它确定了不同类别数据之间的分隔线，下面我们将这个分割线画出来，代码如下：
'''
画分类的分割线
'''
def plotBestFit(dataList, labelList, weights):
dataArr = np.array(dataList)
m = np.shape(dataArr)[0] # 获取样本的数量
xcord0 = []
ycord0 = []
xcord1 = []
ycord1 = []
# 对不同标记的数据进行分类
for i in range(m):
if int(labelList[i]) == 0:
xcord0.append(dataArr[i, 1])
ycord0.append(dataArr[i, 2])
else:
xcord1.append(dataArr[i, 1])
ycord1.append(dataArr[i, 2])
# 画散点图
plt.scatter(xcord0, ycord0, s=30, c='red', marker='s')
plt.scatter(xcord1, ycord1, s=30, c='green')
# 画分割线
x1 = np.arange(-3.0, 3.0, 0.1)
x2 = (-weights[0] - weights[1] * x1) / weights[2]
plt.plot(x1, x2)
plt.xlabel('X1')
plt.ylabel('X2')
plt.show()上述代码中，比较令人不解的是x2 = (-weights[0] - weights[1] * x1) / weights[2]这个是怎么来的？因为对于样本数据来说
是两个分类（类别1和类别0）的分界处，根据该式代入x_1即可求x_2从而这条线可以在图中表示出来。运行如下代码：
dataList, labelList = loadDataSet("testSet.txt")
weights = gradAscent(dataList, labelList)
plotBestFit(dataList, labelList, weights.getA())输出结果如下图所示，从图中可以看到分类效果相当不错，从图上看值只分错了几个点。
然而跟线性回归中使用梯度下降算法求最小值类似，这种方法需要大量的计算，每次更新回归系数都需要遍历整个数据集，假如数据量一多，该方法的计算复杂度就很高了。所以我们对上述的批量梯度下降算法进行改变，改进方法是一次仅用一个样本点来更新回归系数，该方法称为随机梯度下降算法。下面是随机梯度下降算法的实现代码
'''
随机梯度下降算法
'''
def stocGradAscent0(dataList, labelList):
dataArr = np.array(dataList) # x数据转化为矩阵
m, n = np.shape(dataArr)
alpha = 0.01
weights = np.ones(n)
for i in range(m):
h = sigmoid(np.sum(dataArr[i]*weights))
error = (h-labelList[i])
weights = weights-alpha*error*dataArr[i]
return weights运行如下代码，得到了该方法的分类效果图：
dataList, labelList = loadDataSet("testSet.txt")
weights = stocGradAscent0(dataList, labelList)
plotBestFit(dataList, labelList, weights)
相比批量梯度下降算法而言，效果较差。直接比较效果是有点不公平的，因为两者使用的数据集次数不同。而判断一个优化算法优劣的可靠方法是看它是否收敛，也就是说参数是否达到稳定值。通过下面这段代码，来绘制三个回归系数的变化情况：
'''
记录随机梯度下降过程中的历史权重值
'''
def stocGradAscent0WeightHistory(dataList, labelList):
dataMat = np.array(dataList) # x数据转化为矩阵
m, n = np.shape(dataMat)
alpha = 0.01
times = 100
weights = np.ones(n)
weightsHistory = np.zeros((times*m, n) )
for j in range(times):
for i in range(m):
h = sigmoid(np.sum(dataMat[i] * weights))
error = (h - labelList[i])
weights = weights - alpha * error * dataMat[i]
weightsHistory[j*m+i, :] = weights
return weightsHistory
'''
绘制历史权值图
'''
def plotSDGError(weightsHistory):
plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
fig = plt.figure()
ax = fig.add_subplot(311)
ax.plot(weightsHistory[:,0])
plt.ylabel('X0')
ax = fig.add_subplot(312)
ax.plot(weightsHistory[:,1])
plt.ylabel('X1')
ax = fig.add_subplot(313)
ax.plot(weightsHistory[:,2])
plt.ylabel('X2')
plt.xlabel(u'迭代次数')
plt.show()
def main():
dataList, labelList = loadDataSet("testSet.txt")
weightsHistory = stocGradAscent0WeightHistory(dataList, labelList)
plotSDGError(weightsHistory)
if __name__ == '__main__':
main()最终绘制的三个回归系数的变化情况如下：
从图中可以看到在大的波动之后,有些系数还会存在一些小的周期性波动。这种现象的原因是因为存在一些不能正确分类的样本点，在每次迭代时会引发系数的剧烈改变。我们希望算法能避免来回波动，从而收敛于某个值，同时是加快收敛速度，因此我们将上述的随机梯度下降算法进行更改。更改为如下：
'''
改进的随机梯度下降算法
'''
def stocGradAscent1(dataList, labelList, numIter=150):
dataArr = np.array(dataList)
m,n = np.shape(dataArr)
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.01 # 步长为动态变化
rand = int(np.random.uniform(0, len(dataIndex)))
choseIndex = dataIndex[rand]
h = sigmoid(np.sum(dataArr[choseIndex]*weights))
error = h-labelList[choseIndex]
weights = weights-alpha*error*dataArr[choseIndex]
del(dataIndex[rand])
return weights相对第一个的随机梯度下降算法，上述算法一方面使alpha在每次迭代的时候都会进行调整，缓解了数据波动或者高频波动。虽然alpha会随着迭代次数不断减小，但是不会减少至0，而这样子做也保证了在多次迭代之后新数据仍然具有一定的影响。另一方面，通过随机选取样本来更新回归系数，从而减少周期性的波动。同时每次随机从列表中选出一个值之后，就从列表中删除该值。
上述算法对机器学习实战中的使用的算法有进行更改，作者个人认为书上在选出一个值并删除该值之后的实现方式不对。按照书上实现的算法，比如我们dataIndex序列是[0, 1, 2, 3, 4, 5]，那么假如第一次随机产生的数是0，那么使用的则是dataArr[0]这个样本，之后使用del删除了0下标这个数。那么dataIndex序列变成了[1, 2, 3, 4, 5]。之后再随机产生一个数，比如还是0，那么使用的还是dataArr[0]这个样本，之后del删除了下标为0的这个数，那么序列就变成了[2, 3, 4, 5]了。那么这种实现方式跟书的作者意思就有点不一致了。书中作者的意思应该是从[0, 1, 2, 3, 4, 5]这个列表中随机选取一个值，比如选取了0，然后使用第0个样本，之后删除0这个数，那么下次就变成从[1, 2, 3, 4, 5]这个列表中随机选取一个值，比如选取了1，然后使用第1个样本，之后删掉1这个数。
运行如下代码，得到的效果如下图所示，该效果达到了跟批量梯度下降算法差不多的效果，但是所使用的计算机更少
dataList, labelList = loadDataSet("testSet.txt")
weights = stocGradAscent1(dataList, labelList)
plotBestFit(dataList, labelList, weights)
使用更改之后的算法，来显示每次迭代时各个回归系数的变化情况，实现的代码如下
'''
记录改进的随机梯度下降过程中的历史权重值
'''
def stocGradAscent1WeightHistory(dataList, labelList, numIter=150):
dataArr = np.array(dataList)
m, n = np.shape(dataArr)
weights = np.ones(n)
weightsHistory = np.zeros((numIter*m, n))
for j in range(numIter):
dataIndex = list(range(m))
for i in range(m):
alpha = 4 / (1.0 + j + i) + 0.01 # 步长为动态变化
rand = int(np.random.uniform(0, len(dataIndex)))
choseIndex = dataIndex[rand]
h = sigmoid(np.sum(dataArr[choseIndex] * weights))
error = h - labelList[choseIndex]
weights = weights - alpha * error * dataArr[choseIndex]
weightsHistory[j*m+i, :] = weights
del (dataIndex[rand])
return weightsHistory运行如下代码，得到的效果图如图所示。从图中可以看到没有出现周期性的波动，这归功于样本随机选择机制，同时收敛的速度也更快。
dataList, labelList = loadDataSet("testSet.txt")
weightsHistory = stocGradAscent1WeightHistory(dataList, labelList, 40)
plotSDGError(weightsHistory)
