逻辑回归

课程与数据集

本文对应 Andrew Ng（吴恩达）《Machine Learning》逻辑回归部分。代码使用的 ex2data1.txt 为课程官方编程练习数据：两名考试分数与是否被录取（0/1）。

数据与参考代码可在黄海广教授整理的中文笔记仓库下载：fengdu78/Coursera-ML-AndrewNg-Notes（见 code/ex2/ 目录）。课程主页：Coursera Machine Learning。

---

一、什么是逻辑回归？

逻辑回归是监督学习的一种。与线性回归的目的是拟合直线或曲线不同，逻辑回归的目的是分类。

拿最经典简单的逻辑回归为例：给你许多苹果的数据，告诉你它们的大小，以及它们分别是普通苹果还是红富士苹果，然后让你判断一个新苹果是哪一种。

• 普通苹果，记为 0
• 红富士苹果，记为 1

逻辑回归和线性回归的区别就显而易见了：逻辑回归的真值往往就是有限个简单的数字。上面的例子只要求分为两类，我们自然而然 label 为 0 和 1（这样标记在后面会看到其实大有用处）。当自变量有两个的时候，我们也不必建立 z 轴显示 0 和 1，而是用颜色表示。现在的任务就是利用已有的数据去预测新的苹果是普通苹果，还是红富士苹果。

说明：苹果例子用于理解"二分类"思想；本笔记代码实际使用的是吴恩达课程中的 ex2data1.txt（两门考试分数 → 是否录取），见下文数据可视化部分。

我们不妨先把数据可视化为散点图，好直观地理解为什么线性回归不能用了，以及数据是否是可分类的。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

path = 'ex2data1.txt'
data = pd.read_csv(path, header=None, names=['exam1', 'exam2', 'admitted'])
print(data.head())

# 接下来我们取最后一列（只有 0 和 1 的那一列），将 0 和 1 的行分别记为 negative, positive
positive = data.loc[data['admitted'].isin([1]), :]
negative = data.loc[data['admitted'].isin([0]), :]  # 用 isin() 方法筛选出真值为 0 和 1 的行，将赋值 bool 值
print(positive.head())
print(negative.head())

这里的 isin 方法经常和 loc、iloc 配合。在列表里输入 1 时，它会给值为 1 的元素赋为 True，不符合条件的记为 False。当把 data['admitted'].isin([1]) print 出来，将会得到：

0    False
1    False
2    False
3     True
4     True
Name: admitted, dtype: bool

loc 会将显示为 True 的行或者列留下。

现在我们绘制散点图：

fig, ax = plt.subplots(figsize=(12, 8))
ax.scatter(positive['exam1'], positive['exam2'], s=50, c='b', marker='o', label='Admitted')
ax.scatter(negative['exam1'], negative['exam2'], s=50, c='r', marker='x', label='Not Admitted')
ax.legend()
ax.set_xlabel('Exam1 Score')
ax.set_ylabel('Exam2 Score')
plt.show()

---

二、逻辑回归的数学原理

2.1 数学模型——sigmoid 函数

对于这种分类算法，依靠线性回归的数学模型显然不能满足我们的需求。下面我们将讨论 sigmoid 函数。学过概统的，现在应该已经知道为什么我们会选择它。下面让我们看看 sigmoid 函数的性质，理解为什么它适合逻辑回归。

sigmoid 函数：

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

这个函数高中生在脑子里就可以绘制出来了。我们也可以用 matplotlib 绘制一下：

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

nums = np.arange(-10, 10, step=0.5)
fig, ax = plt.subplots(figsize=(12, 8))
ax.plot(nums, sigmoid(nums), 'r')
plt.show()

显然在红富士苹果的例子中，苹果越大，越有可能是红富士苹果。我们只需要很简单的 sigmoid 函数，尺寸越大，值越接近 1，尺寸越小，值越接近 0。这样子看来它给出的值确实非常适合评估某个苹果是否是红富士苹果，我们可以把它的函数值当成数学期望，好像确实是很好的评估工具。

但是很显然，ex2data1.txt 里面的数据可不是给个大小让你去看它是红富士苹果的期望大不大。它给出了两个数据，而且并不是越大越可能是 1。这个时候，我们就得对式子里的 $z$ 做文章了。

怎样强化 sigmoid 函数才能使得我们可以去调控多个变量的系数，并且可以用梯度下降使得它可以朝我们需要的方向发展呢？这听起来是线性回归干的活。于是我们做如下处理：

