Logistic 回归 概述
Logistic 回归 或者叫逻辑回归 虽然名字有回归，但是它是用来做分类的。其主要思想是: 根据现有数据对分类边界线(Decision Boundary)建立回归公式，以此进行分类。
须知概念
Sigmoid 函数
回归 概念
假设现在有一些数据点，我们用一条直线对这些点进行拟合（这条直线称为最佳拟合直线），这个拟合的过程就叫做回归。进而可以得到对这些点的拟合直线方程，那么我们根据这个回归方程，怎么进行分类呢？请看下面。
二值型输出分类函数
我们想要的函数应该是: 能接受所有的输入然后预测出类别。例如，在两个类的情况下，上述函数输出 0 或 1.或许你之前接触过具有这种性质的函数，该函数称为 海维塞得阶跃函数(Heaviside step function)，或者直接称为 单位阶跃函数。然而，海维塞得阶跃函数的问题在于: 该函数在跳跃点上从 0 瞬间跳跃到 1，这个瞬间跳跃过程有时很难处理。幸好，另一个函数也有类似的性质（可以输出 0 或者 1 的性质），且数学上更易处理，这就是 Sigmoid 函数。 Sigmoid 函数具体的计算公式如下:
下图给出了 Sigmoid 函数在不同坐标尺度下的两条曲线图。当 x 为 0 时，Sigmoid 函数值为 0.5 。随着 x 的增大，对应的 Sigmoid 值将逼近于 1 ; 而随着 x 的减小， Sigmoid 值将逼近于 0 。如果横坐标刻度足够大， Sigmoid 函数看起来很像一个阶跃函数。
因此，为了实现 Logistic 回归分类器，我们可以在每个特征上都乘以一个回归系数（如下公式所示），然后把所有结果值相加，将这个总和代入 Sigmoid 函数中，进而得到一个范围在 0~1 之间的数值。任何大于 0.5 的数据被分入 1 类，小于 0.5 即被归入 0 类。所以，Logistic 回归也是一种概率估计，比如这里Sigmoid 函数得出的值为0.5，可以理解为给定数据和参数，数据被分入 1 类的概率为0.5。想对Sigmoid 函数有更多了解，可以点开此链接跟此函数互动。
基于最优化方法的回归系数确定
Sigmoid 函数的输入记为 z ，由下面公式得到:
如果采用向量的写法，上述公式可以写成 
，它表示将这两个数值向量对应元素相乘然后全部加起来即得到 z 值。其中的向量 x 是分类器的输入数据，向量 w 也就是我们要找到的最佳参数（系数），从而使得分类器尽可能地精确。为了寻找该最佳参数，需要用到最优化理论的一些知识。我们这里使用的是——梯度上升法（Gradient Ascent）。
梯度上升法
梯度上升法的思想
要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。如果梯度记为 ▽ ，则函数 f(x, y) 的梯度由下式表示:
这个梯度意味着要沿 x 的方向移动 
，沿 y 的方向移动 
。其中，函数f(x, y) 必须要在待计算的点上有定义并且可微。下图是一个具体的例子。
上图展示的，梯度上升算法到达每个点后都会重新估计移动的方向。从 P0 开始，计算完该点的梯度，函数就根据梯度移动到下一点 P1。在 P1 点，梯度再次被重新计算，并沿着新的梯度方向移动到 P2 。如此循环迭代，直到满足停止条件。迭代过程中，梯度算子总是保证我们能选取到最佳的移动方向。
上图中的梯度上升算法沿梯度方向移动了一步。可以看到，梯度算子总是指向函数值增长最快的方向。这里所说的是移动方向，而未提到移动量的大小。该量值称为步长，记作 α 。用向量来表示的话，梯度上升算法的迭代公式如下:
该公式将一直被迭代执行，直至达到某个停止条件为止，比如迭代次数达到某个指定值或者算法达到某个可以允许的误差范围。
介绍一下几个相关的概念：
例如：y = w0 + w1x1 + w2x2 + ... + wnxn
梯度：参考上图的例子，二维图像，x方向代表第一个系数，也就是 w1，y方向代表第二个系数也就是 w2，这样的向量就是梯度。
α：上面的梯度算法的迭代公式中的阿尔法，这个代表的是移动步长（step length）。移动步长会影响最终结果的拟合程度，最好的方法就是随着迭代次数更改移动步长。
步长通俗的理解，100米，如果我一步走10米，我需要走10步；如果一步走20米，我只需要走5步。这里的一步走多少米就是步长的意思。
▽f(w)：代表沿着梯度变化的方向。
问：有人会好奇为什么有些书籍上说的是梯度下降法（Gradient Decent）?
答： 其实这个两个方法在此情况下本质上是相同的。关键在于代价函数（cost function）或者叫目标函数（objective function）。如果目标函数是损失函数，那就是最小化损失函数来求函数的最小值，就用梯度下降。 如果目标函数是似然函数（Likelihood function），就是要最大化似然函数来求函数的最大值，那就用梯度上升。在逻辑回归中， 损失函数和似然函数无非就是互为正负关系。
只需要在迭代公式中的加法变成减法。因此，对应的公式可以写成
局部最优现象 （Local Optima）
上图表示参数 θ 与误差函数 J(θ) 的关系图 (这里的误差函数是损失函数，所以我们要最小化损失函数)，红色的部分是表示 J(θ) 有着比较高的取值，我们需要的是，能够让 J(θ) 的值尽量的低。也就是深蓝色的部分。θ0，θ1 表示 θ 向量的两个维度（此处的θ0，θ1是x0和x1的系数，也对应的是上文w0和w1）。
可能梯度下降的最终点并非是全局最小点，可能是一个局部最小点，如我们上图中的右边的梯度下降曲线，描述的是最终到达一个局部最小点，这是我们重新选择了一个初始点得到的。
看来我们这个算法将会在很大的程度上被初始点的选择影响而陷入局部最小点。
Logistic 回归 原理
Logistic 回归 工作原理
每个回归系数初始化为 1
重复 R 次:
计算整个数据集的梯度
使用 步长 x 梯度 更新回归系数的向量
返回回归系数
Logistic 回归 开发流程
收集数据: 采用任意方法收集数据
准备数据: 由于需要进行距离计算，因此要求数据类型为数值型。另外，结构化数据格式则最佳。
分析数据: 采用任意方法对数据进行分析。
训练算法: 大部分时间将用于训练，训练的目的是为了找到最佳的分类回归系数。
测试算法: 一旦训练步骤完成，分类将会很快。
使用算法: 首先，我们需要输入一些数据，并将其转换成对应的结构化数值；接着，基于训练好的回归系数就可以对这些数值进行简单的回归计算，判定它们属于哪个类别；在这之后，我们就可以在输出的类别上做一些其他分析工作。
Logistic 回归 算法特点
优点: 计算代价不高，易于理解和实现。
缺点: 容易欠拟合，分类精度可能不高。
适用数据类型: 数值型和标称型数据。
附加 方向导数与梯度
Logistic 回归 项目案例
项目概述
在一个简单的数据集上，采用梯度上升法找到 Logistic 回归分类器在此数据集上的最佳回归系数
开发流程
收集数据: 可以使用任何方法
准备数据: 由于需要进行距离计算，因此要求数据类型为数值型。另外，结构化数据格式则最佳
分析数据: 画出决策边界
训练算法: 使用梯度上升找到最佳参数
测试算法: 使用 Logistic 回归进行分类
使用算法: 对简单数据集中数据进行分类
收集数据: 可以使用任何方法
我们采用存储在 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
绘制在图中，如下图所示:
准备数据: 由于需要进行距离计算，因此要求数据类型为数值型。另外，结构化数据格式则最佳
import numpy as np data_arr = [] label_arr = [] f = open('D:\mlInAction\data\5.Logistic\TestSet.txt', 'r') for line in f.readlines(): line_arr = line.strip().split() # 为了方便计算，我们将 X0 的值设为 1.0 ，也就是在每一行的开头添加一个 1.0 作为 X0 data_arr.append([1.0, np.float(line_arr[0]), np.float(line_arr[1])]) label_arr.append(int(line_arr[2]))
分析数据: 采用任意方法对数据进行分析，此处不需要
训练算法: 使用梯度上升找到最佳参数
定义sigmoid阶跃函数
def sigmoid(x): # 这里其实非常有必要解释一下，会出现的错误 RuntimeWarning: overflow encountered in exp # 这个错误在学习阶段虽然可以忽略，但是我们至少应该知道为什么 # 这里是因为我们输入的有的 x 实在是太小了，比如 -6000之类的，那么计算一个数字 np.exp(6000)这个结果太大了，没法表示，所以就溢出了 # 如果是计算 np.exp（-6000），这样虽然也会溢出，但是这是下溢，就是表示成零 # 去网上搜了很多方法，比如 使用bigfloat这个库（我竟然没有安装成功，就不尝试了，反正应该是有用的 return 1.0 / (1 + np.exp(-x))
Logistic 回归梯度上升优化算法
def grad_ascent(data_arr, class_labels): """ 梯度上升法，其实就是因为使用了极大似然估计，这个大家有必要去看推导，只看代码感觉不太够 :param data_arr: 传入的就是一个普通的数组，当然你传入一个二维的ndarray也行 :param class_labels: class_labels 是类别标签，它是一个 1*100 的行向量。 为了便于矩阵计算，需要将该行向量转换为列向量，做法是将原向量转置，再将它赋值给label_mat :return: """ # 注意一下，我把原来 data_mat_in 改成data_arr,因为传进来的是一个数组，用这个比较不容易搞混 # turn the data_arr to numpy matrix data_mat = np.mat(data_arr) # 变成矩阵之后进行转置 label_mat = np.mat(class_labels).transpose() # m->数据量，样本数 n->特征数 m, n = np.shape(data_mat) # 学习率，learning rate alpha = 0.001 # 最大迭代次数，假装迭代这么多次就能收敛2333 max_cycles = 500 # 生成一个长度和特征数相同的矩阵，此处n为3 -> [[1],[1],[1]] # weights 代表回归系数， 此处的 ones((n,1)) 创建一个长度和特征数相同的矩阵，其中的数全部都是 1 weights = np.ones((n, 1)) for k in range(max_cycles): # 这里是点乘 m x 3 dot 3 x 1 h = sigmoid(data_mat * weights) error = label_mat - h # 这里比较建议看一下推导，为什么这么做可以，这里已经是求导之后的 weights = weights + alpha * data_mat.transpose() * error return weights
大家看到这儿可能会有一些疑惑，就是，我们在迭代中更新我们的回归系数，后边的部分是怎么计算出来的？为什么会是 alpha * dataMatrix.transpose() * error ?因为这就是我们所求的梯度，也就是对 f(w) 对 w 求一阶导数。具体推导如下:
可参考http://blog.csdn.net/achuo/article/details/51160101
画出数据集和 Logistic 回归最佳拟合直线的函数
import matplotlib.pyplot as plt def plot_best_fit(data_mat, label_mat, weights): """ 可视化 :param weights: :return: """ data_arr = np.array(data_mat) n = np.shape(data_mat)[0] x_cord1 = [] y_cord1 = [] x_cord2 = [] y_cord2 = [] for i in range(n): if int(label_mat[i]) == 1: x_cord1.append(data_arr[i, 1]) y_cord1.append(data_arr[i, 2]) else: x_cord2.append(data_arr[i, 1]) y_cord2.append(data_arr[i, 2]) fig = plt.figure() ax = fig.add_subplot(111) ax.scatter(x_cord1, y_cord1, s=30, color='k', marker='^') ax.scatter(x_cord2, y_cord2, s=30, color='red', marker='s') x = np.arange(-3.0, 3.0, 0.1) # print(x) y = (-weights[0] - weights[1] * x) / weights[2] # type(y) y = np.ravel(y) # y原来是一个二维，需要转化为1维 """ y的由来，卧槽，是不是没看懂？ 首先理论上是这个样子的。 dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])]) w0*x0+w1*x1+w2*x2=f(x) x0最开始就设置为1叻， x2就是我们画图的y值，而f(x)被我们磨合误差给算到w0,w1,w2身上去了 所以： w0+w1*x+w2*y=0 => y = (-w0-w1*x)/w2 """ ax.plot(x, y) plt.xlabel('x1') plt.ylabel('y1') plt.show()
注意
梯度上升算法在每次更新回归系数时都需要遍历整个数据集，该方法在处理 100 个左右的数据集时尚可，但如果有数十亿样本和成千上万的特征，那么该方法的计算复杂度就太高了。一种改进方法是一次仅用一个样本点来更新回归系数，该方法称为 随机梯度上升算法。由于可以在新样本到来时对分类器进行增量式更新，因而随机梯度上升算法是一个在线学习(online learning)算法。与 “在线学习” 相对应，一次处理所有数据被称作是 “批处理” （batch） 。
随机梯度上升算法可以写成如下的伪代码:
所有回归系数初始化为 1 对数据集中每个样本 计算该样本的梯度 使用 alpha * gradient 更新回归系数值 返回回归系数值
以下是随机梯度上升算法的实现代码:
def stoc_grad_ascent0(data_mat, class_labels): """ 随机梯度上升，只使用一个样本点来更新回归系数 :param data_mat: 输入数据的数据特征（除去最后一列）,ndarray :param class_labels: 输入数据的类别标签（最后一列数据） :return: 得到的最佳回归系数 """ m, n = np.shape(data_mat) alpha = 0.01 weights = np.ones(n) for i in range(m): # sum(data_mat[i]*weights)为了求 f(x)的值， f(x)=a1*x1+b2*x2+..+nn*xn, # 此处求出的 h 是一个具体的数值，而不是一个矩阵 h = sigmoid(sum(data_mat[i] * weights)) error = class_labels[i] - h # 还是和上面一样，这个先去看推导，再写程序 weights = weights + alpha * error * data_mat[i] return weights
可以看到，随机梯度上升算法与梯度上升算法在代码上很相似，但也有一些区别: 第一，后者的变量 h 和误差 error 都是向量，而前者则全是数值；第二，前者没有矩阵的转换过程，所有变量的数据类型都是 NumPy 数组。
判断优化算法优劣的可靠方法是看它是否收敛，也就是说参数是否达到了稳定值，是否还会不断地变化？下图展示了随机梯度上升算法在 200 次迭代过程中回归系数的变化情况。其中的系数2，也就是 X2 只经过了 50 次迭代就达到了稳定值，但系数 1 和 0 则需要更多次的迭代。如下图所示:
针对这个问题，我们改进了之前的随机梯度上升算法，如下:
def stoc_grad_ascent1(data_mat, class_labels, num_iter=150): """ 改进版的随机梯度上升，使用随机的一个样本来更新回归系数 :param data_mat: 输入数据的数据特征（除去最后一列）,ndarray :param class_labels: 输入数据的类别标签（最后一列数据 :param num_iter: 迭代次数 :return: 得到的最佳回归系数 """ m, n = np.shape(data_mat) weights = np.ones(n) for j in range(num_iter): # 这里必须要用list，不然后面的del没法使用 data_index = list(range(m)) for i in range(m): # i和j的不断增大，导致alpha的值不断减少，但是不为0 alpha = 4 / (1.0 + j + i) + 0.01 # 随机产生一个 0～len()之间的一个值 # random.uniform(x, y) 方法将随机生成下一个实数，它在[x,y]范围内,x是这个范围内的最小值，y是这个范围内的最大值。 rand_index = int(np.random.uniform(0, len(data_index))) h = sigmoid(np.sum(data_mat[data_index[rand_index]] * weights)) error = class_labels[data_index[rand_index]] - h weights = weights + alpha * error * data_mat[data_index[rand_index]] del(data_index[rand_index]) return weights
import numpy as np
data_arr = []
label_arr = []
f = open('D:\mlInAction\data\5.Logistic\TestSet.txt', 'r')
for line in f.readlines():
line_arr = line.strip().split()
# 为了方便计算，我们将 X0 的值设为 1.0 ，也就是在每一行的开头添加一个 1.0 作为 X0
data_arr.append([1.0, np.float(line_arr[0]), np.float(line_arr[1])])
label_arr.append(int(line_arr[2]))
data_arr
label_arr
def sigmoid(x):
# 这里其实非常有必要解释一下，会出现的错误 RuntimeWarning: overflow encountered in exp
# 这个错误在学习阶段虽然可以忽略，但是我们至少应该知道为什么
# 这里是因为我们输入的有的 x 实在是太小了，比如 -6000之类的，那么计算一个数字 np.exp(6000)这个结果太大了，没法表示，所以就溢出了
# 如果是计算 np.exp（-6000），这样虽然也会溢出，但是这是下溢，就是表示成零
# 去网上搜了很多方法，比如 使用bigfloat这个库（我竟然没有安装成功，就不尝试了，反正应该是有用的
return 1.0 / (1 + np.exp(-x))
def grad_ascent(data_arr, class_labels):
"""
梯度上升法，其实就是因为使用了极大似然估计，这个大家有必要去看推导，只看代码感觉不太够
:param data_arr: 传入的就是一个普通的数组，当然你传入一个二维的ndarray也行
:param class_labels: class_labels 是类别标签，它是一个 1*100 的行向量。
为了便于矩阵计算，需要将该行向量转换为列向量，做法是将原向量转置，再将它赋值给label_mat
:return:
"""
# 注意一下，我把原来 data_mat_in 改成data_arr,因为传进来的是一个数组，用这个比较不容易搞混
# turn the data_arr to numpy matrix
data_mat = np.mat(data_arr)
# 变成矩阵之后进行转置
label_mat = np.mat(class_labels).transpose()
# m->数据量，样本数 n->特征数
m, n = np.shape(data_mat)
# 学习率，learning rate
alpha = 0.001
# 最大迭代次数，假装迭代这么多次就能收敛2333
max_cycles = 500
# 生成一个长度和特征数相同的矩阵，此处n为3 -> [[1],[1],[1]]
# weights 代表回归系数， 此处的 ones((n,1)) 创建一个长度和特征数相同的矩阵，其中的数全部都是 1
weights = np.ones((n, 1))
for k in range(max_cycles):
# 这里是点乘 m x 3 dot 3 x 1
h = sigmoid(data_mat * weights)
error = label_mat - h
# 这里比较建议看一下推导，为什么这么做可以，这里已经是求导之后的
weights = weights + alpha * data_mat.transpose() * error
return weights
weights = grad_ascent(data_arr, label_arr)
weights
import matplotlib.pyplot as plt
def plot_best_fit(data_mat, label_mat, weights):
"""
可视化
:param weights:
:return:
"""
data_arr = np.array(data_mat)
n = np.shape(data_mat)[0]
x_cord1 = []
y_cord1 = []
x_cord2 = []
y_cord2 = []
for i in range(n):
if int(label_mat[i]) == 1:
x_cord1.append(data_arr[i, 1])
y_cord1.append(data_arr[i, 2])
else:
x_cord2.append(data_arr[i, 1])
y_cord2.append(data_arr[i, 2])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(x_cord1, y_cord1, s=30, color='k', marker='^')
ax.scatter(x_cord2, y_cord2, s=30, color='red', marker='s')
x = np.arange(-3.0, 3.0, 0.1)
# print(x)
y = (-weights[0] - weights[1] * x) / weights[2]
# type(y)
y = np.ravel(y) # y原来是一个二维，需要转化为1维
"""
y的由来，卧槽，是不是没看懂？
首先理论上是这个样子的。
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
w0*x0+w1*x1+w2*x2=f(x)
x0最开始就设置为1叻， x2就是我们画图的y值，而f(x)被我们磨合误差给算到w0,w1,w2身上去了
所以： w0+w1*x+w2*y=0 => y = (-w0-w1*x)/w2
"""
ax.plot(x, y)
plt.xlabel('x1')
plt.ylabel('y1')
plt.show()
plot_best_fit(data_arr, label_arr, weights)
def stoc_grad_ascent0(data_mat, class_labels):
"""
随机梯度上升，只使用一个样本点来更新回归系数
:param data_mat: 输入数据的数据特征（除去最后一列）,ndarray
:param class_labels: 输入数据的类别标签（最后一列数据）
:return: 得到的最佳回归系数
"""
m, n = np.shape(data_mat)
alpha = 0.01
weights = np.ones(n)
for i in range(m):
# sum(data_mat[i]*weights)为了求 f(x)的值， f(x)=a1*x1+b2*x2+..+nn*xn,
# 此处求出的 h 是一个具体的数值，而不是一个矩阵
h = sigmoid(sum(data_mat[i] * weights))
error = class_labels[i] - h
# 还是和上面一样，这个先去看推导，再写程序
weights = weights + alpha * error * data_mat[i]
return weights
def stoc_grad_ascent1(data_mat, class_labels, num_iter=150):
"""
改进版的随机梯度上升，使用随机的一个样本来更新回归系数
:param data_mat: 输入数据的数据特征（除去最后一列）,ndarray
:param class_labels: 输入数据的类别标签（最后一列数据
:param num_iter: 迭代次数
:return: 得到的最佳回归系数
"""
m, n = np.shape(data_mat)
weights = np.ones(n)
for j in range(num_iter):
# 这里必须要用list，不然后面的del没法使用
data_index = list(range(m))
for i in range(m):
# i和j的不断增大，导致alpha的值不断减少，但是不为0
alpha = 4 / (1.0 + j + i) + 0.01
# 随机产生一个 0～len()之间的一个值
# random.uniform(x, y) 方法将随机生成下一个实数，它在[x,y]范围内,x是这个范围内的最小值，y是这个范围内的最大值。
rand_index = int(np.random.uniform(0, len(data_index)))
h = sigmoid(np.sum(data_mat[data_index[rand_index]] * weights))
error = class_labels[data_index[rand_index]] - h
weights = weights + alpha * error * data_mat[data_index[rand_index]]
del(data_index[rand_index])
return weights
weights1 = stoc_grad_ascent1(np.array(data_arr), np.array(label_arr))
weights1
plot_best_fit(data_arr, label_arr, weights1)
