机器学习笔记(六)——Logistic回归

1 模型描述

逻辑斯谛回归（Logistic regression）是机器学习中的一种经典的 分类 算法，虽然名字中带有 回归 ，但其并不是 回归 算法，这是由于历史上命名原因造成的。

1.1 模型假设

给定一个包含两种类别的数据集：

$T = \{(x_1, y_1),(x_2, y_2),\cdots, (x_m, y_m)\}$

在线性可分的二分类问题中，要将两个类别区分开来，一般是寻找一个由 $w$ （权重）和 $b$ （偏置）决定的最优超平面（决策边界）：$w\cdot x + b=0$，两种类别的样本分别处在该超平面的两侧。

在感知机模型中，使用的是符号函数确定样本类别（作为激活函数），将两种样本类别分别记为 $-1,+1$ 两类，对于超平面 $w$ 的样本 $(x_i,y_i)$，认为其为反例样本(-1)，而 $w$ 的样本 $(x_i,y_i)$，认为其为正例样本(+1)。

而在 $Logistic$ 回归中，我们则使用 sigmoid函数确定样本类别（作为激活函数），将两种类别分别记为 $0,1$ 两类，而 sigmoid函数将 $w$ 的结果映射到区间 $(0, 1)$，当 $w$ 时，$sigmoid(w\cdot x_i+b) \in(0,0.5)$ ；而当$w$ 时，$sigmoid(w\cdot x_i+b) \in[0.5,1)$ 。

sigmoid函数映射的数值表示其为类别 1 的概率，越接近 1 ，为类别 1 的概率越大，则判定其为类别 1；反之，结果越接近 0， 则说明为类别1 的概率越小，可判断其为类别 0。则可以认为，当映射数值在 $(0, 0.5)$ 则其为 0 类，而在 $(0.5, 1)$ 则判定其为 1 类。

为方便与课程推导对应，将权值向量 $w$ 替换为 $\theta$ ，则可以得出其 $Logsitic$ 回归的二分类模型的假设函数为：

$h_{\theta}(x_i) = g(\theta^{T}x_i+b)$

参数说明如下：

1. $m$ 和 $n$ 代表样本数量和特征种类数量；
2. $x_i \in R^n$，代表第 $i$ 个样本的特征空间，$x_{i,j}$ 代表第 $i$ 个样本的 第 $j$ 个特征，其中 $i=1,2,\cdots,m.\ j=1,2,\cdots,n$；
3. $y_i\in\{0, 1\}$ 代表第 $i$ 个样本的类别；
4. $g$ 为 $sigmoid$ 函数，即为：

$g(z)=\frac{1}{1 + e^{-z}}$

5. $\theta\in R^n$ 为超平面的法向量，同时也为 $Logistic$ 回归模型的参数，代表权值向量(weight vector)，$\theta^{T}$ 代表 $\theta$ 的转置，$\theta_i,i=1,2,\cdots,n$ 代表第 $i$ 个权值参数的值；
6. $b\in R$ 也为超平面的截距，同时也为 $Logistic$ 回归模型的参数，代表偏置(bias)。

为了方便简化计算，我们给每个样本的特征加上 $x_{i,0} = 1$ 这个特征，让偏置 $b$ 加到权值向量 $\theta$ 中成为 $\theta_0$ ，让假设函数变为如下简化形式：

$\begin{aligned} h_\theta(x_i)&=g(\theta_0x_{i,0}+\theta_1x_{i,1}+\theta_2x_{i,2}+\cdots+\theta_nx_{i,n})\\ &=g(\theta^Tx_i)\\ &=\frac{1}{1+e^{-\theta^Tx_i}} \end{aligned}$

1.2 代价函数

在线性回归中，我们选取的是平方误差函数作为代价函数，但若应用在 $Logistic$ 回归的假设函数中作为代价函数，可以发现其并不是一个凸函数，存在局部最优解，在使用梯度下降的过程中可能会在局部最小值处停止，而无法达到全局最小值。因此不能使用其作为代价函数。

合理的 $Logistic$ 回归的代价函数如下：

$J(\theta) = \frac{1}{m}\sum_{i=1}^{m}Cost(h_\theta(x_i), y_i)$

$Cost(x_i,y_i) = \begin{cases} -log(h_\theta(x_i)), &y_i=1\\ -log(1-h_\theta(x_i)), &y_i=0 \end{cases}$

将上式简化得到：

$\begin{aligned} J(\theta) &= \frac{1}{m}\sum_{i=1}^{m}(-y_i\log(h_\theta(x_i))-(1-y_i)\log(1-h_\theta(x_i)))\\ &= -\frac{1}{m}\sum_{i=1}^{m}(y_i\log(h_\theta(x_i))+(1-y_i)\log(1-h_\theta(x_i))) \end{aligned}$

1.3 梯度计算

在 1.2 节 中我们已经我们已经得到了 $Logistic$ 回归模型的代价函数，我们的目标时找到一组最优的参数 $\theta$ ，使代价函数达到最小化，即：

$arg\,\min_{\theta}J(\theta)$

则我们使用梯度下降进行寻优，其在第 $j$ 个参数的梯度为：

