Machine Learning in Action Chapter05 Logistic regression

Machine Learning in Action Chapter05 Logistic regression

一般过程:

1
2
3
4
5
6
(1) 收集数据:采用任意方法收集数据。  
(2) 准备数据:由于需要进行距离计算,因此要求数据类型为数值型。另外,结构化数据 格式则最佳。
(3) 分析数据:采用任意方法对数据进行分析。
(4) 训练算法:大部分时间将用于训练,训练的目的是为了找到最佳的分类回归系数。
(5) 测试算法:一旦训练步骤完成,分类将会很快。
(6) 使用算法:首先,我们需要输入一些数据,并将其转换成对应的结构化数值; 接着,基于训练好的回归系数就可以对这些数值进行简单的回归计算,判定它们属于 哪个类别;在这之后,我们就可以在输出的类别上做一些其他分析工作。

基于 Logistic 回归和 Sigmoid 函数的分类

  • 优点:计算代价不高,易于理解和实现
  • 缺点:容易欠拟合,分类精度可能不高
  • 使用数据类型:数值型和标称型数据

Sigmoid 函数:\(\sigma (z)=\frac{1}{1+e^{-z}}\)

为了实现 Logistic 回归分类器,可以在每一个特征上都乘以一个回归系数,然后所有的结果相加,将这个总和带入到 Sigmoid 函数中。进而得到一个范围在 0~1 之间的数值。任何大于 0.5 的数据被分入 1 类,小于 0.5 即被归入 0 类。所以,Logistic 回归也可以被看成是一种概 率估计。

基于最优化方法的最佳回归系数确定

Sigmoid 函数的输入记为 z,由下面公式得出:

\[z=w_0x_0+w_1x_1+w_2x_2+\dots+w_nx_n\]

采用向量写法,即 $z=w^Tx $

梯度上升法

要找到某函数的最大值,最好的方法是沿着该函数的梯度方向探寻。如果梯度记为 ▽ ,则函数 f(x, y) 的梯度由下式表示:

这个梯度意味着,要沿 x 的方向移动 $ \(,沿 y 的方向移动\) $。其中,函数 f(x,y) 保持在待计算的点上有定义且可微,具体例子如下:

梯度上升算法的迭代公式如下:

\[w := w+\alpha \nabla_wf(w)\]

直到达到某个停止条件为止,比如迭代次数达到某个指定值或算法达到某个可以允许的误差范围。

Note: 我们常听到的是梯度下降算法,它与这里的梯度上升算法是一样的,只是公式中的加法需要变成减法。因此,对应的公式可以写成:

\[w := w-\alpha \nabla_wf(w)\]

一个 Logistic 回归分类器的应用例子:

5.2.2 训练算法:使用梯度上升找到最佳参数

梯度上升法的伪代码如下:

1
2
3
4
5
每个回归系数初始化为 1
重复 R 次:
计算整个数据集的梯度
使用 alpha X gradient 更新回归系数的向量
返回回归系数

