逻辑回归实战练习——根据学生成绩预测是否被录取

逻辑回归实战练习——根据学生成绩预测是否被录取

前言:

在学习了梯度下降和逻辑回归的基本算法后,选取此案例来进行实践练习,本次练习主要通过python中的三大块pandas、numpy和matplotlib来实现,基本不涉及到sklearn库的调用,一方面是自己编写公式时可以加深对推导公式的理解,另一方面也是为了复习numpy中的一些基本运算命令。

OK,Let's think !!!

本次练习所涉及到的算法理论及分析思路均参考自吴恩达老师的《机器学习》和唐宇迪老师的《机器学习实践》

研究背景:

根据两次考试的结果来决定每个申请人的录取机会。现有申请人的历史数据,我们可以用它作为逻辑回归的训练集。我们将建立一个分类模型,根据考试成绩来估计入学概率。

一、数据预处理

数据预处理的目的在于对样本数据形成初步认识,从而决定数据清晰的内容

  1. 数据预览
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

df=pd.read_csv('LogiReg_data.txt',header=None,names=['Exam 1', 'Exam 2', 'Admitted'])
df.info()
df.head()

从表中可以看出样本中有三行,分别是

Exam1(第一次考是成绩)
Exam2(第二次考是成绩)
Admitted(是否呗录取) 0表示不录取,1表示录取

2. 可视化预览

逻辑回归是典型的分类问题,从表中并不能直观的看出Admitted是否有明显的边界,因此可视化是个很好的方法。在可视化之前,我们先对Admitted定义新标签positive/negative

这里我是将positive和negative分别提取出来,然后画在同一个图中,当然也可以先groupby然后直接用Pandas画图:

#设定标签 
positive=df[df['Admitted']==1]    #'Admitted'==1 为正向类
negative=df[df['Admitted']==0]    #'Admitted'==0 为负向类
#可视化数据预览
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(1,1,1)
plt.scatter(x=positive['Exam 1'],y= positive['Exam 2'],marker='x',s=50,label='positive')
plt.scatter(x=negative['Exam 1'],y= negative['Exam 2'],marker='o',s=50,label='negative')
plt.legend()
ax.set_xlabel('Exam 1')
ax.set_ylabel('Exam 2')

从图中可以清晰的观察到positive和negative之间存在这明显的分界线,可以通过逻辑回归来拟合出一条回归线。

二、建立模型

1. 逻辑回归模型的假设函数(hypothesis)如下:

#sigmoid函数
def sigmoid(z):
    return 1/(1+np.exp(-z))
#假设函数h(x)
def model(X,theta):
    return sigmoid(np.dot(X,theta.T))

这里我们需要添加一列x0(x0=1)的值,并把特征矩阵X和观测值y提取出来

#添加x0的特征值为1
df.insert(0,'ones',1)
#转化为矩阵matrix
data=df.as_matrix()
cols=data.shape[1]  #计算data的列数
X=data[:,:cols-1]
y=data[:,cols-1:cols]    #这里必须设置cols-1:cols,否则y的shape不对
#建立theta的矩阵
theta=np.zeros([1,3])

建立好预测函数后我们可以试着去运行一下,这样可以加深对矩阵运算法则的印象

model(X,theta)

得到的是一个100X1的矩阵,值全部为0.5,这是因为我们给了theta一个初始值[0,0,0],这样z值全为0,在sigmoid函数中,当z=0是,其函数值(概率)为0.5

2. 代价函数和线性回归的代价函数不同:

#代价函数cost
def cost(X,y,theta):
    left=np.multiply(-y, np.log(model(X,theta)))
    right=np.multiply(1-y, np.log(1-model(X,theta)))
    return np.sum(left-right)/ (len(X))
cost(X,y,theta)

再来检查一下运算结果:

cost(X,y,theta)

我个人对这个值的理解是当theta=[0,0,0]时,所有样本的平均误差为0.6931471805599453

3. 梯度下降的推导式与线性回归时一样的:

#计算梯度方向
def gradient(X,y,theta):
    grad=np.zeros_like(theta)
    error=(model(X,theta)-y).ravel()
    for j in range(len(theta.ravel())):
        terms=(np.multiply(error,X[:,j]))
        grad[0,j]=np.sum(terms)/(len(X))
    return grad