数学模型：假设有 $n$ 个变量构成 $\overrightarrow{x}$，则定义：

$$f(\overrightarrow{x})=g(h(\overrightarrow{x}))$$

其中：$g(z)={\displaystyle \frac{1}{1+e^{-z}}}$，$h(x)=\overrightarrow{\theta }\cdot \overrightarrow{x}$

即：

$$f(\overrightarrow{x})=\frac{1}{1+e^{-\overrightarrow{\theta }\cdot \overrightarrow{x}}}$$

这样就可以利用线性回归拟合的功能，调整系数使之符合我们的方向。

2.2 代价函数

对于一组数据的代价函数我们定义为：

条件	损失 $L(f(\overrightarrow{x^{(i)}}),y^{(i)})$
$y^{(i)}=1$	$-\log (f(\overrightarrow{x^{(i)}}))$
$y^{(i)}=0$	$-\log (1-f(\overrightarrow{x^{(i)}}))$

此处 log 没有底数表示任何底数都能用于误差的评估，后面求导的时候默认底数为自然对数 $e$。可以直接取 $e$ 为底，我们一旦做出图像就会发现，当 $f(\overrightarrow{x^{(i)}})$ 真值为 1 时，越接近 1，代价函数越小，越接近 0，代价函数越大，这符合我们的要求。这个时候，已经有人在写 if-else 语句了，而我们有更好的选择：

$$L(f(\overrightarrow{x^{(i)}}),y^{(i)})=-y^{(i)}\log (f(\overrightarrow{x^{(i)}}))-(1-y^{(i)})\log (1-f(\overrightarrow{x^{(i)}}))$$

把 $y^{(i)}$ 和 $1-y^{(i)}$ 代入可以发现二者等价。倘若有 $m$ 组数据：

$$L(\overrightarrow{\theta })=\frac{1}{m}\sum _{i=1}^m[-y^{(i)}\log (f(\overrightarrow{x^{(i)}}))-(1-y^{(i)})\log (1-f(\overrightarrow{x^{(i)}}))]$$

2.3 代价函数求导

为了实现梯度下降，我们还要知道代价函数对 ${\theta }_j$ 求偏导后的形式（此处为了简便，偏导默认为对第 $j$ 个系数进行）：

$$f(\overrightarrow{x})=\frac{1}{1+e^{-\overrightarrow{\theta }\cdot \overrightarrow{x}}}\to f^{\prime }(\overrightarrow{x})=\frac{x_j{e}^{-\overrightarrow{\theta }\cdot \overrightarrow{x}}}{(1+e^{-\overrightarrow{\theta }\cdot \overrightarrow{x}}{)}^2}$$

注意到：

$$\frac{f^{\prime }(\overrightarrow{x})}{f(\overrightarrow{x})}=x_j(1-f(\overrightarrow{x}))$$$$\frac{-f^{\prime }(\overrightarrow{x})}{1-f(\overrightarrow{x})}=-x_{j}f(\overrightarrow{x})$$

使用链式法则对代价函数求偏导：

$$\frac{\partial L}{\partial {\theta }_j}=-y^{(i)}\frac{f^{\prime }(\overrightarrow{x})}{f(\overrightarrow{x})}-(1-y^{(i)})\frac{-f^{\prime }(\overrightarrow{x})}{1-f(\overrightarrow{x})}$$

代入我们经典的 $f(\overrightarrow{x})$，注意到 $y^{(i)}$ 和 $1-y^{(i)}$：

$$\frac{\partial L}{\partial {\theta }_j}=-y^{(i)}{x}_j(1-f(\overrightarrow{x}))-(1-y^{(i)})(-x_{j}f(\overrightarrow{x}))$$

最终代价函数的和式表达式为：

$$\frac{\partial L}{\partial {\theta }_j}=\frac{1}{m}\sum _{i=1}^m(f(\overrightarrow{x^{(i)}})-y^{(i)})x_j^{(i)}$$

很显然对于常数项：

$$\frac{\partial L}{\partial {\theta }_0}=\frac{1}{m}\sum _{i=1}^m(f(\overrightarrow{x^{(i)}})-y^{(i)})$$

有了函数模型和代价函数及其偏导数，就可以开始写代码了。

---

三、代码编写

3.1 导入库

首先我们导入必要的库：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.optimize as opt

你会发现我们导入了一个新的库，这个库的 opt 功能我们之后会讲解要注意的事项。

