Re:从零开始的机器学习 - Machine Learning(二) 逻辑回归LR

从我对整个职业生涯的规划出发，我不仅想做一些高质量的应用(软件工程的角度),还想做一些激动人心的应用，所以我希望能在机器学习的方向走，尽管我在大学粗浅的学了些皮毛，但如果要把机器学习作为职业发展的话这些还差得远，所以我开始写了这个系列的文章。

我希望通过这个系列能对机器学习领域的知识点做一个总结，所以对于Machine Learning这个部分，我的目标是写出能让高中生看得懂。

前言

继上一章 Re:从零开始的机器学习 - Machine Learning(一) 线性回归，我还是继续按照斯坦福的路线，这一章来说说逻辑回归（Logistic Regression）

分类(Classification)

和回归(Regression)一样，分类(Classification)问题也是机器学习里面很大的一块。

分类问题是机器学习非常重要的一个组成部分，它的目标是根据已知样本的某些特征，判断一个新的样本属于哪种已知的样本类。

其实常见的例子很多，判断一个邮件是否是垃圾邮件之类的，预测一个用户是否对我的商品感兴趣，以及图像处理里面对图像进行的分类。

分类问题有简单的二分类也有多分类。

正文

这篇博客主要讲的是逻辑回归(Logistic regression)。

逻辑回归LR(Logistic Regression)

看到名字的时候你可能会有一些奇怪，为什么明明叫逻辑“回归”却用在分类问题上。虽然这个名字似乎指示着什么，但是逻辑回归实际上是分类算法。我猜它之所以这样命名是因为在它的学习方法上和线性回归相似，但是损失(loss)和梯度(gradient)函数的表达不同。特别的，逻辑回归使用 S型函数(sigmoid)而不是线性回归的连续输出。当我们深入到实现中去，我们会了解到更多。

首先我们先把逻辑回归放到一边，之前也说了逻辑回归是用来解决分类问题的，对于分类问题，我们实际上是希望得到一个分类器(Classifier)，当输入数据之后，这个分类器能给我预测这个数据属于某一类的概率，也就是说我们需要的是一个概率。

上一节我们介绍的线性回归，其输出的是预测值，其假设函数(Hypothesis Function)也就是输出预测值的函数是这样的。

而逻辑回归则是预测属于某一类的 概率，所以我们让其假设函数是下面这个

这个函数的意义其实是 当输入为x时输出y=1的概率，其实就是二分类问题里面是某个东西的概率。读者现在可能会对这个函数有所疑问，比如为什么是这个函数，这个留在后面会讨论。

这里出现了条件概率，实际上就是指事件A在另外一个事件B已经发生条件下的发生概率。高中生可以补这部分也可以暂时不补

我们把这个函数g(z)叫做sigmoid函数，很明显这个函数的值域是0到1的开区间。接下来我们会详细介绍一下这个函数。

Sigmoid

Sigmoid函数的函数表达式如下

sigmoid表达式

sigmoid函数图像

为什么是Sigmoid

损失函数(Loss Function)

上一小节我们也说过，为了修正参数Θ我们需要有个手段来衡量当前参数Θ的优秀程度。损失函数(Loss Function)就是用来衡量假设函数(hypothesis function)的准确性。

对于逻辑回归来说，我们希望的是当预测概率约接近实际情况(0或1)的时候误差最小，而且不希望曲线是一条直线，而是对于越接近的地方变化越小，约远离的地方变化越大的函数。

下面就是逻辑回归的损失函数。

逻辑回归的损失函数

我们可以将函数合并一下，毕竟这种分段函数处理起来不是很舒服，其实就是下图这样，也很好理解，毕竟二分类训练数据y只有0和1两个值。

这样我们就可以算出在一个训练集中基于当前参数Θ得到结果的误差了。

矢量化编程

其实上一小节的实战中我们已经用到了矢量化编程。

# 计算损失，用了矢量化编程而不是for循环
def computeLoss(X, y, theta):  
    inner = np.power(((X * theta.T) - y), 2)
    return np.sum(inner) / (2 * len(X))
复制代码

矢量化编程是提高算法速度的一种有效方法。为了提升特定数值运算操作（如矩阵相乘、矩阵相加、矩阵-向量乘法等）的速度，数值计算和并行计算的研究人员已经努力了几十年。矢量化编程的思想就是尽量使用这些被高度优化的数值运算操作来实现我们的学习算法。

