机器学习实战系列---Logistic回归
Logistic回归或者叫逻辑回归，虽然名字里有回归二字，但它是用来做分类的。其主要思想为：根据现有数据对分类界线建立回归公式，以此进行分类。
logistic回归是一种分类方法，常用于两分类问题。为概率型非线性回归模型，是研究二分类观察结果
相关概念
回归概念
假设现在有一些数据点，我们用一条直线对这些点进行拟合（这条直线称为最佳拟合直线），这个拟合的过程就叫做回归。进而可以得到对这些点的拟合直线方程，那么我们根据这个回归方程，怎么进行分类呢？请看下面。
二值输出分类函数
我们想要的函数应该是: 能接受所有的输入然后预测出类别。例如，在两个类的情况下，上述函数输出 0 或 1.或许你之前接触过具有这种性质的函数，该函数称为 海维塞得阶跃函数(Heaviside step function)，或者直接称为 单位阶跃函数。然而，海维塞得阶跃函数的问题在于: 该函数在跳跃点上从 0 瞬间跳跃到 1，这个瞬间跳跃过程有时很难处理。幸好，另一个函数也有类似的性质（可以输出 0 或者 1 的性质），且数学上更易处理，这就是 Sigmoid 函数。 Sigmoid 函数具体的计算公式如下:
因此，为了实现 Logistic 回归分类器，我们可以在每个特征上都乘以一个回归系数（如下公式所示），然后把所有结果值相加，将这个总和代入 Sigmoid 函数中，进而得到一个范围在 0~1 之间的数值。任何大于 0.5 的数据被分入 1 类，小于 0.5 即被归入 0 类。所以，Logistic 回归也是一种概率估计，比如这里Sigmoid 函数得出的值为0.5，可以理解为给定数据和参数，数据被分入 1 类的概率为0.5。
基于最优化方法的回归系数确定
设Sigmoid函数的输入为z，则z=w1x1+w2x2+w3x3+...+wnxn=WTX，其中向量X是分类器的输入数据，向量W就是我们要找到的最佳回归系数。为了寻找该最佳参数，需要用到最优化理论的一些知识。我们这里使用的是——梯度上升法（Gradient Ascent）。
梯度上升法与梯度下降法
要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。如果梯度记为 ▽ ，则函数 f(x, y) 的梯度由下式表示:
这个梯度意味着要沿 x 的方向移动 
，沿 y 的方向移动 
。其中，函数f(x, y) 必须要在待计算的点上有定义并且可微。
梯度算子总是指向函数值增长最快的方向。这里所说的是移动方向，而未提到移动量的大小。该量值称为步长，记作 α 。用向量来表示的话，梯度上升算法的迭代公式如下:
该公式将一直被迭代执行，直至达到某个停止条件为止，比如迭代次数达到某个指定值或者算法达到某个可以允许的误差范围。
问：有人会好奇为什么有些书籍上说的是梯度下降法（Gradient Decent）?
答： 其实这个两个方法在此情况下本质上是相同的。关键在于代价函数（cost function）或者叫目标函数（objective function）。如果目标函数是损失函数，那就是最小化损失函数来求函数的最小值，就用梯度下降。 如果目标函数是似然函数（Likelihood function），就是要最大化似然函数来求函数的最大值，那就用梯度上升。在逻辑回归中， 损失函数和似然函数无非就是互为正负关系。
只需要在迭代公式中的加法变成减法。因此，对应的公式可以写成
Logistic回归 代价函数
定义逻辑回归的代价函数为：
，其中
hΘ(x)与Cost(hΘ(x),y)之间的关系如下图：
这样构建的Cost(hΘ(x),y)的函数特点是：当实际的y=1且hΘ(x)也为1时误差为0，当y=1但hΘ(x)不为1时误差随着hΘ(x)变小而变大；当实际的y=0且hΘ(x)也为0时误差为0，当y=0但hΘ(x)不为0时误差随着hΘ(x)的变大而变大。
将构建的Cost(hΘ(x),y)简化如下：
带入代价函数得：
Logistic回归 工作原理
每个回归系数初始化为1
重复R次：
计算整个数据集的梯度
使用 步长*梯度 更新回归系数向量
返回回归系数
Logistic回归 算法特点
优点：计算代价不高，易于理解和实现。
缺点：容易欠拟合，分类精度不高。
使用数据类型：数值型和标称型数据。
有关回归系数W更新的公式推导
http://blog.csdn.net/achuo/article/details/51160101
项目案例
我们采用存储在 TestSet.txt 文本文件中的数据，存储格式如下:
-0.017612 14.053064 0
-1.395634 4.662541 1
-0.752157 6.538620 0
-1.322371 7.152853 0
0.423363 11.054677 0
绘制在图中，如图所示：
代码实现
from math import exp
from numpy import *
from matplotlib import pyplot as plt
from matplotlib.font_manager import FontProperties
def loadDataSet(filename):
dataMat = []
labelMat = []
f = open(filename)
for line in f.readlines():
lineArr = line.strip().split()
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
labelMat.append(int(lineArr[2]))
return dataMat, labelMat
def sigmoid(inX):
return 1.0/(1+exp(-inX))
def gradAscent(dataMatIn, classLabels):
'''
梯度上升算法
:param dataMatIn: 数据集，一个二维Numpy数组，每列分别代表每个不同的特征，每行则代表每个训练样本。
:param classLabels: 类别标签，是一个（1，100）的行向量。
:return: 返回回归系数
'''
dataMatrix = mat(dataMatIn)
labelMat = mat(classLabels).transpose()
# m->样本数，n->特征数 在本例中m=100,n=3
m, n = shape(dataMatrix)
# alpha代表向目标移动的步长
alpha = 0.001
# 迭代次数
maxCycles = 500
# 权重值初始全为1
weights = ones((n, 1))
for k in range(maxCycles):
h = sigmoid(dataMatrix * weights) # shape = (100,1)
error = labelMat - h # shape = (100,1)
weights += alpha * dataMatrix.transpose() * error
return weights
def stocGradAscent(dataMatIn, classLabels, numIter=150):
'''
随机梯度上升法
:param dataMatIn: 数据集，一个二维Numpy数组，每列分别代表每个不同的特征，每行则代表每个训练样本。
:param classLabels: 类别标签，是一个（1，100）的行向量。
:param numIter: 外循环次数
:return: 返回回归系数
'''
dataMatrix = mat(dataMatIn)
m, n = shape(dataMatrix)
weights = ones((n, 1))
# 随机梯度, 循环150,观察是否收敛
for j in range(numIter):
# [0, 1, 2 .. m-1]
dataIndex = list(range(m))
for i in range(m):
# i和j的不断增大，导致alpha的值不断减少，但是不为0
alpha = 4/(1.0+j+i)+0.0001 # alpha 会随着迭代不断减小，但永远不会减小到0，因为后边还有一个常数项0.0001
# 随机产生一个 0～len(dataIndex)之间的一个值
# random.uniform(x, y) 方法将随机生成下一个实数，它在[x,y]范围内,x是这个范围内的最小值，y是这个范围内的最大值。
randIndex = int(random.uniform(0, len(dataIndex)))
choose = dataIndex[randIndex]
# sum(dataMatrix[i]*weights)为了求 f(x)的值， f(x)=w1*x1+w2*x2+..+wn*xn
h = sigmoid(sum(dataMatrix[choose]*weights))
error = classLabels[choose] - h
weights += alpha * error * dataMatrix[choose].transpose()
del(dataIndex[randIndex])
return weights
def showData(dataArr, labelMat):
'''
展示数据集分布情况
:param dataArr:
:param labelMat:
:return: None
'''
n = shape(dataArr)[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(figsize=(8, 6))
ax = fig.add_subplot(111)
p1 = ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
p2 = ax.scatter(xcord2, ycord2, s=30, c='green')
plt.legend([p1, p2], ['Class 1', 'Class 0'], loc='lower right', scatterpoints=1)
plt.show()
def plotBestFit(dataArr, labelMat, weights1, weights2):
'''
将我们得到的数据可视化展示出来
:param dataArr:样本数据的特征
:param labelMat:样本数据的类别标签，即目标变量
:param weights:回归系数
:return:None
'''
font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
n = shape(dataArr)[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(figsize=(8, 6))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
ax1.scatter(xcord1, ycord1, s=30, c='red', marker='s')
ax1.scatter(xcord2, ycord2, s=30, c='green')
ax2.scatter(xcord1, ycord1, s=30, c='red', marker='s')
ax2.scatter(xcord2, ycord2, s=30, c='green')
x1 = arange(-3.0, 3.0, 0.1)
y1 = (-weights1[0] - weights1[1] * x1) / weights1[2]
x2 = arange(-3.0, 3.0, 0.1)
y2 = (-weights2[0] - weights2[1] * x2) / weights2[2]
ax1.plot(x1, y1)
ax2.plot(x2, y2)
ax1_title_text = ax1.set_title(u'梯度上升算法', FontProperties=font)
ax2_title_text = ax2.set_title(u'随机梯度上升算法', FontProperties=font)
plt.xlabel('X')
plt.ylabel('Y')
#plt.setp(ax1_title_text, size=20, weight='bold', color='black')
#plt.setp(ax2_title_text, size=20, weight='bold', color='black')
plt.show()
def testLR():
dataMat, classLabels = loadDataSet('data/TestSet.txt')
dataArr = array(dataMat)
weights1= gradAscent(dataArr, classLabels)
weights2 = stocGradAscent(dataArr, classLabels)
test(dataArr, classLabels)
plotBestFit(dataArr, classLabels, weights1, weights2)
if __name__ == '__main__':
testLR()
结果图：
