文章目录
在逻辑回归+正则化的分类问题中,使用决策边界划分0和1,现在将完成多类分类(多个logistic回归)
作业题目:
在本练习中,您将实现一个one-vs-all逻辑回归和神经网络来识别手写数字(0-9)。
下图为本人的大致笔记思路
1.Dataset
在之前的学习中数据均为文本数据可以直接打开看,此次作业的数据是matlab的本机格式,需要使用scipy工具(loadmat)来加载数据。
import numpy as np
import pandas as pd
from scipy.io import loadmat
import matplotlib.pyplot as plt
data = loadmat('ex3data1.mat')
'''
{'__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)}
'''
可以看出y是标签,标签从10-1,分别代表数字9-0
同时查看一下X矩阵、y矩阵、标签分类的数目
def load_data(path):
data = loadmat(path)
X = data['X']
y = data['y']
return X,y
X, y = load_data('ex3data1.mat')
X.shape,y.shape,np.unique(y)
'''
((5000, 400),
(5000, 1),
array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dtype=uint8))
'''
表示有5000个样本,每个样本是20*20像素的数字的灰度图像。每个像素代表一个浮点数,表示该位置的灰度强度。20×20的像素网格被展开成一个400维的向量。在我们的数据矩阵X中,每一个样本都变成了一行,这给了我们一个5000×400矩阵X,每一行都是一个手写数字图像的训练样本,手写数字的像素点有值为255,空白像素点的位置为0
提示:
单个像素点通过8位(二进制数)的灰度值(0-255)来表示。例如一幅500*500像素的单通道灰度图是由500X500=250000个不同灰度的像素点组成。
2.Visualizing the data
随机打印一个数字
def plot_an_image(X):
pick_one = np.random.randint(0, 5000)
#从0-5000取出一个随机数
image = X[pick_one, :]
#取出该行的所有列(即展开的灰度值)
fig, ax = plt.subplots(figsize=(1, 1))
ax.matshow(image.reshape((20, 20)), cmap='gray_r')
plt.xticks([]) # 去除刻度,美观
plt.yticks([])
plt.show()
print('this should be {}'.format(y[pick_one]))
plot_an_image(X)
显示效果如下:
随机打印100个数字
def plot_100_image(X):
sample_idx = np.random.choice(np.arange(X.shape[0]), 100)
# 从X的所有row里面随机选100个样本即为sample_idx数列
sample_images = X[sample_idx, :] # (100,400)
fig, ax_array = plt.subplots(nrows=10, ncols=10, sharey=True, sharex=True, figsize=(8, 8))
#10×10
for row in range(10):
for column in range(10):
ax_array[row, column].matshow(sample_images[10 * row + column].reshape((20, 20)),
cmap='gray_r')
plt.xticks([])
plt.yticks([])
plt.show()
plot_100_image(X)
效果如下:
3.Vectorizing Logistic Regression
我们将使用多个one-vs-all(一对多)logistic回归模型来构建一个多类分类器。由于有10个类,需要训练10个独立的分类器。为了提高训练效率,重要的是向量化。所以,我们将实现一个不使用任何for循环的向量化的logistic回归版本。
正则化的logistic regression的代价函数:
对于每个样本i都要求sigmoid和h(x),事实上我们可以对所有的样本用矩阵乘法来快速的计算。让我们如下来定义 X 和 θ :
然后通过计算矩阵积X * θ ,我们可以得到:
在最后一个等式中,我们用到了一个定理,如果 a 和 b 都是向量,那么 a.T ✖ b = b.T ✖ a 这样我们就可以用一行代码计算出所有的样本。
同正则化逻辑回归一样
def sigmoid(z):
return 1 /(1 + np.exp(-z))
def regularized_cost(theta, X, y, l):
thetaReg = theta[1:]
first = (-y*np.log(sigmoid(X@theta))) + (y-1)*np.log(1-sigmoid(X@theta))
reg = (thetaReg@thetaReg)*l / (2*len(X))
return np.mean(first) + reg
4.Vectorizing the gradient
同正则化逻辑回归一样的
def regularized_gradient(theta, X, y, l):
thetaReg = theta[1:]
first = ( X.T @ (sigmoid(X @ theta) - y)) / len(X)
# 这里人为插入一维0,使得对theta_0不惩罚,方便计算
reg = np.concatenate([np.array([0]), (l / len(X)) * thetaReg])
return first + reg
4.One-vs-all Classification
这部分我们将实现一对多分类通过训练多个正则化logistic回归分类器。
对于这个任务,我们有10个可能的类,并且由于logistic回归只能一次在2个类之间进行分类,每个分类器在“类别 i”和“不是 i”之间决定。
from scipy.optimize import minimize
def one_vs_all(X, y, l, num_labels):#num_labels = 10
all_theta = np.zeros((num_labels, X.shape[1]))
# 大θ矩阵包含θ1-theta10维度为(10, 401)
for i in range(1, num_labels+1):
theta = np.zeros(X.shape[1])#求θi
y_i = np.array([1 if label == i else 0 for label in y])
ret = minimize(fun=regularized_cost, x0=theta, args=(X, y_i, l), method='TNC',
jac=regularized_gradient, options={'disp': True})
all_theta[i-1,:] = ret.x #将求好的θi赋值给大θ矩阵
return all_theta
此时的all_theta为10✖401维,每一行对应一个分类器的θi如(判别是不是数字1的θ1,判别是不是数字2的θ2)
初始化X,y,求出all_theta
raw_X, raw_y = load_data('ex3data1.mat')
X = np.insert(raw_X, 0, 1, axis=1) # 在第0列,插入数字1,axis控制维度为列(5000, 401)
y = raw_y.flatten() # 这里消除了一个维度,方便后面的计算 or .reshape(-1) (5000,)
all_theta = one_vs_all(X, y, 1, 10)
all_theta # 每一行是一个分类器的一组参数
'''
array([[-2.38187334e+00, 0.00000000e+00, 0.00000000e+00, ...,
1.30433279e-03, -7.29580949e-10, 0.00000000e+00],
[-3.18303389e+00, 0.00000000e+00, 0.00000000e+00, ...,
4.46340729e-03, -5.08870029e-04, 0.00000000e+00],
[-4.79638233e+00, 0.00000000e+00, 0.00000000e+00, ...,
-2.87468695e-05, -2.47395863e-07, 0.00000000e+00],
...,
[-7.98700752e+00, 0.00000000e+00, 0.00000000e+00, ...,
-8.94576566e-05, 7.21256372e-06, 0.00000000e+00],
[-4.57358931e+00, 0.00000000e+00, 0.00000000e+00, ...,
-1.33390955e-03, 9.96868542e-05, 0.00000000e+00],
[-5.40542751e+00, 0.00000000e+00, 0.00000000e+00, ...,
-1.16613537e-04, 7.88124085e-06, 0.00000000e+00]])
'''
接下来通过theta来求预估的hθ(x),便于后面与真实值做对比求模型的准确度
def predict_all(X, all_theta):
h = sigmoid(X @ all_theta.T)
# 注意的这里的all_theta需要转置
h_argmax = np.argmax(h, axis=1)
#找到每行的最大索引(0-9)
h_argmax = h_argmax + 1
# +1后变为每行的分类器 (1-10)
return h_argmax
求出模型的准确度
y_pred = predict_all(X, all_theta)
accuracy = np.mean(y_pred == y)
print ('accuracy = {0}%'.format(accuracy * 100))
#accuracy = 94.46%