python逻辑模型

逻辑模型"/>

python逻辑模型

目录

提示：以下是本篇文章正文内容，下面案例可供参考

logistic回归(Logistic Regression)是一种广义线性回归(Generalized Linear Model),在机器学习中是最常见的一种用于二分类的算法模型，由于数学原理简单，方便理解，作用高效，其实际运用相当广泛。为了通过自变量的线性组合来预测类别因变量的取值，logistic回归模型应运而生。logistic回归的因变量可以是二分类的，也可以是多分类的，但是二分类的更为常用，也更加容易解释，多类可以使用softmax方法进行处理。实际中最为常用的就是二分类的logistic回归。虽然带有回归二字，但实则是分类模型，下面从最基础的logit变换开始理解。对于研究因变量(结果)与引发其变化的因素自变量(因素)的关系时，我们想到最基础的方法就是建立因变量与自变量的多元线性关系：

其中

 为模型的参数，可视为该一变量的权重。如果因变量是数值型的话，可以解释为不同变量

变化导致结果Y因此而变化多少。如果因变量

是用来刻画结果是否发生或者更一般的来刻画某特定结果发生的概率的情况，一个自变量

的改变可能会让Y结果改变的微乎其微，有时候甚至忽略不计。然而实际生活中，我们知道某些关键因素会直接导致某一结果的发生，例如最典型的算法AHP中影响去不同的景点的因素权重不同导致选取方向发生改变几率大不相同。我们需要让不显著的线性关系变得显著，使得模型能够很好解释随因素的变化，结果也会发生较显著的变化，这时候，人们想到了logit变换，下图是对数函数图像：

从对数函数的图像来看，其在

之间的因变量的变化是很迅速的，也就是说自变量的微小变化会导致因变量的巨大变化，这就符合了之前想要的效果。于是，对因变量进行对数变换，右边依然保持线性关系，有下面式子:

虽然上式解决了因变量随自变量变化的敏感性问题，同时也约束了y的取值范围为( 0 , + ∞ )。我们知道概率是用来描述某件事发生的可能性，一件事情发生与否，更应该是调和对称的，也就是说该事件发生与不发生有对立性，结果可以走向必然发生(概率为1)，也可以走向必然不发生(概率为0)，概率的取值范围为( 0 , 1 ) ，而等式左边 y 的取值范围是( 0 , + ∞ )，所以需要进一步压缩，又引进了几率。

三、几率

几率和机率是两个不同的词，前者是指事件发生的概率与不发生的概率之比，而后者只是指事情发生的概率。

假设事件 A 发生的概率为

，不发生的概率为

,那么事件 A 的几率为

几率恰好反应了某一事件两个对立面，具有很好的对称性，下面我们再来看一下概率和几率的关系:

首先，我们看到概率从0.01不断增大到 0.99，几率也从0.01随之不断变大到99，两者具有很好的正相关系，我们再对 p  向两端取极限有:

于是，几率的取值范围就在( 0 ， + ∞ )，这符合我们之前的因变量取值范围的假设。

四、Logistic模型

从上述推到可以发现，我们可以利用几率来代替Y结果，二者范围都为(0，+∞)。且几率与概率密切对等。这样既能满足结果对特定因素的敏感性，又能满足对称性 。于是便可推导出：

现在对这个公式进行化简：让等式左边对数变成自然对数

,等式右边改为向量乘积的形式，便有

其中

,

,我们解得

其中

是自然常数，保留5位小数是2.71828。这就是我们常见的logistic模型表达式，作出其函数图像如下

import matplotlib.pyplot as plt

import math

import numpy as np

def sigmoid(z):

y=[]

for i in z:

y.append(1/(1+math.exp(-i)))

return y

if __name__ == '__main__':

x=np.arange(-10,10)

y=sigmoid(x)

ax=plt.gca()

ax.spines['right'].set_color('none')

ax.spines['top'].set_color('none')

