logistic回归模型python实现_逻辑回归模型原理推导及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')