梯度下降法求解逻辑回归_python实现

逻辑回归是一种经典的二分类算法（当然也可以用作多分类），也是我们在进行数据分析和挖掘过程中最常使用的算法之一。
通常我们在选择机器学习算法的时候，往往都是先选取简单的模型，如逻辑回归模型，因为它的原理简单、可解释性强，我们将简单的模型和复杂的模型取得的结果进行对比，如果它们效果差不多，那我们也没有必要选择复杂的模型了，逻辑回归就是这样一种有效、易用的算法。
关于逻辑回归的推导过程，学习资料有很多，这里不再进行赘述。在这里 我主要给大家梳理一下：①逻辑回归算法的计算逻辑是怎么样的；②以及如何用python来模块化实现这个流程。

逻辑回归的计算逻辑

逻辑回归的计算逻辑可以用以上一张图来表示。举例来说，假设我们现在拿到一组样本数据：exam1和exam2属性代表两门考试的成绩，addimitted表示最终是否被录取（1-录取，0-未录取），即每一行代表一个学生的两门考试成绩以及最终是否被录取，现在要使用逻辑回归算法来预测某学生最终是否能被录取？

exam1	exam2	addimitted
34.623660	78.024693	0
30.286711	43.894998	0
35.847409	72.902198	0
60.182599	86.308552	1
79.032736	75.344376	1
…（100个样本）	…	…

①首先我们要将原始数据分成X和y，用于后面的训练，X包含exam1和exam2数据，y包含addimitt数据
②指定一个θ初始值以及学习率α，如果可以的话，还可以自己设定一个batchsize来选择梯度下降的方法。
③计算损失函数和梯度值。
④判断是否达到停止条件（如损失值的变化很小了、梯度值基本不再变化了、又或者达到了指定的迭代次数），如果达到了就返回θ值，否则继续迭代更新θ值。
以上就是逻辑回归算法实现的基本思路。

使用python实现

以下是使用python实现以上算法流程时主要要完成的功能模块，计算损失值、梯度值、进行参数更新等，除了这些主要模块，还有停止策略等模块没有列出。
代码如下：

#加载相关的库
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import numpy.random
import time
%matplotlib inline

#载入文件
import os
path = 'data' + os.sep +'LogiReg_data.txt'
pdData = pd.read_csv(path, header = None, names = ['Exam 1','Exam 2', 'Admitted'])

#一、对数据预处理
#这里主要是在数据最前面插入了数值全为1的列
pdData.insert(0, 'Ones', 1)
orig_data = pdData.as_matrix()
cols = orig_data.shape[1]
X = orig_data[:, 0:cols - 1]
y = orig_data[:, cols-1:cols]
theta = np.zeros([1,3])

#二、实现各个功能模块
#①sigmoid函数转换
def sigmoid(z):
    return 1/(1 + np.exp(-z))

#②输入X和θ得到预测值
def model(X, theta):
    return sigmoid(np.dot(X, theta.T))

#③计算损失函数
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)

#④计算梯度
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

#⑤设定停止策略
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 np.linalg.norm(value) < threshold

#⑥对数据进行洗牌
def shuffleData(data):
    np.random.shuffle(data)
    cols = data.shape[1]
    X = data[:, 0:cols-1]
    y = data[:, cols-1:]
    return X,y
   
#三、使用以上模块定义梯度下降求解的函数
def descent(data, theta, batchSize, stopType, thresh, alpha):
    #梯度下降求解
    #batchSize参数指定了使用批量、随机还是小批量梯度下降
    #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 #去batch数量个数倍
        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):
            break
    return theta, i-1, costs, grad, time.time() - init_time
 
 #四、定义运行的函数，这里主要是为了让结果更加直观，主要部分我们已经在上述步骤完成了
 def runExpe(data, theta, batchSize, stopType, thresh, alpha):
    theta, iter, costs, grad, dur = descent(data, theta, batchSize, stopType, thresh, 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(thresh)
    elif stopType == STOP_COST: strStop = 'costs change < {}'.format(thresh)
    else:strStop = 'gradient_norm < {}'.format(thresh)
    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('Iterations')
    ax.set_ylabel('Cost')
    ax.set_title(name.upper() + '-Error vs. Iteration')
    return theta

运行结果

①采用随机梯度下降，设置学习率为0.001，迭代次数为5000次
可以看到损失函数没有收敛的迹象，原因可能是学习率太大，或者迭代次数太少，因此调整学习率和迭代次数。
②采用随机梯度下降，改变学习率为0.000002，迭代次数为15000次
可以看到损失函数存在波动，但是最终能够收敛。
随机梯度下降胜在速度快，但是稳定性差

③采用小批量梯度下降，批量为16，学习率为0.001，迭代次数为15000次
可以看损失函数波动很大，这里可以尝试对数据进行标准化。
④采用小批量梯度下降，批量为16，学习率为0.001，迭代次数为15000次，但是对数据进行标准化
可以看损失函数能够收敛了。所以，有一个原则：当我们得到的结果不太好，如发生浮动的时候，我们应该先从数据下手，看是否能对数据进行一系列改变之后能不能使得结果更好，先改数据、后改模型，这是基本思路！【数据预处理很重要】