ax.spines['bottom'].set_position(('data', 0))

ax.spines['left'].set_position(('data', 0))

plt.plot(x,y)

plt.show()

我们看到logistic函数图像是一条S型曲线，又名sigmoid曲线，以 ( 0 , 0.5 ) 为对称中心，随着自变量 x 不断增大，其函数值不断增大接近1，随自变量 x 不断减小，其函数值不断降低接近0，函数的取值范围在 ( 0 , 1 ) 之间，且函数曲线在中心位置变化速度最快，在两端的变化速率较慢。

从上面的操作，我们可以看到逻辑回归模型从最初的线性回归模型基础上对因变量进行 logit 变换，使得因变量对自变量显著，同时约束因变量取值范围为0到正无穷大，然后用概率表示几率，最后求出概率关于自变量的表达式，把线性回归的结果压缩在 ( 0 , 1 ) 范围内，这样最后计算出的结果是一个0到1之间的概率值，表示某事件发生的可能性大小，可以做概率建模，这也是为什么逻辑回归叫逻辑回归，而不叫逻辑分类了。

五、基于最优化方法的最佳回归系数确定

根据上述推导我们已知

，其中的向量x是分类器的输入数据，向量

也就是我们要找的最佳参数系数，从而使得分类器尽可能地精确。为了寻找最佳参数，需要用到最优化理论的一些知识。

5.1梯度上升算法

5.1.1梯度

梯度上升算法基于的思想是：要找到函数的最大值，最好的方法是沿着该函数的梯度方向探寻。对此我们需要明白代表函数变化快慢的导数以及偏导数，在此基础上我们不仅要知道函数在坐标轴正方向上的变化率(即偏导数)，而且还要设法求得函数在其他特定方向上的变化率。而方向导数就是函数在其他特定方向上的变化率。梯度与方向导数有一定的关联性，在微积分里面，对多元函数的参数求

偏导数，把求得的各个参数的偏导数以向量的形式写出来，就是梯度。

例如我们对函数

求得它的梯度：

对其求偏导可知：

所以

。

函数在某一点的梯度是这样一个向量，它的方向与取得最大方向导数的方向一致，而它的模为方向导数的最大值。

5.1.2使用梯度上升找到最佳参数

可以假设为爬山运动，我们总是往向着山顶的方向攀爬，当爬到一定角度以后也会驻足停留下观察自身角度是否是朝着山顶的角度上攀爬。并且我们需要总是指向攀爬速度最快的方向爬。

关于梯度上升的几个概念：

)步长(learning rate)：步长决定了在梯度下降迭代过程中，每一步沿梯度负方向前进的长度

2)特征(feature)：指的是样本中输入部门，比如样本(x0，y0)，(x1，y1)，则样本特征为x，样本输出为y

3)假设函数(hypothesis function)：在监督学习中，为了拟合输入样本，而使用的假设函数，记为

。比如对于样本

(

,

)(i=1,2,...n),可以采用拟合函数如下：

4)损失函数(loss function)：为了评估模型拟合的好坏，通常用损失函数来度量拟合的程度。损失函数极小化，意味着拟合程度最好，对应的模型参数即为最优参数。在线性回归中，损失函数通常为样本输出和假设函数的差取平方。比如对于样本(

,

)(i=1,2,...n),采用线性回归，损失函数为：

其中

​表示样本特征x的第i个元素，

表示样本输出y的第i个元素，

为假设函数。

梯度上升算法的基本思想是：要找到某函数的最大值，最好的方法就是沿着该函数的梯度方向搜寻。我们假设步长为

，用向量来表示的话，梯度上升算法的迭代公式如下：

。该公式停止的条件是迭代次数达到某个指定值或者算法达到某个允许的误差范围。

关于梯度算法将直接用代码演示，这里不作详细解释，在之后会详细说明这一算法，先直接进入实战：

