编者注：本文用logistic回归来识别多分类问题的手写数字，是之前logisitic回归二分类问题的延续，该篇文章关于其思想以及编程原理见本人之前文章，在这里只注重识别及其编程过程。
01数据准备
数据为Matlab加载格式(mat)，包含y、X变量，数据来源为(大家可以去这下载)：
/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%。
写作不易，特别是技术类的写作，请大家多多支持，关注、点赞、转发等等，也欢迎大家关注.
