机器学习练习二：用Python实现逻辑回归数据介绍：在本部分练习中，您将构建一个逻辑回归模型来预测学生是否被大学录取。

1.逻辑回归

数据介绍：在本部分练习中，您将构建一个逻辑回归模型来预测学生是否被大学录取。假设您是一个大学部门的管理员，您希望根据两次考试的结果来确定每个申请人的录取机会。您可以使用以前申请者的历史数据作为逻辑回归的训练集。

1.1 画出原始数据图

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
# 展示原始数据
data = pd.read_csv('ex2data1.txt',names=['score1','score2','admitted'])
# 布尔值索引data['Admitted'].isin([1])：False,True,False...
positive = data[data['admitted'].isin([1])]
negative = data[data['admitted'].isin([0])]
# 子画布
fig, ax = plt.subplots(figsize=(12,8))
ax.scatter(positive['score1'], positive['score2'], s=50, c='g', marker='o', label='Admitted')
ax.scatter(negative['score1'], negative['score2'], s=50, c='r', marker='x', label='Not Admitted')
ax.legend()  # 标签
ax.set_xlabel('Score1')
ax.set_ylabel('Score2')
plt.show()

其中data['admitted'].isin([1])能够将admitted列中满足值为1的位置变为True，不满足的位置变为False，即生成这样的Series:

然后通过布尔值取data中的数据，生成两个Dataframe：positive、negative。

1.2 数据处理

插入全1列的目的是匹配偏置 $\theta_0$

# 插入全1列x0
data.insert(0, 'Ones', 1)
# set X (training data) and y (target variable)
cols = data.shape[1]
X = data.iloc[:,0:cols-1]  # X是所有行，去掉最后一列
y = data.iloc[:,cols-1:cols]  # X是所有行，最后一列
theta = np.array([0,0,0])

1.3 sigmoid函数

逻辑回归需要使用sigmoid函数来当做激活函数，其中 $h_\theta(x)$ 的值表示样本为1类别的概率。

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

# 画出sigmoid
nums = np.arange(-10,10,step=1)
plt.figure(figsize=(20,8),dpi=100)
plt.plot(nums,sigmoid(nums))
plt.show()

1.4 代价函数和梯度下降

代价函数： $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)
    X = np.matrix(X)
    y = np.matrix(y)
    # X(100,3),y(100,1),theta(1,3)
    first = np.multiply(-y, np.log(sigmoid(X * theta.T)))
    second = np.multiply((1 - y), np.log(1 - sigmoid(X * theta.T)))
    return np.sum(first - second) / (len(X))

梯度下降：

代码实现：

这里梯度下降函数只实现了 $\frac{\partial J\left( \theta \right)}{\partial {{\theta }_{j}}}=\frac{1}{m}\sum\limits_{i=1}^{m}{({{h}_{\theta }}\left( {{x}^{(i)}} \right)-{{y}^{(i)}})x_{_{j}}^{(i)}}$ 部分且为向量化实现。

# 构造梯度下降函数(只做了J(θ)偏导数部分)
def gradientDescent(theta, X, y):
    X = np.matrix(X)
    y = np.matrix(y)
    theta = np.matrix(theta)
    return (1 / len(X)) * X.T @ (sigmoid(X @ theta.T) - y)

1.5 进行梯度下降并计算准确率

1.5.1 梯度下降计算

这里使用了scipy库中的方法，关于该方法参数的意义可以使用help查看。

# 使用 scipy.optimize.minimize 去寻找参数
import scipy.optimize as opt
res = opt.minimize(fun=cost, x0=theta, args=(X, y), method='TNC', jac=gradientDescent)
res

输出如下：

其中fun参数是最终的损失值，x参数是最终的 $\theta$ 值。

1.5.2 准确度计算

获取训练集的预测结果

# 获取训练集预测结果
def predict(theta, X):
    theta = np.mat(theta)
    y_predict = sigmoid(X @ theta.T)
    return [1 if x >= 0.5 else 0 for x in y_predict]

判断准确度

# 判断准确度
y_pre = np.mat(predict(res.x, X))  # 预测值矩阵
y_true = np.mat(y).ravel()  # 真实值矩阵
# 矩阵进行比较返回各元素比对布尔值矩阵，列表进行比较返回整个列表的比对布尔值
accuracy = np.mean(y_pre == y_true)
print('accuracy = {}%'.format(accuracy * 100))

先将y_predict和y_true转换为一维矩阵，一维矩阵的布尔运算返回的是矩阵中各位置元素的布尔值，形如：[True,False,False...]。

而np.mean()求得矩阵中True元素占总元素数量的比例，即预测准确度。

1.6 绘制决策边界

当 $\theta^Tx=0$ 时 $h_\theta(x)=0.5$ ，即所有概率为类别1的样本的集合，亦即决策边界。

显然本次逻辑回归的决策边界是线性的，即直线 $x_2=\frac{-\theta_0-\theta_1x_1}{\theta_2}$ 。

代码实现：

