使用逻辑回归来识别手写数字(0到9)
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.io import loadmat
data = loadmat('ex3data1.mat')
data
{'__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)}
data['X'].shape, data['y'].shape
((5000, 400), (5000, 1))
图像在martix X中表示为400维向量(其中有5,000个)。 400维“特征”是原始20 x 20图像中每个像素的灰度强度。类标签在向量y中作为表示图像中数字的数字类。
第一个任务是将我们的逻辑回归实现修改为完全向量化(即没有“for”循环)。这是因为向量化代码除了简洁外,还能够利用线性代数优化,并且通常比迭代代码快得多。但是,如果从练习2中看到我们的代价函数已经完全向量化实现了,所以我们可以在这里重复使用相同的实现。
sigmoid 函数
g 代表一个常用的逻辑函数(logistic function)为S形函数(Sigmoid function),公式为:
g
(
z
)
=
1
1
+
e
−
z
g\left( z \right)=\frac{1}{1+{{e}^{-z}}}
g(z)=1+e−z1
合起来,我们得到逻辑回归模型的假设函数:
h
θ
(
x
)
=
1
1
+
e
−
θ
T
X
{{h}_{\theta }}\left( x \right)=\frac{1}{1+{{e}^{-{{\theta }^{T}}X}}}
hθ(x)=1+e−θTX1
def sigmoid(z):
return 1 / (1 + np.exp(-z))
代价函数:
J ( θ ) = 1 m ∑ i = 1 m [ − y ( i ) log ( h θ ( x ( i ) ) ) − ( 1 − y ( i ) ) log ( 1 − h θ ( x ( i ) ) ) ] + λ 2 m ∑ j = 1 n θ j 2 J\left( \theta \right)=\frac{1}{m}\sum\limits_{i=1}^{m}{[-{{y}^{(i)}}\log \left( {{h}_{\theta }}\left( {{x}^{(i)}} \right) \right)-\left( 1-{{y}^{(i)}} \right)\log \left( 1-{{h}_{\theta }}\left( {{x}^{(i)}} \right) \right)]}+\frac{\lambda }{2m}\sum\limits_{j=1}^{n}{\theta _{j}^{2}} J(θ)=m1i=1∑m[−y(i)log(hθ(x(i)))−(1−y(i))log(1−hθ(x(i)))]+2mλj=1∑nθj2
Tip:可以通过np.matrix()函数将一个变量转换为numpy型矩阵
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)
h = sigmoid(theta @ X.T) # (1 x 5000)
m = len(y)
cross_cost = (np.log(h) @ (-y) - np.log(1 - h) @ (1 - y)) / m
cross_cost = cross_cost[0,0]
# STEP3:根据公式计算损失函数中的正则化部分
# your code here (appro ~ 1 lines)
theta_1 = theta[:, 1:theta.shape[1]]
reg = learningRate / (2 * m) * np.sum(np.power(theta_1, 2))
# STEP4:把上两步当中的结果加起来得到整体损失函数
# your code here (appro ~ 1 lines)
whole_cost = cross_cost + reg
return whole_cost
如果我们要使用梯度下降法令这个代价函数最小化,因为我们未对
θ
0
{{\theta }_{0}}
θ0 进行正则化,所以梯度下降算法将分两种情形:
R
e
p
e
a
t
u
n
t
i
l
c
o
n
v
e
r
g
e
n
c
e
{
θ
0
:
=
θ
0
−
a
1
m
∑
i
=
1
m
[
h
θ
(
x
(
i
)
)
−
y
(
i
)
]
x
0
(
i
)
θ
j
:
=
θ
j
−
a
1
m
∑
i
=
1
m
[
h
θ
(
x
(
i
)
)
−
y
(
i
)
]
x
j
(
i
)
+
λ
m
θ
j
}
R
e
p
e
a
t
Repeat\text{ }until\text{ }convergence\text{ }\!\!\{\!\!\text{ } \\ \text{ }{{\theta }_{0}}:={{\theta }_{0}}-a\frac{1}{m}\sum\limits_{i=1}^{m}{[{{h}_{\theta }}\left( {{x}^{(i)}} \right)-{{y}^{(i)}}]x_{_{0}}^{(i)}} \\ \text{ }{{\theta }_{j}}:={{\theta }_{j}}-a\frac{1}{m}\sum\limits_{i=1}^{m}{[{{h}_{\theta }}\left( {{x}^{(i)}} \right)-{{y}^{(i)}}]x_{j}^{(i)}}+\frac{\lambda }{m}{{\theta }_{j}} \\ \text{ }\!\!\}\!\!\text{ } \\ Repeat
Repeat until convergence { θ0:=θ0−am1i=1∑m[hθ(x(i))−y(i)]x0(i) θj:=θj−am1i=1∑m[hθ(x(i))−y(i)]xj(i)+mλθj } Repeat
梯度函数
向量化的梯度函数
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) # 5000 x 401
y = np.matrix(y)
# STEP2:将theta矩阵拉直(转换为一个向量)
# your code here (appro ~ 1 lines)
#parameters = int(theta.ravel().shape[1])
# STEP3:计算预测的误差
# your code here (appro ~ 1 lines)
error = sigmoid(X * theta.T) - y # 5000 x 1
# STEP4:根据上面的公式计算梯度
# your code here (appro ~ 1 lines)
h = sigmoid(X @ theta.T) # 5000 x 1
m = len(y)
grad = (X.T @ error).T / m + learningRate / m * theta
# STEP5:由于j=0时不需要正则化,所以这里重置一下
# your code here (appro ~ 1 lines)
grad[0, 0] = error.T @ X[:,0] / m
return np.array(grad).ravel()
多分类 one_vs_all
现在我们已经定义了代价函数和梯度函数,现在是构建分类器的时候了。 对于这个任务,我们有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))
# minimize the objective function
fmin = minimize(fun=cost, x0=theta, args=(X, y_i, learning_rate), method='TNC', jac=gradient)
all_theta[i-1,:] = fmin.x
return all_theta
这里需要注意的几点:首先,我们为theta添加了一个额外的参数(与训练数据一列),以计算截距项(常数项)。 其次,我们将y从类标签转换为每个分类器的二进制值(要么是类i,要么不是类i)。 最后,我们使用SciPy的较新优化API来最小化每个分类器的代价函数。 如果指定的话,API将采用目标函数,初始参数集,优化方法和jacobian(渐变)函数。 然后将优化程序找到的参数分配给参数数组。
实现向量化代码的一个更具挑战性的部分是正确地写入所有的矩阵,保证维度正确。
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, y_0.shape, theta.shape, all_theta.shape
((5000, 401), (5000, 1), (401,), (10, 401))
注意,theta是一维数组,因此当它被转换为计算梯度的代码中的矩阵时,它变为(1×401)矩阵。 我们还检查y中的类标签,以确保它们看起来像我们想象的一致。
np.unique(data['y'])#看下有几类标签
array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dtype=uint8)
让我们确保我们的训练函数正确运行,并且得到合理的输出。
all_theta = one_vs_all(data['X'], data['y'], 10, 1)
all_theta
array([[-2.38291641e+00, 0.00000000e+00, 0.00000000e+00, ...,
1.30399244e-03, -6.56319761e-10, 0.00000000e+00],
[-3.18263827e+00, 0.00000000e+00, 0.00000000e+00, ...,
4.46199339e-03, -5.08731233e-04, 0.00000000e+00],
[-4.79748863e+00, 0.00000000e+00, 0.00000000e+00, ...,
-2.86334542e-05, -2.46854791e-07, 0.00000000e+00],
...,
[-7.98858976e+00, 0.00000000e+00, 0.00000000e+00, ...,
-8.95444499e-05, 7.21593525e-06, 0.00000000e+00],
[-4.57292023e+00, 0.00000000e+00, 0.00000000e+00, ...,
-1.33526616e-03, 9.98755415e-05, 0.00000000e+00],
[-5.40304059e+00, 0.00000000e+00, 0.00000000e+00, ...,
-1.16560035e-04, 7.87934267e-06, 0.00000000e+00]])
预测
我们现在准备好最后一步 - 使用训练完毕的分类器预测每个图像的标签。 对于这一步,我们将计算每个类的类概率,对于每个训练样本(使用当然的向量化代码),并将输出类标签为具有最高概率的类。
Tip:可以使用np.argmax()函数找到矩阵中指定维度的最大值
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(all_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))
accuracy = 94.46%