4. 吴恩达机器学习课程-作业4-神经网络学习

javascript:void(0)/article/details/104573555
https://blog.csdn.net/Cowry5/article/details/80399350

def plot_100_images(self, X):
        """随机画100个数字"""
        index = np.random.choice(range(5000), 100)
        images = X[index]
        fig, ax_array = plt.subplots(10, 10, sharey=True, sharex=True, figsize=(8, 8))
        for r in range(10):
            for c in range(10):
                ax_array[r, c].matshow(images[r * 10 + c].reshape(20, 20), cmap='gray_r')
        plt.xticks([])
        plt.yticks([])
        plt.show()

    def main(self):
        pass

    def __init__(self):
        data = sio.loadmat('ex4data1.mat')  # return a dict
        X = np.array(data['X'])
        y = np.array(data['y']).flatten()
        print(X.shape)
        print(y.shape)
        # 随机画出100个图片
        self.plot_100_images(X)
        pass

(5000, 400)
(5000,)

def serialize(self, a, b):
        """
        展开参数，扁平化参数
        :param a:
        :param b:
        :return: 拼接后的矩阵
        """
        # np.r_是按列连接两个矩阵，就是把两矩阵上下相加，要求列数相等。
        # 拼接矩阵的意思
        return np.r_[a.flatten(), b.flatten()]

theta = serialize(t1, t2)  # 扁平化参数，25*401+10*26=10285
theta.shape  # (10285,)

def deserialize(self, seq):
        """
        提取参数,扁平化的参数，还原theta1,theta2
        :return:  theta1,theta2
        """
        return seq[:25 * 401].reshape(25, 401), seq[25 * 401:].reshape(10, 26)

def feed_forward(self, theta, X):
        """
        得到每层的输入和输出
        前置输出
        :param theta: 扁平化后的theta
        :param X:
        :return:
        """
        t1, t2 = self.deserialize(theta)
        # 前面已经插入过偏置单元，这里就不用插入了
        a1 = X
        z2 = a1.dot(t1.T)
        a2 = np.insert(self.sigmoid(z2), 0, 1, axis=1)
        z3 = a2.dot(t2.T)
        a3 = self.sigmoid(z3)

        return a1, z2, a2, z3, a3
      
def sigmoid(self, z):
        """
        sigmoid函数实现了将结果从R转化到0-1，用来表示概率。实现sigmoid函数
        :param z: z值
        :return:
        """
        return 1 / (1 + np.exp(-z))

def cost(self, theta, X, y):
        """
        计算神经网络的代价函数（不带正则化项）
        :param theta: 扁平化的参数
        :param X:
        :param y:
        :return: 代价值
        """
        a1, z2, a2, z3, h = self.feed_forward(theta, X)
        J = 0
        for i in range(len(X)):
            first = - y[i] * np.log(h[i])
            second = (1 - y[i]) * np.log(1 - h[i])
            J = J + np.sum(first - second)
        J = J / len(X)
        return J

def regularized_cost(self, theta, X, y, l=1):
        """
        正则化时忽略每层的偏置项，也就是参数矩阵的第一列'
        :param theta: 扁平化后的参数
        :param X:
        :param y:
        :param l: 惩罚力度
        :return: 正则化后的代价函数
        """
        t1, t2 = self.deserialize(theta)
        reg = np.sum(t1[:, 1:] ** 2) + np.sum(t2[:, 1:] ** 2)  # or use np.power(a, 2)
        return l / (2 * len(X)) * reg + self.cost(theta, X, y)

def sigmoid_gradient(self,z):
        """
        这里可以手动推导，并不难。
        :param z:
        :return:
        """
        return self.sigmoid(z) * (1 - self.sigmoid(z))

def random_init(self,size):
        """
        从服从的均匀分布的范围中随机返回size大小的值
        :param size: 随机取出的数量
        :return: 
        """
        return np.random.uniform(-0.12, 0.12, size)

def gradient(theta, X, y):
    '''
    unregularized gradient, notice no d1 since the input layer has no error 
    return 所有参数theta的梯度，故梯度D(i)和参数theta(i)同shape，重要。
    '''
    t1, t2 = deserialize(theta)
    a1, z2, a2, z3, h = feed_forward(theta, X)
    d3 = h - y # (5000, 10)
    d2 = d3 @ t2[:,1:] * sigmoid_gradient(z2)  # (5000, 25)
    D2 = d3.T @ a2  # (10, 26)
    D1 = d2.T @ a1 # (25, 401)
    D = (1 / len(X)) * serialize(D1, D2)  # (10285,)
    
    return D

