**前言：**
生活中，人们经常会遇到各种最优化问题，比如如何在最短时间从一个地点到另外一个地点？如何在投入最少的资金而却能得到最高的受益？如何设计一款芯片使其功耗最低而性能最好？这一节就要学习一种最优化算法——Logistic回归，设计最优化算法的目的依然是用于分类。在这里，**Logistic回归的主要思想是根据现有的数据对分类边界线建立回归公式，达到分类的目的。**假设我们有一堆数据，需要划一条线(最佳直线)对其分类，这就是Logistic回归的目的。
而“Logistic回归”中的“回归”又是代表什么？数学认为，回归是一个拟合过程，回归分析本质上就是一个函数估计的问题，就是找出因变量和自变量之间的因果关系。具体到例子，假设我们有一些数据点，现在使用一条直线对这些点进行拟合，使得这条线尽可能地表示数据点的分布，这个拟合过程就称作“回归”。
在机器学习任务中，训练分类器时就是寻找最佳拟合曲线的过程，因此接下来将使用最优化算法。在实现算法之前，先总结Logistic回归的一些属性：
* 优点：计算代价不高，易于理解和实现
* 缺点：容易欠拟合，分类精度可能不高
* 适用的数据类型：数值型和标称型数据
**目录：**
一、基于Sigmoid函数的Logistic回归
1. Sigmoid函数
2. 基于最优化方法的最佳回归系数确定
3. 梯度上升法的实现
4. 随机梯度上升算法的实现
二、一个实例：从疝气病症预测病马的死亡率
1. 准备数据
2. 测试算法，使用Logistic回归进行分类
三、总结
===============================================
**一、基于Sigmoid函数的Logistic回归**
**1.Sigmoid函数**
Logistic回归想要得到的函数是，能接受所有的输入然后返回预测的类别，比如，在两类情况下函数应输出类别0或1。sigmoid函数可是胜任这一工作，它像是一个阶跃函数。其公式如下：
其中：
向量***w***称为回归系数，也就是我们要找到的最佳参数， ***x***是`n`维的特征向量，是分类器的输入数据。
下面附上书本上的图，该图是函数在不同坐标尺度下的曲线图：
为了实现Logistic回归分类器，我们可以在每个特征上乘以一个回归系数，然后把所有的结果值相加，将这个总和结果代入sigmoid函数中，进而得到一个范围在`0~1`之间的数值。任何大于`0.5`的数据被分入`1`类，小于`0.5`的数据被归入`0`类。所以，Logistic回归也可以被看成是一种概率估计。
**2.基于最优化方法的最佳回归系数确定**
上面提到Sigmoid函数里的一部分：
其中向量***w***称为回归系数，也就是我们要找到的最佳参数， ***x***是`n`维的特征向量，是分类器的输入数据。接下来将介绍几种需找最佳参数的方法：
* 梯度上升法：基于的思想是要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。
* 梯度下降法：与梯度上升的思想相似，但方向相反，即如果要找到某函数的最小值，最好的方法是沿着该函数的梯度方向的反方向寻找。
这里使用梯度上升法，对于一个函数`f(x,y)`，其梯度表示方法如下：
该梯度意味着要沿`x`和`y`的方向分别移动一定的距离，这其实是**确立了算法到达每个点后下一步移动的方向**。其中，函数`f(x,y)` 必须要在待计算的点上有定义并且可微。
**移动方向**确定了，这里我们定义移动的大小为步长，用`α`表示，用向量来表示的话，梯度上升算法的迭代公式如下：
该公式表明，最佳参数***w***的结果收到梯度变化和步长`α`的影响，这个公式将一直被迭代执行，直到满足某个停止条件。
**3.梯度上升法的实现**
以下使用Python实现了梯度上升法，并找到了最佳参数，其中使用的样本数据有`100`个，每个数据包含两个数值型的特征：`X1`和`X2` ，这些数据点可以被绘制成二维散点图。在找出最佳参数***w***后，再编写函数画出不同类别数据之间的分割线，观察最优化算法的效果。
~~~
# -*- coding: utf-8 -*-
"""
Created on Sat Sep 12 22:53:07 2015
@author: Herbert
"""
from numpy import *
sys.setrecursionlimit(1500)
def loadData():
dataMat = []; labelMat = []
fr = open('testSet.txt')
for line in fr.readlines():
lineSplit = line.strip().split()
dataMat.append([1.0, float(lineSplit[0]), float(lineSplit[1])])
labelMat.append(int(lineSplit[2]))
return dataMat, labelMat
def sigmoid(x):
return 1.0 / (1 + exp(-x))
def gradAscent(data, label):
dataMat = mat(data)
labelMat = mat(label).transpose()
m, n = shape(dataMat)
alpha = 0.001
maxCycles = 500
w = ones((n, 1))
for k in range(maxCycles):
p = sigmoid(dataMat * w)
error = labelMat - p
w = w + alpha * dataMat.transpose() * error
return dataMat, labelMat, w
data, label = loadData()
dataMat, labelMat, w = gradAscent(data, label)
print labelMat
def plotBestSplit(w):
import matplotlib.pyplot as plt
dataMat, labelMat = loadData()
dataArr = array(dataMat)
n = shape(dataArr)[0]
xcord1 = []; ycord1 = []
xcord2 = []; ycord2 = []
for i in range(n):
if int(labelMat[i]) == 1:
xcord1.append(dataArr[i, 1]) # array only
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 = 30, c = 'red', marker = 's')
ax.scatter(xcord2, ycord2, s = 30, c = 'blue')
x = arange(-3.0, 3.0, 0.1)
y = (-w[0] - w[1] * x) / w[2]
ax.plot(x, y)
plt.xlabel('x1'); plt.ylabel('x2');
plt.show()
plotBestSplit(w.getA())
~~~
结果如下，可看到使用了Logistic回归的分类结果还是挺不错的，尽管有三四个样本点被错分：
**4.随机梯度上升算法的实现**
梯度上升算法在处理`100`个左右的数据集时尚可，但如果有数十亿样本和成千上万的特征，那么该方法的计算复杂度将变得很恐怖。改进方法为随机梯度上升算法，该方法一次仅用一个样本点来更新回归系数。它占用更少的计算资源，是一种在线算法，可以在数据到来时就完成参数的更新，而不需要重新读取整个数据集来进行批处理运算。一次处理所有的数据被称为批处理。以下代码实现了随机梯度上升算法：
~~~
# 随机梯度上升算法
def stocGradAscent(data, label):
m,n = shape(data)
alpha = 0.01
w = ones(n)
for i in range(m):
h = sigmoid(sum(data[i]) * w)
error = label[i] - h
w = w + alpha * error *data[i]
return w
wNew = stocGradAscent(data, label)
plotBestSplit(wNew)
~~~
结果：
由上图可以看出，随机梯度上升算法分错了很多样本点，其分类效果并没有原先实现的普通梯度上升算法好。
书中给出的解释是，普通梯度上升算法是在整个数据集上迭代`500`次得到的结果，而随机梯度上升算法只是迭代了`100`次。一个判断算法优劣的可靠方法是看它是否收敛，也就是说求解的参数是否达到了稳定值，是否还会不断变化。
书中给出了一个例子，让随机梯度上升算法在整个数据集上运行`200`次，迭代过程中 **w向量** 的3个参数的变化如下图：
图中，**w**的第二个参数最先达到稳定，而其他两个参数则还需要更多的迭代次数来达到稳定。而我们也发现，在整个迭代过程中，**w**的三个参数都有不同程度的波动。产生这种现象的原因是存在一些被错分的样本点，这些样本在参与每次迭代的过程中都会引起参数的剧烈变化。为了避免这种剧烈波动，并改善算法的性能，采用以下策略对随机梯度上升算法做了改进：
* 在每次迭代时更新步长alpha的值，使得alpha的值不断减小，但是不能减少为`0`，这样做的原因是为了保证在多次迭代后新数据仍然具有一定的影响。
* 通过随机采取样本来更新回归参数，这种方法将减少周期性的波动，**每次随机从列表中取一个值，然后该值会从列表中删除。**
实现代码：
~~~
def stocGradAscentAdvanced(data, label, numIter = 150):
m,n = shape(data)
w = ones(n)
for i in range(numIter):
dataIndex = 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(data[randIndex] * w))
error = label[randIndex] - h
w = w + alpha * error* array(data[randIndex])
del(dataIndex[randIndex])
return w
wNewAdv = stocGradAscentAdvanced(data, label, numIter = 150)
plotBestSplit(wNewAdv)
~~~
在使用了优化策略后，各个回归系数参数的变化情况大幅改善，参数的收敛得更快了：
使用改进后的随机梯度上升算法对样本点进行划分，效果与普通梯度上升算法相当，但其中所用的计算量更少了：
**二、一个实例：从疝气病症预测病马的死亡率**
**1.准备数据**
该实例使用Logistic回归来预测患有疝病的马的存活问题。这里的数据来自2010年1月11日的UCI机器学习数据库，其中包含`368`个样本和**28**个特征。**这里的数据集是有`30%`的数据缺失的**。-.-
在实现算法之前，先了解一些关于处理数据中缺失值的方法：
* 使用可用特征的均值来填补缺失值；
* 使用特殊值来填补缺失值，如-1；
* 忽略有缺失值的样本；
* 使用相似样本的均值来填补缺失值；
* 使用另外的机器学习算法预测缺失值；
* 对于类别标签丢失的数据，只能将该数据丢弃。
这里我们使用实数`0`来替换所有缺失的值，这种方法恰好适用于Logistic回归；而对于类别标签丢失的数据，则须将该数据丢弃。
**2.测试算法，使用Logistic回归进行分类**
基于以上工作，以下代码把测试集上的每个特征向量乘以最优化算法得来的回归系数w，再将该乘积结果求和，最后输入到Sigmoid函数，如果对应的Sigmoid函数值大于`0.5`，则将该样本的类别判定为`1`，否则判定为`0`；最后，统计判定结果与实际结果的误差，由于误差具有不确定性，程序最后使用了`10`次分类结果求平均的方法，得出算法的平均分类错误率。
~~~
def classifyVector(x, w):
prob = sigmoid(sum(x * w))
if prob > 0.5: return 1.0
else: return 0.0
def horseColicTest():
frTrain = open('horseColicTraining.txt')
frTest = open('horseColicTest.txt')
trainingData = []; trainingLabel = []
for line in frTrain.readlines():
currLine = line.strip().split('\t')
lineArr = []
for i in range(21):
lineArr.append(float(currLine[i]))
trainingData.append(lineArr)
trainingLabel.append(float(currLine[21]))
w = stocGradAscentAdvanced(array(trainingData), trainingLabel, numIter = 500)
errorCount = 0.0; numOfTest = 0.0
for line in frTest.readlines():
numOfTest += 1.0
currLine = line.strip().split('\t')
lineArr = []
for i in range(21):
lineArr.append(float(currLine[i]))
if int(classifyVector(array(lineArr), w)) != int(currLine[21]):
errorCount += 1
errorRate = float(errorCount) / numOfTest
print "The error rate of this test is: %f" % errorRate
return errorRate
def finalTest():
numOfTest = 10; errorSum = 0.0
for i in range(numOfTest):
errorSum += horseColicTest()
print "After %d iterations the average error rate is: %f" \
% (numOfTest, errorSum / float(numOfTest))
finalTest()
~~~
分类错误率是`39.4%`，乍一看是挺高的。但如果你知道所使用的数据集是有`30%`的数据缺失的，你就不会这么想了…-.-
**三、总结**
Logistic回归的目的是寻找一个非线性函数sigmoid的最佳拟合参数，求解过程可以由最优化算法来完成。在最优化算法中，最常用的就是梯度上升算法，而梯度上升算法又可以简化为随机梯度上升算法。
随机梯度上升算法和梯度上升算法的效果相当，但占用更少的计算资源。随机梯度上升算法是一种在线算法，可以在数据到来时就完成参数的更新，而不需要重新读取整个数据集来进行批处理运算。
