机器学习系列（二十六）——Logistic回归

2019-07-20 本文已影响23人 Ice_spring

用回归方法做分类

Logistic回归听名字像是回归算法，其实它解决的是分类问题，为什么要叫回归呢，这是因为Logistic回归是用回归的方法来解决分类问题。那么问题来了，回归的办法如何解决分类问题呢？回归得到的结果是连续的值如何对应到离散的类别呢？一个简单的方式就是选定阀值，以二分类问题为例，利用回归我们可以通过每个样本的特征得到对应的预测数值，对于这些预测数值我们设定一个阀值，如果预测数值大于阀值则判定为类别1，否则判定为类别2，这样一来便可达到分类的目的。
Logostc的输出和概率联系在一起，它预测的是样本发生的概率，于是我们对概率设定阀值，便得到了分类输出。
$\widehat{p}=\sigma(y),\widehat{y}=\left\{ \begin{aligned} &1,&\widehat{p}\geq0.5\\ &0,&\widehat{p}\leq0.5\\ \end{aligned} \right.$

因此逻辑回归既可以看作回归算法（回归求概率 $\widehat{p}$ ），也可以看作分类算法（根据 $\widehat{p}$ 计算 $\widehat{y}$ ）。作为分类算法用时，只可以解决二分类问题，对于多分类问题Logistic回归本身是不支持的，不过可以经过改进使之支持多分类。
逻辑回归是基于线性回归的，线性回归计算的是 $y=\theta^T.x_b$ ，它的值域是 $(-\infty,+\infty)$ ，为了能让线性回归的输出转变为概率进而支持分类，需要一个函数 $\sigma(t)$ ，使得 $\sigma(y)$ 的输出在(0,1)之间，这里
$\sigma(t)=\frac{1}{1+e^{-t}}$

叫做Sigmoid函数，该函数是机器学习中非常常用的函数，我们绘制Sigmoid函数图像，来看一下该函数的形状：

import numpy as np
import matplotlib.pyplot as plt

def sigmoid(t):
    return 1/(1+np.exp(-t))
'''绘图sigmoid'''
x = np.linspace(-10,10,500)
y = sigmoid(x)
plt.plot(x,y)
plt.show()

sigmoid函数

对于Sigmoid函数，它的值域为(0,1)，恰好能对应地表达概率。线性回归的结果经过Sigmoid函数后就被映射到(0,1)上成为输出概率，再对概率进行阀值设置，就可以得到二分类输出。一般阀值设定为0.5，这也很好理解，当样本有50%以上的概率属于某一类别时就将其判断为该类别。不过对于具体问题还可以调整阀值，以增强模型性能。
逻辑回归的损失函数为：
$J(\theta)=-\frac{1}{m}\sum_{i=1}^{m}y^{(i)}log(p^{(i)})-(1-y^{(i)})log(1-p^{(i)})$

于是一个二分类问题变为：对于给定的样本数据集，我们需要找到参数 $\theta$ ，使得用这样的方式，可以最大程度获得样本数据集X对应的分类输出y。这样定义的损失函数是很好的，也很容易解释，对于一个本身是0（1）的样本，如果将其预测为了1（0），损失函数将给出非常大的惩罚。
该损失函数是没有数学解析解的，因此只能通过梯度下降法求解最小值，可以通过简单推导（这里不再给出过程，读者自行推导）得到其梯度表达式的向量形式为：
$\nabla{J(\theta)}=\frac{1}{m}.X_{b}^{T}.(\sigma(X_{b}\theta)-y)$

有心的读者可以发现，这里逻辑回归损失函数的向量式和多元线性回归损失函数的梯度向量式在形式上非常接近。

Logistic Regression的实现

由于Logstic回归的梯度形式与多元线性回归损失函数的梯度相似，于是可以稍微修改多元线性回归损失函数的梯度下降法来实现Logistic回归的梯度下降求解，在play_Ml中新建LogisticRegression.py：

import numpy as np
from .metrics import accuracy_score

class LogisticRegression:
    def __init__(self):
        self.coef_ = None
        self.interception_ =None
        self._theta = None
    def _sigmoid(self,t):
        return 1./(1.+np.exp(-t))
    def fit(self,X_train,y_train,eta=0.01, n_iters=1e4):
        assert X_train.shape[0] == y_train.shape[0],"must be the same!"

        def J(theta, X_b, y):
            y_hat = self._sigmoid(X_b.dot(theta))
            try:
                return np.sum(y*np.log(y_hat)+(1-y)*log(1-y_hat))/len(y)
            except:
                return float('inf')
    
        def dJ(theta, X_b, y):
            
            #向量化的运算
            return X_b.T.dot(self._sigmoid(X_b.dot(theta))-y)/len(X_b)
            
        def gradient_descent(X_b, y, initial_theta, eta, n_iters = 1e4, epsilon = 1e-5):
            theta = initial_theta
            i_iter = 0
            while i_iter < n_iters:
                gradient = dJ(theta, X_b, y)
                last_theta = theta
                theta = theta - eta * gradient#参数更新

                if (abs(J(theta, X_b, y)-J(last_theta, X_b, y))<epsilon):
                    break
                i_iter += 1
            return theta
        X_b = np.hstack([np.ones((len(X_train),1)),X_train])
        initial_theta = np.zeros(X_b.shape[1])
        self._theta = gradient_descent(X_b,y_train,initial_theta,eta,n_iters)
        self.interception_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

    def predict_proba(self,X_predict):
        '''返回预测结果概率向量'''
        assert self.interception_ is not None and self.coef_ is not None,"must fit first!"
        assert X_predict.shape[1] == len(self.coef_),"feature must be the same!"

        X_b = np.hstack([np.ones((len(X_predict),1)),X_predict])
        return self._sigmoid(X_b.dot(self._theta))

    def predict(self,X_predict):
        assert self.interception_ is not None and self.coef_ is not None,"must fit first!"
        assert X_predict.shape[1] == len(self.coef_),"feature must be the same!"

        proba = self.predict_proba(X_predict)
        return np.array(proba>=0.5,dtype=int)

    def score(self,X_test,y_test):
        y_predict = self.predict(X_test)
        return accuracy_score(y_test,y_predict)
        
    def __repr__(self):
        return "LogisticRegression()"

定义好逻辑回归类后，使用鸢尾花数据集来测试我们的逻辑回归。由于鸢尾花数据集是三个类别，而简单逻辑回归解决的是二分类，故只取鸢尾花数据集的前两个类别，而且为了可视化方便，只取数据集的前两个特征：

from sklearn import datasets
iris = datasets.load_iris()
X = iris.data
y = iris.target
'''逻辑回归只能解决二分类，而鸢尾花有三类，将其变为两类,只取标签为0和1的，而且为了可视化的方便，只取前两个特征'''
X = X[y<2,:2]
y = y[y<2]
X.shape
plt.scatter(X[y==0,0],X[y==0,1],color='r')
plt.scatter(X[y==1,0],X[y==1,1],color='b')
plt.show()

数据集图示

接下来使用我们的逻辑回归：

'''使用我们的逻辑回归'''
from play_Ml.model_selection import train_test_split
X_train,X_test,y_train,y_test = train_test_split(X,y,seed=666)
from play_Ml.LogisticRegression import LogisticRegression
log_reg = LogisticRegression()
log_reg.fit(X_train,y_train)
log_reg.score(X_test,y_test)

准确率

可以看到，对于我们从鸢尾花数据集分离出来的数据，逻辑回归在测试集上给出了100%正确的分类。相应地我们可以查看一下各个特征的预测概率、特征对应的真值以及预测值：

预测结果

与我们预想的一致，对于概率大于0.5的预测为了类别1，否则预测为类别0。

机器学习系列（二十六）——Logistic回归

用回归方法做分类

Logistic Regression的实现

猜你喜欢

热点阅读