【从零开始的机器学习】-07 分类问题与逻辑回归

导语

如我们在第二章中提到的，有监督学习主要分为回归问题和分类问题。之前的章节我们已经介绍过一元线性回归问题，多元线性回归问题，从本章开始我们将进入另一个方向——分类问题 (Classification)。

1. 什么是分类问题？

分类问题主要针对“是不是”和“有没有”的问题，大致分为：

• 二分类问题：比如猫狗识别，判断一张图片中是猫还是狗（是不是）
• 多分类问题：比如阿拉伯数字识别，判断一张图片中的数字是几（是不是）
• 多标签二分类问题：比如一个包含多个物品的图片，判断是否存在物品A，B，C，…（有没有）
• 多标签多分类问题：图片降噪，每个像素都有0~255，总共256个类别，因此为多标签256分类问题

2. 是否可以用线性回归来解决？

举个例子，垃圾邮件过滤。假如我们用一些手段获得了一系列邮件样本X (其中包含若干个特征）以及真实标签Y（表示每个样本是否是垃圾邮件），然后设置了一个线性模型h(x)，通过样本X和真实标签Y进行学习，优化参数，那我们就会得到一个拟合X的模型。在这个过程中，真实标签Y其实只有两种离散的结果：0表示非垃圾邮件，1表示垃圾邮件。为我们的线性模型h(x)输出的是连续值，所以我们需要设定一个阈值(threshold) 来判断h(x)的输出应该对应0还是1，一般的，我们习惯在h(x)<阈值时认定为0，h(x)>阈值的时候认定为1。这里为了简单，我们只用1个特征，绘制图像如下：

我们假设当h(x)大于等于0.5时，就认定预测标签为1，当h(x)小于0.5时，认定预测标签为0。如此划分貌似很合理，我们很正确的区分了垃圾邮件和正常邮件，纵线表示的是我们的阈值（大约3.5），从此x大于等于这个阈值时就是垃圾邮件，小于这个值时就是正常邮件。

但是，假如我们再加一个x很大的，y=1的点，会发生什么呢？答：阈值会向右移，导致误判。

因此，使用线性回归来处理分类问题是不合适的。从另一个角度讲，我们的标签只有0和1，但h(x)的值域却可以是负无穷到正无穷，二者并不统一。

3. 二分类

3.1 假设函数

在上面的例子中，我们使用$h_{\theta }(x)={\theta }^{T}x$作为模型，但效果不好。这时我们进行转换，设$h_{\theta }(x)=g({\theta }^{T}x)$，$z={\theta }^{T}x$，而$g(z)=\frac{1}{1+e^{-z}}$，则：

$$h_{\theta }(x)=\frac{1}{1+e^{-{\theta }^{T}x}}$$

这就是sigmoid函数（也称逻辑函数，Logistic Function，逻辑回归问题也因此得名，虽然叫回归，但其实是分类问题）

sigmoid函数的图形（例如，theta=4时）如下：

sigmoid函数有几个特性：

1. 它的输出是0~1，满足我们对预测标签范围的要求
2. 它可以表示待测样本是正样本的概率，即 $h(x)=P(y=1|x;\theta )$，h(x,θ)=0.7意味着特征为x，参数为θ的情况下，y=1的可能性是70%，而y=0的可能是就是30%，二者相加一定等于100%
3. 它关于x=0，y=0.5中心对称，可以用y=0.5做为分水岭

至于为什么要用sigmoid函数，涉及到最大熵，在此不做展开

3.2 决定边界

我们使用sigmoid函数，并规定$h_{\theta }(x)\ge 0.5$时，我们预测结果为1；$h_{\theta }(x)<0.5$时，我们预测结果为0。

通过观察sigmoid函数的图像，我们发现当$z\ge 0$时，$h_{\theta }(x)\ge 0.5$，而$z<0$时，$h_{\theta }(x)<0.5$。

又因为$z={\theta }^{T}x$，所以上面的观察又可以替换成：当${\theta }^{T}x\ge 0$时，$h_{\theta }(x)\ge 0.5$；当${\theta }^{T}x<0$时，$h_{\theta }(x)<0.5$，

那么，${\theta }^{T}x\ge 0$ 和 ${\theta }^{T}x<0$又意味着什么呢？

假设我们的模型函数是：$h_{\theta }(x)=g({\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2)$，那么为了使$h_{\theta }(x)\ge 0.5$，应该有 $$\begin{gathered} {\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2\ge 0 \\ {x}_2\ge \frac{-{\theta }_1{x}_1-{\theta }_0}{{\theta }_2} \\ {x}_2\ge -\frac{{\theta }_1}{{\theta }_2}{x}_1-\frac{{\theta }_0}{{\theta }_2} \end{gathered}$$

我们将$x_2$视为一个关于$x_1$的因变量的话，就可以画出一条直线（形如：y=kx+b），例如下面这种情况，圆点表示负样本（y=0)，×点表示正样本（y=1），红色虚线就是决定边界（Decision Boundary）。决定边界就是y=0和y=1的分界线，下图中边界以下的点我们会预测为0，边界以上的点我们会预测为1

而决定边界不一定是一条直线，还可能是圆，椭圆，三角，曲线，各种线，图形。比如：如果有两个特征x1和x2，满足决定边界关系为$x_1^2+x_2^2\ge 1$，那么决定边界将是一个半径为1，以(0,0)为中心的圆，圆内的样本我们会预测为0，圆外的样本我们会预测为1，随着特征的量和指数的增加，决定边界可以变成各种各种的形状。

3.3 代价函数

线性回归的代价函数，我们使用平方误差函数，$J(\theta )=\frac{1}{2m}{\sum }_{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)}{)}^2$，但如果我们直接沿用这个函数，会得到一个非凸函数（代价会上下波动），导致我们无法通过求导和梯度下降法最优化参数。

