本次编程作业的实现环境是Python3、Anaconda3(64-bit)、Jupyter Notebook。是在深度之眼“机器学习训练营”作业基础上完成的,个别代码有修改,供交流学习之用。
编程作业 3 - 多类分类¶
对于此练习,我们将使用逻辑回归来识别手写数字(0到9)。 我们将扩展我们在练习2中写的逻辑回归的实现,并将其应用于一对一的分类。 让我们开始加载数据集。 它是在MATLAB的本机格式,所以要加载它在Python,我们需要使用一个SciPy工具。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.io import loadmat
data = loadmat('ex3data1.mat')
data
Out[84]:
{'__header__': b'MATLAB 5.0 MAT-file, Platform: GLNXA64, Created on: Sun Oct 16 13:09:09 2011',
'__version__': '1.0',
'__globals__': [],
'X': array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]]),
'y': array([[10],
[10],
[10],
...,
[ 9],
[ 9],
[ 9]], dtype=uint8)}
In [85]:
data['X'].shape, type(data['X']), data['y'].shape, type(data['y'])
#len(data['X'])
Out[85]:
((5000, 400), numpy.ndarray, (5000, 1), numpy.ndarray)
sigmoid 函数
def sigmoid(z):
return 1 / (1 + np.exp(-z))
代价函数
def cost(theta, X, y, learningRate):
# INPUT:参数值theta,数据X,标签y,学习率
# OUTPUT:当前参数值下的交叉熵损失
# TODO:根据参数和输入的数据计算交叉熵损失函数
# STEP1:将theta, X, y转换为numpy类型的矩阵
# your code here (appro ~ 3 lines)
theta = np.matrix(theta)
X = np.matrix(X)
y = np.matrix(y)
# STEP2:根据公式计算损失函数(不含正则化)
# your code here (appro ~ 2 lines)
cross_cost = -1/len(X) * (y.T * np.log(sigmoid(X * theta.T)) + (1 - y.T) * np.log(1-sigmoid(X * theta.T)))
#cross_cost = -1/len(X) * np.sum(np.multiply(y, np.log(sigmoid(X * theta.T))) + np.multiply((1 - y), np.log(1 - sigmoid(X * theta.T))))
# STEP3:根据公式计算损失函数中的正则化部分
# your code here (appro ~ 1 lines),应使theta0为0(theta从1开始计算?)
reg = learningRate / (2*len(X)) * (np.power(theta[1:], 2).sum())
#reg = (learningRate / (2 * len(X))) * np.sum(np.power(theta[1:], 2))
# STEP4:把上两步当中的结果加起来得到整体损失函数
# your code here (appro ~ 1 lines)
whole_cost = cross_cost + reg
return whole_cost
向量化的梯度函数
#X = np.insert(data['X'], 0, values=np.ones(5000), axis=1)
#y = data['y']
#theta = np.zeros(401)
#learningRate =1
def gradient(theta, X, y, learningRate):
# INPUT:参数值theta,数据X,标签y,学习率
# OUTPUT:当前参数值下的梯度
# TODO:根据参数和输入的数据计算梯度
# STEP1:将theta, X, y转换为numpy类型的矩阵
# your code here (appro ~ 3 lines)
theta = np.matrix(theta)
X = np.matrix(X)
y = np.matrix(y)
# STEP2:将theta矩阵拉直(转换为一个向量)
# your code here (appro ~ 1 lines)
#parameters = theta.ravel()
# STEP3:计算预测的误差
# your code here (appro ~ 1 lines)
error = sigmoid(X @ theta.T) - y
# STEP4:根据上面的公式计算梯度
# your code here (appro ~ 1 lines)
grad1 = (1/len(X)) * X.T @ error
#grad = ((X.T * error) / len(X)).T + ((learningRate / len(X)) * theta)
# STEP5:由于j=0时不需要正则化,所以这里重置一下
# your code here (appro ~ 1 lines)
reg = learningRate / len(X) * theta[0,1:].T
grad1[1:,0] += reg
grad = grad1
#grad[0, 0] = np.sum(np.multiply(error, X[:,0])) / len(X)
return np.array(grad).ravel()
构建分类器:
有10个可能的类,并且由于逻辑回归只能一次在2个类之间进行分类,我们需要多类分类的策略。 在本练习中,实现一对一全分类方法,其中具有k个不同类的标签就有k个分类器,每个分类器在“类别 i”和“不是 i”之间决定。 我们将把分类器训练包含在一个函数中,该函数计算10个分类器中的每个分类器的最终权重,并将权重返回为k X(n + 1)数组,其中n是参数数量。
from scipy.optimize import minimize
def one_vs_all(X, y, num_labels, learning_rate):
rows = X.shape[0]
params = X.shape[1]
# k X (n + 1) array for the parameters of each of the k classifiers
all_theta = np.zeros((num_labels, params + 1))
# insert a column of ones at the beginning for the intercept term
X = np.insert(X, 0, values=np.ones(rows), axis=1)
# labels are 1-indexed instead of 0-indexed
for i in range(1, num_labels + 1):
theta = np.zeros(params + 1)
y_i = np.array([1 if label == i else 0 for label in y])
y_i = np.reshape(y_i, (rows, 1)) #变为rows行1列的形式
# minimize the objective function
fmin = minimize(fun=cost, x0=theta, args=(X, y_i, learning_rate), method='TNC', jac=gradient, options={'disp': True})
all_theta[i-1,:] = fmin.x
return all_theta
实现向量化代码的一个更具挑战性的部分是正确地写入所有的矩阵,保证维度正确。
rows = data['X'].shape[0]
params = data['X'].shape[1]
all_theta = np.zeros((10, params + 1))
X = np.insert(data['X'], 0, values=np.ones(rows), axis=1)
theta = np.zeros(params + 1)
y_0 = np.array([1 if label == 0 else 0 for label in data['y']])
y_0 = np.reshape(y_0, (rows, 1))
X.shape, type(X), y_0.shape, type(y_0), theta.shape, type(theta), all_theta.shape, type(all_theta)
Out[90]:
((5000, 401),
numpy.ndarray,
(5000, 1),
numpy.ndarray,
(401,),
numpy.ndarray,
(10, 401),
numpy.ndarray)
In [91]:
np.unique(data['y'])#看下有几类标签
Out[91]:
array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dtype=uint8)
让我们确保我们的训练函数正确运行,并且得到合理的输出。
data['X'].shape, type(data['X']), data['y'].shape, type(data['y'])
Out[92]:
((5000, 400), numpy.ndarray, (5000, 1), numpy.ndarray)
In [93]:
all_theta = one_vs_all(data['X'], data['y'], 10, 1)
print('theta:',all_theta.shape,)
最后一步 - 使用训练完毕的分类器预测每个图像的标签。 对于这一步,我们将计算每个类的类概率,对于每个训练样本(使用当然的向量化代码),并将输出类标签为具有最高概率的类。
def predict_all(X, all_theta):
# INPUT:参数值theta,测试数据X
# OUTPUT:预测值
# TODO:对测试数据进行预测
# STEP1:获取矩阵的维度信息
rows = X.shape[0]
#params = X.shape[1]
#num_labels = all_theta.shape[0]
# STEP2:把矩阵X加入一列壹元素
# your code here (appro ~ 1 lines)
X = np.insert(X, 0, values=np.ones(rows), axis=1)
# STEP3:把矩阵X和all_theta转换为numpy型矩阵
# your code here (appro ~ 2 lines)
#X = np.matrix(X)
#all_theta = np.matrix(theta)
# STEP4:计算样本属于每一类的概率
# your code here (appro ~ 1 lines)
h = sigmoid(X @ all_theta.T)
# STEP5:找到每个样本中预测概率最大的值
# your code here (appro ~ 1 lines)
h_argmax = np.argmax(h, axis=1)
# STEP6:因为我们的数组是零索引的,所以我们需要为真正的标签+1
h_argmax = h_argmax + 1
return h_argmax
现在我们可以使用predict_all函数为每个实例生成类预测,看看我们的分类器是如何工作的。
y_pred = predict_all(data['X'], all_theta)
correct = [1 if a == b else 0 for (a, b) in zip(y_pred, data['y'])]
accuracy = (sum(map(int, correct)) / float(len(correct)))
print ('accuracy = {0}%'.format(accuracy * 100))
输出的准确率应为94.46%。