编者注:本文用logistic回归来识别多分类问题的手写数字,是之前logisitic回归二分类问题的延续,该篇文章关于其思想以及编程原理见本人之前文章,在这里只注重识别及其编程过程。
01数据准备
数据为Matlab加载格式(mat),包含y、X变量,数据来源为(大家可以去这下载):
https://github.com/jdwittenauer/ipython-notebooks/blob/master/data/ex3data1.mat ,至于为什么要用Matlab格式,因为X是个矩阵不方便用excel之类软件来保存。
image.png其中,X为5000X400维矩阵,5000代表5000个训练样本,400则代表每个训练样本(也就是手写数字图像)20 X20像素的灰度强度,随机选择其中100个样本,在Matlab可视化的结果如下:
image.pngy 是一个5000行1列的列向量,取值包括(1,2,3,4,5,6,7,8,9,10)T,注意,由于Matlab下标是从1开始的,故用 10 表示数字 0。在python中使用scipy导入mat数据,代码如下:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.io import loadmat
data = loadmat('E:\studypythonml\ex3data1.mat')
data
data['X'].shape, data['y'].shape
02 构建相关函数
在python中分别构建sigmoid、梯度下降、成本函数等。
def sigmoid(z):
return 1 / (1 + np.exp(-z))
def cost(theta, X, y, learningRate):
theta = np.matrix(theta)
X = np.matrix(X)
y = np.matrix(y)
first = np.multiply(-y, np.log(sigmoid(X * theta.T)))
second = np.multiply((1 - y), np.log(1 - sigmoid(X * theta.T)))
reg = (learningRate / 2 * len(X)) * np.sum(np.power(theta[:,1:theta.shape[1]], 2))
return np.sum(first - second) / (len(X)) + reg
def gradient(theta, X, y, learningRate):
theta = np.matrix(theta)
X = np.matrix(X)
y = np.matrix(y)
parameters = int(theta.ravel().shape[1])
error = sigmoid(X * theta.T) - y
grad = ((X.T * error) / len(X)).T + ((learningRate / len(X)) * theta)
# intercept gradient is not regularized
grad[0, 0] = np.sum(np.multiply(error, X[:,0])) / len(X)
return np.array(grad).ravel()
这其中的编程思想在之前写的相关随笔:梯度下降法求解线性回归、logisitic回归等出现过,大家可以参考之前文章来了解构建函数的过程。特别注意的是,从最开始的求解一元线性回归时,都注重了将X变量的向量化,对于这里更复杂的X变量(矩阵形式)起到了作用,避免采用更复杂的code,如for循环来引入。
03 构建one vs all函数
对于logisitic多分类问题,在stanford Coursera公开课中有过专门介绍,所谓多分类问题,是指分类的结果为三类以上。其分类思想与二分类类似,对于k个类别,在预测某个类i的时候,将其作为一类,剩下的作为另一类,这样的话,就会有k个logisitic分类器,如下图:
image.png
就需要N个假设函数(预测模型),也即需要N组模型参数θ,对于每一个类 i 训练一个逻辑回归模型的分类器h(i)θ(x),并且预测 y = i时的概率;对于一个新的输入变量x, 分别对每一个类进行预测,取概率最大的那个类作为分类结果 。因此,在这里构建一个函数将分类器训练全部包含进去,这个函数计算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] #X的行数
params = X.shape[1] #X的列数
# 对于k个分类器,构建k*(n+1)维向量组
all_theta = np.zeros((num_labels, params + 1))
# 在X第一列之前插入一列全为1的列向量作为常数项
X = np.insert(X, 0, values=np.ones(rows), axis=1)
# 对于y若将某个类别i拿出来之后剩下的类别构成一类
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))
# 采用梯度下降法最小化目标函数(cost)
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
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
np.unique(data['y'])
all_theta = one_vs_all(data['X'], data['y'], 10, 1)
all_theta
得到10个分类器的最优参数为:
image.png
04 得到预测值并与原值比较计算准确率
根据训练生成的all_theta参数构建预测函数,得到X对应的预测值
def predict_all(X, all_theta):
rows = X.shape[0]
params = X.shape[1]
num_labels = all_theta.shape[0]
# 与前相同,插入一列全部为1的列向量
X = np.insert(X, 0, values=np.ones(rows), axis=1)
# 转换为矩阵
X = np.matrix(X)
all_theta = np.matrix(all_theta)
# 计算每个训练实例上每个类的类概率
h = sigmoid(X * all_theta.T)
# 选取最高的那个概率为该实例的预测数字标签并构建数组
h_argmax = np.argmax(h, axis=1)
# 由于该数组在训练时是基于图片的0-9而预测的,所以要+1以匹配y
h_argmax = h_argmax + 1
return h_argmax
代入,得到预测精度
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)
首先构建预测值函数,然后将该值与原始类别值0,1比较,计算其正确的精度,结果为97.46%。
写作不易,特别是技术类的写作,请大家多多支持,关注、点赞、转发等等,也欢迎大家关注.