吴恩达机器学习第三章作业答案 part1：多类别逻辑回归

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.io import loadmat
from sklearn.metrics import classification_report
from scipy.optimize import minimize


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


#  向量化计算代价函数
def cost(theta, X, y, learnRate):
    # print('here 1')
    theta = np.mat(theta)
    X = np.mat(X)
    y = np.mat(y)

    first = np.multiply(-y, np.log(sigmoid(X * theta.T)))

    # 当theta = np.ones((1, X.shape[1]))的时候，要采用下列代码，即加上1e-5
    # 这里是因为sigmoid（X*theta.T）得到的接近于1，则如果不带1e-5的话，log得到的结果接近于无穷，造成计算错误
    # second = np.multiply(1 - y, np.log(1 - sigmoid(X * theta.T) + 1e-5))

    second = np.multiply(1 - y, np.log(1 - sigmoid(X * theta.T)))
    reg = np.sum(np.power(theta[:, 1:], 2)) * (learnRate / (2 * X.shape[0]))
    return np.sum(first - second) / X.shape[0] + reg


#  向量化计算梯度
def gradient(theta, X, y, learnRate):
    # theta, X, y均是np.mat
    # print('here 2')
    # 这里theta必须强制转换为np.mat不管是行向量还是列向量，因为minimize里的x0是（n，）的ndarray行向量
    theta = np.mat(theta)
    X = np.mat(X)
    y = np.mat(y)

    g = (X.T * (sigmoid(X * theta.T) - y) + (learnRate * theta.T)) / X.shape[0]
    g[0, 0] = np.sum(sigmoid(X * theta.T) - y) / X.shape[0]
    # 返回的梯度向量是一个行向量或者列向量都可以正常运行
    # 猜测minimize函数内部对向量进行加减法的话是对应位置元素进行，所以不管是行向量还是列向量都是一样的
    return g.T


#  计算每一个类别的theta
def one_vs_all(X, y, num_labels, learnRate):
    all_theta = np.zeros((num_labels, X.shape[1]))
    # print(type(all_theta))
    for i in range(1, num_labels + 1):
        y_i = np.mat([1 if a == i else 0 for a in y]).T

        # theta如果采用全1的矩阵，可能会出现log计算得到无穷的情况，需要对cost函数进行一定的修正
        # 因此一般采用全零作为theta的初始值
        # theta = np.ones((1, X.shape[1]))
        theta = np.zeros((1, X.shape[1]))

        # 如果有X * theta行列数不对应的错误是因为minimize函数里的x0的类型是np.ndarray而且其shape是(n, )，
        # 所以在自己的cost和gradient函数中必须对theta进行矩阵的转换
        all_theta[i - 1] = minimize(fun=cost, x0=theta, args=(X, y_i, learnRate), method='TNC', jac=gradient).x
    return all_theta


def predict_all(X, all_theta):
    X = np.mat(X)
    all_theta = np.mat(all_theta)
    # 这里其实有没有激活函数sigmoid都是一样的，因为激活函数时单调增的，找最大值所在的序号其实有没有激活函数没区别
    predict = sigmoid(X * all_theta.T)
    # predict = X * all_theta.T
    h_argmax = np.argmax(predict, axis=1)
    return h_argmax + 1


data = loadmat('ex3data1.mat')

# the type of X and y is np.array
X = np.mat(np.insert(data['X'], 0, 1, axis=1))
y = np.mat(data['y'])
# print(y)
# print(X[0])
# 从X的行数中随机取得100个行号，进行后续展示, replace=False代表无放回选取
sample_index = np.random.choice(np.arange(X.shape[0]), 100, replace=False)
# 获得100个随机得到的样本的参数
sample = X[sample_index, :]
# print(sample.shape)

# 生成10*10的图片矩阵
# fig, ax = plt.subplots(nrows=10, ncols=10, sharex=True, sharey=True, figsize=(12,12))
# 用样本写入到图片矩阵里的每个小图片
# for i in range(10):
#     for j in range(10):
#         ax[i, j].matshow(sample[10*i+j].reshape(20, 20).T, cmap=plt.cm.binary)
#         plt.xticks(np.array([]))
#         plt.yticks(np.array([]))
# plt.show()
all_theta = one_vs_all(X, y, 10, 1)
predict = predict_all(X, all_theta)

print(classification_report(y, predict))

代码是参考黄博士的答案写的，还有两个自己遇到的问题也进行了查找整理，如下所示。

1. 进行矩阵乘法的两个矩阵行列数不对应

报错信息如下：

ValueError: shapes (5000,401) and (1,401) not aligned: 401 (dim 1) != 1 (dim 0)

这个是因为minimize函数里的x0的类型是np.ndarray且其shape是(n, )。而且其转置的shape也是(n, )，这就会导致X * theta.T得不到我们希望的结果，所以在自己的cost和gradient函数中对theta进行类型转换即可

theta = np.mat(theta)

2. 对theta进行初始化的时候采用np.ones()会产生溢出

报错信息如下：

RuntimeWarning: divide by zero encountered in log

这个产生的原因是在求解代价函数的时候：

second = np.multiply(1 - y, np.log(1 - sigmoid(X * theta.T)))

这里sigmoid(X * theta.T)会得到一个约等于1的值，导致产生np.log(0)的计算，这会产生一个无穷的值，即溢出。解决办法就是给他一个精度就可以了，类似于下面：

second = np.multiply(1 - y, np.log(1 - sigmoid(X * theta.T) + 1e-5 ))