这里我们可以看看数据分布，制作散点图：

import matplotlib.pyplot as plt

from numpy import *

def loadData():

dataMat=[]

labelMat=[]

fr = open('sample1.txt')

for line in fr.readlines():

lineArr = line.strip().split()

dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])

labelMat.append(int(lineArr[2]))

# print(dataMat)

# print(labelMat)

return dataMat, labelMat

def plotDot():

dataMat,labelMat=loadData()

dataArr = array(dataMat)

n = shape(dataArr)[0]

x0=[]

x1=[]

y0=[]

y1=[]

for i in range(n):

if labelMat[i]==1:

x1.append(dataMat[i][1])

y1.append(dataMat[i][2])

else:

x0.append(dataMat[i][1])

y0.append(dataMat[i][2])

print(x0)

print(x1)

fig = plt.figure()

ax = fig.add_subplot(111)#分割为1行1列第1块

ax.scatter(x0, y0, s=30, c='red', marker='s')

ax.scatter(x1, y1, s=30, c='green')

plt.show()

plotDot()

数据集中一共有100个点，每个点包含两个数值型特征：X1和X2。因此可以将数据在一个二维平面上展示出来。我们可以将第一列数据(X1)看作x轴上的值，第二列数据(X2)看作y轴上的值。而最后一列数据即为分类标签。根据标签的不同，对这些点进行分类。

梯度上升法的伪代码：

每个回归系数初始化为1

重复R次：

计算整个数据集的梯度

使用alpha * gradient更新回归系数的向量

返回回归系数

代码实现：

from numpy import *

#预处理数据

def loadDataSet():

#创建两个列表

dataMat=[];labelMat=[]

#打开文本数据集

fr=open('sample1.txt')

#遍历文本的每一行

for line in fr.readlines():

#对当前行去除首尾空格，并按空格进行分离

lineArr=line.strip().split()

#将每一行的两个特征x1，x2，加上x0=1,组成列表并添加到数据集列表中

dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])

#将当前行标签添加到标签列表

labelMat.append(int(lineArr[2]))

#返回数据列表，标签列表

return dataMat,labelMat

#定义sigmoid函数

def sigmoid(inx):

return 1.0/(1+exp(-inx))

#梯度上升法更新最优拟合参数

#@dataMatIn：数据集

#@classLabels：数据标签

def gradAscent(dataMatIn,classLabels):

#将数据集列表转为Numpy矩阵

dataMatrix=mat(dataMatIn)

#将数据集标签列表转为Numpy矩阵，并转置

labelMat=mat(classLabels).transpose()

#获取数据集矩阵的行数和列数

m,n=shape(dataMatrix)

#学习步长

alpha=0.001

#最大迭代次数

maxCycles=500

#初始化权值参数向量每个维度均为1.0

weights=ones((n,1))

print(weights)

#循环迭代次数

for k in range(maxCycles):

#求当前的sigmoid函数预测概率

h=sigmoid(dataMatrix*weights)

#***********************************************

#此处计算真实类别和预测类别的差值

#对logistic回归函数的对数释然函数的参数项求偏导

error=(labelMat-h)

#更新权值参数

weights=weights+alpha*dataMatrix.transpose()*error

#***********************************************

return weights

dataMat,labelMat=loadDataSet()

m=gradAscent(dataMat,labelMat)

print(m)

我们知道对回归系数进行更新的公式为：

,其中

是对参数w求偏导数。则我们可以通过求导验证logistic回归函数对参数w的梯度为(

-sigmoid(

))*

。

5.2、画出决策边界

def plotBestFit(wei):

import matplotlib.pyplot as plt

weights = wei.getA()

dataMat, labelMat = loadDataSet()

dataArr = array(dataMat)

n = shape(dataArr)[0]

xcord1 = []; ycord1 = []

xcord2 = []; ycord2 = []

for i in range(n):

if int(labelMat[i]) == 1:

xcord1.append(dataArr[i,1]); ycord1.append(dataArr[i,2])

else: xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2])

fig = plt.figure()

ax = fig.add_subplot(111)

ax.scatter(xcord1, ycord1, s = 30, c = 'red', marker='s')

ax.scatter(xcord2, ycord2, s = 30, c = 'green')

x = arange(-3.0, 3.0, 0.1)

y = (-weights[0]- weights[1]*x)/weights[2]

ax.plot(x, y)

plt.xlabel('X1');

plt.ylabel('X2');

plt.show()

plotBestFit(m)

在拟合分类线中，我们设定sigmoid函数为0。0是两个分类的分解处，因此我们设定

,然后根据函数制定x坐标表示为x1，y坐标表示为x2，画出分类线。尽管分类结果不错，但是这个方法却需要大量计算。

5.3、  随机梯度上升

梯度上升算法在每次更新回归系数时需要遍历整个数据集，该方法在处理100个左右的数据集时尚可，但如果有数十亿样本和成千上万的特征，那么该方法就很费劲了。所有我们思考了优化，一次仅用一个样本点来更新回归系数，该方法称为随机梯度上升法。对应的更新公式是：

由于可以在新的样本到来时对分类器进行增量式更新，因而随机梯度上升算法是一个在线学习算法。与之对应的是批处理算法。

随机梯度上升伪代码：

所有回归系数初始化为1

对数据集中每个样本

计算该样本的梯度

使用alpha*gradient更新回归系数值

返回回归系数值

代码：

def stocGradAscent0(dataMatrix,classLabels):

m,n=shape(dataMatrix)

alpha=0.01

weights=ones(n)

dataMatrix=array(dataMat)

for i in range(m):

h=sigmoid(sum([dataMatrix[i]*weights]))

error = classLabels[i]-h

weights=weights+alpha*error*dataMatrix[i]

print(type(dataMatrix[i]))

return weights

随机梯度上升法和批量梯度上升法是两个极端，一个采用所有数据来梯度上升，一个用一个样本来梯度上升。自然各自的优缺点都非常突出。对于训练速度来说，随机梯度上升法由于每次仅仅采用一个样本来迭代，训练速度很快，而批量梯度上升法在样本量很大的时候，训练速度不能让人满意。对于准确度来说，随机梯度上升法用于仅仅用一个样本决定梯度方向，导致解很有可能不是最优。对于收敛速度来说，由于随机梯度上升法一次迭代一个样本，导致迭代方向变化很大，不能很快的收敛到局部最优解。 但值得一提的是，随机梯度上升法在处理非凸函数优化的过程当中有非常好的表现，由于其上升方向具有一定随机性，因此能很好的绕开局部最优解，从而逼近全局最优解。

一个判断优化算法的可靠办法是看它是否收敛，也就是说参数是否达到了稳定值，是否还会变化。对此观察随机梯度上升算法的3个参数的变化：

def stocGradAscent0(dataMatrix,classLabels):

m,n=shape(dataMatrix)

print(n)

x1 = arange(0, 142)

y1=[]

y2=[]

y3=[]

alpha=0.01

weights=ones(n)

dataMatrix=array(dataMat)

for i in range(142):

h=sigmoid(sum([dataMatrix[i]*weights]))

error = classLabels[i]-h

weights=weights+alpha*error*dataMatrix[i]

y1.append(weights[0])

y2.append(weights[1])

y3.append(weights[2])

fig=plt.figure(figsize=(14,14),dpi=100)

plt.subplot(3, 1, 1)

plt.plot(x1,y1)

plt.subplot(3, 1, 2)

plt.plot(x1, y2)

plt.subplot(3, 1, 3)

plt.plot(x1, y3)

# ax = fig.add_subplot(111)

# bx = fig.add_subplot(211)

# cx = fig.add_subplot(311)

# ax.plot(x1, y1)

