吴恩达《机器学习》课后测试Ex2：逻辑回归（详细Python代码注解）

最新推荐文章于 2022-01-11 23:24:15 发布

Mengo_x

最新推荐文章于 2022-01-11 23:24:15 发布

阅读量1k

点赞数 2

分类专栏：机器学习 Python与AI 文章标签： python 机器学习人工智能

本文链接：https://blog.csdn.net/qq_44577070/article/details/120644061

版权

Python与AI 同时被 2 个专栏收录

20 篇文章 0 订阅

订阅专栏

机器学习

16 篇文章 5 订阅

订阅专栏

基于吴恩达《机器学习》课程
参考黄海广的笔记

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.optimize as opt    # 用于高级优化

Part1

在第一部分，我们将要构建一个逻辑回归模型来预测，某个学生是否被大学录取。

设想你是大学相关部分的管理者，想通过申请学生两次测试的评分，来决定他们是否被录取。现在你拥有之前申请学生的可以用于训练逻辑回归的训练样本集。对于每一个训练样本，你有他们两次测试的评分和最后是被录取的结果。为了完成这个预测任务，我们准备构建一个可以基于两次测试评分来评估录取可能性的分类模型。

sigmoid 函数

$h_\theta \left( x \right)=\frac{1}{1+{{e}^{-\theta^{T}X}}}$

# 实现sigmoid函数
def sigmoid(z):
    return 1 / (1 + np.exp(-z))

做一个快速的检查，来确保它可以工作。

# 仅用于测试sigmoid函数运行
def sig_test():
    nums = np.arange(-10, 10, step=1)  # 以-10为起点，10为终点，步长为1的列表
    fig, ax = plt.subplots(figsize=(12, 8))  # 以其他关键字参数**fig_kw来创建图
    ax.plot(nums, sigmoid(nums), 'r')
    plt.show()

代价函数

$J\left( \theta \right)=-\frac{1}{m}\sum\limits_{i=1}^{m}{[{{y}^{(i)}}\log \left( {h_\theta}\left( {{x}^{(i)}} \right) \right)+\left( 1-{{y}^{(i)}} \right)\log \left( 1-{h_\theta}\left( {{x}^{(i)}} \right) \right)]}$

def cost(theta, X, y):
    theta = np.matrix(theta)  # 将theta转换为矩阵
    X = np.matrix(X)  # 将X转换为矩阵
    y = np.matrix(y)  # 将y转换为矩阵
    first = np.multiply(-y, np.log(sigmoid(X * theta.T)))  # 公式第一项  np.multiply将数组或矩阵对应位置相乘。
    second = np.multiply((1 - y), np.log(1 - sigmoid(X * theta.T)))  # 公式第二项
    return np.sum(first - second) / (len(X))

first和second都是一个矩阵。

计算梯度 $\frac{\partial }{\partial {\theta_j}}J\left( \theta \right)$

注意，我们实际上没有在这个函数中执行梯度下降，我们仅仅在计算一个梯度步长。梯度下降可以参考Ex1中的代码。

这次，我们不自己写代码实现梯度下降，我们会调用一个已有的库。这就是说，我们不用自己定义迭代次数和步长，功能会直接告诉我们最优解。Andrew NG在课程中用的是Octave的“fminunc”函数，由于我们使用Python，我们可以用scipy.optimize.fmin_tnc做同样的事情。

# 实现梯度计算（并没有更新θ），没有梯度下降
def gradient(theta, X, y):
    theta = np.matrix(theta)  # 将theta转换为矩阵
    X = np.matrix(X)  # 将X转换为矩阵
    y = np.matrix(y)  # 将y转换为矩阵
    temp = np.matrix(np.zeros(theta.shape))  # np.zeros(theta.shape)=[0.,0.],然后将temp变为矩阵[0.,0.]
    parameters = int(theta.ravel().shape[1])
    # theta.ravel()：将多维数组theta降为一维，.shape[1]是统计这个一维数组有多少个元素。parameters表示参数
    grad = np.zeros(parameters)
    error = sigmoid(X * theta.T) - y
    for i in range(parameters):
        term = np.multiply(error, X[:, i])  # 将误差与训练数据相乘，得到偏导数
        grad[i] = np.sum(term) / len(X)
    return grad