换句话说就是尽量避免使用for循环，毕竟矩阵相乘这种场景非常适合并行计算，在巨量的数据面前性能收益非常明显。

如果刚刚的损失函数用矢量化编程的思想来表示的话

如果一时不理解的话我先解释一下，我们先假设共m个数据，而这个模型中变量有n个。则矩阵h就是(m, n) X (n, 1)也就是(m,1)矩阵，矩阵h的意义就是这m个数据的预测值。

损失函数中y的转置为(1, m)，相乘后得到(1, 1)也就是一个值，这两个矩阵相乘的意义则是对应的预测值取对数乘以对应的实际值，最后加在一起。

(m,n)表示维度为m行n列的矩阵，如果学过矩阵的乘法应该知道矩阵相乘(m, n) X (n, k)得到的矩阵是(m, k)

凹函数

逻辑回归的梯度下降法(Gradient Descent)

关于梯度下降法我在上一小节里面已经做了比较具体的描述，如果忘记了可以回去翻翻。我们刚刚知道了怎么评价当前参数Θ的好坏，现在我们需要做的是使用梯度下降法来调整参数。

依旧是对损失函数求偏导数，别忘记α是学习速率的意思。

矢量化表示为

损失函数偏导数求解过程(选修)

过拟合问题

对于一个训练数据集，可视化后如下图所示。

对于三个不同的分类器划分出的边界的三种情况，我们对其有不同的称呼

第一种，分类非常不准，这种我们叫 欠拟合(underfitting)

第二种，分类得恰到好处，这种其实没有特别的称呼。

第三种，分类太过于契合训练数据了。这种我们称为过拟合(overfitting)

过拟合所产生的问题也很明显，它实在太过于契合训练集了，对于我们来说，第二个曲线才是我们想要的，过拟合的结果太过于契合训练数据，实用性可想而知的低。

而解决过拟合的方法主要有两种

减少特征的数量，这个很好理解，更多的特征意味着划分出来的函数曲线可以越复杂。这个可以扩展到以后会讲的特征工程(Feature Engineering)

使用正则化项， 保持所有的特征，但是保证参数θj不会变得巨大。正则化项非常适合在我们拥有很多稍微有点用的特征的时候。

正则化项(regularizer)

正则化项其实也叫惩罚项(penalty term)，其作用是减缓过拟合问题，其实就是在损失函数后面加一个含有各个Θ的项，这样做的目的是让Θ也参与损失函数的计算，这样由于我们需要求的是损失函数的最小值，这个项就会限制Θ的大小。

这个正则化项的目的其实是一个权衡，我们即希望参数Θ能在训练数据集上表现得比较好，又不希望参数Θ训练出来的值非常大而产生一些奇葩的划分曲线，就像下图这样的。

实践

在这篇文章中我们将目标从预测连续值的回归问题转到将结果进行分类的分类问题。

数据都可以在我的GitHub库上下载到。

环境

如果你不想被配环境烦死的话，我真的推荐装Anaconda，除此之外要说的就是我用的都是Python3.x。

背景

本篇的实战主要是利用逻辑回归来帮助招生办筛选申请人。假设你想要基于两次考试的结果预测每个申请人是否会被录取。你有之前历史申请人的历史数据，包括两次考试的分数以及最后他们是否被录取。为了完此目的，我们将使用逻辑回归建立一个基于考试分数的分类模型（classification model ）以估计录取的可能性。

代码及注释

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

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

# 读入训练数据
# windows用户路径可能需要修改下，后期有时间可能会做统一
def loadData(path):
    trainingData = pd.read_csv(path, header=None, names=['Exam 1', 'Exam 2', 'Admitted'])

    trainingData.head()

    trainingData.describe()

    positive = trainingData[trainingData['Admitted'].isin([1])]
    negative = trainingData[trainingData['Admitted'].isin([0])]

    fig, ax = plt.subplots(figsize=(12,8))
    ax.scatter(positive['Exam 1'], positive['Exam 2'], s=50, c='b', marker='o', label='Admitted')
    ax.scatter(negative['Exam 1'], negative['Exam 2'], s=50, c='r', marker='x', label='Not Admitted')
    ax.legend()
    ax.set_xlabel('Exam 1 Score')
    ax.set_ylabel('Exam 2 Score')
    # plt.show()
    return trainingData

