吴恩达机器学习作业3：多类分类（Python实现）

机器学习练习 3 - 多类分类

在本练习中，您将实现一对一的逻辑回归和神经网络来识别手写的数字。在开始编程练习之前，我们强烈建议您观看视频讲座，并完成相关主题的复习问题。要开始这个练习，您需要下载启动代码并将其内容解压缩到您希望完成这个练习的目录中。自动手写数字识别在今天被广泛使用——从识别邮件信封上的邮政编码（邮政编码）到识别银行支票上所写的金额。本练习将展示您所学习到的方法如何用于此分类任务。在练习的第一部分中，您将扩展您以前的逻辑回归实现，并将其应用于一对一的分类。

1 Multi-class Classification

1.1 Dataset

在ex3data1.mat中有一个数据集，其中包含5000个手写数字的训练示例。mat格式意味着数据被保存为矩阵格式，而不是像csv文件那样的文本(ASCII)格式。这些矩阵可以通过使用loadmat命令直接读取到您的程序中。加载后，正确的维度和值的矩阵将出现在程序的内存中。该矩阵将已经被命名，因此您不需要为它们指定名称。

先导入需要的函数库：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.io import loadmat

将ex3data1.mat加载到data中，创建一个load_data函数，传入参数：path
这里的数据为MATLAB的格式，所以要使用SciPy.io的loadmat函数。

def load_data(path):
    data=loadmat(path)
    X=data["X"]
    y=data["y"]
    return X,y

传入参数，调用load_data函数

X, y = load_data('ex3data1.mat')

#查看y的标签数量
print(np.unique(y))

[ 1  2  3  4  5  6  7  8  9 10]

#查看X,y的形状
X.shape,y.shape

((5000, 400), (5000, 1))

可以看到一共有5000个训练样本，每个样本是$20\ast 20$像素的数字灰度图像，每个像素用一个浮点数表示该位置的灰度强度。20×20的像素网格被展开成一个400维的向量。而矩阵X：
$X=\left[\begin{array}{c}\cdot \cdot \cdot (X^{(1)}{)}^T\cdot \cdot \cdot \\ \cdot \cdot \cdot (X^{(2)}{)}^T\cdot \cdot \cdot \\ .\\ .\\ .\\ \cdot \cdot \cdot (X^{(m)}{)}^T\cdot \cdot \cdot \end{array}\right]$

1.2 Visualizing the data

第一个任务:将逻辑回归实现修改为完全向量化,简洁高效并且可以利用线性代数优化。

先编写绘制一个图的函数：plot_an_image,参数是矩阵X。

def plot_an_image(X):
    pick_one=np.random.randint(0,5000) 
    image=X[pick_one,:]
    #image (400:)
    fix,ax=plt.subplots(figsize=(1,1))
    ax.matshow(image.reshape(20,20),cmap='gray_r')#重塑为20*20的形状,并且设置为黑体白框(gray是白体黑框)
    #去掉x,y轴刻度
    plt.xticks([])
    plt.yticks([])
    plt.show()
    print("It's %d"%y[pick_one])

plot_an_image(X)

​

​

It's 9

编写绘制100个图的函数：plot_100_image,参数是矩阵X。

def plot_an_image(X):
    #从[0,5000)中随机选取100个
    sample_index=np.random.choice(np.arange(X.shape[0]),100)
    sample_imagex=X[sample_index,:]
    #image (100:400)
    fix,ax_array=plt.subplots(nrows=10,ncols=10,figsize=(8,8),sharey=True,sharex=True)#行列下标从0开始
    for row in range(10):
        for col in range(10):
            ax_array[row,col].matshow(sample_imagex[row*10+col].reshape(20,20),cmap='gray_r')#重塑第row*10+col个图形为20*20的形状,并且设置为黑体白框(gray是白体黑框)
    #去掉x,y轴刻度
    plt.xticks([])
    plt.yticks([])
    plt.show()

plot_an_image(X)

​

​

1.3 Vectorizing Logistic Regression

将使用多个一对多logistic回归模型来构建一个多类分类器。因为有10个类，所以您需要训练10个单独的logistic回归分类器。为了使这种训练高效，确保代码得到良好的向量化是很重要的。在本节中，您将实现一个逻辑回归的逻辑回归，它不对循环使用任何版本。您可以使用最后一个练习中的代码作为本练习的起点。