3.2 可视化数据

可视化数据，这个开头做过了：

path = 'ex2data1.txt'
data = pd.read_csv(path, header=None, names=['exam1', 'exam2', 'admitted'])
print(data.head())
# 接下来我们取最后一列（只有 0 和 1 的那一列），将 0 和 1 的行分别记为 negative, positive
positive = data.loc[data['admitted'].isin([1]), :]
negative = data.loc[data['admitted'].isin([0]), :]  # 用 isin() 方法筛选出真值为 0 和 1 的行，将赋值 bool 值
print(positive.head())
print(negative.head())

fig, ax = plt.subplots(figsize=(12, 8))
ax.scatter(positive['exam1'], positive['exam2'], s=50, c='b', marker='o', label='Admitted')
ax.scatter(negative['exam1'], negative['exam2'], s=50, c='r', marker='x', label='Not Admitted')
ax.legend()
ax.set_xlabel('Exam1 Score')
ax.set_ylabel('Exam2 Score')
# plt.show()

3.3 数据预处理

数据处理和之前的一样，这里要注意 Y 数组要严格一维数组，否则广播的时候会出错：

data.insert(0, 'Ones', 1)
cols = data.shape[1]
X = data.iloc[:, 0:cols - 1]  # X 是所有行去掉最后一列
Y = data.iloc[:, cols - 1:cols]
X = np.array(X.values)
Y = np.array(Y.values).reshape(-1, 1)
theta_begin = np.zeros(3)

3.4 代价函数

代价函数：参数的位置不能乱放，cost_function 和 gradient 的第一个参数为 theta，这是 opt.fmin_tnc() 的要求。

这里的 theta 需要先做一步 reshape(1, -1) 看起来没什么变化，但是 scipy 的 opt 只能接受一维数组，全局的 theta_begin 必须是 (3,)，在函数里要变成 (1, 3) 防止广播问题。

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

# 代价函数
def cost_function(theta, x, y):
    theta = theta.reshape(1, -1)
    f = sigmoid(x @ theta.T)
    term1 = y * np.log(f)
    term2 = (1 - y) * np.log(1 - f)
    cost = -np.mean(term1 + term2)
    return cost

p = cost_function(theta_begin, X, Y)
print(p)

运行结果：0.6931471805599453

3.5 梯度函数

现在我们只写一个可以计算梯度下降第一步步长的函数，之后的工作交给 scipy。这里建议不要用循环去写，直接 x.T @ error 即可得到我们想要的 sum 数组，然后除以数量即可。（我用循环写出来的数据不一样，问题还没有搞清楚）

def gradient(theta, x, y):
    theta = theta.reshape(-1, 1)  # (3, 1)
    h = sigmoid(x @ theta)  # (m, 1)
    error = h - y.reshape(-1, 1)  # (m, 1)
    grad = (x.T @ error) / len(y)  # (3, 1)
    return grad.ravel()  # (3,)

grad1 = gradient(theta_begin, X, Y)
print(grad1)

运行结果：[-0.1 -12.00921659 -11.26284221]

3.6 使用 scipy 优化

使用 opt.fmin_tnc() 进行优化：

# 使用优化算法求解最优参数
result = opt.fmin_tnc(func=cost_function, x0=theta_begin, fprime=gradient, args=(X, Y))
optimal_theta = result[0]
print(f'Optimal theta: {optimal_theta}')
# 计算最优参数下的代价
final_cost = cost_function(optimal_theta, X, Y)
print(f'Final cost: {final_cost}')

运行结果：

Optimal theta: [-25.16131846  0.20623159  0.20147148]
Final cost: 0.20349770158947528

3.7 预测与准确率

我们接下来看看我们的模型到底准不准，我们把所有的数据重新按得出的预测模型预测一遍，将结果与真值比较：

# 预测函数的值
def predict(theta, x):
    probabilities = sigmoid(x @ theta.T)
    return (probabilities >= 0.5).astype(int)

显然函数 predict 第一个参数是 result[0]，第二个参数就是数据 X：

# 查看最终结果预测效果
theta_min = np.array(result[0])
predictions = predict(theta_min, X)
accuracy = np.mean(predictions == Y.flatten()) * 100
print(f'accuracy = {accuracy}%')

结果的准确率高达：

accuracy = 89.0%

---

下一篇：当线性/logistic 边界不够用时，见《正则化》笔记（数据集 ex2data2.txt，同样来自上述仓库）。