【机器学习】logistic回归原理分析及python实现
1.sigmoid函数和logistic回归分类器
2.梯度上升最优化算法
3.数据中的缺失项处理
4.logistic实现马疝气病预测
首先阐述logistic回归的定义，然后介绍一些最优化算法，其中包括基本的梯度上升算法和改进的随机梯度上升算法，这些最优化算法用于分类器的训练，最后给出logistic回归实例，预测一匹有疝气病的马是否被治愈（二分类）。
一.sigmoid函数和logistic回归分类器
1.什么是回归？
假设现在有一些数据点，我们用一条直线来对这些点进行拟合（该线称为最佳拟合直线），这个拟合过程就称做回归。
2.sigmoid函数
我们想要的函数应该是，接受所有的输入，然后预测出类别。例如，在二类情况下，输出0和1。像单位阶跃函数就可以实现，但该函数在跳跃点上从0->1是瞬间跳跃，这个瞬间跳跃过程有时很难处理。于是另一个有类似性质且数学上更易处理的函数出现了-----sigmoid函数。
sigmoid函数表达式：
实现：
"""
函数：sigmoid函数
"""
def sigmoid(z):
return 1.0/(1+np.exp(-z) )3.sigmoid函数如何用于二分类？
为了实现logistic回归分类器，可以在每个特征上都乘以一个回归系数w，然后把所有的结果值相加得到z值，将这个z值带入sigmoid函数中，会输出一个在【0，1】内的数值。分类：z>0.5，输出1；z<0.5，输出0。
1）输入样本： X=(x0,x1……xn)
2）如何将样本值转化为sigmoid的输入？x-> z
相应的回归系数W=(w0,w1……wn)，样本特征值与相应系数相乘求和：
3)带入sigmoid函数：
4.logistic 回归的优点与缺点？
优点：计算代价不高，易于理解与实现（简单）。
缺点：容易欠拟合，分类精度可能不高。
适用数据类型：数值型和标称型数据。
5.logistic回归的一般过程？
1、收集数据：任何方式
2、准备数据：由于要计算距离，因此要求数据都是数值型的，另外结构化数据格式最佳。
3、分析数据：采用任一方是对数据进行分析
4、训练算法：大部分时间将用于训练，训练的目的为了找到最佳的分类回归系数
5、测试算法：一旦训练步骤完成，分类将会很快
6、使用算法：首先，我们需要输入一些数据，并将其转化成对应的结构化数值；接着基于训练好的回归系数就可以对这些数值进行简单的回归计算，判定它们属于哪一类别；在这之后，我们就可以在输出的类别上做一些其他的分析工作。
二.梯度上升最优化算法
1.如何确定最佳回归系数？
确定了分类器的函数形式之后，现在的问题变成了：最佳回归系数W=(w0,w1……wn)是多少？如何确定它们的大小？
为了寻找最佳参数，需要用到最优化理论的知识。
1.1梯度上升法
梯度上升法的基本思想是：要找到某函数的最大值，最好的方法是沿着该函数的梯度方向去探寻。
函数f(x,y)的梯度记为：
这个梯度意味着f沿着x方向移动
梯度上升算法公式：
梯度上升算法用来求函数的最大值，梯度下降算法（将上式的+改为 - 号）用来求函数的最小值。
何时停止迭代？
该公式一直被迭代，直到达到某个停止条件为止，比如迭代次数达到某指定值或者误差达到允许的范围内。
梯度上升算法的python实现：
"""
函数：梯度上升算法
"""
def gradAscent(dataMat,labelMat):
dataSet=np.mat(dataMat) # m*n
labelSet=np.mat(labelMat).transpose() # 1*m->m*1
m,n=np.shape(dataSet) # m*n: m个样本，n个特征
alpha=0.001 # 学习步长
maxCycles=500 # 最大迭代次数
weights=np.ones( (n,1) )
for i in range(maxCycles):
y=sigmoid(dataSet * weights) # 预测值
error=labelSet -y
weights=weights+ alpha *dataSet.transpose()*error
#print(type(weights))
return weights.getA(),weights ##getA():将Mat转化为ndarray,因为mat不能用index画出决策边界：
"""
函数：画出决策边界
"""
def plotBestFit(weight):
dataMat, labelMat = loadDataSet()
dataArr=np.array(dataMat)
m,n=np.shape(dataArr)
x1=[] #x1,y1:类别为1的特征
x2=[] #x2,y2:类别为2的特征
y1=[]
y2=[]
for i in range(m):
if (labelMat[i])==1:
x1.append(dataArr[i,1])
y1.append(dataArr[i,2])
else:
x2.append(dataArr[i,1])
y2.append(dataArr[i,2])
fig=plt.figure()
ax=fig.add_subplot(111)
ax.scatter(x1,y1,s=30,c='red',marker='s')
ax.scatter(x2,y2,s=30,c='green')
#画出拟合直线
x=np.arange(-3.0, 3.0, 0.1)
y=(-weights[0]-weights[1]*x)/weights[2] #直线满足关系：0=w0*1.0+w1*x1+w2*x2
ax.plot(x,y)
plt.xlabel('X1')
plt.ylabel('X2')
plt.show()梯度上升算法的实际效果：
分析：
梯度上升算法的分类效果挺好，但这个方法需要大量的计算。因为每次更新回归系数时都需要遍历整个数据集，该方法在处理小样本、小数据集时尚可，但如果有数十亿样本和成千上万的特征，那么该方法的计算复杂度就太高了。下面将介绍一种改进方法----随机梯度上升。
1.2一种改进方法----随机梯度上升
随机梯度上升，一次仅用一个样本点来更新回归系数。
实现：
"""
函数：随机梯度上升算法0.0
改进：每次用一个样本来更新回归系数
"""
def stocGradAscent0(dataMat,labelMat):
m, n = np.shape(dataMat) # m*n: m个样本，n个特征
alpha = 0.001 # 学习步长
maxCycles=500
weights = np.ones(n)
for cycle in range(maxCycles):
for i in range(m):
y = sigmoid(sum(dataMat[i] * weights) ) # 预测值
error = labelMat[i] - y
weights = weights + alpha * error* dataMat[i]
# print(type(weights))
return weights
梯度上升算法的实际效果：
分析：
一个判断优化算法优劣的可靠方法是看它是否收敛，即参数是否达到了稳定值。在随机上升算法中，在500次迭代后差不多达到稳定，即收敛。
1.3改进的随机梯度上升法
def stocGradAscent1(dataMat,labelMat):
m, n = np.shape(dataMat) # m*n: m个样本，n个特征
maxCycles = 150
weights = np.ones(n)
for cycle in range(maxCycles):
dataIndex=list( range(m))
for i in range(m):
alpha = 4 / (1.0 + cycle + i) + 0.01 # 学习步长
randIndex=int(np.random.uniform(0,len(dataIndex) )) #随机选取样本
y = sigmoid(sum(dataMat[randIndex] * weights )) # 预测值
error = labelMat[randIndex] - y
weights = weights + alpha * error * dataMat[i]
del(dataIndex[randIndex])
# print(type(weights))
return weights
默认150次迭代，结果如下：
2.可视化：迭代次数与回归系数的关系
此部分参考：https://blog.csdn.net/c406495762/article/details/77851973（感谢好帖！）
分别展示梯度上升法、随机梯度上升法、改进的随机梯度上升法的迭代次数与回归次数的关系如下：
疑问：实现的迭代过程中发现前两种方法未实现收敛，但从分类结果来看，是挺好的，这是为什么呢？若有人复现出与我不同的结果，请留言或私信告知，感谢！
代码如下：
# -*- coding:utf-8 -*-
"""
@author:Lisa
@file:visual_iteration_process.py
@note:可视化迭代次数与回归系数的关系
@time:2018/7/12 0012下午 8:08
"""
import numpy as np
import matplotlib.pyplot as plt
"""
函数：加载数据集
"""
def loadDataSet():
dataMat=[] #列表list
labelMat=[]
txt=open('testSet.txt')
for line in txt.readlines():
lineArr=line.strip().split() #strip():返回一个带前导和尾随空格的字符串的副本
#split():默认以空格为分隔符，空字符串从结果中删除
dataMat.append( [1.0, float(lineArr[0]), float(lineArr[1]) ]) #将二维特征扩展到三维，第一维都设置为1.0
labelMat.append(int(lineArr[2]) )
return dataMat,labelMat
"""
函数：sigmoid函数
"""
def sigmoid(z):
return 1.0/(1+np.exp(-z) )
"""
函数：梯度上升算法
"""
def gradAscent(dataMat,labelMat):
dataSet=np.mat(dataMat) # m*n
labelSet=np.mat(labelMat).transpose() # 1*m->m*1
m,n=np.shape(dataSet) # m*n: m个样本，n个特征
alpha=0.001 # 学习步长
maxCycles=500 # 最大迭代次数
weights=np.ones((n,1))
weights_array=np.array([])
for i in range(maxCycles):
y=sigmoid(dataSet * weights) # 预测值
error=labelSet -y
weights=weights+ alpha *dataSet.transpose()*error
weights_array=np.append(weights_array,weights)
weights_array = weights_array.reshape(maxCycles, n)
#print(weights_array)
return weights.getA(),weights_array ##getA():将Mat转化为ndarray,因为mat不能用index
"""
函数：随机梯度上升算法0.0
改进：每次用一个样本来更新回归系数
"""
def stocGradAscent0(dataMat,labelMat):
m, n = np.shape(dataMat) # m*n: m个样本，n个特征
alpha = 0.001 # 学习步长
maxCycles=500
weights = np.ones(n)
weights_array = np.array([])
for cycle in range(maxCycles):
for i in range(m):
y = sigmoid(sum(dataMat[i] * weights) ) # 预测值
error = labelMat[i] - y
weights = weights + alpha * error* dataMat[i]
weights_array = np.append(weights_array, weights)
# print(type(weights))
weights_array = weights_array.reshape(maxCycles*m, n)
return weights,weights_array
"""
函数：改进的随机梯度上升法1.0
改进：1.alpha随着迭代次数不断减小，但永远不会减小到0
2.通过随机选取样本来更新回归系数
"""
def stocGradAscent1(dataMat,labelMat):
m, n = np.shape(dataMat) # m*n: m个样本，n个特征
maxCycles = 500
weights = np.ones(n)
weights_array=np.array([])
for cycle in range(maxCycles):
dataIndex=list(range(m))
for i in range(m):
alpha = 4 / (1.0 + cycle + i) + 0.01 # 学习步长
randIndex=int(np.random.uniform(0,len(dataIndex) )) #随机选取样本
y = sigmoid(sum(dataMat[randIndex] * weights )) # 预测值
error = labelMat[randIndex] - y
weights = weights + alpha * error * dataMat[randIndex]
del(dataIndex[randIndex])
# print(type(weights))
weights_array = np.append(weights_array, weights)
weights_array = weights_array.reshape(maxCycles*m, n)
return weights,weights_array
"""
函数：可视化迭代过程
"""
def visual(weights_array1,weights_array2,weights_array3):
import matplotlib as mpl
# 指定默认字体
mpl.rcParams['font.sans-serif'] = ['SimHei'] #显示汉字
mpl.rcParams['axes.unicode_minus'] = False # 能够显示符号（负号）
fig,axs=plt.subplots(3,3,figsize=(15,10))
x1=np.arange(0,len(weights_array1),1)
#绘制w0和迭代次数的关系
axs[0][0].plot(x1,weights_array1[:,0])
axs[0][0].set_title('梯度上升算法：回归系数与迭代次数的关系')
axs[0][0].set_ylabel('w0',)
# 绘制w1和迭代次数的关系
axs[1][0].plot(x1, weights_array1[:, 1])
axs[1][0].set_ylabel('w1')
# 绘制w2和迭代次数的关系
axs[2][0].plot(x1, weights_array1[:, 2])
axs[2][0].set_ylabel('w2')
axs[2][0].set_xlabel('迭代次数')
x2 = np.arange(0, len(weights_array2), 1)
# 绘制w0和迭代次数的关系
axs[0][1].plot(x2, weights_array2[:, 0])
axs[0][1].set_title('随机梯度上升算法')
axs[0][1].set_ylabel('w0')
# 绘制w1和迭代次数的关系
axs[1][1].plot(x2, weights_array2[:, 1])
axs[1][1].set_ylabel('w1')
# 绘制w2和迭代次数的关系
axs[2][1].plot(x2, weights_array2[:, 2])
axs[2][1].set_ylabel('w2')
axs[2][1].set_xlabel('迭代次数')
x3 = np.arange(0, len(weights_array3), 1)
# 绘制w0和迭代次数的关系
axs[0][2].plot(x3, weights_array3[:, 0])
axs[0][2].set_title('改进的随机梯度上升算法')
axs[0][2].set_ylabel('w0')
# 绘制w1和迭代次数的关系
axs[1][2].plot(x3, weights_array3[:, 1])
axs[1][2].set_ylabel('w1')
# 绘制w2和迭代次数的关系
axs[2][2].plot(x3, weights_array3[:, 2])
axs[2][2].set_ylabel('w2')
axs[2][2].set_xlabel('迭代次数')
plt.show()
if __name__ == '__main__':
dataMat, labelMat = loadDataSet()
weights1, weights_array1 = gradAscent(dataMat, labelMat)
weights2, weights_array2 = stocGradAscent0(np.array(dataMat), labelMat)
weights3, weights_array3 = stocGradAscent1(np.array(dataMat), labelMat)
visual(weights_array1, weights_array2,weights_array3)三.数据中的缺失项处理
数据中的缺失值是一个非常棘手的问题，很多文献都致力于解决这个问题。那么，数据缺失究竟带来了什么问题？假设有100个样本和20个特征，这些数据都是机器收集回来的。若机器上的某个传感器损坏导致一个特征无效时该怎么办？它们是否还可用？答案是肯定的。因为有时候数据相当昂贵，扔掉和重新获取都是不可取的，所以必须采用一些方法来解决这个问题。下面给出了一些可选的做法：
- 使用可用特征的均值来填补缺失值；
- 使用特殊值来填补缺失值，如-1；
- 忽略有缺失值的样本；
- 使用相似样本的均值添补缺失值；
- 使用另外的机器学习算法预测缺失值。
四.logistic实现马疝气病预测
1.预处理
由于该数据集存在30%的缺失，那么首先必须对数据集进行预处理。
预处理数据做两件事：
- 如果测试集中一条数据的特征值已经确实，那么我们选择实数0来替换所有缺失值，因为本文使用Logistic回归。因此这样做不会影响回归系数的值。sigmoid(0)=0.5，即它对结果的预测不具有任何倾向性。
- 如果测试集中一条数据的类别标签已经缺失，那么我们将该类别数据丢弃，因为类别标签与特征不同，很难确定采用某个合适的值来替换。
所有处理后的数据集可见：（感谢Jack-Cui大神处理的数据集）
https://github.com/LisaPig/ML_LogisticRegression
2.完整代码
# -*- coding:utf-8 -*-
"""
@author:Lisa
@file:horse_colic.py
@note:疝气马能否痊愈预测
@time:2018/7/12 0012下午 10:11
"""
import numpy as np
import logisticRegression as LR
"""
函数：sigmoid回归分类
"""
def classifyVector(dataIn,weights):
h=LR.sigmoid(sum(dataIn*weights))
if h>0.5:
return 1.0
else:
return 0.0
"""
函数：疝气预测
"""
def colicTest():
trainData=open('data\horseColicTraining.txt')
testData = open('data\horseColicTest.txt')
trainSet=[]
trainLabel=[]
for line in trainData.readlines():
curLine=line.strip().split('\t')
lineArr=[]
for i in range(21):
lineArr.append(float (curLine[i]))
trainSet.append(lineArr)
trainLabel.append(float(curLine[21]))
trainWeights=LR.stocGradAscent1(np.array(trainSet),trainLabel)
errorCount=0
numTestVec=0
testSet = []
testLabel = []
for line in testData.readlines():
numTestVec+=1
curLine=line.strip().split('\t')
lineArr=[]
for i in range(21):
lineArr.append(float (curLine[i]))
if int(classifyVector(np.array(lineArr),trainWeights))!=int(curLine[21]):
errorCount+=1
errorRate=float(errorCount/numTestVec)
print("the error rate is: %f" % errorRate)
return errorRate
def multiTest():
numTests=10
errorSum=0
for k in range(numTests):
errorSum+=colicTest()
print("after %d iterations the average error rate is: %f"
%(numTests,errorSum/float(numTests)))
if __name__ == '__main__':
multiTest()运行结果：
——————————————————— end ———————————————————————
参考：
1.《机器学习实战》
2. https://blog.csdn.net/c406495762/article/details/77851973