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')