这里的error就是 $\delta=a-y=\frac{\partial C}{\partial {z}}$ ， $a = g^{'} (z)$ ，X[:, j]则是 $\frac{\partial z}{\partial {w}}=a$ ，所以term为 $\frac{\partial C}{\partial {w}}=\frac{\partial z}{\partial {w}}\frac{\partial C}{\partial {z}}$ ，即为代价函数偏导数。可参考笔记反向传播详解。

评价逻辑回归模型的预测

predictA：预测学生是否录取。如果一个学生exam1得分45，exam2得分85，那么他录取的概率应为0.776

predictB：看模型在训练集上的正确率怎样。给出数据以及参数后，predictB的函数会返回“1”或者“0”。然后再把这个predict函数用于训练集上，看准确率怎样。

def predictA(theta, X):
    return sigmoid(X @ theta)

def predictB(theta, X):
    probability = sigmoid(X * theta.T)
    return [1 if x >= 0.5 else 0 for x in probability]  
    #大于0.5的数取值为1，小于0.5的为0

这里注意 * 和@的区别，* 星乘表示矩阵内各对应位置相乘，@等价于.dot()点乘表示求矩阵内积。

逻辑回归代码

此处代码是连续的，中途输出和一些说明单独拿出来了。

def ex2_1():
    # 从检查数据开始
    data1 = pd.read_csv('ex2data1.txt', names=['text1', 'text2', 'admitted'])
    print(data1.head())
    print(data1.describe())    # 查看数据的基本情况，包括:count非空值数、mean平均值、std 标准差、max、min、(25%、50%、75%)分位数
    
    # 让我们创建两个分数的散点图，并使用颜色编码来可视化，正的（被接纳）或负的（未被接纳）
    positive = data1[data1['admitted'].isin([1])]  # isin()为筛选函数，令positive为数据集中为1的数组
    negative = data1[data1['admitted'].isin([0])]  # isin（）为筛选函数，令positive为数据集中为0的数组
    fig, ax = plt.subplots(figsize=(12, 8))  # 以其他关键字参数**fig_kw来创建图
    # figsize=(a,b):figsize 设置图形的大小,b为图形的宽,b为图形的高,单位为英寸
    ax.scatter(positive['text1'], positive['text2'], s=50, c='b', marker='o', label='admitted')
    # x为Exam 1，y为Exam 2，散点大小s设置为50,颜色参数c为蓝色，散点形状参数marker为圈，以关键字参数**kwargs来设置标记参数labele是Admitted
    ax.scatter(negative['text1'], negative['text2'], s=50, c='r', marker='x', label='not admitted')
    # x为Exam 1，y为Exam 2，散点大小s设置为50,颜色参数c为红色，散点形状参数marker为叉，以关键字参数**kwargs来设置标记参数labele是Not Admitted
    ax.legend()  # legend为显示图例函数
    ax.set_xlabel('text 1 Score')  # 设置x轴变量
    ax.set_ylabel('text 2 Score')  # 设置y轴变量
    plt.show()

此处输出如下：

       text1      text2  admitted
0  34.623660  78.024693         0
1  30.286711  43.894998         0
2  35.847409  72.902198         0
3  60.182599  86.308552         1
4  79.032736  75.344376         1

            text1       text2    admitted
count  100.000000  100.000000  100.000000
mean    65.644274   66.221998    0.600000
std     19.458222   18.582783    0.492366
min     30.058822   30.603263    0.000000
25%     50.919511   48.179205    0.000000
50%     67.032988   67.682381    1.000000
75%     80.212529   79.360605    1.000000
max     99.827858   98.869436    1.000000

在这里插入图片描述

# 代码接上
    # 初始化X，y，theta
    data1.insert(0, 'Ones', 1)  # 在第0列插入表头为“ONE”的列，数值为1
    cols = data1.shape[1]  # 获取表格df的列数
    X = data1.iloc[:, 0:cols - 1]  # 除最后一列外，取其他列的所有行，即X为O和人口组成的列表
    y = data1.iloc[:, cols - 1:cols]  # 取最后一列的所有行，即y为利润
    X = np.array(X.values)  # 将X的各个值转换为数组
    y = np.array(y.values)  # 将y的各个值转换为数组
    theta = np.zeros(3)  # 初始化theta数组为(0,0,0)
    # 检查矩阵的维度
    print(X.shape, theta.shape, y.shape)
    # 用初始θ计算代价
    print(cost(theta, X, y))
    # 看看梯度计算结果
    print(gradient(theta, X, y))

输出如下：

(100, 3) (3,) (100, 1)
0.6931471805599453
[ -0.1        -12.00921659 -11.26284221]

# 代码接上
    # 用SciPy's truncated newton（TNC）实现寻找最优参数。
    result = opt.fmin_tnc(func=cost, x0=theta, fprime=gradient, args=(X, y))
    print(result)
    # 看看在这个结论下代价函数计算结果是什么个样子
    print(cost(result[0], X, y))  # result[0]:取数组result中的第一个元素

result返回值说明：数组：返回优化问题目标值；nfeval整数：function evaluations的数目；rc。

在进行优化的时候,每当目标优化函数被调用一次,就算一个function evaluation。在一次迭代过程中会有多次function evaluation，这个参数不等同于迭代次数，而往往大于迭代次数。

输出结果如下：

(array([-25.16131863,   0.20623159,   0.20147149]), 36, 0)
0.20349770158947458

# 代码接上
    # 画出决策曲线
    plotting_x1 = np.linspace(30, 100, 100)
    plotting_h1 = (- result[0][0] - result[0][1] * plotting_x1) / result[0][2]
    # theta^T*x=0为判定边界，x0theta0+x1theta1+x2theta2=0可得：x2=(-x0theta0-x1theta1)/theta2
    fig, ax = plt.subplots(figsize=(12, 8))
    ax.plot(plotting_x1, plotting_h1, 'y', label='Prediction')
    ax.scatter(positive['text1'], positive['text2'], s=50, c='b', marker='o', label='Admitted')
    ax.scatter(negative['text1'], negative['text2'], s=50, c='r', marker='x', label='Not Admitted')
    ax.legend()
    ax.set_xlabel('Exam 1 Score')
    ax.set_ylabel('Exam 2 Score')
    plt.show()

${\theta^{T}}x=0$ 这条线便是判定边界(decision boundary)。

输出图像：
在这里插入图片描述

# 代码接上
    # predictA预测学生是否录取。
    print(predictA(result[0], [1, 45, 85]))
    
    # predictB预测准确率。
    theta_min = np.matrix(result[0])  # 将最优化最优参数中优化问题目标值数组x转化为矩阵
    predictions = predictB(theta_min, X)
    correct = [1 if ((a == 1 and b == 1) or (a == 0 and b == 0)) else 0 for (a, b) in zip(predictions, y)]
    # predictions为预测值,y为实际值.zip为打包对象函数，将(predictions, y)打包成元组，并赋值为（a,b），然后判断正确的则返回1，判断错误的返回0
    accuracy = (sum(map(int, correct)) % len(correct))
    # map(int, correct):将crrect列表内容的类型映射成int型，原来应该布尔型.sum为求和函数，%为求模运算，计算除法的余数。
    # 实际是求1占数据总数的比例
    print('accuracy = {0}%'.format(accuracy))  # format为格式化字符串的函数,{0}为设置指定位置.

predictB这是训练集的准确性。我们没有保持住了设置或使用交叉验证得到的真实逼近，所以这个数字有可能高于其真实值。

输出结果如下：

0.7762906240463825
accuracy = 89%

Part2

在第二部分，我们将实现加入正则项提升逻辑回归算法。
设想你是工厂的生产主管，你有一些芯片在两次测试中的测试结果，测试结果决定是否芯片要被接受或抛弃。你有一些历史数据，帮助你构建一个逻辑回归模型。

正则化的代价函数

$J\left( \theta \right)=\frac{1}{m}\sum\limits_{i=1}^{m}{[-{{y}^{(i)}}\log \left( {h_\theta}\left( {{x}^{(i)}} \right) \right)-\left( 1-{{y}^{(i)}} \right)\log \left( 1-{h_\theta}\left( {{x}^{(i)}} \right) \right)]}+\frac{\lambda }{2m}\sum\limits_{j=1}^{n}{\theta _{j}^{2}}$

# 正则化的代价函数
def costreg(theta, X, y, learningRate):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    first = np.multiply(-y, np.log(sigmoid(X * theta.T)))
    second = np.multiply((1 - y), np.log(1 - sigmoid(X * theta.T)))
    reg = (learningRate / (2 * len(X)) * np.sum(np.power(theta[:, 1:theta.shape[1]], 2)))
    return np.sum(first - second) / len(X) + reg

learningRate就是正则化参数 $\lambda$ 。

正则化的梯度计算

我们未对 $\theta_0$ 进行正则化，所以梯度下降算法应分两种情形：

$R e p e a t$ $u n t i l$ $c o n v e r g e n c e$ {

${\theta_0}:={\theta_0}-a\frac{1}{m}\sum\limits_{i=1}^{m}{(({h_\theta}({{x}^{(i)}})-{{y}^{(i)}})x_{0}^{(i)}})$

${\theta_j}:={\theta_j}-a[\frac{1}{m}\sum\limits_{i=1}^{m}{(({h_\theta}({{x}^{(i)}})-{{y}^{(i)}})x_{j}^{\left( i \right)}}+\frac{\lambda }{m}{\theta_j}]$ （ $f o r$ $j = 1, 2, . . . n$ ）

}

# 正则化的梯度计算，没有梯度下降
def gradientReg(theta, X, y, learningRate):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    parameters = int(theta.ravel().shape[1])
    #theta.ravel():将多维数组转换为一维数组，并不会产生源数据的副本. shape[1]:获取列数
    grad = np.zeros(parameters) #创建零数组
    error = sigmoid(X * theta.T) - y
    for i in range(parameters):
        term = np.multiply(error, X[:,i])  # term即偏导数
        if (i == 0):  # 不对第0个参数正则化
            grad[i] = np.sum(term) / len(X)
        else:
            grad[i] = (np.sum(term) / len(X)) + ((learningRate / len(X)) * theta[:,i])
    return grad

正则化的逻辑回归代码

此处代码是连续的，中途输出和一些说明单独拿出来了。

def ex2_2():
    # 第一部分很像，从数据可视化开始
    data2 = pd.read_csv('ex2data2.txt', names=['text1', 'text2', 'accepted'])
    print(data2.head())
    positive = data2[data2['accepted'].isin([1])]  # isin（）为筛选函数，令positive为数据集中为1的数组
    negative = data2[data2['accepted'].isin([0])]  # isin（）为筛选函数，令positive为数据集中为0的数组
    fig, ax = plt.subplots(figsize=(12, 8))  # 以其他关键字参数**fig_kw来创建图
    ax.scatter(positive['text1'], positive['text2'], s=50, c='b', marker='o', label='accepted')
    # x为Exam 1，y为Exam 2，散点大小s设置为50,颜色参数c为蓝色，散点形状参数marker为圈，以关键字参数**kwargs来设置标记参数labele是Admitted
    ax.scatter(negative['text1'], negative['text2'], s=50, c='r', marker='x', label='not accepted')
    # x为Exam 1，y为Exam 2，散点大小s设置为50,颜色参数c为红色，散点形状参数marker为叉，以关键字参数**kwargs来设置标记参数labele是Not Admitted
    ax.legend()  # 显示图例
    ax.set_xlabel('text1 Score')  # 设置x轴变量
    ax.set_ylabel('text2 Score')  # 设置y轴变量
    plt.show()

输出如下：

      text1    text2  accepted
0  0.051267  0.69956         1
1 -0.092742  0.68494         1
2 -0.213710  0.69225         1
3 -0.375000  0.50219         1
4 -0.513250  0.46564         1

在这里插入图片描述

以上图片显示，这个数据集不能像之前一样使用直线将两部分分割。而逻辑回归只适用于线性的分割，所以，这个数据集不适合直接使用逻辑回归。

一种更好的使用数据集的方式是为每组数据创造更多的特征（特征映射）：

# 代码接上
    degree = 5
    x1 = data2['text1']
    x2 = data2['text2']
    data2.insert(3, 'Ones', 1)  # 在第4列插入标题为One,各行值为1的列
    # 这里是把输入x1,x2重新映射成5-1=4阶了，x0=1,x1,x1^2,x1x2,x1^3,x1^2x2,x1x2^2...
    for i in range(1, degree):
        for j in range(0, i):
            data2['F' + str(i) + str(j)] = np.power(x1, i - j) * np.power(x2, j)
            # np.power(A,B)   ## 对A中的每个元素求B次方
    data2.drop('text1', axis=1, inplace=True)  # 删除表头为Test 1的列，axis=1表示默认删除行或列，inplace=True表示原数组被data2替换.
    data2.drop('text2', axis=1, inplace=True)  # 删除表头为Test 2的列，axis=1表示默认删除行或列，inplace=True表示原数组被data2替换.
    print(data2.head())

输出如下：

   accepted  Ones       F10       F20  ...       F40       F41       F42       F43
0         1     1  0.051267  0.002628  ...  0.000007  0.000094  0.001286  0.017551
1         1     1 -0.092742  0.008601  ...  0.000074 -0.000546  0.004035 -0.029801
2         1     1 -0.213710  0.045672  ...  0.002086 -0.006757  0.021886 -0.070895
3         1     1 -0.375000  0.140625  ...  0.019775 -0.026483  0.035465 -0.047494
4         1     1 -0.513250  0.263426  ...  0.069393 -0.062956  0.057116 -0.051818
[5 rows x 12 columns]

接下来就像在第一部分中做的一样，初始化变量。

# 代码接上
    cols = data2.shape[1]  # 获取data2的列数
    X2 = data2.iloc[:, 1:cols]  # 取所有列每一行的数据
    y2 = data2.iloc[:, 0:1]  # 取第一和第二列每一行的数据
    X2 = np.array(X2.values)  # 将数值转换为数组
    y2 = np.array(y2.values)
    theta2 = np.zeros(11)  # 创建长度为11的零数组
    # 正则化参数取合理值，我们可以之后再折腾这个。
    learningRate = 1
    print(costreg(theta2, X2, y2, learningRate))  # 用初始参数计算代价
    print(gradientReg(theta2, X2, y2, learningRate))  # 看看梯度计算结果

输出如下：

0.6931471805599454
[0.00847458 0.01878809 0.05034464 0.01150133 0.01835599 0.00732393
 0.00819244 0.03934862 0.00223924 0.01286005 0.00309594]

# 代码接上
# 使用和第一部分相同的优化函数来计算优化后的结果。
    result2 = opt.fmin_tnc(func=costreg, x0=theta2, fprime=gradientReg, args=(X2, y2, learningRate))
    print(result2)

# 使用第1部分中的预测函数predictB来查看我们的方案在训练数据上的准确度。
    theta_min = np.matrix(result2[0])
    predictions = predictB(theta_min, X2)
    correct = [1 if ((a == 1 and b == 1) or (a == 0 and b == 0)) else 0 for (a, b) in zip(predictions, y2)]
    # predictions为预测值,y为实际值.zip为打包对象函数，将(predictions, y)打包成元组，并赋值为（a,b），然后判断正确的则返回1，判断错误的返回0
    accuracy = (sum(map(int, correct)) % len(correct))
    # map(int, correct):将crrect列表内容的类型映射成int型，原来应该布尔型.sum为求和函数，%为求模运算，计算除法的余数。
    # 实际是求1占数据总数的比例
    print('accuracy = {0}%'.format(accuracy))

输出如下：

(array([ 0.53010247,  0.29075567, -1.60725764, -0.58213819,  0.01781027,
       -0.21329508, -0.40024142, -1.3714414 ,  0.02264304, -0.9503358 ,
        0.0344085 ]), 22, 1)
accuracy = 78%