【机器学习】吴恩达机器学习视频作业-逻辑回归多分类

阅读量：2019 次

发布时间：2019-04-28

本文共 5509 字，大约阅读时间需要 18 分钟。

使用逻辑回归进行二分类案例：

使用逻辑回归进行多分类

对于此练习，将使用逻辑回归算法手写数字（0到9）识别。扩展上个练习中逻辑回归的实现，并将其应用于一对多的分类。源数据是在MATLAB的数据格式，需要使用一个SciPy读取。

主要思路：one vs all ，假如需要去将数据样本分为5类，就会创建5个逻辑回归模型。每个逻辑回归模型中，一种分类设置标签为1，即为正样本，其它类别标签设置为0，即为负样本，依此类推构建5个分类模型。预测时，将一条数据在5个模型中运行，得到在各个模型中分类的概率结果，再判断哪个数值最大，那么就认为是哪一类。

1 加载数据与数据可视化

import numpy as npimport pandas as pdimport matplotlib.pyplot as pltfrom scipy.io import loadmat   # 加载matlab数据%matplotlib inline

加载数据，数据内容如下。

data = loadmat('ex3data1.mat')   # 加载数据data

实验数据

查看样本数据，一共5000个样本，每个样本为400列，也就是说将一张图片（20*20）被压平了。

data['X'].shape, data['y'].shape"""((5000, 400), (5000, 1))"""

查看分类标签有哪些值，其中数值0使用标签10代替。

np.unique(np.array(data['y']))   # 查看标签有哪些值  """array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10], dtype=uint8)"""

可视化数据，数据分布：0-499为数字0，标签为10，500-999为数字1，标签为1，以此类推。

plt.imshow(data['X'][0].reshape(20, 20))   # 查看数据图形data['y'][0]                               # 对应标签"""array([10], dtype=uint8)"""

数字0

plt.imshow(data['X'][500].reshape(20, 20))data['y'][500]"""array([1], dtype=uint8)"""

数字1

plt.imshow(data['X'][1000].reshape(20, 20))data['y'][1000]"""array([2], dtype=uint8)"""

数字2

2 构建sigmoid函数以及线性回归函数

模型函数创建，特征对应的线性函数如下：

Expected node of symbol group type, but got node of type ordgroup

逻辑回归模型函数如下：

h_\theta(x) = g(\theta^T x)

其中g(z)表达式如下：

\frac{1}{1+e^{-z}}

模型使用，我们预测：

当 $h_\theta(x) \ge 0.5$ 时，预测 $y = 1$

当 $h_\theta(x) < 0.5$ 时，预测 $y = 0$

在多分类中，选择模型模型结果最大的那个类别。

# 定义一个sigmoid函数def sigmoid(z):    # z可以是一个数，也可以是np中的矩阵    return 1.0/(1+np.exp(-z))def model(X, theta):    # 单个数据样本维一个行向量 n个特征， theta为长度为n的行向量    # 当X为m行n列的数据，返回结果为m行1列的数据    theta = theta.reshape(1, -1)  # 转成符合矩阵计算    return sigmoid(X.dot(theta.T))

3 正则化损失函数(Regularize cost)

$J(\theta)=\frac{1}{m} \sum_{i=1}^{m}\left[-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)\right]+\frac{\lambda}{2 m} \sum_{j=1}^{n} \theta_{j}^{2}$

特点：比普通的损失函数多了一项

\frac{\lambda}{2 m} \sum_{j=1}^{n} \theta_{j}^{2}

，注：此时j从1开始，在之前的损失函数基础上进行改进即可

# 损失函数def cost(theta, X, y, l=1):    # X为特征数据集 单个特征数据为行向量形状要为mxn    # y为标签数据集 为一个列向量形状要为mx1    # theta为一个行向量 形状要为1xn    # l 正则化参数，这里默认设置为1    # 返回一个标量    left = -y*np.log(model(X, theta))    right = (1-y)*np.log(1 - model(X, theta))    return np.sum(left - right)/len(X) + l*np.power(theta[1:], 2).sum()/(2*X.shape[0])