1.3.1 Vectorizing the cost function

我们将从编写代价函数的一个向量化版本开始。回想一下，在（非正则化的）逻辑回归中，代价函数为:

$J(\theta )=\frac{1}{m}{\sum }_{i=1}^m[-y^{(i)}log(h_{\theta }(x^{(i)}))-(1-y^{(i)})log(1-h_{\theta }(x^{(i)}))]+\frac{\lambda }{2m}{\sum }_{j=1}^n{\theta }_j^2$

为了计算总和中的每个元素，我们必须计算每个例子的$h_{\theta }(x^{(i)}),h_{\theta }(x^{(i)})=g({\theta }^T{x}^{(i)}),g(z)=\frac{1}{1+e^{-z}}$

事实上可以对所有的样本用矩阵乘法来快速的计算:

$X=\left[\begin{array}{c}\cdot \cdot \cdot (X^{(1)}{)}^T\cdot \cdot \cdot \\ \cdot \cdot \cdot (X^{(2)}{)}^T\cdot \cdot \cdot \\ .\\ .\\ .\\ \cdot \cdot \cdot (X^{(m)}{)}^T\cdot \cdot \cdot \end{array}\right]\theta =\left[\begin{array}{c}{\theta }_0\\ {\theta }_1\\ \cdot \\ \cdot \\ \cdot \\ {\theta }_n\end{array}\right]$

通过计算$X\theta$，可以得到：$X\theta =\left[\begin{array}{c}\cdot \cdot \cdot (X^{(1)}{)}^T\theta \cdot \cdot \cdot \\ \cdot \cdot \cdot (X^{(2)}{)}^T\theta \cdot \cdot \cdot \\ .\\ .\\ .\\ \cdot \cdot \cdot (X^{(m)}{)}^T\theta \cdot \cdot \cdot \end{array}\right]$=$\left[\begin{array}{c}\cdot \cdot \cdot {\theta }^T(X^{(1)})\cdot \cdot \cdot \\ \cdot \cdot \cdot {\theta }^T(X^{(2)})\cdot \cdot \cdot \\ .\\ .\\ .\\ \cdot \cdot \cdot {\theta }^T(X^{(m)})\cdot \cdot \cdot \end{array}\right]$

在最后一个等式中，我们使用了$a^{T}b=b^{T}a$，如果a和b是向量。这允许我们在一行代码中计算我们所有训练样本的${\theta }^{T}x$。

def sigmoid(z):
    return 1/(1+np.exp(-z))

def regularized_cost(Theta,X,y,l):
    #不需要惩罚Theta[0]
    ThetaReg=Theta[1:]
    cost=(-y*np.log(sigmoid(X@Theta)))-(1-y)*np.log(1-sigmoid((X@Theta)))
    reg=(ThetaReg@ThetaReg)*l/(2*len(X)) #@在矩阵中是矩阵相乘，在向量中是内积
    return np.mean(cost)+reg

1.3.2 Vectorizing the gradient

回想一下，（非正则化的）逻辑回归成本的梯度是一个向量，第$j$个元素被定义为$\frac{\partial J}{\partial {\theta }_j}=\frac{1}{m}{\sum }_{i=1}^m((h_{\theta }(x^{(i)})-y^i)x_j^{(i)})$

为了在数据集上向量化这个操作，我们首先显式地写出所有${\theta }_j$的偏导数，

其中：

注意$x^{(i)}$是向量，而$(h_{\theta }(x^{(i)})-y^{(i)})$是单一的数。

为了理解推导的最后一步，令${\beta }_i=(h_{\theta }(x^{(i)})-y^{(i)})$，并且注意：

上面的表达式允许我们计算所有的偏导数，而没有任何循环。如果你掌握了线性代数，鼓励你通过上面的矩阵乘法。

1.3.3 Vectorizing regularized logistic regression

在实现了逻辑回归的向量化后，现在将向成本函数中添加正则化。回想一下，对于正则化的逻辑回归，成本函数被定义为:

请注意，不应该规范化偏差项${\theta }_0$。

相应地，对${\theta }_j$的正则化逻辑回归代价的偏导数被定义为:

