1 importnumpy as np2 importscipy.io as scio3 importmatplotlib.pyplot as plt4 importscipy.optimize as op5

6 #X:5000*400

7 #Y:5000*10

8 #a1:5000*401(后5000*400)

9 #z2:5000*25

10 #a2:5000*26(后5000*25)

11 #z3:5000*10

12 #a3:5000*10

13 #Theta1:25*401

14 #Theta2:10*26

15 #delta3:5000*10

16 #delta2:5000*25

17 #bigDelta1:25*401

18 #bigDelta2:10*26

19 #Theta1_grad:25*401

20 #Theta2_grad:10*26

21

22

23 #显示图片数据

24 defdisplayData(X):25 m = np.size(X, 0) #X的行数，即样本数量

26 n = np.size(X, 1) #X的列数，即单个样本大小

27 example_width = int(np.round(np.sqrt(n))) #单张图片宽度

28 example_height = int(np.floor(n / example_width)) #单张图片高度

29 display_rows = int(np.floor(np.sqrt(m))) #显示图中，一行多少张图

30 display_cols = int(np.ceil(m / display_rows)) #显示图中，一列多少张图片

31 pad = 1 #图片间的间隔

32 display_array = - np.ones((pad + display_rows * (example_height +pad),33 pad + display_cols * (example_width + pad))) #初始化图片矩阵

34 curr_ex = 0 #当前的图片计数

35 #将每张小图插入图片数组中

36 for j inrange(0, display_rows):37 for i inrange(0, display_cols):38 if curr_ex >=m:39 break

40 max_val =np.max(abs(X[curr_ex, :]))41 jstart = pad + j * (example_height +pad)42 istart = pad + i * (example_width +pad)43 display_array[jstart: (jstart + example_height), istart: (istart + example_width)] =\44 np.array(X[curr_ex, :]).reshape(example_height, example_width) /max_val45 curr_ex = curr_ex + 1

46 if curr_ex >=m:47 break

48 display_array =display_array.T49 plt.imshow(display_array,cmap=plt.cm.gray)50 plt.axis('off')51 plt.show()52

53

54 #计算hθ(z)

55 defsigmoid(z):56 g = 1.0 / (1.0 + np.exp(-z))57 returng58

59

60 #初始化Θ，保持在[-ε，ε]

61 defrandInitializeWeights(sizeList):62 epsilon_init = 0.12

63 theta1_lx = sizeList['theta1_lx']64 theta1_ly = sizeList['theta1_ly']65 theta2_lx = sizeList['theta2_lx']66 theta2_ly = sizeList['theta2_ly']67 theta_size = theta1_lx * theta1_ly + theta2_lx *theta2_ly68 W = np.random.uniform(-epsilon_init, epsilon_init, theta_size)69 returnW70

71

72 #把一维的矩阵改写为多维

73 defchangeForm(theta_vector, theta1_lx, theta1_ly, theta2_lx, theta2_ly):74 theta1 = np.array(theta_vector[0: theta1_lx *theta1_ly]).reshape(theta1_lx, theta1_ly)75 theta2 = np.array(theta_vector[theta1_lx * theta1_ly: theta1_lx * theta1_ly + theta2_lx *theta2_ly])\76 .reshape(theta2_lx, theta2_ly)77 theta = {'Theta1': theta1, 'Theta2': theta2}78 returntheta79

80

81 #计算正向激励的参数a

82 defcomputeA(nn_params, X):83 theta1 = nn_params['Theta1']84 theta2 = nn_params['Theta2']85 m =np.size(X, 0)86

87 #第二层计算

88 one =np.ones(m)89 a1 = np.insert(X, 0, values=one, axis=1)90 a2 =sigmoid(np.dot(a1, theta1.T))91 #第三层计算

92 one =np.ones(np.size(a2, 0))93 a2 = np.insert(a2, 0, values=one, axis=1)94 a3 =sigmoid(np.dot(a2, theta2.T))95 a_res = {'a1': a1, 'a2': a2, 'a3': a3}96 returna_res97

