import numpy as npimport pandas as pdimport matplotlib.pyplot as pltfrom scipy.io import loadmat
data = loadmat(\'ex3data1.mat\')data
data[\'X\'].shape, data[\'y\'].shape
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
#以下是原始代码是使用for循环的梯度函数：def gradient_with_loop(theta, X, y, learningRate):theta = np.matrix(theta)X = np.matrix(X)y = np.matrix(y)parameters = int(theta.ravel().shape[1])grad = np.zeros(parameters)error = sigmoid(X * theta.T) - yfor i in range(parameters):term = np.multiply(error, X[:,i])if (i == 0):grad[i] = np.sum(term) / len(X)else:grad[i] = (np.sum(term) / len(X)) + ((learningRate / len(X)) * theta[:,i])return grad
#向量化的梯度函数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) - ygrad = ((X.T * error) / len(X)).T + ((learningRate / len(X)) * theta)# intercept gradient is not regularizedgrad[0, 0] = np.sum(np.multiply(error, X[:,0])) / len(X)return np.array(grad).ravel()
#现在是构建分类器的时候了。 对于这个任务，我们有10个可能的类，并且由于逻辑回归只能一次在2个类之间进行分类，我们需要多类分类的策略。 在本练习中，我们的任务是实现一对一全分类方法，其中具有k个不同类的标签就有k个分类器，每个分类器在“类别 i”和“不是 i”之间决定。 我们将把分类器训练包含在一个函数中，该函数计算10个分类器中的每个分类器的最终权重，并将权重返回为k X（n + 1）数组，其中n是参数数量。from scipy.optimize import minimizedef one_vs_all(X, y, num_labels, learning_rate):rows = X.shape[0]params = X.shape[1]# k X (n + 1) array for the parameters of each of the k classifiersall_theta = np.zeros((num_labels, params + 1))# insert a column of ones at the beginning for the intercept termX = np.insert(X, 0, values=np.ones(rows), axis=1)# labels are 1-indexed instead of 0-indexedfor 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))# minimize the objective functionfmin = minimize(fun=cost, x0=theta, args=(X, y_i, learning_rate), method=\'TNC\', jac=gradient)all_theta[i-1,:] = fmin.xreturn all_theta
#这里需要注意的几点：首先，我们为theta添加了一个额外的参数（与训练数据一列），以计算截距项（常数项）。 其次，我们将y从类标签转换为每个分类器的二进制值（要么是类i，要么不是类i）。 最后，我们使用SciPy的较新优化API来最小化每个分类器的代价函数。 如果指定的话，API将采用目标函数，初始参数集，优化方法和jacobian（渐变）函数。 然后将优化程序找到的参数分配给参数数组。#实现向量化代码的一个更具挑战性的部分是正确地写入所有的矩阵，保证维度正确。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
#注意，theta是一维数组，因此当它被转换为计算梯度的代码中的矩阵时，它变为（1×401）矩阵。 我们还检查y中的类标签，以确保它们看起来像我们想象的一致。np.unique(data[\'y\'])#看下有几类标签
#让我们确保我们的训练函数正确运行，并且得到合理的输出all_theta = one_vs_all(data[\'X\'], data[\'y\'], 10, 1)all_theta
#我们现在准备好最后一步 - 使用训练完毕的分类器预测每个图像的标签。 对于这一步，我们将计算每个类的类概率，对于每个训练样本（使用当然的向量化代码），并将输出类标签为具有最高概率的类。def predict_all(X, all_theta):rows = X.shape[0]params = X.shape[1]num_labels = all_theta.shape[0]# same as before, insert ones to match the shapeX = np.insert(X, 0, values=np.ones(rows), axis=1)# convert to matricesX = np.matrix(X)all_theta = np.matrix(all_theta)# compute the class probability for each class on each training instanceh = sigmoid(X * all_theta.T)# create array of the index with the maximum probabilityh_argmax = np.argmax(h, axis=1)# because our array was zero-indexed we need to add one for the true label predictionh_argmax = h_argmax + 1return h_argmax
#现在我们可以使用predict_all函数为每个实例生成类预测，看看我们的分类器是如何工作的。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))
# 神经网络模型图示
爱站程序员基地
