Machine Learning(Andrew Ng）ex3.neural network

import matplotlib.pyplot as plt
import numpy as np
import scipy.io as sio
import scipy.optimize as opt
import matplotlib
from sklearn.metrics import classification_report

def load_data(path,transpose=True):
    data=sio.loadmat(path)
    X=data.get('X')#初始数据为m*400
    y=data.get('y')
    y=y.reshape(y.shape[0])
    if transpose:
        X=np.array([im.reshape(20,20).T for im in X] )
        X = np.array([im.reshape(400) for im in X])
    return X,y
        

X,y=load_data('ex3data1.mat')
print(X.shape)
print(y.shape)#ncarry数据类型当是一维向量的时候，shape属性 （5000，）代表是行向量，（5000，1代表是列向量

(5000, 400)
(5000,)

def plot_an_image(image):
    fig,ax=plt.subplots(figsize=(1,1))
    ax.matshow(image.reshape(20,20),cmap=matplotlib.cm.binary)
    plt.xticks(np.array([]))  # 清空坐标轴
    plt.yticks(np.array([])) 

pick_one=4999#np.random.randint(0,5000)
plot_an_image(X[pick_one,:])
plt.show()

def plot_100_image(X):
    sample_idx=np.random.choice(np.arange(100,5000),100)
    sample_images=X[sample_idx,:]
    print(sample_images.shape)
    fig,ax_array=plt.subplots(nrows=10,ncols=10,figsize=(8,8))
    for i in range(10):
        for b in range(10):
            ax_array[i,b].matshow(sample_images[10*i+b].reshape(20,20))
        

plot_100_image(X)
plt.show()

(100, 400)

raw_X, raw_y = load_data('ex3data1.mat')

X = np.insert(raw_X, 0, values=np.ones(raw_X.shape[0]), axis=1)#插入了第一列（全部为1）
raw_y.shape#0是位置，axis区分行和列

(5000,)

y_matrix=[]
for k in range(1,11):
    y_matrix.append((raw_y==k).astype(int))#y标签替换为矩阵，把bool值变换为0/1 raw==k会返回一个（5000，0）的向量，一共返回10个

y_matrix  =[y_matrix[-1]]+y_matrix[:-1] #把代表0的10换到前面来,10->0,同时也
y=np.array(y_matrix)
y.shape

(10, 5000)

training one vs one classification

def sigmoid(z):
    return 1/(1+np.exp(-z))

def cost(theta, X, y):
    return np.mean(-y * np.log(sigmoid(X @ theta)) - (1 - y) * np.log(1 - sigmoid(X @ theta)))

def gradient_decent(theta, x, y):
    return (1 / len(x)) * (((sigmoid(x @ theta) - y).T) @ x)

def regularized_cost(theta,x,y,l=1):
    theta_new=theta[1:]
    regularized_term=(l/(2*len(X)))*np.sum(np.power(theta_new,2))
    return cost(theta,x,y)+regularized_term 

def regularized_gradient(theta,X,y,l=1):
    theta_new=theta[1:]
    regularized_term=np.concatenate([np.array([0]),(l/len(X))*theta_new])
    return np.array(gradient_decent(theta,X,y)+regularized_term)

def logistic_regression(x,y,l=1):
    theta = np.zeros(X.shape[1])

    # train it
    res = opt.minimize(fun=regularized_cost,
                       x0=theta,
                       args=(X, y, l),
                       method='TNC',
                       jac=regularized_gradient,
                       options={'disp': True})
    # get trained parameters
    final_theta = res.x

    return final_theta

def predict(x, theta):
    prob = sigmoid(x @ theta)
    return (prob >= 0.5).astype(int)

t0 = logistic_regression(X, y[0,:])

np.random.rand(4,3)

array([[0.4869344 , 0.03665572, 0.14567694],
       [0.49515102, 0.53962993, 0.80768169],
       [0.99563649, 0.22991872, 0.61401559],
       [0.510981  , 0.93994985, 0.25907377]])

print(t0.shape)
y_pred = predict(X, t0)

(401,)

k_theta=np.array([logistic_regression(X,y[k]) for k in range(10)])
print(k_theta.shape)

(10, 401)

np.arange(5)

array([0, 1, 2, 3, 4])

np.array([5])

array([5])

数据集验证(进行验证）

result=sigmoid(X@k_theta.T)#5000*400 和400*10 最后结果为5000*10
print(result)

[[9.95771804e-01 2.39667000e-09 5.34994023e-04 ... 6.46684260e-05
  3.91525846e-05 1.72080362e-03]
 [9.98344908e-01 9.75193339e-08 5.60550184e-05 ... 9.66506095e-05
  2.90261239e-06 8.48960099e-05]
 [9.91398298e-01 2.59112979e-10 5.68000376e-04 ... 6.54784184e-06
  2.65501671e-02 1.97470938e-03]
 ...
 [6.78400891e-07 4.13968655e-02 3.20831682e-03 ... 1.27221531e-04
  2.97190666e-03 7.07609131e-01]
 [1.84236062e-05 1.32313489e-07 8.86612554e-08 ... 1.64811410e-03
  6.81769126e-02 8.61127936e-01]
 [2.88053490e-02 1.49936381e-10 1.29731954e-04 ... 3.66153840e-01
  4.97620855e-03 1.48243836e-01]]

y_max=np.argmax(result,axis=1)
y_max

array([0, 0, 0, ..., 9, 9, 7], dtype=int64)

y_answer = raw_y.copy()
y_answer[y_answer==10] = 0
z=0
bupt=np.array(y_answer==y_max)
print(bupt)
for i in bupt:
    if(bool(1-i)):
        z+=1
print(z)

[ True  True  True ...  True  True False]
277

precision=(5000-z)/5000
print(precision)#模型精确度为94.46%

0.9446

神经网络

def load_weight(path):
    data=sio.loadmat(path)
    return data['Theta1'],data['Theta2']

$a_i^{(j)}$代表第? 层的第 ? 个激活单元。 ?(?)代表从第 ? 层映射到第? + 1 层时的权重的矩
阵，例如?(1)代表从第一层映射到第二层的权重的矩阵。其尺寸为：以第 ? + 1层的激活单
元数量为行数，以第 ? 层的激活单元数加一为列数的矩阵。例如：隐藏层为3,输入层也为3的神经网络
?(1)的尺寸为 3*4
$$\begin{array}{l}{a}_1^{(2)}=g\left({\Theta }_{10}^{(1)}{x}_0+{\Theta }_{11}^{(1)}{x}_1+{\Theta }_{12}^{(1)}{x}_2+{\Theta }_{13}^{(1)}{x}_3\right)\\ a_2^{(2)}=g\left({\Theta }_{20}^{(1)}{x}_0+{\Theta }_{21}^{(1)}{x}_1+{\Theta }_{22}^{(1)}{x}_2+{\Theta }_{23}^{(1)}{x}_3\right)\\ a_3^{(2)}=g\left({\Theta }_{30}^{(1)}{x}_0+{\Theta }_{31}^{(1)}{x}_1+{\Theta }_{32}^{(1)}{x}_2+{\Theta }_{33}^{(1)}{x}_3\right)\\ h_{\Theta }(x)=g\left({\Theta }_{10}^{(2)}{a}_0^{(2)}+{\Theta }_{11}^{(2)}{a}_1^{(2)}+{\Theta }_{12}^{(2)}{a}_2^{(2)}+{\Theta }_{13}^{(2)}{a}_3^{(2)}\right)\end{array}$$
这里$\theta$下面的的两个坐标可以分别理解为,前一个为生成第二层unit的序号，和前一层unit的序号，这样就描绘了unit之间的箭头

theta1,theta2=load_weight('ex3weights.mat')#维度由上解释
theta1.shape,theta2.shape#第二层有25个单元，最后一层10个单元，

((25, 401), (10, 26))

X, y = load_data('ex3data1.mat',transpose=False)
y[y==10]=0
X = np.insert(X, 0, values=np.ones(X.shape[0]), axis=1)  # intercept

X.shape, y.shape

((5000, 401), (5000,))

feedforward prediction（前向传播）

a1=X

z2 = a1 @ theta1.T # (5000, 401) @ (25,401).T = (5000, 25)
z2 = np.insert(z2, 0, values=np.ones(z2.shape[0]), axis=1)
a2 = sigmoid(z2)
z3 = a2 @ theta2.T
a3 = sigmoid(z3)

y_pred = np.argmax(a3, axis=1)+1
y_pred[y_pred==10]=0#每一行取最大值再加1，matlab从1开始
#y_pred.shape,y.shape
print(y_pred)

print(classification_report(y, y_pred))

              precision    recall  f1-score   support

           0       0.98      0.99      0.99       500
           1       0.97      0.98      0.97       500
           2       0.98      0.97      0.97       500
           3       0.98      0.96      0.97       500
           4       0.97      0.97      0.97       500
           5       0.98      0.98      0.98       500
           6       0.97      0.99      0.98       500
           7       0.98      0.97      0.97       500
           8       0.98      0.98      0.98       500
           9       0.97      0.96      0.96       500

    accuracy                           0.98      5000
   macro avg       0.98      0.98      0.98      5000
weighted avg       0.98      0.98      0.98      5000

因为matlab是从1开始的，但python从0开始，所以存在两种索引方式的转换，使用这行y_pred = np.argmax(a3, axis=1)+1代码的原因，因为前向传播算法提供的矩阵参数（theta）的第一行（也就是python中的0）是用来分类数字1的，而theta的第10行（python中的9），是用来分类数字0的。所以最后结果才要把$h_{\theta }(x)$函数最大值（索引值+1）然后对python中的9（也就是参数对应第10行准换成0）用来对比。