# 计算损失，用了矢量化编程而不是for循环，公式在博客中有详细描述和证明。
def computeLoss(X, y, theta):  
    theta = np.copy(theta)
    X = np.copy(X)
    y = np.copy(y)
    m = X.shape[0]
    h = sigmoid(np.matmul(X, theta.T))
    first = np.matmul(-(y.T), np.log(h))
    second = np.matmul((1 - y).T, np.log(1 - h))
    return np.sum(first - second) / m

# 梯度下降部分
def gradientDescent(X, y, theta, alpha, iters):  
    m = X.shape[0] # 数据项数m
    temp = np.matrix(np.zeros(theta.shape))
    # parameters = 1
    cost = np.zeros(iters)

    for i in range(iters):
        error = sigmoid(np.matmul(X, theta.T)) - y


        theta = theta - ((alpha/m) * np.matmul(X.T, error)).T
        cost[i] = computeLoss(X, y, theta)

    return theta, cost

def predict(theta, X):  
    probability = sigmoid(np.matmul(X, theta.T))
    return [1 if x >= 0.5 else 0 for x in probability]

trainingData = loadData(os.getcwd() + '/../../data/ex2data1.txt')

# 插入常数项
trainingData.insert(0, 'Ones', 1)

cols = trainingData.shape[1]  
X = trainingData.iloc[:,0:cols-1]
y = trainingData.iloc[:,cols-1:cols]

# 初始化X、Y以及theta矩阵
X = np.matrix(X.values)
y = np.matrix(y.values)
theta = np.matrix(np.zeros(3))

# 计算训练前的损失值
computeLoss(X, y, theta)

# 使用梯度下降得到模型参数
alpha = 0.001
iters = 20000
theta_fin, loss = gradientDescent(X, y, theta, alpha, iters)


# 计算训练后的参数的损失值 (不优化)
computeLoss(X, y, theta_fin)  # 


# 损失随着迭代次数的变化 (不优化)
# fig, ax = plt.subplots(figsize=(12,8))  
# ax.plot(np.arange(iters), loss, 'r')  
# ax.set_xlabel('Iterations')  
# ax.set_ylabel('Cost')  
# ax.set_title('Error vs. Training Epoch') 
# plt.show()

# 不理解为什么不优化的会这么低，是学习速率没动态变化么？
predictions = predict(theta_fin, X)  
correct = [1 if ((a == 1 and b == 1) or (a == 0 and b == 0)) else 0 for (a, b) in zip(predictions, y)]  
accuracy = (sum(map(int, correct)) % len(correct))  
print('accuracy 1 = {0}%'.format(accuracy)) # 60%


# 使用scipy的optimize来做优化
import scipy.optimize as opt
# 换了下参数位置让其符合fmin_tnc
def gradient(theta, X, y):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    
    parameters = int(theta.ravel().shape[1])
    grad = np.zeros(parameters)
    
    error = sigmoid(X * theta.T) - y
    
    for i in range(parameters):
        term = np.multiply(error, X[:,i])
        grad[i] = np.sum(term) / len(X)
    
    return grad
# 换了下参数位置让其符合fmin_tnc
def computeLoss2(theta, X, y):  
    theta = np.copy(theta)
    X = np.copy(X)
    y = np.copy(y)
    m = X.shape[0]
    h = sigmoid(np.matmul(X, theta.T))
    first = np.matmul(-(y.T), np.log(h))
    second = np.matmul((1 - y).T, np.log(1 - h))
    return np.sum(first - second) / m
result = opt.fmin_tnc(func=computeLoss2, x0=theta, fprime=gradient, args=(X, y))

predictions = predict(np.matrix(result[0]), X)  
correct = [1 if ((a == 1 and b == 1) or (a == 0 and b == 0)) else 0 for (a, b) in zip(predictions, y)]  
accuracy = (sum(map(int, correct)) % len(correct))  
print('accuracy 2 = {0}%'.format(accuracy)) # 89%
复制代码

解释其实注释里面都比较清楚，就不赘述了。

结果

数据可视化

训练集中正确率

疑问

代码中使用了两个方法，一个是和上一章一样的手动梯度下降更新theta值，另一个是使用scipy的优化方法fmin_tnc。最后的准确率差别很大不知道为什么。

本文章来源于 - 梁王(lwio、lwyj123)

主要参考斯坦福ML课程 johnwittenauer