为了解决这个问题，分类问题我们采用另一种代价函数：
$$\begin{gathered} J(\theta )=\frac{1}{m}\sum _{i=1}^{m}Cost(h_{\theta }(x^{(i)},y^{(i)}) \\ Cost(h_{\theta }(x),y)=-log(h_{\theta }(x))ify=1 \\ Cost(h_{\theta }(x),y)=-log(1-h_{\theta }(x))ify=0 \end{gathered}$$

这个代价函数的图像如下（蓝色为y=1的代价函数，红色为y=0的代价函数）：

横坐标为$h_{\theta }(x)$，范围总是0~1，纵坐标为代价$J(\theta )$，当y=1时，$h_{\theta }(x)$越接近1，代价越小，当二者相等时，代价为0，而$h_{\theta }(x)$越接近0，代价越大，当$h_{\theta }(x)=0$时，代价为$+\infty$。类似地，当y=0时，$h_{\theta }(x)$越接近0，代价越小，越靠近1，代价越大。

可是，分段函数还是不太方便，介于我们是针对二分类问题，我们可以将cost函数合并成1个：
$$Cost(h_{\theta }(x),y)=-ylog(h_{\theta }(x))-(1-y)log(1-h_{\theta }(x))$$

这个函数和上面的分段函数其实是等价的，因为当y=1时，1-y=0，后一项会变成0；当y=0时，前一项会变成0，但如此操作之后，我们可以更加方便的进行计算和代码实现了。将$Cost$函数代入代价函数$J(\theta )$，得到：
$$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)}))]$$

如果我们用矩阵表示的话，代价函数可以化简为：
$$J(\theta )=\frac{1}{m}\cdot (-y^{T}log(h)-(1-y{)}^{T}log(1-h))，h=g(X\theta )$$

3.4 梯度下降法

现在我们知道了模型，知道了代价函数，接下来该最优化参数了。我们同样使用梯度下降法，则每一次迭代得到的新参数为（实际求导和化简的过程省略）：
$${\theta }_j:={\theta }_j-\frac{\alpha }{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})x_j^{(i)}$$

看上去和我们在线性回归中得到的学习算法一致，但线性回归中的模型$h(x)={\theta }^{T}x$，而这里的模型为$h_{\theta }(x)=\frac{1}{1+e^{-{\theta }^{T}x}}$，所以其实是不一样的。

如果用矩阵表示的话，上面的学习算法可以化简为：
$$\theta :=\theta -\frac{\alpha }{m}{X}^T(g(X\theta )-\vec{y}))$$

3.5 代码实现

import numpy as np
from matplotlib import pyplot as plt
import csv
from sklearn.linear_model import LogisticRegression


def legend_without_duplicate_labels(ax):
    handles, labels = ax.get_legend_handles_labels()
    unique = [(h, l) for i, (h, l) in enumerate(zip(handles, labels)) if l not in labels[:i]]
    ax.legend(*zip(*unique))


def plotData(X, y, m):
    pos = np.argwhere(y==1)
    neg = np.argwhere(y==0)
    fig, ax = plt.subplots(facecolor="w")
    ax.plot(X[pos, 0], X[pos, 1], '+', color='g', label="Admitted")
    ax.plot(X[neg, 0], X[neg, 1], 'o',color='r', label="Not Admitted")
    plt.xlabel('Exam 1 Score')
    plt.ylabel('Exam 2 Score')
    legend_without_duplicate_labels(ax)
    # plt.show()

def gradientDescent(X,y,theta,alpha,iterations):
    m = len(y)
    J_history = []

    for iter in range(1,iterations+1):
        J,grad = computeCost(X,y,theta)
        theta = theta-grad*alpha
        J_history.append([iter,J])
    return theta,J_history