值得注意的是,在编写此公式时需要用到ravel()——扁平化,例如:

model(X,theta)-y
结果是100X1的二维数组
(model(X,theta)-y).ravel()
结果是一个一维的数组

因此这里需要将二维转化为一维,从shape上也可以进一步确认。如果没有扁平化,那么error应该是一个(100,1)的二维数组,X[:,j]是一个一维数组,如果通过multiply对应相乘,那么得到的结果将会是一个100X100的二维数组,显然不是我们想要的结果。

4. 梯度下降停止的策略

梯度下降停止的策略有三种:

STOP_ITER = 0 根据迭代次数
STOP_COST = 1 根据代价函数的变化值
STOP_GRAD = 2 根据梯度变化

这里学到了一个小技巧是将三种策略当作不同的变量来使用而并非字符串

'''
type:停止策略
value:针对不同的策略具有不同含义的变量
threshold:阈值
'''
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:
        #梯度下降的方向小于阈值即可
        #np.linalg.norm求范数,默认情况下=各值的平方和开根号
        return np.linalg.norm(value) < threshold

停止策略由stopcriterion函数来完成,既然时停止,所以它只能一个判断条件二嫔妃值,所以在这里stopcriterion函数返回的时是一个布尔值,用于梯度下降求解时的循环停止的条件

5. 梯度下降方式:

  • 全批量下降:每次选择所有样本进行迭代,速度慢
  • 随机下降:每次选择一个样本迭代,不能保证收敛
  • 小批量下降:选择部分样本进行迭代,实用

梯度下降方式的选择主要与选择迭代的样本数量(batchsize)有关,因此可以通过这一参数来选择下降方式

6. 特征缩放:

目的:将所有特征的尺度都缩放在-1 到1之间,这样可以帮助梯度下降算法更快的收敛

方法:用特征值减去其均值,然后除以方差即可,这样所有数据的均值都为0,方差为1,但是对与X0这一特征值不需要进行缩放

#调用sklearn 会很方便
from sklearn import preprocessing as pp  
scaled_data=data.copy()
scaled_data[:,1:3]=pp.scale(data[:,1:3])  #对x0不进行缩放

7. 洗牌

在后续梯度下降的切结过程中会使用到小批量下降的方式,因此需要抽取一定量的样本,既然涉及到抽样,那么洗牌是必不可少的工作,同时对于全批量下降来说也能打乱顺序,减小偶然性。

import numpy.random
def shuffleData(data):
    np.random.shuffle(data)     #洗牌
    cols=data.shape[1]          #提取data的列数
    X=data[:,:cols-1]           #提取特征矩阵X
    y=data[:,cols-1:cols]       #提取观测值y
    return X, y

8. 梯度下降求解:

将上述已经建立好的函数与梯度下降迭代公式进行揉和即可,这里的难点是建一个复杂函数,既要包含下降方式,又要包含停止策略,对新手来说不态友好,自己也是研究了很久,其核心思想还是通过循环和判断来完成:

data:样本数据
theta:先要给定一个初始化的theta向量,然后进行迭代
batchsize:选择迭代的样本数量
stoptype:停止策略
thresh:阈值
alpha:学习率
n=100  #样本总数为100
cost1=[]
def descent(data,theta,batchsize,stoptype,thresh,alpha):
    i=0  #迭代次数
    k=0  #batch
    X,y=shuffleData(data)  #洗牌
    grad=np.zeros_like(theta)  #梯度方向
    costs=[cost(X,y,theta)]  #计算代价函数的损失值
    
    while True:
        #1.计算梯度方向grad,batchsize始值选择样本的数量
        grad=gradient(X[k:k+batchsize],y[k:k+batchsize],theta)  
        #设置循环,每次迭代都选择新的一组[k,k+batchsize]的样本进行迭代
        k= k + batchsize
        #如果循环过程中k超过了样本总量,则让它回归至0,重新洗牌,再次循环
        if k >=n:
            k=0
            X,y=shuffleData(data)
        #2.根据梯度下降公式求解theta值
        theta= theta - alpha*grad
        #同时将每次迭代中代价函数的损失值计算出来
        costs.append(cost(X,y,theta)) 
        cost1.append(cost(X,y,theta))
        #3.选择下降停止策略:
        #循环计算次数i
        i= i+1  
        #如果选择根据迭代次数停止,则value=迭代次数
        if stoptype==STOP_ITER:
           value=i
        #如果根据代价损失来停止,则value=代价函数
        elif stoptype==STOP_COST:
           value=costs
        #若根据梯度方向来停止,则value=梯度下降方向
        elif  stoptype==STOP_GRAD:
           value=grad
        #直到所选择的value都满足阈值thresh,则停止循环
        if stopcriterion(stoptype,value,thresh):
            break
        #4.返回最后更新的theta值,迭代次数,损失值和梯度方向
        
    return theta,i-1,grad