$\begin{aligned} \nabla_{\theta_j} &= \frac{\partial J(\theta)}{\partial \theta_j}\\ &=\frac{1}{m}\sum_{i=1}^{m}(h_\theta(x_i) - y_i)x_{i,j} \end{aligned}$

则每次更新梯度的公式为：

$\theta_j := \theta_j - \frac{\alpha}{m}\sum_{i=1}^{m}(h_\theta(x_i) - y_i)x_{i,j}$

其中 $\alpha$ 为学习率。

2 案例实现

采用的数据集为吴恩达课后作业 ex2 中的 ex2data1.txt，下载地址。

数据有两种特征，分别为课程1的分数和课程2的分数，类别有两种，类别0 代表被录取，类别1 代表未被录取。

需要我们判断两门课程分别为多少时，其能被录取。

2.1 代码实现

import numpy as np
import matplotlib.pyplot as plt
import time


class Model:
    """
    二分类逻辑斯特回归
    使用梯度下降法
    """

    def sigmoid(self, z):
        """
        假设函数 h
        :param x: 传入特征
        :return: 特征 x 在此假设函数下的映射
        """
        return 1. / (1 + np.exp(-z))

    def __J(self):
        """
        代价函数
        :return:
        """
        eps = 1e-6
        hx = self.sigmoid(self.x.dot(self.theta))  # h_theta(x) 的值
        J = -(np.inner(np.log(hx + eps), self.y) + np.inner(np.log(1 - hx + eps), 1 - self.y)) / self.m
        return J

    def __init__(self, alpha=0.01, eps=1e-5, iters=10000):
        """
        构造函数
        :param alpha: 学习步长
        :param eps: 精度
        """
        self.alpha = alpha
        self.eps = eps
        self.iters = iters

    def fit(self, x, y):
        """
        进行训练
        :param x: 特征向量
        :param y: 标签
        """
        self.m, self.n = x.shape[0], x.shape[1] + 1
        self.y = y
        c = np.ones(shape=(self.m, 1))
        self.x = np.concatenate((c, x), axis=1)

        # 初始化参数 theta 为0
        self.theta = np.ones(self.n)

        self.j_list = [self.__J()]
        for i in range(self.iters): # 进行迭代
            hx = self.sigmoid(self.x.dot(self.theta)).flatten()
            sub = hx - self.y
            self.theta -= self.alpha * (np.dot(self.x.T, sub).flatten()) / self.m
            tj = self.__J()
            self.j_list.append(tj)
            if abs(tj - 0) <= 1e-6:
                print("Well done in step %i! ", i)
                return
        print("Not completed!")

    def get_param(self):
        """
        获取拟合最优的参数
        :return: theta
        """
        return self.theta

    def predict(self, x):
        """
        预测
        :param x: 要预测的特征
        :return: 预测的结果
        """
        c = np.ones(shape=(x.shape[0], 1))
        x = np.concatenate((c, x), axis=1)  # 将特征的第一维添加1，方便计算
        return self.sigmoid(x).flatten()

    def plot_cost_change(self):
        """绘制损失函数的变化趋势"""
        plt.figure()
        plt.plot(self.j_list, linewidth=".4")
        plt.xlabel("iter")
        plt.ylabel("cost")
        plt.title("代价函数值随迭代次数增加的变化趋势")
        plt.show()


# 导入数据集
data = np.loadtxt("./AndrewNg/ex2/ex2data1.txt", delimiter=',')
# print(data)

# 数据提取
X = data[:, :2]  # 取前100样本的前两种特征
y = data[:, 2]


model = Model(alpha=0.003, iters=2000000)

st = time.time()
model.fit(X, y)
ed = time.time()
print("耗时：%.3fs" % (ed - st))

model.plot_cost_change()

theta = model.theta
print(theta)
xx = np.arange(X[:, 0].min(), X[:, 0].max(), 1)
yy = -(theta[1] * xx + theta[0]) / theta[2]

# 绘制结果
plt.figure(2)
plt.plot(xx, yy, c='blue')
plt.scatter(X[y == 0, 0], X[y == 0, 1], label='Not admitted', c="none", marker='o', edgecolors='r')
plt.scatter(X[y == 1, 0], X[y == 1, 1], label='Admitted', c="black", marker='+')
plt.legend()
plt.xlabel("Exam 1 score")
plt.ylabel("Exam 2 score")
plt.title("两种类别的划分结果")
plt.show()

2.2 结果输出

Not completed!
耗时：56.725s
[-24.05743837   0.19740346   0.19253832]

其代价函数的变化为：

其超平面的分割图像如下：

可以看到，超平面大致能分割出两种类别，但却不是最优。使用scipy.optimize.minimize 优化函数进行寻优的话运算速度极快，但也只能同梯度下降的划分结果相近 $\theta = [-25.16130062, 0.20623142, 0.20147143]$，与梯度下降的结果差别不大。

想要更合理的划分两种类别，可以考虑添加高次项的特征如 $x_1^2, x_2^2, x_1x_2,\cdots$ ，再引入 正则化项 进行梯度下降找到更好的超平面。