逻辑是:
1.定义sigmoid函数,目的是将函数值映射到概率;
2.定义model函数,将线性方程值带入sigmoid函数;
3.定义cost损失函数(似然函数变换得到),目的是为了梯度下降求参;
4.定义梯度下降的不同迭代停止原则;
5.对比不同的梯度下降方法;
6.查看模型精度
1. 导入数据
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')
#import os
#os.chdir('C:\\Users\\Administrator.PC-201912302019\\code')
file='LogiReg_data.txt'
pdData=pd.read_csv(file,header=None,
names=['Exam1','Exam2','Admit'])
print(pdData.shape)
pdData.head(2)
2. 探索数据关系
pos=pdData[pdData.Admit==1]
neg=pdData[pdData.Admit==0]
fig,ax=plt.subplots(1,1,figsize=(10,5))
ax.scatter('Exam1','Exam2',data=pos,s=30,c='b',marker='o',label='Adimitted')
ax.scatter('Exam1','Exam2',data=neg,s=30,c='r',marker='x',label='Not Adimitted')
ax.legend()
ax.set_xlabel('Exam1 score')
ax.set_ylabel('Exam2 score')
3. 逻辑回归主体
目标:建立分类器,求解三个参数
设定阈值,将预测的概率值转换为0-1,判断是否录取
要完成的模块
- sigmiod:映射到概率的函数
- model:返回预测结果值
- cost:根据参数计算损失函数
- gradient:计算每个参数的梯度方向
- descent:进行参数更新
- accuracy:计算精度
3.1 映射到概率的函数sigmoid()
#定义函数
def sigmoid(z):
return 1/ (1+np.exp(-z))
#设置小例子查看函数图像
nums=np.arange(-10,10,step=1) #步长为1,生成20个数字
print(len(nums))
print(nums,'\n',sigmoid(nums))
fig,ax=plt.subplots(figsize=(12,4))
ax.plot(nums,sigmoid(nums),'r')
3.2 预测结果函数model()
def model(X,theta):
return sigmoid(np.dot(X,theta.T)) #点乘
#函数结果返回一列概率
pdData.insert(0,'ones',1) #数据框首列插入一列,作为函数式的常数项
orig_data=pdData.as_matrix() #数据框转化为array
cols=orig_data.shape[1]
X=orig_data[:,0:cols-1]
y=orig_data[:,cols-1:cols]
theta=np.zeros([1,3]) #参数先进行占位
#自变量查看
X[:5]
3.3 损失函数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))
3.4 梯度函数gradient()
对theta求偏导
def gradient(X,y,theta):
grad=np.zeros(theta.shape)
error=(model(X,theta)-y).ravel() #转变为数组形式
for j in range(len(theta.ravel())):
term=np.multiply(error,X[:,j])
grad[0,j]=np.sum(term)/len(X)
return grad
3.5 定义不同的迭代停止原则
- 根据迭代次数(迭代到一定次数停止)
- 根据损失函数(损失函数达到阈值停止)
- 根据梯度(梯度变化很小时停止)
STOP_ITER=0
STOP_COST=1
STOP_GRAD=2
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 numpy.linalg.norm(value) < threshold
3.6 比较三种不同梯度下降方法
- 批量梯度下降:考虑所有样本(易得最优解,但速度慢)
- 随机梯度下降:随机找一个样本(不一定朝着收敛的方向)
- 小批量梯度下降:每次更新一小部分数据来算
#先定义一个洗牌函数
import numpy.random
def shuffleData(data):
np.random.shuffle(data)
cols=data.shape[1]
X=data[:,0:cols-1]
y=data[:,cols-1:]
return X,y
import time
def descent(data,theta,batchSize,stopType,threshold,alpha):
init_time=time.time()
i=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 #取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,threshold): break
return theta,i-1,costs,grad,time.time()-init_time
3.7 定义运行函数
def runExpe(data,theta,batchSize,
stopType,threshold,alpha):
#(数据,参数,梯度下降样本量,停止规则,阈值,学习率)
theta,iter,costs,grad,dur=\
descent(data,theta,batchSize,stopType,threshold,alpha)
name='original' if (data[:1]>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(threshold)
elif stopType==STOP_COST: strstop='costs change <{}'.format(threshold)
else: strstop='gradient norm <{}'.format(threshold)
name += strstop
print('***{}\nTheta:{}-Iter:{}-Last cost:{:03.2f}-Duration:{:03.2f}s'.format(
name,theta,iter,costs[-1],dur))
fig,ax=plt.subplots(figsize=(12,4))
ax.plot(np.arange(len(costs)),costs,'r')
ax.set_xlabel('Iteration')
ax.set_ylabel('Cost')
ax.set_title(name.upper()+'-Error vs. Iteration')
return theta
4.对比不同的停止策略
4.1 设定迭代次数
#选择的梯度下降方法是基于所有样本的
n=100
runExpe(orig_data,theta,n,STOP_ITER,threshold=5000,alpha=0.000001)
4.2 根据损失函数值停止
设定阈值1E-6,差不多需要110000次迭代
runExpe(orig_data,theta,n,STOP_COST,threshold=0.000001,alpha=0.001)
4.3 根据梯度变化停止
设定阈值0.05,差不多需要4000次迭代
runExpe(orig_data,theta,n,STOP_GRAD,threshold=0.05,alpha=0.001)
5. 对比不同的梯度下降方法
5.1 Stochastic descent
runExpe(orig_data,theta,1,STOP_ITER,threshold=5000,alpha=0.001)
不稳定,试着把学习率调小一点
runExpe(orig_data,theta,1,STOP_ITER,threshold=5000,alpha=0.000002)
速度快,但稳定性差,需要很小的学习率
5.2 Mini-batch descent
runExpe(orig_data,theta,16,STOP_ITER,threshold=15000,alpha=0.001)
浮动仍然比较大,尝试下对数据进行标准化,将数据按其属性(按列进行)减去其均值,除以方差。最后得到的结果是,对每个属性/每列来说所有数据都聚集在0附近,方差为1
6. 不同参数组合建模
6.1 标准化后用全样本,停止原则用迭代次数
from sklearn import preprocessing as pp
scaled_data=orig_data.copy()
scaled_data[:,1:3]=pp.scale(orig_data[:,1:3])
runExpe(scaled_data,theta,n,STOP_ITER,threshold=5000,alpha=0.001)
较之前好一些,原始数据,只能达到0.61,这里可以达到0.38.所以对数据进行预处理是十分重要的
6.2 标准化后用全样本,停止原则用梯度
runExpe(scaled_data,theta,n,STOP_GRAD,threshold=0.02,alpha=0.001)
更多的迭代次数会使损失下降的更多!
6.3 标准化后用随机梯度 停止原则梯度
theta=runExpe(scaled_data,theta,1,STOP_GRAD,threshold=0.002/5,alpha=0.001)
随机梯度也下降的很快,但需要迭代的次数也更多,用mini-batch更合适!!
6.4 标准化后batch 停止原则为梯度
runExpe(scaled_data,theta,16,STOP_GRAD,threshold=0.002*2,alpha=0.001)
7. 精度
#设定阈值
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]
predictions=predict(scaled_X,theta)
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={0}%'.format(accuracy))
accuracy=89%