逻辑回归模型原理推导及python实现

1 . 模型定义：

逻辑回归模型定义为具有如下概率分布的模型：

其中

是输入（特征空间）
是输出（标记）
为模型参数——权值向量(weight)
为模型参数——偏置(bias)

简洁表达起见，将 扩充入 ，同时 后加入一项全为1的列，使得

其中

,

几率

一个事件发生的几率 指 该事件发生的概率与不发生的概率的比值，
如果某一事件发生的概率为 , 那么它的几率为
那么该事件的对数几率函数即

因此逻辑回归模型的对数几率函数为：

上式意味着， 输出 的 对数几率 是 输入 的线性函数， 即逻辑回归模型。

2. 模型参数估计

我们的目标是从训练数据中学到模型参数， 采用极大似然估计法。
对于给定的训练数据集

为了应用到逻辑回归模型中， 用于估计模型参数 ，需要首先假设

,

则似然函数为

其中

为 训练数据大小

上式不能直接对 w 求导， 转为对数似然函数：

目标是求得 的最大值，对数似然函数 对 求偏导

根据凸优化理论，可以利用梯度上升法求 的最大值， 即

为步长， 也称为学习率

4. Python 实现代码

import time
import numpy as np
'''
数据集: MNIST
训练集大小: 6000
测试集大小: 1000
下载地址：https://pjreddie.com/projects/mnist-in-csv/
-----
运行结果：
start read trainSet
start read testSet
start to train
start to test
the accuracy is : 0.9704
time span: 44.21207141876221
'''

def loadData(fileName):
    '''
    加载数据集
    :param fileName: 文件路径
    :return: dataList，labelList分别为特征集X和标记Y. 均为list

    '''
    dataList = [] 
    labelList = []

    f = open(fileName, 'r')

    for line in f.readlines():
        curline = line.strip().split(',')
        '''
        这里考虑到我用的文件是csv格式，所以用split(',')
        
        Mnsit有0-9十个标记
        文件每行开头第一个数字为该行的label标记
        这里为了简化，限定二分类任务，所以将标记0的作为1(正例)，其余为0(反例)
        '''
        if (int(curline[0]) == 0):
            labelList.append(1)
        else:
            labelList.append(0)

        '''
        加入特征集X
        由于开头第一个数字为标记，故从下标 1 开始
        这里转为int类型
        /255 是为了归一化，有效减少数字爆炸。
        '''
        dataList.append([int(num)/255 for num in curline[1:]])

    #读取完毕关闭文件
    f.close()

    #返回训练集X_train, y_train
    return dataList, labelList


def predict(x, w):
    '''
    预测新数据的标记
    :param x: 用于预测的样本,为matrix[1 * m]
    :param w: 训练后的的到的w ,为matrix[m * 1]
    :return: 预测结果——标记1或0
    '''
    prob = sigmoid(np.sum(x * w)) 
    P1 = sigmoid(prob)
    if P1 > 0.5:
        return 1.0
    return 0.0

def sigmoid(inX):
    '''
    sigmoid函数
    '''
    return 1.0 / (1 + np.exp(-inX))

def LR_train(X_train, y_train, alpha, iter=40):
    '''
    逻辑回归训练
    :param X_train: 训练集X
    :param y_train: 标记Y
    :param iter:迭代次数
    :return: 模型参数 w
    '''

    #这里根据前面提到的合并将(w·x +b ) 变换为(w·x)， 训练集特征X后面增加一列全为1 的项
    for i in range(len(X_train)):
        X_train[i].append(1)

    #使用矩阵可以提高运算速度
    X_train = np.matrix(X_train)
    y_train = np.matrix(y_train).T

    #初始化w = {w1,w2,...,wn}全为0
    w = np.zeros(X_train.shape[1])
    w = np.matrix(w).T #转换后w大小为(n, 1)


    '''
    梯度上升法求解 w 的最大值
    公式在前文部分已经给出，原本该利用xi 计算一点的梯度来更新 wi，
    但这样计算时间大大增加，因为要循环跑完整个数据集，才能更新整个w，
    这样下来 就是双层循环，外层 为迭代次数，内层为数据集的大小，
    本例中的开销大概是 20 * 6000，需要很长时间。 大约需要半小时
    
    这里参考了《机器学习实战》 p86处的矩阵算法
    直接利用整个数据集计算梯度来更新w
    时间开销最终不到1min.
    '''
    for i in range(iter):
        error = y_train - sigmoid(X_train * w)
        w += alpha * X_train.T * error

    #训练完毕，返回w
    return w


def test(X_test, y_test, w):
    '''
    验证测试数据集
    :param X_test: 测试集特征
    :param y_test: 测试集标记
    :param w: 训练后得到的w
    :return: 准确率
    '''

    '''
    这里对测试集特征空间进行扩充处理，同前面训练集一样，增加一列全为1的项
    '''
    for i in range(len(X_test)):
        X_test[i].append(1)

    #错误计数
    errorCount = 0

    #后续发现这里可以优化，采用矩阵乘法
    #考虑到测试集不大，有时间再改
    for i in range(len(X_test)):
        if y_test[i] != int(predict(np.matrix(X_test[i]), w)):
            errorCount += 1

    #返回在测试集上的准确率
    return 1 - errorCount / len(X_test)


if __name__ == '__main__':
    start = time.time()

    print('read trainSet')
    trainData, trainLabel = loadData('D:/PythonLearn/MLA/mnist_train.csv')

    print('read testSet')
    testData, testLabel = loadData('D:/PythonLearn/MLA/mnist_test.csv')

    print('train')
    w = LR_train(trainData, trainLabel, alpha=0.001)

    print('testing')

    accuracy = test(testData, testLabel, w)
    end = time.time()

    print('the accuracy is :', accuracy)
    print('time span:', end - start, 's')