def a_numeric_grad(self, plus, minus, X, y, e):
        """
        对每个参数theta_i计算数值梯度，即理论梯度。
        :param e:
        :param y:
        :param X:
        :param plus:
        :param minus:
        :return:
        """
        return (self.regularized_cost(plus, X, y) - self.regularized_cost(minus, X, y)) / (e * 2)

    def gradient_checking(self, theta, X, y, e):
        """
        荼毒检查
        :param theta:
        :param X:
        :param y:
        :param e:
        :return:
        """
        numeric_grad = []
        for i in range(len(theta)):
            plus = theta.copy()  # deep copy otherwise you will change the raw theta
            minus = theta.copy()
            plus[i] = plus[i] + e
            minus[i] = minus[i] - e
            grad_i = self.a_numeric_grad(plus, minus, X, y, e)
            numeric_grad.append(grad_i)
        numeric_grad = np.array(numeric_grad)
        analytic_grad = self.regularized_gradient(theta, X, y)
        diff = np.linalg.norm(numeric_grad - analytic_grad) / np.linalg.norm(numeric_grad + analytic_grad)

        print(
            'If your backpropagation implementation is correct,\nthe relative difference will be smaller than 10e-9 (assume epsilon=0.0001).\nRelative Difference: {}\n'.format(
                diff))

def regularized_gradient(self, theta, X, y, l):
        """
        反向传播算法
        :param theta: ndarray，一维参数向量
        :param X: ndarray,输入层的输入值
        :param y: ndarray,数据的标签
        :param l: float,惩罚参数
        :return: ndarray，下降后一维参数向量
        """
        a1, z2, a2, z3, a3 = self.feed_forward(theta, X)
        theta1, theta2 = self.deserialize(theta)  # theta1(25,401), theta2(10,26)
        m = X.shape[0]
        d3 = a3 - y  # d3(5000,10)
        d2 = d3 @ theta2[:, 1:] * self.sigmoid_gradient(z2)  # d2(5000,25)
        theta1 = np.insert(np.delete(theta1, 0, axis=1), 0, 0, axis=1)
        theta2 = np.insert(np.delete(theta2, 0, axis=1), 0, 0, axis=1)
        D1 = (1 / m) * d2.T.dot(a1) + (l / m) * theta1  # D1(25,401)
        D2 = (1 / m) * d3.T.dot(a2) + (l / m) * theta2  # D2(10,26)
        return self.serialize((D1, D2))

def nn_training(self, X, y):
        """
        开始训练
        :param X:
        :param y:
        :return:
        """
        # 初始化参数
        print("========开始初始化参数，开始训练============")
        init_theta = self.random_init(10285)  # 25*401 + 10*26
        # 参数优化
        res = opt.minimize(fun=self.regularized_cost,
                           x0=init_theta,
                           args=(X, y, 1),
                           method='TNC',
                           jac=self.regularized_gradient,
                           options={'maxiter': 400})
        return res

X = self.X
        y = self.y
        res = self.nn_training(X, y)
        print("========训练结束，开始预测===========")
        _, _, _, _, h = self.feed_forward(theta=res.x, X=X)
        y_pred = np.argmax(h, axis=1) + 1
        print(y_pred.shape)
        print(y_pred)
        # 开始预测，评价
        # self.raw_y (5000,) y_pred(5000,)
        print("self.raw_y.shape---", self.raw_y.shape)
        print("y_pred.shape", y_pred.shape)
        print("打印算法评估", classification_report(self.raw_y, y_pred))

def plot_hidden(self, theta):
        t1, _ = self.deserialize(theta)
        t1 = t1[:, 1:]
        fig, ax_array = plt.subplots(5, 5, sharex=True, sharey=True, figsize=(6, 6))
        for r in range(5):
            for c in range(5):
                ax_array[r, c].matshow(t1[r * 5 + c].reshape(20, 20), cmap='gray_r')
                plt.xticks([])
                plt.yticks([])
        plt.show()

self.plot_hidden(res.x)

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