98

99 #计算g'(z)

100 defsigmoidGradient(z):101 g = np.multiply(sigmoid(z), 1 -sigmoid(z))102 returng103

104

105 #计算 J

106 defnnCostFunction(nn_params, X, Y, lamb, sizeList):107 theta =changeForm(nn_params,108 sizeList['theta1_lx'], sizeList['theta1_ly'],109 sizeList['theta2_lx'], sizeList['theta2_ly'])110 theta1 = theta['Theta1']111 theta2 = theta['Theta2']112 m =np.size(X, 0)113 a_res =computeA(theta, X)114 a3 = a_res['a3']115 #计算J

116 J = 1 / m * np.sum(-np.multiply(Y, np.log(a3)) - np.multiply((1 - Y), np.log(1 -a3)))117 #规格化

118 theta1_copy = theta1[:, 1:]119 theta2_copy = theta2[:, 1:]120 J = J + lamb / (2 * m) * (np.sum(theta1_copy ** 2) + np.sum(theta2_copy ** 2))121 print(J)122 returnJ123

124

125 #计算 D

126 defnnGradient(nn_params, X, Y, lamb, sizeList):127 theta =changeForm(nn_params,128 sizeList['theta1_lx'], sizeList['theta1_ly'],129 sizeList['theta2_lx'], sizeList['theta2_ly'])130 theta1 = theta['Theta1']131 theta2 = theta['Theta2']132 m =np.size(X, 0)133 a_res =computeA(theta, X)134 a1 = a_res['a1']135 a2 = a_res['a2']136 a3 = a_res['a3']137 theta1_copy = theta1[:, 1:]138 theta2_copy = theta2[:, 1:]139 #计算δ

140 delta3 = a3 -Y141 delta2 =np.multiply(np.dot(delta3, theta2_copy), sigmoidGradient(np.dot(a1, theta1.T)))142 #计算Δ

143 bigDeilta1 =np.dot(delta2.T, a1)144 bigDeilta2 =np.dot(delta3.T, a2)145 #计算D

146 theta1_grad = bigDeilta1 / m + lamb / m *theta1147 theta2_grad = bigDeilta2 / m + lamb / m *theta2148 theta1_grad[:, 0] = bigDeilta1[:, 0] /m149 theta2_grad[:, 0] = bigDeilta2[:, 0] /m150 #当使用高级优化方法来优化神经网络时，需要将多个参数矩阵展开，才能传入优化函数

151 grad =np.r_[theta1_grad.flatten(), theta2_grad.flatten()]152 #print(np.size(grad))

153 returngrad154

155

156 #测试参数的初始化

157 defdebugInitializeWeights(L_out, L_in):158 W = np.arange(1, L_out * (L_in + 1)+1)159 W =np.sin(W)160 W = np.array(W).reshape(L_out, (L_in + 1)) / 10;161 returnW162

163

164 #数值方法计算梯度

165 defcomputeNumericalGradient(theta, X, Y ,lamb, sizeList):166 numgrad =np.zeros(np.size(theta))167 perturb =np.zeros(np.size(theta))168 e = 1e-4

169 for p inrange(0, np.size(theta)):170 perturb[p] =e171 theta_minus = theta -perturb172 theta_plus = theta +perturb173 loss1 =nnCostFunction(theta_minus, X, Y, lamb, sizeList)174 loss2 =nnCostFunction(theta_plus, X, Y, lamb, sizeList)175 numgrad[p] = (loss2 - loss1) / (2 *e)176 perturb[p] =0177 returnnumgrad178

179

180 #梯度检测函数

181 defcheckNNGradients(lamb):182 #设置测试参数

183 input_layer_size = 3;184 hidden_layer_size = 5;185 num_labels = 3;186 lamb = 1

187 m = 5;188 sizeList = {'theta1_lx': hidden_layer_size,189 'theta1_ly': input_layer_size + 1,190 'theta2_lx': num_labels,191 'theta2_ly': hidden_layer_size + 1} #保存θ大小的参数

