Life's monolog

用python实现逻辑回归

Word count: 1,347 / Reading time: 6 min
2018/04/18 Share

基本形式

逻辑回归主要用于解决分类问题。逻辑回归函数的输出是一个处于0-1之间的值,代表输出是1的概率。
\begin{align}& h_\theta(x) = P(y=1 | x ; \theta) = 1 - P(y=0 | x ; \theta) \newline& P(y = 0 | x;\theta) + P(y = 1 | x ; \theta) = 1\end{align}

基本形式如下:
\begin{align}& h_\theta (x) = g ( \theta^T x ) \newline \newline& z = \theta^T x \newline& g(z) = \dfrac{1}{1 + e^{-z}}\end{align}

可以看到逻辑回归其实就是在线性回归上套了一个$g$函数,也称为Sigmoid函数,形式如下图:
sigmoid

可以看到,当$z$的值大于0的时候,$g(z)$是大于0.5的,当$z$小于0时,$g(z)$是小于0.5的。所以如果假设$h_\theta(x) \geq 0.5$时有$y = 1$(即当输出概率大于等于0.5时,判定为类别1),就等价于$\theta^T x \geq 0 \Rightarrow y = 1$。

模型求解

求解模型首先要了解模型的损失函数,损失函数可以用来衡量模型对数据的拟合程度。逻辑回归采用对数损失函数,如下:

当$y=1$时,即这个样本是正类:

  1. 如果此时$h_\theta(x) = 1$,则对于这个样本而言是分类正确的,此时$cost=-log(1)=0$,表示对于分类正确的样本没有惩罚;
  2. 如果此时$h_\theta(x) = 0$,则对于这个样本而言是分类错误的,此时$cost=-log(0) \to \infty$,表示对于分类错误的样本要给予很大的惩罚。

$y=0$时同理,不再赘述。

可以将上述损失函数简化合并,变成如下形式:

考虑所有的样本,就变成了如下的形式:

逻辑回归模型的求解过程就是最小化上面的损失函数。这里还是采用梯度下降法来求解,还是用损失函数对每个参数$\theta_j$求偏导,求导过程不再涉及,参数$\theta_j$的梯度(偏导)为$\frac{1}{m} \sum_{i=1}^m (h_\theta(x^{(i)}) - y^{(i)}) x_j^{(i)}$,梯度下降算法就是对每个参数进行如下更新:

为了便于实现,这里还是给出损失函数和梯度下降的向量化形式,如下:
\begin{align}
& h = g(X\theta)\newline
& J(\theta) = \frac{1}{m} \cdot \left(-y^{T}\log(h)-(1-y)^{T}\log(1-h)\right)\newline
& \theta := \theta - \frac{\alpha}{m} X^{T} (g(X \theta ) - \vec{y})
\end{align}

模型实现

有了上面的公式,再使用python自带的numpy矩阵运算库,可以很方便地实现逻辑回归模型。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
"""逻辑回归"""
import numpy as np


class LogisticRegression:
def sigmoid(self, z):
"""sigmoid函数"""
g = 1 / (1 + np.exp(-z))
return g

def costFunction(self, X, y):
"""损失函数
Args:
X: 训练数据
y: 标签数据

Returns:
J: 模型当前的损失
grad: 损失函数的梯度
"""
n_samples = len(y)
h = self.sigmoid(X.dot(self.theta))
J = (1 / n_samples) * (-1 * y.T.dot(np.log(h)) - (1 - y).T.dot(np.log(1 - h)))
grad = (1 / n_samples) * X.T.dot(h - y)
return J, grad

def fit(self, X, y, n_iterations=1000, learning_rate=0.01):
"""训练模型
Args:
X: 训练数据
y: 标签数据
n_iterations: 梯度下降算法迭代次数
"""
n_features = X.shape[1]
self.theta = np.zeros([n_features, 1])
for i in range(n_iterations):
J, grad = self.costFunction(X, y)
self.theta -= learning_rate * grad

def predict(self, X):
"""用模型进行预测"""
y_pred = np.round(self.sigmoid(X.dot(self.theta))).astype(int)
return y_pred

通过可视化训练过程中损失函数的变化,可以判断梯度下降算法是否正常工作,如下图所示:
lr_gd
可以看到模型的损失随着迭代次数不断下降。下面使用sklearn自带的数据集iris来粗略地测试一下我们的分类器的效果:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# 使用sklearn自带的测试数据集
iris = datasets.load_iris()
X = iris.data[:100]
y = np.array(iris.target[:100]).reshape(100, 1)
# 分割训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=19)

# 用训练集训练模型,用测试集测试模型
model = LogisticRegression()
model.fit(X_train, y_train)
y_pred = model.predict(X_test)

# 计算模型的各项评分
accuracy = accuracy_score(y_test, y_pred)
precision = precision_score(y_test, y_pred)
recall = recall_score(y_test, y_pred)
print(f'accuracy {accuracy}, precision {precision}, recall {recall}')

可以看到最后模型的accuracy、precision、recall都为1。感兴趣的读者还可以使用其他数据集来测试,完整代码在这里

总结

本文用梯度下降算法实现了逻辑回归模型,由于梯度下降过程中考虑了所有的训练样本,所以是批量梯度下降。当数据量很大时,应该采用随机梯度下降或者小批量梯度下降来加快训练速度。或许读者会看到有的博客或者书中求解逻辑回归模型用的是梯度上升算法,这取决于损失函数的定义。如果原始的问题转化为求某个损失函数的极小值,那么就应该采取梯度下降算法,如本文中所示;如果原始的问题转化为求某个损失函数的最大值,那么就应该用梯度上升算法。梯度上升与梯度下降算法的区别仅仅在于更新参数时是加上梯度还是减去梯度。

CATALOG
  1. 1. 基本形式
  2. 2. 模型求解
  3. 3. 模型实现
  4. 4. 总结