吴恩达机器学习python作业之逻辑回归解决多元分类问题

参考链接：
https://blog.csdn.net/weixin_44750583/article/details/88603512
python读取文件——python读取和保存mat文件_行志的博客-CSDN博客_python读取mat数据

注意点1：本来打算进行矩阵计算，于是直接将所有参加矩阵运算的参数转换成matrix类型。最后发现optimize.fmin_ncg中的X0参数类型要求是array类型。
注意点2：不可变类型进行传参时，函数内部会修改参数。

import numpy as np

a = np.array([1,2,3,4,5])

def test(a):
    b = a
    b[0] = 0
    print(a) [1,2,3,4,5]

test(a)
print(a)
[1,2,3,4,5]

#定义正则化的梯度函数
#错误原因：theta传参进去后，temp与theta指向了同一个内存空间
#所以每次重新迭代时，theta_0都会被重新设置为0，这样一来最后的准确率就是37%左右
def gradientRegularzation(theta,x,y,lamda = 0.1):
    temp = theta
    temp[0] = 0
    return gradient(theta,x,y) + lamda / m * theta

import numpy as np
import pandas as pd
from scipy import io
from matplotlib import pyplot as plt
from scipy import optimize as op

#1.使用scipy库中的loadmat函数读取数据并进行格式转换
#数据集分布：ex3data1中使用dict格式存储数据X，一共5000行，每行20*20个像素点roll开后变成400，因此样本数为5000,特征数量为400。
#同时还使用字典格式存储了Y，从0-9，其中0代表10
dt = io.loadmat("E:\机器学习\吴恩达\data_sets\ex3data1.mat")
x = np.array(dt["X"]) #(5000, 400)
y = np.array(dt["y"]) #(5000, 1)
m = x.shape[0] #5000
cols = x.shape[1]+1

#2.查看某一行的数据集
# temp = x[1040,:].reshape(20,20)
# temp = temp.T
# plt.imshow(temp,cmap='gray')
# plt.show()
# print(y[1040]) #[2]

#3.计算代价函数与梯度函数
#定义sigmoid函数
def sigmoid(z):
    return 1/(1+np.exp(-z))

#定义未正则化的代价函数
#实际调用时，theta,x,y都是array类型
def costFunction(theta,x,y):
    h = sigmoid(x@theta) #(5000,)
    temp = y @ np.log(h)+(1-y) @ np.log(1-h)
    return -1 / m * temp


#定义正则化的代价函数
def costFunctionRegularzation(theta,x,y,lamda = 0.1):
    #正则化项中theta_0设置为0
    theta_temp = theta[1:]
    return costFunction(theta,x,y) + lamda / (2*m) * np.sum(theta_temp * theta_temp)

#定义未正则化的梯度函数
def gradient(theta,x,y):
    h = sigmoid(x @ theta) #(5000, 1)
    temp = x.T @ (h - y)  #(401,)
    return 1 / m * temp



#定义正则化的梯度函数
def gradientRegularzation(theta,x,y,lamda = 0.1):
    temp = lamda / m * theta
    temp[0] = 0
    return gradient(theta,x,y) + temp


#多元分类函数
def multiClassfication(x,y):
    real_theta = np.zeros((10,cols)) #(10,401)的矩阵
    #real_theta的每一行表示一个数字对应的特征，一共10行
    #计算每一个数字对应的特征可以直接使用按行来比对
    #每走一次循环，训练好表示一个数字的一行特征出来
    #其中手写数字0-9 对应的标签是10,1-9，因此选择映射的方案为第0行表示标签为1的特征
    for i in range(10):
        real_y = np.array([ 1 if j == i+1 else 0 for j in y ])
        #real_y构造出一个新的(5000,1)的矩阵
        # 如果我们这轮寻找的是对应2的特征，那么real_y是原来y中为2的赋值为1，不为2的赋值为0
        #x0必须是ndarray
        temp_theta = op.fmin_ncg(f=costFunctionRegularzation,fprime=gradientRegularzation,x0 = real_theta[i:i+1,:],args=(x,real_y),maxiter=400)
        real_theta[i:i+1,:] = temp_theta
    return real_theta.T


#4.做好调用函数的数据准备
#进行矩阵乘法时最好还是将array类型转换为matrix类型，因为martrix类型更符合一般线性代数矩阵乘法的原则
#转换为矩阵后，最好使用@实现矩阵乘法
a = np.ones((m,1))
x = np.concatenate((a,x),1)
theta = multiClassfication(x,y)



#5.进行调用函数
# 预测的概率矩阵 (5000, 10)
def probability(x, theta):
    return sigmoid(x @ theta)


# 预测的y值
def predict(prob):
    y_predict = np.zeros((prob.shape[0], 1))
    for i in range(prob.shape[0]):
        # 查找第i行的最大值并返回它所在的位置,再加1就是对应的类别
        y_predict[i] = np.unravel_index(np.argmax(prob[i, :]), prob[i, :].shape)[0] + 1
    return y_predict


# 精度
def accuracy(y_predict, y=y):
    m = y.size
    count = 0
    for i in range(y.shape[0]):
        if y_predict[i] == y[i]:
            j = 1
        else:
            j = 0
        count = j + count  # 计数预测值和期望值相等的项
    return count / m


prob = probability(x, theta)
y_predict = predict(prob)
accuracy(y_predict)
print('accuracy = {0}%'.format(accuracy(y_predict) * 100))