192 theta1 =debugInitializeWeights(hidden_layer_size, input_layer_size)193 theta2 =debugInitializeWeights(num_labels, hidden_layer_size)194 theta =np.r_[theta1.flatten(), theta2.flatten()]195 X = debugInitializeWeights(m, input_layer_size - 1)196 y = np.random.randint(0, num_labels, (m, 1))197 #对y进行改写，改为 m*num_labels 规格的矩阵

198 Y =np.zeros((m, num_labels))199 for i inrange(0, m):200 Y[i, y[i, 0]] = 1

201 grad =nnGradient(theta, X, Y, lamb, sizeList)202 numGrad =computeNumericalGradient(theta, X, Y, lamb, sizeList)203 diff = np.linalg.norm(numGrad - grad) / np.linalg.norm(numGrad +grad)204 print('check NN Gradient: diff =', diff)205

206

207 #使用模型进行预测

208 defpredict(theta1, theta2, X):209 m =np.size(X,0)210 p = np.zeros((np.size(X, 0), 1))211 #第二层计算

212 one =np.ones(m)213 X = np.insert(X, 0, values=one, axis=1)214 a2 =sigmoid(np.dot(X, theta1.T))215 #第三层计算

216 one =np.ones(np.size(a2,0))217 a2 = np.insert(a2, 0, values=one, axis=1)218 a3 =sigmoid(np.dot(a2, theta2.T))219 p = a3.argmax(axis=1) + 1 #y的值为1-10，所以此处0-9要加1

220 returnp.flatten()221

222

223 #——————————————主函数————————————————————

224 #初始化数据

225 input_layer_size = 400

226 hidden_layer_size = 25

227 num_labels = 10

228 sizeList = {'theta1_lx': hidden_layer_size,229 'theta1_ly': input_layer_size + 1,230 'theta2_lx': num_labels,231 'theta2_ly': hidden_layer_size + 1} #保存θ大小的参数

232 lamb = 1

233

234 #加载数据文件

235 data = scio.loadmat('ex4data1.mat')236 X = data['X']237 m =np.size(X, 0)238 y = data['y']239 #对y进行改写，改为5000*10规格的矩阵，第0-9个位置分别表示1，2，...，9，0

240 Y =np.zeros((m, num_labels))241 for i inrange(0, m):242 Y[i, y[i, 0] - 1] = 1

243 rand_indices = np.random.randint(0, m, 100)244 sel =X[rand_indices, :]245 displayData(sel)246

247 #测试数据θ

248 theta = scio.loadmat('ex4weights.mat')249 theta1 = theta['Theta1']250 theta2 = theta['Theta2']251 nn_theta =np.r_[theta1.flatten(), theta2.flatten()]252

253 #测试nnCostFunction

254 #J = nnCostFunction(nn_theta, X, Y, 3, sizeList)

255 #print(J)

256

257 #测试nnGradient

258 print(nnGradient(nn_theta, X, Y, lamb, sizeList))259

260 #初始化参数

261 nn_params =randInitializeWeights(sizeList)262

263 #梯度检测

264 #checkNNGradients(lamb)

265

266 #训练模型

267 res = op.minimize(fun=nnCostFunction,268 x0=nn_params,269 args=(X, Y, lamb, sizeList),270 method='TNC',271 jac=nnGradient,272 options={'maxiter': 100})273 print(res)274

275 #计算准确率

276 all_theta = changeForm(res.x, sizeList['theta1_lx'], sizeList['theta1_ly'],277 sizeList['theta2_lx'], sizeList['theta2_ly'])278 res_theta1 = all_theta['Theta1']279 res_theta2 = all_theta['Theta2']280 pred =predict(res_theta1, res_theta2, X)281 acc = np.mean(pred == y.flatten())*100

282 print('Training Set Accuracy:',acc,'%')283

284 #显示中间隐藏层

285 displayData(res_theta1[:, 1:])