|
| 1 | +import numpy as np |
| 2 | +import pandas as pd |
| 3 | +import matplotlib.pyplot as plt |
| 4 | +from scipy.io import loadmat |
| 5 | +from scipy.optimize import minimize |
| 6 | + |
| 7 | + |
| 8 | +def load_data(path): |
| 9 | + data = loadmat(path) |
| 10 | + x = data['X'] |
| 11 | + y = data['y'] |
| 12 | + return x, y |
| 13 | + |
| 14 | + |
| 15 | +def plot_an_image(X, y): |
| 16 | + pick_one = np.random.randint(0, 5000) |
| 17 | + image = X[pick_one, :] |
| 18 | + plt.matshow(image.reshape((20, 20))) |
| 19 | + plt.show() |
| 20 | + print('this should be {}'.format(y[pick_one])) |
| 21 | + |
| 22 | + |
| 23 | +def plot_100_image(X): |
| 24 | + """ |
| 25 | + 随机画100个数字 |
| 26 | + """ |
| 27 | + sample_idx = np.random.choice(np.arange(X.shape[0]), 100) # 随机选100个样本 |
| 28 | + sample_images = X[sample_idx, :] # (100,400) |
| 29 | + |
| 30 | + fig, ax_array = plt.subplots(nrows=10, ncols=10, sharex=True, sharey=True, figsize=(8, 8)) |
| 31 | + |
| 32 | + for row in range(10): |
| 33 | + for column in range(10): |
| 34 | + ax_array[row, column].matshow(sample_images[10 * row + column].reshape((20, 20))) |
| 35 | + plt.xticks([]) |
| 36 | + plt.yticks([]) |
| 37 | + plt.show() |
| 38 | + |
| 39 | + |
| 40 | +def sigmoid(z): |
| 41 | + return 1 / (1 + np.exp(-z)) |
| 42 | + |
| 43 | + |
| 44 | +def regularized_cost(theta, X, y, l): |
| 45 | + """ |
| 46 | + don't penalize theta_0 |
| 47 | + args: |
| 48 | + X: feature matrix, (m, n+1) # 插入了x0=1 |
| 49 | + y: target vector, (m, ) |
| 50 | + l: lambda constant for regularization |
| 51 | + """ |
| 52 | + thetaReg = theta[1:] |
| 53 | + first = (-y * np.log(sigmoid(X @ theta))) + (y - 1) * np.log(1 - sigmoid(X @ theta)) |
| 54 | + reg = (thetaReg @ thetaReg) * l / (2 * len(X)) |
| 55 | + return np.mean(first) + reg |
| 56 | + |
| 57 | + |
| 58 | +def regularized_gradient(theta, X, y, l): |
| 59 | + """ |
| 60 | + don't penalize theta_0 |
| 61 | + args: |
| 62 | + l: lambda constant |
| 63 | + return: |
| 64 | + a vector of gradient |
| 65 | + """ |
| 66 | + thetaReg = theta[1:] |
| 67 | + first = (1 / len(X)) * X.T @ (sigmoid(X @ theta) - y) |
| 68 | + # 这里人为插入一维0,使得对theta_0不惩罚,方便计算 |
| 69 | + reg = np.concatenate([np.array([0]), (l / len(X)) * thetaReg]) |
| 70 | + return first + reg |
| 71 | + |
| 72 | + |
| 73 | +def one_vs_all(X, y, l, K): |
| 74 | + """generalized logistic regression |
| 75 | + args: |
| 76 | + X: feature matrix, (m, n+1) # with incercept x0=1 |
| 77 | + y: target vector, (m, ) |
| 78 | + l: lambda constant for regularization |
| 79 | + K: numbel of labels |
| 80 | + return: trained parameters |
| 81 | + """ |
| 82 | + all_theta = np.zeros((K, X.shape[1])) # (10, 401) |
| 83 | + |
| 84 | + for i in range(1, K + 1): |
| 85 | + theta = np.zeros(X.shape[1]) # (401,) |
| 86 | + y_i = np.array([1 if label == i else 0 for label in y]) |
| 87 | + |
| 88 | + ret = minimize(fun=regularized_cost, x0=theta, args=(X, y_i, l), method='TNC', |
| 89 | + jac=regularized_gradient, options={'disp': True}) |
| 90 | + all_theta[i - 1, :] = ret.x |
| 91 | + print(all_theta.shape) |
| 92 | + return all_theta |
| 93 | + |
| 94 | + |
| 95 | +def predict_all(X, all_theta): |
| 96 | + # compute the class probability for each class on each training instance |
| 97 | + h = sigmoid(X @ all_theta.T) # 注意的这里的all_theta需要转置 |
| 98 | + # create array of the index with the maximum probability |
| 99 | + # Returns the indices of the maximum values along an axis. |
| 100 | + h_argmax = np.argmax(h, axis=1) |
| 101 | + # because our array was zero-indexed we need to add one for the true label prediction |
| 102 | + h_argmax = h_argmax + 1 |
| 103 | + |
| 104 | + return h_argmax |
| 105 | + |
| 106 | + |
| 107 | +def load_weight(path): |
| 108 | + data = loadmat(path) |
| 109 | + return data['Theta1'], data['Theta2'] |
| 110 | + |
| 111 | + |
| 112 | +def neural_networks(): |
| 113 | + # Neural Networks |
| 114 | + theta1, theta2 = load_weight('ex3weights.mat') |
| 115 | + print(theta1.shape, theta2.shape) |
| 116 | + X, y = load_data('ex3data1.mat') |
| 117 | + y = y.flatten() |
| 118 | + X = np.insert(X, 0, 1, axis=1) # intercept |
| 119 | + print(X.shape) |
| 120 | + a1 = X |
| 121 | + z2 = a1 @ theta1.T |
| 122 | + z2 = np.insert(z2, 0, 1, axis=1) |
| 123 | + a2 = sigmoid(z2) |
| 124 | + z3 = a2 @ theta2.T |
| 125 | + a3 = sigmoid(z3) |
| 126 | + y_pred = np.argmax(a3, axis=1) + 1 |
| 127 | + accuracy = np.mean(y_pred == y) |
| 128 | + print('accuracy = {0}%'.format(accuracy * 100)) # accuracy = 97.52% |
| 129 | + |
| 130 | + |
| 131 | +def main(): |
| 132 | + x, y = load_data('ex3data1.mat') |
| 133 | + # plot_an_image(x, y) |
| 134 | + # plot_100_image(x) |
| 135 | + x = np.insert(x, 0, 1, axis=1) # (5000, 401) |
| 136 | + y = y.flatten() # 这里消除了一个维度,方便后面的计算 or .reshape(-1) (5000,) |
| 137 | + all_theta = one_vs_all(x, y, 1, 10) # 每一行是一个分类器的一组参数 |
| 138 | + y_pred = predict_all(x, all_theta) |
| 139 | + accuracy = np.mean(y_pred == y) |
| 140 | + print('accuracy = {0}%'.format(accuracy * 100)) |
| 141 | + |
| 142 | + |
| 143 | +if __name__ == '__main__': |
| 144 | + # main() |
| 145 | + neural_networks() |
0 commit comments