def featureNormalize(X)->[]:
    mu = np.mean(X,axis=0)
    sigma = np.std(X,axis=0)
    X = (X-mu)/sigma
    return [X,mu,sigma]

def computeCost(X,y,theta):
    m = len(y)
    J = sum(np.log(sigmoid(X.dot(theta))).T.dot(-1 * y) - np.log(1 - sigmoid(X.dot(theta))).T.dot(1 - y)) / m
    grad = X.T.dot(sigmoid(X.dot(theta))-y)/m
    return J,grad

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

# def computeCost_Reg(X,y,theta,la):
#     m = len(y)
#     J = sum(np.log(sigmoid(X.dot(theta))).T.dot(-1 * y) - np.log(1 - sigmoid(X.dot(theta))).T.dot(1 - y)) / m
#     J = J + sum(np.square(theta[1:]))*la/(2*m)
#     grad = X.T.dot(sigmoid(X.dot(theta))-y)/m
#     temp = theta*la/m
#     temp[0,0] = 0
#     grad = grad+temp
#     return J,grad


def predict_norm(X,mu,sigma,theta):
    X = (X - mu) / sigma
    m, n = X.shape
    X = np.hstack((np.ones((m, 1)), X))
    return (sigmoid(X.dot(theta))>=0.5)


def predict(X,theta):
    m = len(X)
    X = np.hstack((np.ones((m, 1)), X))
    return (sigmoid(X.dot(theta))>=0.5)

def main():
    data = []
    with open('data1.txt') as f:
        csv_reader = csv.reader(f,delimiter=',')
        for row in csv_reader:
            data.append(row)
    f.close()

    data = np.mat(data,dtype='float32')
    m,n = data.shape
    X = data[:,0:n-1]

    y = np.mat(data[:,n-1],dtype='int32')

    plotData(X,y,m)

    test_X = X
    X,mu,sigma = featureNormalize(X)

    print(sigmoid(0))

    X = np.hstack((np.ones((m, 1)), X))
    theta = np.zeros((n, 1))

    # cost, grad = computeCost(X, y, theta)
    # print(cost)
    # print(grad)

    test_theta = np.mat([[-25.1613],[0.2062],[0.2015]])
    # cost,grad = computeCost(X,y,test_theta)
    # print(cost)
    # print(grad)

    alpha = 0.3
    iterations = 10000
    theta,history = gradientDescent(X,y,theta,alpha,iterations)
    print(theta)

    print(predict_norm(np.mat([45, 85]), mu, sigma, theta))
    print(predict_norm(np.mat([55, 75]), mu, sigma, theta))
    p = predict_norm(test_X,mu,sigma,theta)
    ans = p-y
    correct = np.argwhere(ans==0)
    print(len(correct)/m*100)

    p2 = predict(test_X,test_theta)
    ans2 = p2 - y
    correct2 = np.argwhere(ans2 == 0)
    print(len(correct2) / m * 100)

    clf = LogisticRegression()
    clf.fit(test_X,y)
    print(clf.score(test_X,y))

    print(computeCost_Reg(X,y,theta,0.1))

if __name__ == '__main__':
    main()

4. 多元分类

上面我们提到了二分类问题，那么对于多分类问题（比如，阿拉伯数字识别），我们可以将其拆分成多个二分类问题。

假设我们有三种类，下面是我们的样本，绿色表示$I$类，红色表示$II$类，黑色表示$III$类。那么我们可以建立3套二分类模型，分别识别“是否是$I$类”，“是否是$II$类”，以及“是否是$III$类”三个问题。

对于$I$类，我们设模型为$h_{\theta }^{(1)}(x)$，这个模型输出的结果为：在特征x和参数$\theta$下，y=1的概率，剩下两类都会被认为是y=0，我们只关心是否是$I$类。因此，从$I$类的视角来看，上面的样本图应该如下所示，青色代表无关样本，绿色代表：$I$类样本，红色虚线表示决定边界。

对于$II$类，我们设模型为$h_{\theta }^{(2)}(x)$，这个模型输出的结果为：在特征x和参数$\theta$下，y=1的概率，剩下两类都会被认为是y=0，我们只关心是否是$II$类。因此，从$II$类的视角来看，上面的样本图应该如下所示，青色代表无关样本，红色代表：$II$类样本，红色虚线表示决定边界。

对于$III$类，我们设模型为$h_{\theta }^{(3)}(x)$，这个模型输出的结果为：在特征x和参数$\theta$下，y=1的概率，剩下两类都会被认为是y=0，我们只关心是否是$III$类。因此，从$III$类的视角来看，上面的样本图应该如下所示，青色代表负样本，黑色代表：$III$类样本，红色虚线表示决定边界。

在分别最优化模型参数之后，对于一个待测样本x，我们分别用三个模型去计算，得到的就是x是$I$类，$II$类，$III$类的概率，然后去概率最大的，就是我们的预测结果。