# bx.plot(x1,y2)

# cx.plot(x1,y3)

plt.show()

return weights

我们发现第三个参数仅仅在50次迭代就达到了稳定值，而其他两个则需要更多次迭代。产生这种现象的原因是存在一些不能正确分类的样本点(数据集并非线性可分)，在每次迭代时会引发系数的剧烈变化。

5.4、优化随机梯度上升

我们对随机梯度上升算法做一些调整：

def stocGradAscent0(dataMatrix,classLabels,numIter=150):

m,n=shape(dataMatrix)

alpha=0.01

weights=ones(n)

dataMatrix=array(dataMat)

for j in range(150):

dataIndex = range(m)

dataIndex=list(dataIndex)

for i in range (m):

#alpha每次迭代时调整

alpha=4/(1.0+j+i)+0.1

#随机选取更新

randIndex=int(random.uniform(0,len(dataIndex)))

h=sigmoid(sum([dataMatrix[randIndex]*weights]))

error = classLabels[randIndex]-h

weights=weights+alpha*error*dataMatrix[randIndex]

del(dataIndex[randIndex])

return weights

该分割线达到了与gradAscent()差不多的效果，而且计算量明显少于批量梯度上升算法。

六、实战：从疝气病症预测病马的死亡率

我们将使用Logistic回归来预测患有疝气病的马存活问题。此数据集包含368个样本和28个特征。另外需要说明的是这些数据集是有部分缺失，我们需要先进行数据预处理再利用回归预测。

对于数据集中缺失数据我们可以采取以下办法进行处理：

使用可用特征的均值来填补缺失值；

使用特殊值来填补缺失值，比如-1；

忽略有缺失值的样本；

使用相似样本的均值添补缺失值；

‘使用另外的机器学习算法预测缺失值；(如KNN)

而对于该数据集我们选择用0代替缺失值，因为根据回归系数的更新公式：

当一数据集特征为0时，

，将不作更新。

另外由于sigmoid(0)=0.5，它对结果的预测不具有任何倾向性，因此上述做法也不会对误差项造成任何影响。

现在我们根据上述所做的工作，现在只需要将特征向量乘以最优化系数得到的值全部相加，作为输入给sigmoid函数中。如果sigmoid值大于0.5则给定标签为1，否则为0。

def colicTest():

frTrain = open('horseColicTraining.txt', encoding='utf8');

frTest = open('horseColicTest.txt', encoding='utf8')

trainingSet = []

trainingLabels = []

print(frTest)

print(frTrain)

for line in frTrain.readlines():

currLine = line.strip().split('\t')

lineArr = []

for i in range(21):

lineArr.append(float(currLine[i]))

trainingSet.append(lineArr)

trainingLabels.append(float(currLine[21]))

trainWeights = stocGradAscent0(array(trainingSet), trainingLabels, 1000)

errorCount = 0

numTestVec = 0.0

for line in frTest.readlines():

numTestVec += 1.0

currLine = line.strip().split('\t')

lineArr = []

for i in range(21):

lineArr.append(float(currLine[i]))

if int(classifyVector(array(lineArr), trainWeights)) != int(currLine[21]):

errorCount += 1

errorRate = (float(errorCount) / numTestVec)

print("该模型错误率为: %f" % errorRate)

return errorRate

由于每次是随机取数，所以我们可将colicTest()迭代次数更多一些：

def multiTest():

numTests = 10

errorSum = 0.0

for k in range(numTests):

errorSum += colicTest()

print("该模型迭代%d次平均错误率为: %f",numTests,(numTests, errorSum / float(numTests)))

运行输出：

可见错误率还是很低，因为数据存在30%的丢失。

Logistic回归利用Sigmoid函数作为最佳拟合函数，利用最优化算法求出参数，再把数据集特征向量矩阵与对应训练好的参数相乘再相加，作为Sigmoid的输入从而得到预测。

参阅：