机器学习练习二：用Python实现逻辑回归

时间  2019-12-02
标签 机器 学习 练习 python 实现 逻辑 回归 栏目 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)]} 代码实现：bash

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
复制代码

输出以下：ui

其中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...]。spa

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

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()
复制代码

绘制结果如图：

2.6 探究\lambda超参数大小对决策边界的影响

2.6.1 当\lambda值过大时

如设置当前\lambda=100，观察决策边界。

能够看到此时的结果明显是欠拟合的，这是由于损失函数中 \frac{\lambda }{2m}\sum\limits_{j=1}^{n}{\theta _{j}^{2}}这一部分对总体的影响过大，使得对真正的损失优化不足。

2.6.2 当\lambda值太小时

如设置当前\lambda=0.0001，观察决策边界。

能够看到此时的结果明显是过拟合的，由于 \lambda=0.0001太小至关于并无进行正则化操做，而过多的特征值会致使过拟合。
相关文章
相关标签/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公众号
   欢迎关注本站公众号,获取更多信息
联系我们 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程