class BP:

    def regularized_gradient(self, theta, X, y, l=1):
        """
        反向传播算法
        :param theta: ndarray，一维参数向量
        :param X: ndarray,输入层的输入值
        :param y: ndarray,数据的标签
        :param l: float,惩罚参数
        :return: ndarray，下降后一维参数向量
        """

        # a1, z2, a2, z3, h = self.feed_forward(theta, X)
        D1, D2 = self.deserialize(self.gradient(theta, X, y))
        t1, t2 = self.deserialize(theta)
        t1[:, 0] = 0
        t2[:, 0] = 0
        reg_D1 = D1 + (l / len(X)) * t1
        reg_D2 = D2 + (l / len(X)) * t2
        return self.serialize(reg_D1, reg_D2)

    def gradient(self, theta, X, y):
        """
        获取整个网络代价函数的梯度。以便在优化算法中求解。
        不规则梯度，注意没有d1，因为输入层没有输入
        :param theta: 扁平化的参数，即权重w1,w2
        :param X:
        :param y:
        :return: 所有参数theta的梯度，故梯度D(i)和参数theta(i)同shape，重要。
        """
        # t1=theta1,t2=theta2
        t1, t2 = self.deserialize(theta)
        # 使用前向传播算法计算出所有的数值
        a1, z2, a2, z3, a3 = self.feed_forward(theta, X)
        # 均方误差的梯度(5000, 10)
        d3 = a3 - y  # (5000, 10)
        # 输出层误差，t2为w2的意思，是隐藏层到输入层的权值
        d2 = d3 @ t2[:, 1:] * self.sigmoid_gradient(z2)  # (5000, 25)

        # a2为隐藏层输出,D2为梯度
        D2 = d3.T @ a2  # (10, 26)
        # D1为损失函数值
        D1 = d2.T @ a1  # (25, 401)

        D = (1 / len(X)) * self.serialize(D1, D2)  # (10285,)
        return D

    def a_numeric_grad(self, plus, minus, X, y, e):
        """
        对每个参数theta_i计算数值梯度，即理论梯度。
        :param e:
        :param y:
        :param X:
        :param plus:
        :param minus:
        :return:
        """
        return (self.regularized_cost(plus, X, y) - self.regularized_cost(minus, X, y)) / (e * 2)

    def gradient_checking(self, theta, X, y, e):
        """
        荼毒检查
        :param theta:
        :param X:
        :param y:
        :param e:
        :return:
        """
        numeric_grad = []
        for i in range(len(theta)):
            plus = theta.copy()  # deep copy otherwise you will change the raw theta
            minus = theta.copy()
            plus[i] = plus[i] + e
            minus[i] = minus[i] - e
            grad_i = self.a_numeric_grad(plus, minus, X, y, e)
            numeric_grad.append(grad_i)
        numeric_grad = np.array(numeric_grad)
        analytic_grad = self.regularized_gradient(theta, X, y)
        diff = np.linalg.norm(numeric_grad - analytic_grad) / np.linalg.norm(numeric_grad + analytic_grad)

        print(
            'If your backpropagation implementation is correct,\nthe relative difference will be smaller than 10e-9 (assume epsilon=0.0001).\nRelative Difference: {}\n'.format(
                diff))

    def sigmoid_gradient(self, z):
        """
        这里可以手动推导，并不难。
        :param z:
        :return: sigmoid(x) 导数值
        """
        return self.sigmoid(z) * (1 - self.sigmoid(z))

    def random_init(self, size):
        """
        从服从的均匀分布的范围中随机返回size大小的值
        :param size: 随机取出的数量
        :return:
        """
        return np.random.uniform(-0.12, 0.12, size)

    def expand_y(self, y):
        """
        首先我们要将标签值（1，2，3，4，…，10）转化成非线性相关的向量，
        向量对应位置（y[i-1]）上的值等于1，例如y[0]=6转化为y[0]=[0,0,0,0,0,1,0,0,0,0]。
        :param y:
        :return:
        """
        result = []
        # 把y中每个类别转化为一个向量，对应的lable值在向量对应位置上置为1
        for i in y:
            y_array = np.zeros(10)
            y_array[i - 1] = 1
            result.append(y_array)

        return np.array(result)

    def serialize(self, a, b):
        """
        展开参数，扁平化参数
        :param a:
        :param b:
        :return: 拼接后的矩阵
        """
        # np.r_是按列连接两个矩阵，就是把两矩阵上下相加，要求列数相等。
        # 拼接矩阵的意思
        return np.r_[a.flatten(), b.flatten()]

    def deserialize(self, seq):
        """
        提取参数,扁平化的参数，还原theta1,theta2
        :return:  theta1,theta2
        """
        return seq[:25 * 401].reshape(25, 401), seq[25 * 401:].reshape(10, 26)

    def feed_forward(self, theta, X):
        """
        得到每层的输入和输出
        前置输出
        :param theta: 扁平化后的theta
        :param X:
        :return:
        """
        t1, t2 = self.deserialize(theta)
        # 前面已经插入过偏置单元，这里就不用插入了
        # X (5000,401) t1 (25,401)
        a1 = X
        # z2 (5000,25)
        z2 = a1.dot(t1.T)
        # a2 (5000,26)
        a2 = np.insert(self.sigmoid(z2), 0, 1, axis=1)
        # a2 (5000,26) .t2(10,26) z3(5000,10)
        z3 = a2.dot(t2.T)
        # a3 (5000,10)
        a3 = self.sigmoid(z3)

        return a1, z2, a2, z3, a3

    def cost(self, theta, X, y):
        """
        计算神经网络的代价函数（不带正则化项）
        :param theta: 扁平化的参数
        :param X:
        :param y:
        :return: 代价值
        """
        a1, z2, a2, z3, h = self.feed_forward(theta, X)
        J = 0
        for i in range(len(X)):
            first = - y[i] * np.log(h[i])
            second = (1 - y[i]) * np.log(1 - h[i])
            J = J + np.sum(first - second)
        J = J / len(X)
        return J

    def regularized_cost(self, theta, X, y, l=1):
        """
        正则化时忽略每层的偏置项，也就是参数矩阵的第一列'
        :param theta: 扁平化后的参数
        :param X:
        :param y:
        :param l: 惩罚力度
        :return: 正则化后的代价函数
        """
        t1, t2 = self.deserialize(theta)
        reg = np.sum(t1[:, 1:] ** 2) + np.sum(t2[:, 1:] ** 2)  # or use np.power(a, 2)
        return l / (2 * len(X)) * reg + self.cost(theta, X, y)

    def sigmoid(self, z):
        """
        sigmoid函数实现了将结果从R转化到0-1，用来表示概率。实现sigmoid函数
        :param z: z值
        :return:
        """
        return 1 / (1 + np.exp(-z))

    def plot_100_images(self, X):
        """
        随机画100个数字
        :param X:
        :return:
        """
        index = np.random.choice(range(5000), 100)
        images = X[index]
        fig, ax_array = plt.subplots(10, 10, sharey=True, sharex=True, figsize=(8, 8))
        for r in range(10):
            for c in range(10):
                ax_array[r, c].matshow(images[r * 10 + c].reshape(20, 20), cmap='gray_r')
        plt.xticks([])
        plt.yticks([])
        plt.show()

    def plot_hidden(self, theta):
        t1, _ = self.deserialize(theta)
        t1 = t1[:, 1:]
        fig, ax_array = plt.subplots(5, 5, sharex=True, sharey=True, figsize=(6, 6))
        for r in range(5):
            for c in range(5):
                ax_array[r, c].matshow(t1[r * 5 + c].reshape(20, 20), cmap='gray_r')
                plt.xticks([])
                plt.yticks([])
        plt.show()

    def nn_training(self, X, y):
        """
        开始训练
        :param X:
        :param y:
        :return:
        """
        # 初始化参数
        print("========开始初始化参数，开始训练============")
        init_theta = self.random_init(10285)  # 25*401 + 10*26
        # 参数优化
        res = opt.minimize(fun=self.regularized_cost,
                           x0=init_theta,
                           args=(X, y, 1),
                           method='TNC',
                           jac=self.regularized_gradient,
                           options={'maxiter': 400})
        return res

    def main(self):
        X = self.X
        y = self.y
        res = self.nn_training(X, y)
        print("========训练结束，开始预测===========")
        _, _, _, _, h = self.feed_forward(theta=res.x, X=X)
        y_pred = np.argmax(h, axis=1) + 1
        print(y_pred.shape)
        print(y_pred)
        # 开始预测，评价
        # self.raw_y (5000,) y_pred(5000,)
        print("self.raw_y.shape---", self.raw_y.shape)
        print("y_pred.shape", y_pred.shape)
        print("打印算法评估", classification_report(self.raw_y, y_pred))
        # 可视化参数
        self.plot_hidden(res.x)

    def __init__(self):
        self.data = sio.loadmat('ex4data1.mat')  # return a dict
        theta = sio.loadmat("ex4weights.mat")
        self.theta1 = theta["Theta1"]  # (25,401)
        self.theta2 = theta["Theta2"]  # (10,26)
        X = np.array(self.data["X"])
        self.raw_y = np.array(self.data["y"]).flatten()
        # 随机画出100个图片
        self.plot_100_images(X)
        y = self.expand_y(self.data["y"])
        self.y = np.array(y)
        self.X = np.insert(X, 0, 1, axis=1)  # 添加一列1，方便计算b

        # 扁平化参数，25*401+10*26=10285
        self.theta = self.serialize(self.theta1, self.theta2)


if __name__ == '__main__':
    begin_time = time()
    print("========程序开始============")
    obj = BP()
    obj.main()
    print("========程序结束============")
    end_time = time()
    run_time = end_time - begin_time
    print('该循环程序运行时间：', run_time)  # 该循环程序运行时间： 1.4201874732