Logistic回归
我的微信公众号： s406205391; 欢迎大家一起学习，一起进步！！！
假设现在有一些数据点，我们用一条直线对这些点进行拟合（该线称为最佳拟合直线），这个拟合过程就称作回归。利用Logistic回归进行分类的主要思想是：根据现有数据对分类界线建立回归公式，以此进行分类。训练分类器时的做法就是寻找最佳拟合参数，使用的是最优化算法。
Logistic回归的一般过程：
- 收集数据：采用任意方法收集数据。
- 准备数据：由于需要进行距离计算，因此要求数据类型为数值型。另外，结构化数据格式则最佳。
- 分析数据：采用任意方法对数据进行分析。
- 训练算法：大部分时间将用于训练，训练的目的是为了找到最佳的分类回归系数。
- 测试算法：一旦训练完成，分类将会很快。
- 使用算法：首先，我们需要输入一些数据，并将其转换成对应的结构化数字；接着，基于训练好的回归系数就可以对这些数值进行简单的回归计算，判定他们属于哪个类别；在这之后，我们就可以在输出的类别上做一些其他分析工作。
基于Logistic回归和Sigmoid函数的分类
Logistic回归的优先在于计算代价不高，易于理解和实现。缺点在于容易欠拟合，分类精度可能不高。适用数据类型为数值型和标称型数据。
Logistic回归中，使用了Sigmoid阶跃函数，该函数的特点是：该函数在跳跃点上可以从0瞬间跳跃到1。
基于这个函数，我们可以将回归方程，带入这个函数。进而得到一个范围在0~1之间的数字。任何大于0.5的数据被分为1类，小于0.5即被归入0类。这样便可以很好的完成分类。公式如下：
Sigmoid： <font size=2, color='grey'>z表示回归方程</font>
回归方程：z = , 也可以写成
接下来，我们要做的事情，便是找到w的最佳值，一般有梯度下降法、梯度上升法、随机梯度下降法、随机梯度上升法等。其最终推导公式如下：
<font size=2, color='grey'>g(x)即为Sigmoid函数，α为随机常数，</font>
具体推到过程，可以参考这篇文章。
梯度上升/下降法
梯度上升和梯度下降其实差不多，只需要把上面的公式的加号变成减号就行了。
本文所用的数据集(提取码：ogge)。下面是一个梯度上升算法的具体实现。load_data_set()函数，从文本中提取了一组数据，作为我们的训练的数据集。sigmoid()函数则是Sigmoid的计算。grad_ascent()函数则是梯度上升求参数最佳值的方法。其实主要就是运用了上面的推导公式。
import matplotlib.pyplot as plt
import numpy as np
def load_data_set():
"""读取测试数据集"""
data_mat = []
label_mat = []
for line in open('data\\Ch05\\testSet.txt'):
line = line.strip().split()
data_mat.append([1.0, float(line[0]), float(line[1])])
label_mat.append(int(line[2]))
return data_mat, label_mat
def sigmoid(in_x):
"""sigmoid阶跃函数"""
return 1 / (1 + np.exp(-in_x))
def grad_ascent(data_mat_in, class_labels):
"""使用梯度上升法，求得系数"""
data_matrix = np.mat(data_mat_in)
label_mat = np.mat(class_labels).transpose() # 转置矩阵
m, n = np.shape(data_matrix)
alpha = 0.001
max_cycles = 500
weights = np.ones((n, 1))
for k in range(max_cycles):
h = sigmoid(data_matrix * weights)
error = (label_mat - h)
weights = weights + alpha * data_matrix.transpose() * error
return weights
def plot_best_fit(weights):
data_mat, label_mat = load_data_set()
data_array = np.array(data_mat)
n = np.shape(data_array)[0]
xcord1 = []
ycord1 = []
xcord2 = []
ycord2 = []
for i in range(n):
if int(label_mat[i]) == 1:
xcord1.append(data_array[i, 1])
ycord1.append(data_array[i, 2])
else:
xcord2.append(data_array[i, 1])
ycord2.append(data_array[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='green')
x = np.arange(-3.0, 3.0, 0.1)
y = (-weights[0] - weights[1] * x) / weights[2]
ax.plot(x, y)
plt.xlabel('x1')
plt.ylabel('x2')
plt.show()
if __name__ == '__main__':
d, l = load_data_set()
w = grad_ascent(d, l)
plot_best_fit(w)
下图是梯度上升法，迭代500次的分类结果。从结果上看，分类的还是比较好的。
随机梯度上升法
梯度上升算法在每次更新回归系数时都需要遍历整个数据集，该方法在处理100个左右的数据集时尚可，但是如果有数十亿样本和成千上万个特征，那么该方法计算复杂度就太高了。一种改进的方法是一次仅用一个样本点来更新回归系数，该方法称为随机梯度上升算法。由于可以在新的样本到来时对分类器进行增量式更新，因而随机梯度上升算法是一个在线学习算法。下面是随机梯度算法的代码实现。
def sto_grad_ascent(data_matrix, class_labels):
"""随机梯度上升法"""
m, n = np.shape(data_matrix)
alpha = 0.001
weights = np.ones(n)
for i in range(m):
h = sigmoid(data_matrix[i] * weights)
error = (class_labels[i] - h)
weights = weights + alpha * data_matrix[i] * error
return weights
if __name__ == '__main__':
d, l = load_data_set()
w = sto_grad_ascent(np.array(d), l)
plot_best_fit(w)
下图是随机梯度上升法，拟合的结果。不如梯度上升法完美，但是也需要考虑，后者是迭代了500次才得到的。一个判断优化算法优劣的可靠方法是看它是否收敛，也就是说参数是否达到了稳定值，是否还会不断地变化。
改进的随机梯度上升算法
对随机梯度算法进行一些改进，得到下面的代码。其改进之处在于。一方面，alpha在每次迭代的时候都会调整，这会缓解随机梯度上升算法的数据波动。另外，虽然alpha会随着迭代次数不断减小，但永远不会减小到0，这是因为公式中存在常数项，必须保证在多次迭代后，该常数项仍然具有一定的影响。另一个改进之处在于，这里通过随机选取样本来更新回归系数。这种方法将减少周期性波动。
def stoc_grad_ascent(data_matrix, class_labels, num_iter=500):
"""改进的随机梯度上升算法"""
m, n = np.shape(data_matrix)
weights = np.ones(n)
for j in range(num_iter):
data_index = list(range(m))
for i in range(m):
alpha = 4 / (1 + j + i) + 0.01
rand_index = int(np.random.uniform(0, len(data_index)))
h = sigmoid(sum(data_matrix[rand_index] * weights))
error = class_labels[rand_index] - h
weights = weights + alpha * error * data_matrix[rand_index]
del data_index[rand_index]
return weights
改进之后的拟合结果，要明显优于改进之前得结果。
示例：从疝气病症预测病马病死率
接下来将使用Logistic回归来预测患有疝病的马的存活问题。这里的数据包含368个样本和28个特征。 我们在日常数据处理中，常会遇到数据集包含缺失值的情况，那么在Logistic回归算法中，缺失值是否可以用0来代替呢？
答案是可以的。我们观察回归系数的更新公式：
weights = weights + alpha * error * dataMatrix[randIndex]
如果dataMatrix的某个特征对应值为0，那么该特征的系数将不做更新，即：
weights = weights
另外，由于sigmoid(0)=0.5, 即它对结果的预测不具有任何倾向性。基于上述原因，将缺失值用0代替既可以保留现有数据，也不需要对优化的算法进行修改。
下面的代码中，classify_vector()构建了一个Logistic分类器，只需要将训练你的到的回归系数带入进去，就可以对样本进行分类了。colic_test()函数则是示例的具体代码
def classify_vector(in_x, weithts):
prob = sigmoid(sum(in_x * weithts))
if prob > 0.5:
return 1.0
else:
return 0.0
def colic_test():
fr_train = open('data\\Ch05\\horseColicTraining.txt')
fr_test = open('data\\Ch05\\horseColicTest.txt')
training_set = []
training_labels = []
for line in fr_train.readlines():
line = line.strip().split('\t')
line_arr = list(map(float, line[0: 21]))
training_set.append(line_arr)
training_labels.append(float(line[21]))
train_weights = stoc_grad_ascent(np.array(training_set), training_labels, 500)
error_count = 0
num_test_vec = 0
for line in fr_test.readlines():
num_test_vec += 1
line = line.strip().split('\t')
line_arr = list(map(float, line[0: 21]))
if int(classify_vector(np.array(line_arr), train_weights)) != int(line[21]):
error_count += 1
error_rate = error_count / num_test_vec
print(f"the error rate of this test is: {error_rate}")
return error_rate
if __name__ == '__main__':
colic_test()