接下来实现梯度函数：

def regularized_gradient(Theta,X,y,l):
    ThetaReg=Theta[1:]
    cost=(X.T@(sigmoid(X@Theta)-y))*(1/len(X))
    reg=np.concatenate([np.array([0]),(l/len(X))*ThetaReg])
    return cost+reg

1.4 One-vs-all Classification

在本部分的练习中，将通过训练多个正则化logistic回归分类器实现一对多分类，每个对应数据集中K类中的一个。

对于这个任务，有10个可能的类，并且由于logistic回归只能一次在2个类之间进行分类，每个分类器在"类别i"和"不是i"之间决定。

将把分类器训练包含在一个函数中，该函数计算10个分类器中的每个分类器的最终权重，并将权重返回$shape$为$(k,(n+1))$数组，其中$n$是参数数量。

首先，为X添加了一列常数项X0=1。其次，将y从类标签转换为每个分类器的二进制值（要么是类i，要么不是类i）。最后，使用SciPy的较新优化API来最小化每个分类器的代价函数。（API将采用目标函数、初始参数集、优化方法和jacobian（渐变）函数。 将优化程序找到的参数分配给参数数组。）

from scipy.optimize import minimize

def one_vs_all(X,y,l,K):
    '''
    X:feature matrix, (m, n+1)
    y:target vector, (m, )
    l: lambda constant for regularization
    K: numbel of labels
    '''
    All_Theta=np.zeros((K,X.shape[1]))#(10,401)
    for i in range(1,K+1):
        Theta=np.zeros(X.shape[1])
        y_i=np.array([1 if labal==i else 0 for labal 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
    return All_Theta

1.4.1 One-vs-all Prediction

在训练了一对多的分类器之后，现在使用它来预测给定图像中包含的数字。对于每个输入，应该使用训练过的逻辑回归分类器来计算它属于每个类的"概率"。一对多预测函数将选择相应的逻辑回归分类器输出最高概率的类，并返回类标签$(1,2,...,或K)$作为输入示例的预测。
你应该看到训练集的准确率约为94.9%（也就是说，它正确地分类了训练集中94.9%的例子）。

$h=(5000,10)$，每行代表一个样本，每列是预测对应数字的概率。取概率最大对应的index加1就是我们分类器最终预测出来的类别。$h\_argmax$是一个$array$，包含$5000$个样本对应的预测值。

def predict(X,All_Theta):
    h=sigmoid(X@All_Theta.T)
    h_argmax=np.argmax(h,axis=1)
    h_argmax=h_argmax+1
    return h_argmax

X=np.insert(X,0,1,axis=1)
y=y.flatten()
print(X,y)
All_Theta=one_vs_all(X, y, 1, 10)
All_Theta

[[1. 0. 0. ... 0. 0. 0.]
 [1. 0. 0. ... 0. 0. 0.]
 [1. 0. 0. ... 0. 0. 0.]
 ...
 [1. 0. 0. ... 0. 0. 0.]
 [1. 0. 0. ... 0. 0. 0.]
 [1. 0. 0. ... 0. 0. 0.]] [10 10 10 ...  9  9  9]





array([[-2.38326459e+00,  0.00000000e+00,  0.00000000e+00, ...,
         1.30473915e-03, -8.24448604e-10,  0.00000000e+00],
       [-3.18289654e+00,  0.00000000e+00,  0.00000000e+00, ...,
         4.45802725e-03, -5.08266637e-04,  0.00000000e+00],
       [-4.79725924e+00,  0.00000000e+00,  0.00000000e+00, ...,
        -2.86397823e-05, -2.46983419e-07,  0.00000000e+00],
       ...,
       [-7.98538798e+00,  0.00000000e+00,  0.00000000e+00, ...,
        -8.95502218e-05,  7.22536632e-06,  0.00000000e+00],
       [-4.57416889e+00,  0.00000000e+00,  0.00000000e+00, ...,
        -1.33600190e-03,  9.99444523e-05,  0.00000000e+00],
       [-5.40408300e+00,  0.00000000e+00,  0.00000000e+00, ...,
        -1.16720634e-04,  7.88817209e-06,  0.00000000e+00]])

y_predict=predict(X, All_Theta)
accuracy=np.mean(y_predict==y)
print("accuracy=%.2f%%"%(accuracy*100))
# print("accuracy={:.2f}%".format(accuracy*100))

accuracy=94.46%

2 Neural Networks

在该部分的前一部分中，实现了多类逻辑回归来识别手写数字。然而，逻辑回归不能形成更复杂的假设，因为它只是一个线性分类器。

在该部分练习中，将使用相同的训练集实现神经网络来识别手写数字。神经网络将能够代表形成非线性假设的复杂模型。该部分练习将使用已经训练过的神经网络中的参数。利用已经训练好了的权重进行预测。

2.1 Model representation

神经网络如图所示:.它有3层:一个输入层，一个隐藏层和一个输出层。其中输入是数字图像的像素值。

def load_weight(path):
    data=loadmat(path)
    return data['Theta1'],data['Theta2']

Theta1,Theta2=load_weight("ex3weights.mat")
Theta1.shape,Theta2.shape

((25, 401), (10, 26))

可以发现数据加载函数中：原始数据做了转置。但是，转置的数据与给定的参数不兼容，因为这些参数是由原始数据训练的。所以为了应用给定的参数，需要使用原始数据（不进行转置）。

X, y = load_data('ex3data1.mat')
X=np.insert(X,0,1,axis=1)
y=y.flatten()
X.shape,y.shape

((5000, 401), (5000,))

2.2 Feedforward Propagation and Prediction

现在将为神经网络实现前馈传播。为每一个样例i实现前馈计算$h_{\theta }(x^{(i)})$并且返回相关联的预测。类似于一对多的分类策略，来自神经网络的预测将是可能性最大的输出$(h_{\theta }(x^{(i)}){)}_k$

a1=X
a1.shape,Theta1.shape

((5000, 401), (25, 401))

z2=a1@Theta1.T
z2.shape

(5000, 25)

z2=np.insert(z2,0,1,axis=1)

a2=sigmoid(z2)
# print(a2)
a2.shape

(5000, 26)

z3=a2@Theta2.T
z3.shape

(5000, 10)

a3=sigmoid(z3)
a3

array([[1.38245045e-04, 2.05540079e-03, 3.04012453e-03, ...,
        4.91017499e-04, 7.74325818e-03, 9.96229459e-01],
       [5.87756717e-04, 2.85026516e-03, 4.14687943e-03, ...,
        2.92311247e-03, 2.35616705e-03, 9.96196668e-01],
       [1.08683616e-04, 3.82659802e-03, 3.05855129e-02, ...,
        7.51453949e-02, 6.57039547e-03, 9.35862781e-01],
       ...,
       [6.27824726e-02, 4.50406476e-03, 3.54510925e-02, ...,
        2.63669734e-03, 6.89448164e-01, 2.74369466e-05],
       [1.01908736e-03, 7.34360211e-04, 3.78558700e-04, ...,
        1.45616578e-02, 9.75989758e-01, 2.33374461e-04],
       [5.90807037e-05, 5.41717668e-04, 2.58968308e-05, ...,
        7.00508308e-03, 7.32814653e-01, 9.16696059e-02]])

y_predict=np.argmax(a3,axis=1)+1#返回每行的最大值的索引
y_predict.shape

(5000,)

准确度的计算：
虽然人工神经网络是非常强大的模型，但训练数据的准确性并不能完美预测实际数据，在这里很容易过拟合。

accuracy=np.mean(y_predict==y)
print("accuracy=%.2f%%"%(accuracy*100))

accuracy=97.52%

from sklearn.metrics import classification_report#这个包是评价报告
print(classification_report(y, y_predict))

              precision    recall  f1-score   support

           1       0.97      0.98      0.97       500
           2       0.98      0.97      0.97       500
           3       0.98      0.96      0.97       500
           4       0.97      0.97      0.97       500
           5       0.98      0.98      0.98       500
           6       0.97      0.99      0.98       500
           7       0.98      0.97      0.97       500
           8       0.98      0.98      0.98       500
           9       0.97      0.96      0.96       500
          10       0.98      0.99      0.99       500

    accuracy                           0.98      5000
   macro avg       0.98      0.98      0.98      5000
weighted avg       0.98      0.98      0.98      5000

​