机器学习（吴恩达）--手写数字识别（逻辑回归实现分类）

学习内容：逻辑回归练习

学习时间：2020/10/7

教学视频、数据集等：

https://www.bilibili.com/video/BV1mt411p7kG?p=1

代码实现：

https://github.com/boxes757/Machine-learning-exercises/blob/master/机器学习2——逻辑回归练习（分类）

学习笔记：

练习题目

1. 手写数字识别：根据图片识别0~9这10个数字，其中数字0的标签值为10

相关工具（包）

numpy（计算，处理多维数组）、matplotlib.pyplot（绘图框架）、scipy.io（读取.mat文件）、scipy.optimize.minimize（一个或多个变量的函数的最小化）

参数

1. sigmoid函数
            $g(z)=\frac{1}{1+e^{-z}}$

   对于$h_{\theta }(x)=g({\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2+...+{\theta }_n{x}_n)$来说，我们计算出的$h_{\theta }(x)$都在(0,1)这个区间内。
   在二分类中，我们将$h_{\theta }(x)$输出的结果理解为对于特征集x来说，y=1的概率。比如$h_{\theta }(x)=0.7$，我们理解为$y=1$的概率为0.7。
   在多分类中，对于第 $i$ 分类器来说，都有一组$[{\theta }_0^{(i)},{\theta }_1^{(i)},...,{\theta }_n^{(i)}]$，我们通过这组值计算出相应的$h_{\theta }^{(i)}(x)$。对于一个样本来说，比较其在各分类器中的$h_{\theta }^{(i)}(x)$值，选取最大的一个，则代表它属于这个分类的概率值最大。

2. 损失函数（带正则项）
   损失函数与之前类似，化简后的损失函数为：
   

3. 梯度下降函数（带正则项）
   （$\theta =\theta -\frac{\partial J(\theta )}{\partial {\theta }_j}$，注意 $j=0$ 时与之前的梯度下降函数一致）
   

读取数据

读取数据我们用到了scipy.io这个包，我们调用其里面的loadmat来读取.mat文件

import scipy.io as sio
data = sio.loadmat('03-neural network/ex3data1.mat')

预览数据集（matplotlib.pyplot）

• numpy.random.choice(a, size=None, replace=True, p=None)
  a : 如果是一维数组，就表示从这个一维数组中随机采样； 如果是int型，就表示从0到a-1这个序列中随机采样。
  size：取样的数量
  replace:True时允许重复；反之不允许
  p : 一维数组，代表a中每个元素采样的概率，默认时a中每个元素被采样的概率相同。

def print_100_picture(raw_X):   ##随机画一百张图
    sample_index = np.random.choice(5000,100) ##
    images = raw_X[sample_index,:]

    fig,ax = plt.subplots(10,10,figsize=(10,10),sharex=True,sharey=True)  ##产生10*10个子图，整个图的大小是10*10,共享x,y的坐标

    for i in range(0,10):
        for j in range(0,10):
            R = images[10*i + j].reshape(20,20) ##第3行第4列是ax[2,3]，对应第24张图，也就是10*2+3=23行的数据
            ax[i,j].imshow(R.T,cmap='gray_r') ##reshape的作用是把该行一维数组变成（20*20）的矩阵，才能形成可视化图片，转置是因为图片视角不对，gray_r 会将其显示为白底黑字
    plt.xticks([])  ##横纵坐标赋空列表，则不显示数值
    plt.yticks([])
    plt.show()

print_100_picture(raw_X)

优化损失函数（scipy.optimize.minimize）

该方法的具体文档可以参考以下链接：
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html (官方英文文档)
https://blog.csdn.net/cazenove/article/details/107820322（用例）
https://www.cnblogs.com/onenoteone/p/12441668.html（用例）
注：

• theta为一维向量，在之前例子中为（n,1）二维向量，但用minimize方法时，x0（初始预测值）需为一维向量！！

• 之前例子中我们使用迭代方法手动进行梯度下降，此处我们只需要计算梯度向量返回给minimize方法中的jac参数即可，我们也不需要手动进行梯度下降时设置的 $\alpha$ 值。

• 带正则化项的梯度向量如下所示：
  

• 一维数组与矩阵的乘法
  左乘作为行向量运算，右乘作为列向量进行运算，结果得到一维数组
  如(2,) $\ast$ (2,4) = (4,)    (3,5) $\ast$ (5,) = (3,)

• np.zeros
  生成由n个0元素组成的一维数组：np.zeros(5)或np.zeros(5,)
  生成（2,3）的由0组成的二维数组：np.zeros((2,3))    注意两层括号

def CostFunction(theta,X,y,lamda):  ##theta为一维向量！！
    R=sigmoid(X@theta)   ##theta一维向量，作矩阵乘法时。左乘作为行向量运算，右乘作为列向量进行运算，计算结果也是一维数组
    A=y * np.log(R)  ##X为（5000,401）矩阵，theta为（401,）一维向量，R是(5000,)一维向量
    B=(1-y) * np.log(1-R)
    ##reg = np.sum(np.power(theta[1:],2)) * (lamda/ （2 * len(X)）)
    reg = theta[1:] @ theta[1:] * (lamda/ (2 * len(X)))  ##两种写法结果一致
    return -np.sum(A+B)/len(X)+reg

def gradient_reg(theta,X,y,lamda):   ##带正则化项的梯度向量
    A = sigmoid(X @ theta)
    R = theta[1:] * (lamda / len(X))
    R = np.insert(R, 0, values=0, axis=0)
    grad1 = (X.T @ (A-y)) / len(X)
    return grad1+R

X = np.insert(raw_X,0,values=1,axis=1)
print(X.shape)

y = raw_y.flatten()  ##损失函数中涉及到y与R的运算，R为（5000，）的一维向量，需要将y也转化为一维向量
print(y.shape)

from scipy.optimize import minimize

def one_vs_all(X,y,lamda,K):  ##K是分类器值
    n = X.shape[1]  ##获取X的第二个维度的值，此处即列数401
    theta_all = np.zeros((K,n))   ##此处为（10,401）矩阵

    for i in range(1,K+1):    ##此处K=10，i取1,2,3,...,10
        theta_i = np.zeros(n,)   ##以0作为初始值，分别为各分类器运行优化函数
        res = minimize(CostFunction,x0=theta_i,args=(X,y == i,lamda),method='TNC',jac=gradient_reg)
        ##args是将额外参数传入目标函数（该处为损失函数），y==i是将该分类器值代入
        ##jac也可以是一个可调用的，返回目标的梯度的函数。在这种情况下，它必须接受与目标函数（该处为损失函数）相同的参数。
        theta_all[i-1,:] = res.x  ##将结果保存到theta_all的相应位置，注意i为分类器的标号，存储时需要-1
    return theta_all

预测及准确性

• numpy的argmax（可参考：https://www.cnblogs.com/touch-skyer/p/8509217.html）

1. 一维数组时，返回的是索引值的值，比如第4个数，返回3；
2. 二维数组时，返回的是由各索引值组成的一维数组
3. axis=0代表在各列中寻找该列某一行的最大值的索引值（寻找各列最大值所在行的索引值）
4. axis=1代表在各行中寻找该行某一列的最大值的索引值（寻找各行最大值所在列的索引值）

def predict(X,theta_final):
    A = sigmoid(X @ theta_final.T) ##X是（5000,401），theta_final是（10,401），A为（5000,10）
    k = np.argmax(A,axis=1)  ##np.argmax返回的是5000行中每一行中最大数所在列的索引值（即该类别概率最大）。
    return k+1               
 ## k+1：当索引值为4时，代表第5个分类器，所以要加1     0是第一列，代表分类器1；1代表分类器2；9代表分类器10，分类器10规定为数字0；
y_predict = predict(X,theta_final)
acc = np.mean(y_predict == y)  ##由于predict方法返回的是分类器的值（索引值加1）组成的一维数组，所以和y比，y也是转化来的一维数组
print(acc)

theta_final及准确值（acc）：

结果展示