逻辑回归 - 简书

逻辑回归
逻辑回归是一种经典的二分类算法，它的决策边界可以是非线性的（包含高阶项）
Sigmoid函数
公式： $g(z)=\frac{1}{1+{{e}^{-z}}}$ ，其中自变量取值为任意实数
解释：将任意的输入映射到了[0,1]区间
比如在线性回归中可以得到一个预测值（ $\hat{y}$ ），再将该值映射到Sigmoid函数，就完成了由值到概率的转换，也就是分类任务
假设预测函数 ${{h}_{\mathbf{\theta }}}(x)=g({{\mathbf{\theta }}^{\mathbf{T}}}\mathbf{X})=\frac{1}{1+{{e}^{-{{\mathbf{\theta }}^{\mathbf{T}}}\mathbf{X}}}}$
分类任务： $\begin{align} & p(y=1|\mathbf{X};\mathbf{\theta })={{h}_{\mathbf{\theta }}}(x) \\ & p(y=0|\mathbf{X};\mathbf{\theta })=1-{{h}_{\mathbf{\theta }}}(x) \\ \end{align}$
相当于0-1分布，整合上述式子得 $p(y|\mathbf{X};\mathbf{\theta })={{({{h}_{\mathbf{\theta }}}(x))}^{y}}{{(1-{{h}_{\mathbf{\theta }}}(x))}^{1-y}}$
似然函数：
$L(\mathbf{\theta })=\prod\limits_{i=1}^{m}{p(y|\mathbf{X};\mathbf{\theta })}=\prod\limits_{i=1}^{m}{{{({{h}_{\mathbf{\theta }}}({{x}_{i}}))}^{{{y}_{i}}}}{{(1-{{h}_{\mathbf{\theta }}}({{x}_{i}}))}^{1-{{y}_{i}}}}}$
取对数似然函数：
$\ln L(\mathbf{\theta })=\sum\limits_{i=1}^{m}{{{y}_{i}}\ln ({{h}_{\theta }}({{x}_{i}}))\text{+(1-}}{{y}_{i}})\ln (1-{{h}_{\theta }}({{x}_{i}}))$
现在要求上述对数似然函数的最大值，如果要用梯度下降法，需要将目标函数转化为求最小值【梯度的方向是值越来越大的方向，求最小值，应该朝着梯度反方向走】，因此引入
$J(\mathbf{\theta })=-\frac{1}{m}\ln L(\mathbf{\theta })=-\frac{1}{m}\sum\limits_{i=1}^{m}{{{y}_{i}}\ln ({{h}_{\theta }}({{x}_{i}}))\text{+(1-}}{{y}_{i}})\ln (1-{{h}_{\theta }}({{x}_{i}}))$
对上式求偏导
$\begin{align} & \frac{\partial J(\mathbf{\theta })}{\partial {{\theta }_{j}}}=-\frac{1}{m}\sum\limits_{i=1}^{m}{\left[ {{y}_{i}}\frac{1}{{{h}_{\mathbf{\theta }}}({{x}_{i}})}\frac{\partial {{h}_{\mathbf{\theta }}}({{x}_{i}})}{\partial {{\theta }_{j}}}-(1-{{y}_{i}})\frac{1}{1-{{h}_{\mathbf{\theta }}}({{x}_{i}})}\frac{\partial {{h}_{\mathbf{\theta }}}({{x}_{i}})}{\partial {{\theta }_{j}}} \right]} \\ & \text{ =}-\frac{1}{m}\sum\limits_{i=1}^{m}{\left[ {{y}_{i}}\frac{1}{{{h}_{\mathbf{\theta }}}({{x}_{i}})}-(1-{{y}_{i}})\frac{1}{1-{{h}_{\mathbf{\theta }}}({{x}_{i}})} \right]}\frac{\partial {{h}_{\mathbf{\theta }}}({{x}_{i}})}{\partial {{\theta }_{j}}} \\ & \text{ =}-\frac{1}{m}\sum\limits_{i=1}^{m}{\left[ {{y}_{i}}\frac{1}{g({{\mathbf{\theta }}^{\mathbf{T}}}{{\mathbf{x}}_{i}})}-(1-{{y}_{i}})\frac{1}{1-g({{\mathbf{\theta }}^{\mathbf{T}}}{{\mathbf{x}}_{i}})} \right]}\frac{\partial g({{\mathbf{\theta }}^{\mathbf{T}}}{{\mathbf{x}}_{i}})}{\partial {{\theta }_{j}}} \\ & \text{ =}-\frac{1}{m}\sum\limits_{i=1}^{m}{\left[ {{y}_{i}}\frac{1}{g({{\mathbf{\theta }}^{\mathbf{T}}}{{\mathbf{x}}_{i}})}-(1-{{y}_{i}})\frac{1}{1-g({{\mathbf{\theta }}^{\mathbf{T}}}{{\mathbf{x}}_{i}})} \right]}\frac{{{e}^{-{{\mathbf{\theta }}^{\mathbf{T}}}{{\mathbf{x}}_{\mathbf{i}}}}}}{1+{{e}^{-{{\mathbf{\theta }}^{\mathbf{T}}}{{\mathbf{x}}_{\mathbf{i}}}}}}x_{i}^{j} \\ & \text{ =}-\frac{1}{m}\sum\limits_{i=1}^{m}{\left[ {{y}_{i}}\frac{1}{g({{\mathbf{\theta }}^{\mathbf{T}}}{{\mathbf{x}}_{i}})}-(1-{{y}_{i}})\frac{1}{1-g({{\mathbf{\theta }}^{\mathbf{T}}}{{\mathbf{x}}_{i}})} \right]}g({{\mathbf{\theta }}^{\mathbf{T}}}{{\mathbf{x}}_{i}})(1-g({{\mathbf{\theta }}^{\mathbf{T}}}{{\mathbf{x}}_{i}}))x_{i}^{j} \\ & \text{ =}-\frac{1}{m}\sum\limits_{i=1}^{m}{\left[ {{y}_{i}}(1-g({{\mathbf{\theta }}^{\mathbf{T}}}{{\mathbf{x}}_{i}}))-(1-{{y}_{i}})g({{\mathbf{\theta }}^{\mathbf{T}}}{{\mathbf{x}}_{i}}) \right]}x_{i}^{j} \\ & \text{ =}-\frac{1}{m}\sum\limits_{i=1}^{m}{\left[ {{y}_{i}}-g({{\mathbf{\theta }}^{\mathbf{T}}}{{\mathbf{x}}_{i}}) \right]}x_{i}^{j} \\ & \text{ =}\frac{1}{m}\sum\limits_{i=1}^{m}{\left[ {{h}_{\mathbf{\theta }}}({{x}_{i}})-{{y}_{i}} \right]}x_{i}^{j} \\ \end{align}$
确定了梯度方向之后，利用上式更新参数，得到
$\theta _{j}^{'}={{\theta }_{j}}-\alpha \frac{1}{m}\sum\limits_{i=1}^{m}{\left[ {{h}_{\mathbf{\theta }}}({{x}_{i}})-{{y}_{i}} \right]}x_{i}^{j}$
按照上式继续迭代
下面是一个关于实现逻辑回归的小例子
目标：建立一个逻辑回归模型来预测一个大学生是否被大学录取。
特征：两次考试成绩
导入数据

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
pdData=pd.read_csv('C:/Users/lenovo/Desktop/LogiReg_data.txt',header=None,names=['Exam 1','Exam 2','Admitted'])#header=None不让数据自己指定列名
pdData.insert(0,'ones',1) #为了考虑常数项
pdData.head()

数据如以下所示

image.png

通常在做数据处理之前，需要做一些数据预处理，比如数据标准化等，这边先不做，看看不做标准化会对结果什么影响。
逻辑回归的关键是计算损失值和梯度方向
为了计算这两个函数，需要引入sigmoid函数【将值映射到概率】

#定义sigmoid函数
def sigmoid_2(z):
    return 1/(1+np.exp(-z))
#预测函数模型
def model(X,theta):
    return sigmoid_2(np.dot(X,theta.T))

切片，分出X和y

orig_data=pdData.as_matrix()# 转换成array类型
cols=orig_data.shape[1]
X=orig_data[:,:cols-1]#前三列,切片012
y=orig_data[:,cols-1:]#和y=orig_data[:,cols-1]结果额维度不一样
theta=np.zeros([1,3]) #theta初始化是0

注：这里要注意一个关于数组切片的问题

image.png

上面画圈的两个语句一个返回一维数组，一个返回三维数组

损失函数

将对数似然函数取相反数
$-\ln L(\mathbf{\theta })=-\sum\limits_{i=1}^{n}{{{y}_{i}}\ln ({{h}_{\theta }}({{x}_{i}}))-\text{(1-}}{{y}_{i}})\ln (1-{{h}_{\theta }}({{x}_{i}}))$
求平均损失
$J(\mathbf{\theta })=\frac{1}{n}\ln L(\mathbf{\theta })=\frac{1}{n}\sum\limits_{i=1}^{n}{\text{-}{{y}_{i}}\ln ({{h}_{\theta }}({{x}_{i}}))\text{-(1-}}{{y}_{i}})\ln (1-{{h}_{\theta }}({{x}_{i}}))$
目的是求平均损失的最小值

def cost(X,y,theta):#对应于J(theta)
    left=np.multiply(-y,np.log(model(X,theta)))#np.multiply对应位置相乘
    right=np.multiply(1-y,np.log(1-model(X,theta)))
    return np.sum(left-right)/len(X)

计算梯度

$\frac{\partial J(\mathbf{\theta })}{\partial {{\theta }_{j}}}\text{=}-\frac{1}{n}\sum\limits_{i=1}^{n}{\left[ {{y}_{i}}-g({{\mathbf{\theta }}^{\mathbf{T}}}{{\mathbf{x}}_{i}}) \right]}x_{i}^{j}\text{=}\frac{1}{n}\sum\limits_{i=1}^{n}{\left[ {{h}_{\mathbf{\theta }}}({{x}_{i}})-{{y}_{i}} \right]}x_{i}^{j}$

def gradient(X,y,theta):
    grad=np.zeros(theta.shape)
    error=(model(X,theta)-y).ravel()
   # print(error.shape)
    for j in range(len(theta.ravel())):#求每个参数的偏导
        term=np.multiply(error,X[:,j])
        grad[0,j]=np.sum(term)/len(X)
    return grad

参数更新

首先设定三种停止策略【迭代次数、损失值差距、梯度】

def stopCriterion(type,value,threshold):
    if type==STOP_ITER:
        return value>threshold #返回逻辑值#迭代次数
    elif type==STOP_COST:
        return abs(value[-1]-value[-2])<threshold#两次迭代结果的损失值差距较小
    elif type==STOP_GRAD:
        return np.linalg.norm(value)<threshold #梯度较小

为了避免原数据存在某种规律【比如先收集女生成绩再收集男生成绩】，将原数据顺序打乱

import numpy.random
#洗牌，将数据顺序打乱
def shuffleData(data):
    np.random.shuffle(data)  #shuffle() 方法将序列的所有元素随机排序
    cols=data.shape[1]
    X=data[:,:cols-1]
    y=data[:,cols-1:]
    return X,y

然后进行参数更新【引入time库是为了比较不同策略的运行时间】

import time
def descent(data,theta,batchSize,stopType,thresh,alpha):
    #batchSize等于总样本数，即批量梯度下降；batchSize等于1，即随机梯度下降；batchSize取1~n之间的数据，即小批量梯度下降
    
    
    #梯度下降求解
    #初始化
    init_time=time.time()  #比较不同梯度下降方法的运行速度
    i=0 #迭代次数，第0次开始
    k=0 #batch
    X,y=shuffleData(data)
    grad=np.zeros(theta.shape)  #计算的梯度
    costs=[cost(X,y,theta)]  #计算损失值
    
    
    
    while True:
        grad=gradient(X[k:k+batchSize],y[k:k+batchSize],theta)
        k += batchSize #取样本数据
        if k >= n:
            k = 0
            X,y=shuffleData(data) #重新洗牌
        theta = theta - alpha*grad  # 参数更新
        costs.append(cost(X,y,theta))
        i += 1
        
        #停止策略
        if stopType==STOP_ITER:
            value = i
        elif stopType==STOP_COST:
            value=costs
        elif stopType==STOP_GRAD:
            value=grad
        if stopCriterion(stopType,value,thresh):#如果if语句为真，就跳出整个循环
            break
        #print(grad,np.linalg.norm(grad))
    return theta,i-1,costs,grad,time.time()-init_time

注：np.linalg.norm([4,3])表示 $\sqrt{{{3}^{2}}+{{4}^{2}}}$

def runExpe(data,theta,batchSize,stopType,thresh,alpha):
    theta,iter,costs,grad,dur=descent(data,theta,batchSize,stopType,thresh,alpha)
#     name='Original' if (data[:,2]>2).sum()>1 else 'Scaled'
#     name += 'data-learning rate: {} -'.format(alpha)
#     if batchSize==n: strDescType='Gradient'
#     elif batchSize==1: strDescType='Stochastic'
#     else:strDescType='mini-batch {} '.format(batchSize)
#     name += strDescType + 'descent-stop: '
#     if stopType==STOP_ITER: strStop='{} iterations'.format(thresh)
#     elif stopType==STOP_COST: strStop='cost change < {}'.format(thresh)
#     else: strStop ='gradient norm < {}'.format(thresh)
#     name += strStop
    fig,ax=plt.subplots(figsize=(12,4))
    ax.plot(np.arange(len(costs)),costs,c='r')
    ax.set_xlabel('Iteration')
    ax.set_ylabel('Cost')
    #ax.set_title(name.upper())
    print('iter: {}, last cost: {:03.2f}, duration: {:03.2f}s'.format(iter,costs[-1],dur))
    return theta

不同的停止策略

以批量梯度下降为例

停止条件为迭代次数

n=100
runExpe(orig_data,theta,n,STOP_ITER,thresh=5000,alpha=0.000001)

image.png

看似损失值已经稳定在最低点0.63

停止条件为损失值
设定阈值为0.000001，需要迭代110000次左右

runExpe(orig_data,theta,n,STOP_COST,thresh=0.000001,alpha=0.001)

image.png

损失值最低为0.38，似乎还可以进一步收敛

停止条件为梯度大小
设定阈值0.05，需要迭代40000次左右

runExpe(orig_data,theta,n,STOP_GRAD,thresh=0.05,alpha=0.001)

image.png

损失值最小为0.49，似乎还可以进一步收敛
综上，基于批量梯度下降方法，上述三种停止条件得到的损失函数值为0.63、0.38和0.49，迭代次数分别为5000次、110000次和40000次，迭代次数越多，损失值越小

对比不同的梯度下降方法

停止策略为迭代次数

随机梯度下降

runExpe(orig_data,theta,1,STOP_ITER,thresh=5000,alpha=0.001)

image.png

波动非常大，迭代过程不稳定，这也是随机梯度下降的主要缺点
尝试降低学习率为0.000001，增加迭代次数为15000

runExpe(orig_data,theta,1,STOP_ITER,thresh=15000,alpha=0.000001)

image.png

效果要好一些，损失值似乎稳定在0.63,根据上面的结果可知，0.63不算是一个特别合理的值

小批量梯度下降

#取样本为16
runExpe(orig_data,theta,16,STOP_ITER,thresh=15000,alpha=0.001)

image.png

上下波动，迭代过程不稳定
尝试调低学习率为0.000001

runExpe(orig_data,theta,16,STOP_ITER,thresh=15000,alpha=0.001)

image.png

降低学习率之后没有效果，迭代过程依旧不稳定
因此，可能不是模型本身的问题，而是数据本身的问题，尝试着对数据做一些变换，此处对数据进行标准化，用标准化后的数据求解

#标准化
from sklearn import preprocessing as pp
scaled_data=orig_data.copy()
scaled_data[:,1:3]=pp.scale(scaled_data[:,1:3])
#用标准化后的数据求解
runExpe(scaled_data,theta,16,STOP_ITER,thresh=15000,alpha=0.001)

image.png

损失值收敛到0.28，比0.63好很多
再尝试一下梯度为停止条件的情况

runExpe(scaled_data,theta,16,STOP_GRAD,thresh=0.004,alpha=0.001)

image.png

迭代次数由15000增加到60000多，损失值由0.28降低到0.22，又改善了一步

计算模型分类结果的精确率

#设定阈值为0.5，大于0.5就可以入学
def predict(X,theta):
    return [1 if x >= 0.5 else 0 for x in model(X,theta)]

scaled_X=scaled_data[:,:3]
y=scaled_data[:,3]
theta = runExpe(scaled_data,theta,16,STOP_GRAD,thresh=0.004,alpha=0.001)
predictions = predict(scaled_X,theta)
correct = [1 if a==b else 0 for (a,b) in zip(predictions,y)]   #真实值与预测值相等，同为1或者同为0
accuracy = sum(correct)/len(correct)

返回accuracy等于0.89
通过这个例子，算是对逻辑回归的基本原理有一个比较清晰的认识了。

sigmoid函数：将值映射到概率的函数
model：返回预测结果值
cost：根据参数计算损失
gradient：计算每个参数的梯度方向
descent：进行每个参数的更新
accuracy：计算精度