这里在函数外添加了一个cost1的列表,其值为每次迭代中代价函数的函数值,主要是用于后续的可视化工作,来观察不同方式下降时代价函数的变化形态。

三、 回归分析

1. 不同的下降策略

1.1 全批量下降:

#全批量,停止方式为当迭代次数为5000次是停止,学习率为0.0000001
descent(scaled_data,theta,100,STOP_ITER,5000,0.0000001)

其返回结果为:

theta=[5.00068745e-05, 1.40634086e-04, 1.25510535e-04]
time.time()-start_time=1.42010498046875
fig1, ax = plt.subplots(figsize=(12,4))
ax.plot(np.arange(len(cost1)), cost1)
横轴为迭代次数,纵轴为代价函数值

从图像可以看出,代价函数呈收敛的趋势,但是下降幅度并不是很大。或许是学习率太小的缘故,接下里我们将学习率设为0.001,不断改变迭代次数来看看代价函数的变化:

可以看出修改学习率后下降效果非常明显,迭代次数从5000变化到40000过程中,代价函数均处于下降的趋势,当迭代此时无讹40000次时,逐渐开始收敛,当迭代次数为10万次时,最低点达到了接近0.2,只是此时耗时也相对较高。

1.2 随机下降

刚开始我尝试用只改变阈值thresh为1,即每次选择一个样本进行迭代,但是发现和上图(全批量)的结果几乎一样,于是这次我设置了学习率为0.001,可以看到如果选择当迭代次数为10万次时,代价函数才会开始收敛但并未达到收敛。

axs2_=[ax1,ax2,ax3,ax4]
fig2=plt.figure(figsize=(20,12))
dictmap={1:5000,2:10000,3:50000,4:100000}
for k,v in dictmap.items():
    cost1=[]
    descent(scaled_data,theta,1,STOP_ITER,v,0.0001)
    axs2_[k-1]=fig2.add_subplot(2,2,k)
    axs2_[k-1].plot(np.arange(len(cost1)), cost1)
fig.subplots_adjust(left=0.125,right=0.9,bottom=0.1,top=0.9,wspace=0.1,hspace = 0.2)

1.3 小批量下降

这次我选择学习率为0.01,通过可视化分析可以看出,当迭代次数为1万次时代价函数几乎已经达到了收敛状态,可视化过程中没有注明每种回归计算的时间,但是从理论上来讲,小批量比全量要快一点。

2. 不同的停止策略

2.1 根据代价函数的损耗来停止

顾名思义,代价函数在迭代过程中不断减小,我们可以将一次更新前后的差值与阈值相比较,如果小于阈值,说明变化已经很小了,可以听下来了。

fig4=plt.figure(figsize=(14,9))
cost1=[]
descent(scaled_data,theta,50,STOP_COST,0.000001,0.001)
ax4_1=fig4.add_subplot(2,1,1)
ax4_1.plot(np.arange(len(cost1)), cost1)
ax4_1.set_title('thresh=0.000001,t=9.27')
cost1=[]
descent(scaled_data,theta,50,STOP_COST,0.0000001,0.001)
ax4_2=fig4.add_subplot(2,1,2)
ax4_2.plot(np.arange(len(cost1)), cost1)
ax4_2.set_title('thresh=0.0000001,t=18.16')
fig4.subplots_adjust(left=0.125,right=0.9,bottom=0.1,top=0.9,wspace=0.1,hspace = 0.2)

