一、前言
logistic回归又称logistic回归分析,是一种广义的线性回归分析模型,常用于数据挖掘,疾病自动诊断,经济预测等领域。例如,探讨引发疾病的危险因素,并根据危险因素预测疾病发生的概率等。以胃癌病情分析为例,选择两组人群,一组是胃癌组,一组是非胃癌组,两组人群必定具有不同的体征与生活方式等。因此因变量就为是否胃癌,值为“是”或“否”,自变量就可以包括很多了,如年龄、性别、饮食习惯、幽门螺杆菌感染等。自变量既可以是连续的,也可以是分类的。然后通过logistic回归分析,可以得到自变量的权重,从而可以大致了解到底哪些因素是胃癌的危险因素。同时根据该权值可以根据危险因素预测一个人患癌症的可能性。
(1)收集数据:采用各种方法收集数据,比如爬虫等;
(2)准备数据:因为需要计算距离,所以数据类型应该是数值型,最好是结构化数据格式;
(3)分析数据:通过业务的角度或者其他的方法分析数据;
(4)训练算法:这是关键的一步,训练的目的是找到最佳的分类回归系数,可以使用随机梯度上升法;
(5)测试算法:训练完成,将数据投入模型进行测试;
(6)使用算法:将需要的数据进行处理成适合模型的结构化数据,输出的是类别,只有0,1两类。
二、Logistic回归函数和Sigmoid函数
事件的几率(odds),是指该事件发生的概率与该事件不发生的概率的比值。如果事件发生的概率是p,那么该事件的几率是p/(1-p)。取该事件发生几率的对数,定义为该事件的对数几率(log odds)或logit函数:
由于事件发生的概率p的取值范围为[0,1],故计算出来的几率只能是非负的。而通过取对数,便可以将输出转换到整个实数范围内。我们就可以将对数几率记为输入特征值的线性表达式 :
其中,p(y =1|x)是条件概率分布,表示当输入为x时,实例被分为1类的概率,依据此概率我们能得到事件发生的对数几率。但是,我们的初衷是做分类器,简单点说就是通过输入特征来判定该实例属于哪一类别或者属于某一类别的概率。所以我们取logit函数的反函数,令
的线性组合为输入,p为输出,经如下推导:
Φ(x)就是一个Sigmoid函数,其图像为:
当x为0时,Sigmoid函数值为0.5。随着x的增大,函数值将逼近于1,随着x的减小,函数值将逼近于0。在x=0点处,Sigmoid函数看起来很像跃阶函数。
在图像中,大于0.5的数据被分为1类,小于0.5即被归入0类。
def sigmoid(x):
return 1.0/(1.0+exp(-x))
三、读入数据
采取的数据来自集美大学acm校队校选,其中数据改变的地方为将时间单位从秒改为分钟。 以下为部分数据。第一列是得分,第二列是用时,第三列是年级(本次实验未使用),第四列是分类,标签为1的同学为进入acm校队,标签为1以外的同学则表示未进入,在此次实验中我们统一用标签0来表示。
在数据读入中,我们需要注意以下几点
1.为了便于后续计算,我将数据集returnMat第一列值设为1.0。
2.第二列和第三列分别表示特征值score和time
3.标签为1的在标签集classLabelVector上还是用1来表示,标签不是1的则用0来表示
#读入数据
def file2matrix(filename):
fr = open(filename) #打开文件
arrayOLines = fr.readlines()#读取所有行
random.shuffle(arrayOLines)#将数据打乱
numberOfLines = len(arrayOLines)
returnMat = zeros((numberOfLines,3))
classLabelVector=[]
index=0
for line in arrayOLines:
line =line.strip() #删除一行中首尾的空白符
listFromLine=line.split(' ')#字符按照空格截取
returnMat[index,0]=1.0
returnMat[index,1]=listFromLine[0]
returnMat[index,2]=listFromLine[1]
if float(listFromLine[-1])==1:
classLabelVector.append(1.0)
else:
classLabelVector.append(0.0)
index+=1
return returnMat,classLabelVector
我们可以对returnMat,classLabelVector进行输出,我们可以发现和我们思考的是一样的。
if __name__=='__main__':
print(file2matrix("data.txt"))
datingDataMat,datingLabels=file2matrix("data.txt")
print(datingDataMat)
print(datingLabels)
returnMat
classLabelVector
四、梯度上升算法求最佳回归系数
要找到某函数的最大值,最好的方法是沿着该函数的梯度方向探寻。如果梯度记为∇ ,则函数f(x,y)的梯度由下式表示:
![](https://i-blog.csdnimg.cn/blog_migrate/321872deedbdc13fa5aa22602ba396c9.png)
这个梯度表示要沿着x的方向移动
,沿着y的方向移动
。需要注意的是,函数f(x, y)在待计算点必须有定义且可微。具体如下图所示:
梯度上升算法到达每个点后都会重新估计移动的方向。从P0开始,计算完该点的梯度,函数就根据梯度移动到下一点P1。在P1点,梯度再次被重新计算,并沿新的梯度方向移动到P2.如此循环迭代,知道满足通知条件。迭代的过程中,梯度算子总是保证我们能选取到最佳的移动方向。
梯度上升算法迭代公式:
梯度下降算法迭代公式:
在这里我具体实现的是梯度上升算法,具体如下:
这个函数有两个参数,第一个参数是数据集45×3的矩阵,一共有45个样本。第二个参数是标签集
#梯度上升算法
def gradAscent(data,classLabels):
dataMat=np.mat(data)#转换为Numpy矩阵数据类型
labelMat=np.mat(classLabels).transpose()#转置矩阵
m,n=np.shape(dataMat)
print("样本数and特征数:", m, n)
alpha=0.001 #向目标移动的步长
maxCycles=500 #迭代次数
weights=np.ones((n,1)) # 表示回归系数,3*1的矩阵,值全为1
for k in range(maxCycles):
h=sigmoid(dataMat*weights) #矩阵相乘
error=(labelMat-h)
weights=weights+alpha*dataMat.transpose()*error
return weights #最佳回归系数
注意:ones函数可以用来创建任意维度和元素个数的数组,且元素值均为1;上面的代码中weights =np. ones((n, 1)) # 表示回归系数,3*1的矩阵,值全为1表示我们自行创建了一个n*1的矩阵(n是特征数,为3,其中后两个是特征,第一个是设置的X_0(值为1)。最后我们通过梯度上升算法迭代公式就能求得最佳回归系数。
if __name__=='__main__':
datingDataMat,datingLabels=file2matrix("data.txt")
a = gradAscent(datingDataMat,datingLabels)
print(a)
#plotDataSet(a)
五、画出决策边界
回顾逻辑回归分类的原理:通过训练的方式求出一个
维权值向量 w ,每当新来一个样本
时,与参数 w 进行点乘,结果带入sigmoid函数,得到的值为该样本发生(即
事件发生)的概率值。如果概率大于0.5,分类为1;否则分类为0。
当z>0时,
,因此
;当z<0时,
,因此![p< 0.5](https://i-blog.csdnimg.cn/blog_migrate/18e0a76f4bbd48f0a91a8b3de1130f32.gif)
也就是,其中有一个边界点
,大于这个边界点,分类为1,小于这个边界点。我们称之为决策边界(decision boundary)。
这是一个直线,可以将数据集分成两类。可以改写成我们熟悉的形式:
具体代码实现如下:
#画出决策边界
def plotDataSet(weights):
dataMat, labelMat =file2matrix("data.txt") # 加载数据集
dataMat = np.array(dataMat) # 转换成numpy的array数组
n = np.shape(dataMat)[0] # 数据个数,即行数
xcord1 = [] ; ycord1 = [] # 正样本
xcord2 = [] ; ycord2 = [] # 负样本
for i in range(n):
if int(labelMat[i]) == 1: #1为正样本
xcord1.append(dataMat[i,1])
ycord1.append(dataMat[i,2])
else: #0为负样本
xcord2.append(dataMat[i,1])
ycord2.append(dataMat[i,2])
fig = plt.figure()
ax = fig.add_subplot(111) #添加subplot
ax.scatter(xcord1,ycord1,s=20,c='red',marker = 's') #绘制正样本
ax.scatter(xcord2,ycord2,s=20,c='green',marker = '^') #绘制负样本
x = arange(20,250,1)
y = (-weights[0] - weights[1] * x) / weights[2] #最佳拟合线
#x = x.reshape(1,230)
y = y.reshape(230,1)
print(x)
print(y)
ax.plot(x,y)
plt.xlabel('score')
plt.ylabel('time') #绘制label
plt.show()
if __name__=='__main__':
datingDataMat,datingLabels=file2matrix("data.txt")
a = gradAscent(datingDataMat,datingLabels)
#print(a)
plotDataSet(a)
我们可以看到在大体上边界划分可以把数据分为两类,但是仍然存在一些不是1类的被划分为1类,不是0类的被划分为1类。为了改进效果,我们可以采取随机梯度上升算法。
五.可能发生的错误
1.X和Y不在同一维度![](https://i-blog.csdnimg.cn/blog_migrate/213ab1f69909ac783283560ddd86e07a.png)
通过对x和y进行输出,我发现x是一个数组,而y是一个1行230列的矩阵(具体可以自己去尝试一下)
y = y.reshape(230,1)
我将y的维度变化为230行1列,错误解决。
2.如果出现了如下错误信息:
return 1.0 / (1 + exp(-inX))
TypeError: only size-1 arrays can be converted to Python scalars
是因为在使用numpy数组的时候,我们最好也使用其中的数学函数处理数据。
解决方法:从numpy模块中导入这个exp
from numpy import *