# 绘制决策边界
x1 = np.linspace(data.score1.min(), data.score1.max(), 100)  # 返回在指定范围内的100个等间隔数
x2 = (- res.x[0] - res.x[1] * x1) / res.x[2]
# 布尔值索引data['Admitted'].isin([1])：False,True,False...
positive = data[data['admitted'].isin([1])]
negative = data[data['admitted'].isin([0])]
# 子画布
fig, ax = plt.subplots(figsize=(12,8))
ax.scatter(positive['score1'], positive['score2'], s=50, c='g', marker='o', label='Admitted')
ax.scatter(negative['score1'], negative['score2'], s=50, c='r', marker='x', label='Not Admitted')
ax.plot(x1, x2, 'r', label='Prediction')
ax.legend()  # 标签
ax.set_xlabel('Score1')
ax.set_ylabel('Score2')
plt.show()

其中x1为在x坐标轴间随机取得的100个点，x2即为满足要求的y值。

2.正则化逻辑回归

数据介绍：在练习的这一部分，您将实现规范的逻辑回归来预测来自制造工厂的微芯片是否通过质量保证(QA)。在QA过程中，每个微晶片都要经过各种测试以确保其正常工作。假设您是工厂的产品经理，您在两个不同的测试中获得了一些微芯片的测试结果。从这两个测试中，您希望确定是否应该接受或拒绝微芯片。为了帮助您做出决定，您有一个关于过去的微芯片的测试结果的数据集，您可以从中获得数据。

2.1 数据展示

data2 = pd.read_csv('ex2data2.txt',names=['test1','test2','pass'])
# 布尔值索引data['Admitted'].isin([1])：False,True,False...
positive = data2[data2['pass'].isin([1])]
negative = data2[data2['pass'].isin([0])]
# 子画布
fig, ax = plt.subplots(figsize=(12,8))
ax.scatter(positive['test1'], positive['test2'], s=50, c='g', marker='o', label='Pass')
ax.scatter(negative['test1'], negative['test2'], s=50, c='r', marker='x', label='Not Pass')
ax.legend()  # 标签
ax.set_xlabel('test1')
ax.set_ylabel('test2')
plt.show()

如图所示，我们的数据集不能通过线性方式划分为正例和反例。因此，简单的逻辑回归应用在这个数据集上的效果并不好，因为逻辑回归只能找到一个线性决策边界。更好地适应数据的一种方法是从每个数据点创建更多的特性，我们将把特征映射到所有的多项式项x1和x2的六次方。

2.2 特征映射与数据分割

多项式展开(polynomial expansion)：

伪代码：

for i in 0..power+1：
    for p in 0..i+1:
        output x^(i-p) * y^p

代码实现：

def poly_expansion(x1,x2,power,dataframe):
    """
    多项式展开
    """
    for i in range(power + 1):
        for j in range(i + 1):
            dataframe['F' + str(i-j) + str(j)] = np.power(x1, i-j) * np.power(x2, j)
    dataframe.drop('test1', axis=1, inplace=True)
    dataframe.drop('test2', axis=1, inplace=True)

映射后data2.shape为(118,29)，即从两个特征拓展为28个特征。

数据分割

# 分割提取数据
# set X (training data) and y (target variable)
cols = data2.shape[1]
X = data2.iloc[:,1:cols]  # X是所有行，去掉第一列
y = data2.iloc[:,0:1]  # y是所有行，第一列，注意取值方法
theta = np.zeros(X.shape[1])

2.3 正则化代价函数

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}}

前半部分与普通逻辑回归一致，后半部分是 $\lambda$ 超参数和 $\theta_j$ 的求和的乘积，其目的是让每一个 $\theta_j$ 的值尽量小，避免过拟合。

代码实现：

def cost(theta, X, y, lambd):
    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)))
    front = np.sum(first - second) / (len(X))
    # 函数后半部分
    last = lambd / (2 * len(X)) * np.sum(np.power(theta, 2))
    return front + last

2.4 梯度下降

原理

因为我们未对 ${{\theta }_{0}}$ 进行正则化，所以梯度下降算法将分两种情形：

$\begin{align} & Repeat\text{ }until\text{ }convergence\text{ }\!\!\{\!\!\text{ } \\ & \text{ }{{\theta }_{0}}:={{\theta }_{0}}-a\frac{1}{m}\sum\limits_{i=1}^{m}{[{{h}_{\theta }}\left( {{x}^{(i)}} \right)-{{y}^{(i)}}]x_{_{0}}^{(i)}} \\ & \text{ }{{\theta }_{j}}:={{\theta }_{j}}-a\frac{1}{m}\sum\limits_{i=1}^{m}{[{{h}_{\theta }}\left( {{x}^{(i)}} \right)-{{y}^{(i)}}]x_{j}^{(i)}}-\frac{\lambda }{m}{{\theta }_{j}} \\ & \text{ }\!\!\}\!\!\text{ } \\ & Repeat \\ \end{align}$

对上面的算法中 j=1,2,...,n 时的更新式子进行调整可得：