在这里我们设定下降策略为批量下降,且batchsize=50,学习率为0.001,观察当阈值thresh分别取0.000001,0.0000001时的梯度下降

可以看出阈值越小所需要的时间也就越长,当一次更新前后的代价值小于0.0000001时,需要迭代4万多次。当然根据上述分析,如果适当的提高学习率会快很多。

2.2 根据梯度方向来停止

理论上说如果代价函数达到收敛,那么他的导函数即梯度下降的方向应该等于0,因此我们可以取一个非常接近于0的阈值,当梯度方向小于这个阈值时,迭代编会停止

刚开始我犯了一个愚蠢的错误,我将阈值设为0,试图去寻找完全收敛的那个点,然后就发现电脑跑死了。这是因为我们在定义stopcriterion函数时,设定的时grad<thresh,而不是grad<=thresh,如果thresh=0,则会一直迭代更新,直到grad为负数,显然这是不可能的。

下面我将设置batchsize为50,学习率为0.001,阈值分别为0.05和0.01来进行考察:

fig5=plt.figure(figsize=(14,9))
cost1=[]
descent(scaled_data,theta,50,STOP_GRAD,0.05,0.001)
ax5_1=fig5.add_subplot(2,1,1)
ax5_1.plot(np.arange(len(cost1)), cost1)
ax5_1.set_title('thresh=0.05,t=3.63')
cost1=[]
descent(scaled_data,theta,50,STOP_GRAD,0.01,0.001)
ax5_2=fig5.add_subplot(2,1,2)
ax5_2.plot(np.arange(len(cost1)), cost1)
ax5_2.set_title('thresh=0.01,t=11.60')
fig5.subplots_adjust(left=0.125,right=0.9,bottom=0.1,top=0.9,wspace=0.1,hspace = 0.2)

阈值为0.01的收敛效果明显比0.05要好,当然还可以取更小的阈值,这将损耗更多的时间来实现。

四、计算模型精度

最后我们来评估一下我们预测得到的theta,设置模型参数:

小批量下降,batchsize=50;停止方式为梯度方向;阈值为0.01;学习率为0.001
descent(scaled_data,theta,50,STOP_GRAD,0.01,0.001)

迭代得到的theta为[0.97353561, 2.42695873, 2.216349 ]

接着我利用Pandas模块来进行计算:

theta_new=np.array([[0.97353561, 2.42695873, 2.216349  ]])
#定义预测函数,带入model计算后,如果概率大于0.5,则为正向类1,否则为负向类0
def predict(X_,theta_):
    return [1 if x>0.5 else 0 for x in model(X_,theta_)]
#选区标准化后的特征矩阵
X_new=scaled_data[:,:3]
data_new =pd.DataFrame(scaled_data,columns=['ones','Exam1','Exam2','Admitted'])
data_new['prediction']=predict(X_new,theta_new)
#判断预测值和观测值是否相等,相等返回True,不相等返回False
accuracy=data_new['prediction']==data_new['Admitted']
#查看True的个数
accuracy.value_counts()

样本总数为100,正确率为89%,说明样本中有89%的观测值符合预测模型

五、 总结

  1. 通过本次练习对于数组的计算有了更加深刻的认识,对于数组的维度掌握额还不够好,也是这次练习中犯错最多的地方
  2. 本次练习没有重点关注学习率这一参数,而学习率又密切影响着代价函数是否会收敛,吴恩达老师指出一般取学习率α=0.01,0.03,0.1,0.3,1等,后面研究中会尝试采用这些学习率进行计算
  3. 对于最后的梯度求解descent函数是参考唐宇迪老师的代码,将所有参数进行揉和,设置了多个判断条件,对于新手来说最好是把每种策略进行拆解,例如对于三种不同的停止策略可以写出三种不同的descent函数,这样难度会小一点。不过将所有函数揉和在一起的好处是可以在后续的研究中直接当一个逻辑回归的模板来套用,非常方便
  4. 利用循环来进行批量画图也算是本次练习的一个小收获,但是代码方面依旧不够优化
  5. 最后计算得到的精度为89%,但是并没有进行评估,而且我是利用pandas来进行计算,斌不是一个通用的方法,在后面的学习中会继续对本次练习改进

发布于 2019-04-05

文章被以下专栏收录