具体代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
def load_data_set():
"""
加载数据集
:return:返回两个数组,普通数组
data_arr -- 原始数据的特征
label_arr -- 原始数据的标签,也就是每条样本对应的类别
"""
data_arr = []
label_arr = []
f = open('./datasets/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]))
return data_arr, label_arr

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

5.2.3 分析数据:画出决策边界

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
def plot_best_fit(weights):
"""
可视化
:param weights:
:return:
"""
import matplotlib.pyplot as plt
data_mat, label_mat = load_data_set()
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)
y = (-weights[0] - weights[1] * x) / weights[2]
"""
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()

训练算法:随机梯度上升

梯度上升算法在每次更新回归系数时都需要遍历整个数据集,在数据集过大时计算成本高。一种改进方法就是一次仅用一个样本点来更新回归系数,该方法称为随机梯度上升法。由于可以在新 样本到来时对分类器进行增量式更新,因而随机梯度上升算法是一个在线学习算法。与“在线学 习”相对应,一次处理所有数据被称作是“批处理”。

伪代码如下:

1
2
3
4
5
所有回归系数初始化为1 
对数据集中每个样本
计算该样本的梯度
使用 alpha × gradient 更新回归系数值
返回回归系数值

实现代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
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

所得线性拟合图效果不错,分错了三分之一的样本。未随机的梯度上升算法在整个数据集上迭代了 500 次才得到只错 2 个的结果。不断调整训练次数,观察三个回归系数的变化情况:

X2 只经过了 50 次迭代就达到了稳定值,但系数 1 和 0 则需要更多次的迭代。另外值 得注意的是,在大的波动停止后,还有一些小的周期性波动。不难理解,产生这种现象的原因是 存在一些不能正确分类的样本点(数据集并非线性可分),在每次迭代时会引发系数的剧烈改变。 我们期望算法能避免来回波动,从而收敛到某个值。另外,收敛速度也需要加快。

改进版的随机梯度下降算法:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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

训练效果如下:

模型效果:

示例:从疝气病预测病马的死亡率

实现流程:

1
2
3
4
5
6
(1) 收集数据:给定数据文件。  
(2) 准备数据:用 Python 解析文本文件并填充缺失值。
(3) 分析数据:可视化并观察数据。
(4) 训练算法:使用优化算法,找到最佳的系数。
(5) 测试算法:为了量化回归的效果,需要观察错误率。根据错误率决定是否回退到训练阶段,通过改变迭代的次数和步长等参数来得到更好的回归系数。
(6) 使用算法:实现一个简单的命令行程序来收集马的症状并输出预测结果并非难事,这可以做为留给读者的一道习题

准备数据:处理数据中的缺失值

  • 使用可用特征的均值来填补缺失值;
  • 使用特殊值来填补缺失值,如 -1;
  • 忽略有缺失值的样本;
  • 使用有相似样本的均值添补缺失值;
  • 使用另外的机器学习算法预测缺失值。

测试算法:用 Logistic 进行分类

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
def classify_vector(in_x, weights):
"""
最终的分类函数,根据回归系数和特征向量来计算 Sigmoid 的值,大于0.5函数返回1,否则返回0
:param in_x: 特征向量,features
:param weights: 根据梯度下降/随机梯度下降 计算得到的回归系数
:return:
"""
# print(np.sum(in_x * weights))
prob = sigmoid(np.sum(in_x * weights))
if prob > 0.5:
return 1.0
return 0.0

def colic_test():
"""
打开测试集和训练集,并对数据进行格式化处理,其实最主要的的部分,比如缺失值的补充(真的需要学会的),人家已经做了
:return:
"""
f_train = open('./datasets/HorseColicTraining.txt', 'r')
f_test = open('./datasets/HorseColicTest.txt', 'r')
training_set = []
training_labels = []
# 解析训练数据集中的数据特征和Labels
# trainingSet 中存储训练数据集的特征,trainingLabels 存储训练数据集的样本对应的分类标签
for line in f_train.readlines():
curr_line = line.strip().split('\t')
if len(curr_line) == 1:
continue # 这里如果就一个空的元素,则跳过本次循环
line_arr = [float(curr_line[i]) for i in range(21)]
training_set.append(line_arr)
training_labels.append(float(curr_line[21]))
# 使用 改进后的 随机梯度下降算法 求得在此数据集上的最佳回归系数 trainWeights
train_weights = stoc_grad_ascent1(np.array(training_set), training_labels, 500)
error_count = 0
num_test_vec = 0.0
# 读取 测试数据集 进行测试,计算分类错误的样本条数和最终的错误率
for line in f_test.readlines():
num_test_vec += 1
curr_line = line.strip().split('\t')
if len(curr_line) == 1:
continue # 这里如果就一个空的元素,则跳过本次循环
line_arr = [float(curr_line[i]) for i in range(21)]
if int(classify_vector(np.array(line_arr), train_weights)) != int(curr_line[21]):
error_count += 1
error_rate = error_count / num_test_vec
print('the error rate is {}'.format(error_rate))
return error_rate

def multi_test():
"""
调用 colicTest() 10次并求结果的平均值
:return: nothing
"""
num_tests = 10
error_sum = 0
for k in range(num_tests):
error_sum += colic_test()
print('after {} iteration the average error rate is {}'.format(num_tests, error_sum / num_tests))

从上面的结果可以看到,10 次迭代之后的平均错误率为 35%。事实上,这个结果并不差,因为有 30%的数据缺失。当然,如果调整 colicTest() 中的迭代次数和 stochGradAscent1() 中的步长,平均错误率可以降到20%左右。

Comments

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×