查看单个模型训练数据状况。

# 初步查看数据X = np.array(data['X'])  X = np.insert(X, 0, values=1, axis=1)    # 插入常数项y = np.array(data['y']).reshape(-1, 1)y_0 = np.array([1 if label == 0 else 0 for label in y]) # 数据样本中类为1设置为1，其它设置为0class_num = 0                            # 获取分类个数all_theta = np.zeros((10, X.shape[1]))   # 生成10*(params + 1)维得参数，因为10个分类，每一个分类对应一组参数thetaX.shape, y.shape, all_theta.shape"""((5000, 401), (5000, 1), (10, 401))"""

模型1初始损失函数

# 模型1的初始损失函数cost(all_theta[0], X, y_0)"""3465.7359027995562"""

4 梯度函数

正则化梯度公式分为两部分：

\frac{\partial J(\theta)}{\partial \theta_{0}}=\frac{1}{m} \sum_{i=1}^{m}\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{j}^{(i)} \quad \text { for } j = 1

$\frac{\partial J(\theta)}{\partial \theta_{j}}=\left(\frac{1}{m} \sum_{i=1}^{m}\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{j}^{(i)}\right)+\frac{\lambda}{m} \theta_{j} \quad \text { for } j \geq 1$

# 参数的梯度def gradient(theta, X, y, l=1):    # X为特征数据集 mxn    # y为标签数据集 mx1    # l 正则化参数，默认为1    # theta为一个行向量 1xn    grad = ((model(X, theta) - y).T@X).flatten()/ len(X)     grad[1:] = grad[1:] + l*theta[1:]/X.shape[0]   # 对非常数项参数加上正则项    # 返回行向量    return grad

查看模型1初始梯度值

gradient(all_theta[0], X, y_0[:, np.newaxis])[:5]"""array([ 4.00000000e-01,  0.00000000e+00,  0.00000000e+00, -7.74530186e-08,        4.91729643e-07])"""

7 使用优化函数计算参数

使用scipy的optimize模块。

import scipy.optimize as optdef one_vs_all(X, y, kinds, l):    # X为特征数据集 mxn(尚未插入常数项1)    # y为标签数据集 mx1    # l 正则化参数，默认为1    # theta为一个行向量 1xn    features = X.shape[1]  # 获取特征个数    all_theta = np.zeros((kinds, features))  # 初始化10个模型的参数    for i in range(1, kinds + 1):        theta = np.zeros(features)  # 当前模型参数值        # one vs all 标签转换        y_i = np.array([ 1 if label == i else 0 for label in y])[:, np.newaxis] # 转成列向量        # 使用TNC 截断牛顿法最小化损失值        fmin = opt.minimize(fun=cost, x0=theta, args=(X, y_i, l), method='TNC', jac=gradient)        all_theta[i-1,:] = fmin.x  # 获取计算的参数值    return all_theta

开始训练，获取各模型的参数。

all_theta = one_vs_all(X, y, 10, 1)all_theta

模型参数

8 使用训练的分类器预测数据

创建获取预测值函数。

def predict_all(X, all_theta):    # X 测试数据    # all_theta 模型参数    # 各模型计算结果后为一个5000x10的矩阵，从每行中选取最大值对应的索引即可    # 由于默认索引从0开始以及实际的标签是从1开始，需要给对应索引值+1    pre = np.argmax(sigmoid(X@all_theta.T), axis=1) + 1     return pre

查看模型预测的准确率。

y_pre = predict_all(X, all_theta)# 比较各数据预测情况correct = [1 if a==b else 0 for (a,b) in zip(y, y_pre)] # 计算精度accuracy = (sum(map(int, correct))/len(X))*100print("accuracy:{}%".format(accuracy))"""accuracy:94.46%"""