${{\theta }_{j}}:={{\theta }_{j}}(1-a\frac{\lambda }{m})-a\frac{1}{m}\sum\limits_{i=1}^{m}{({{h}_{\theta }}\left( {{x}^{(i)}} \right)-{{y}^{(i)}})x_{j}^{(i)}}$

在正则化的梯度下降中，我们对除了 $\theta_0$ 以外的参数进行优化，因为 $\theta_0$ 对特征值的权重没有影响（ $x_0=1$ ）。

仍然只对 $\frac{\partial J\left( \theta \right)}{\partial {{\theta }_{j}}}=\frac{1}{m}\sum\limits_{i=1}^{m}{({{h}_{\theta }}\left( {{x}^{(i)}} \right)-{{y}^{(i)}})x_{_{j}}^{(i)}} + \frac{\lambda}{m}\theta_j$ 这一部分进行代码实现。

代码实现：

# 构造梯度下降函数(只做了J(θ)偏导数部分)
def gradient(theta, X, y, lambd):
    X = np.matrix(X)
    y = np.matrix(y)
    theta = np.matrix(theta)
    first = (1 / len(X)) * (X.T @ (sigmoid(X @ theta.T) - y))
    # theta0不需要优化
    second = (lambd / len(X)) * theta[:,1:].T
    second = np.concatenate((np.array([[0]]),second),axis=0)
    return first + second

second = np.concatenate((np.array([[0]]),second),axis=0)这一句代码是为了在不优化 $\theta_0$ 的同时使得前后两部分矩阵对齐（加了一个0行）。

进行梯度下降优化

和上一章方法一样，不再赘述。

# 使用 scipy.optimize.minimize 去寻找参数
import scipy.optimize as opt
res = opt.minimize(fun=cost, x0=theta, args=(X, y, 100), method='TNC', jac=gradient)
res

准确度计算

# 判断准确度
y_pre = np.mat(predict(res.x, X))  # 预测值矩阵
y_true = np.mat(y).ravel()  # 真实值矩阵
# 矩阵进行比较返回各元素比对布尔值矩阵，列表进行比较返回整个列表的比对布尔值
accuracy = np.mean(y_pre == y_true)
print('accuracy = {}%'.format(accuracy * 100))

2.5 决策边界

与上一章不同的是，现在我们将要进行不规则决策边界的绘制。显然上文所使用的的绘制直线方程的方法不再适用。

从数据图中可以看出，红蓝点的区分界限并不是一条直线，而是一个不规则的形状，这就是不规则决策边界。那么绘制不规则决策边界的方法其实也很简单，就是将特征平面上的每一个点都用我们训练出的模型判断它属于哪一类，然后将判断出的分类颜色绘制出来，就得到了上图所示的效果，那么不规则决策边界自然也就出来了，这个原理类似绘制地形图的等高线，在同一等高范围内的点就是同一类。

引自：程序员说

代码实现:

def plot_decision_boundary(axis,axe):
    # meshgrid函数用两个坐标轴上的点在平面上画格，返回坐标矩阵
    X0, X1 = np.meshgrid(
        # 随机两组数，起始值和密度由坐标轴的起始值决定
        np.linspace(axis[0], axis[1], int((axis[1] - axis[0]) * 100)).reshape(-1, 1),
        np.linspace(axis[2], axis[3], int((axis[3] - axis[2]) * 100)).reshape(-1, 1),
    )
    # ravel()方法将高维数组降为一维数组，c_[]将两个数组以列的形式拼接起来，形成矩阵
    X_grid_matrix = np.c_[X0.ravel(), X1.ravel()]
    # 将特征矩阵转换为DataFrame进行多项式展开
    X_grid_matrix = pd.DataFrame(X_grid_matrix,columns=["test1","test2"])
    x1 = X_grid_matrix.test1
    x2 = X_grid_matrix.test2
    poly_expansion(x1,x2,6,X_grid_matrix)
    # 通过训练好的逻辑回归模型，预测平面上这些点的分类
    y_predict = np.array(predict(res.x, X_grid_matrix))
    # 将预测结果转换为与X0和X1形状相同的矩阵
    y_predict_matrix = y_predict.reshape(X0.shape)
    
    # 设置色彩表
    from matplotlib.colors import ListedColormap
    my_colormap = ListedColormap(['#0000CD', '#40E0D0', '#FFFF00'])
    
    # 绘制等高线，并且填充等高区域的颜色
    ax.contourf(X0, X1, y_predict_matrix, cmap=my_colormap)

绘制边界：

# 布尔值索引data['Admitted'].isin([1])：False,True,False...
positive = data2[data2['pass'].isin([1])]
negative = data2[data2['pass'].isin([0])]
# 子画布
fig, ax = plt.subplots(figsize=(12,8))
plot_decision_boundary(axis=[data2.F10.min(), data2.F10.max(), data2.F01.min(), data2.F01.max()], axe=ax)
ax.scatter(positive['F10'], positive['F01'], s=50, c='g', marker='o', label='Pass')
ax.scatter(negative['F10'], negative['F01'], s=50, c='r', marker='x', label='Not Pass')
ax.legend()  # 标签
ax.set_xlabel('test1')
ax.set_ylabel('test2')
plt.show()