机器学习实战-Logistic回归

1.基于 Logistic 回归和 Sigmoid 函数的分类

逻辑回归适合于01情况的分类就是描述一个问题是或者不是，所以就引入sigmoid函数，因为这个函数可以将所有值变成0-1之间的一个值，这样就方便算概率首先我们可以先看看Sigmoid函数(又叫Logistic函数)将任意的输入映射到了[0,1]区间我们在线性回归中可以得到一个预测值，再将该值映射到sigmoid函数中这样就完成了由值到概率的转换，也就是分类任务，公式如下：

整合成一个公式，就变成了如下公式：

z是一个矩阵，θ是参数列向量(要求解的)，x是样本列向量(给定的数据集)，θ^T表示θ的转置
Sigmoid函数的输入记为z，由下面公式得出:
z=w0x0+w1x1+w2x2+...+wnxn
如果采用向量的写法，上述公式可以写成z = wTx，它表示将这两个数值向量对应元素相乘然后
全部加起来即得到z值。其中的向量x是分类器的输入数据，向量w也就是我们要找到的最佳参数 (系数)，从而使得分类器尽可能地精确。
逻辑回归的简单来说，就是根据数据得到的回归直线方程z=a*x+b方程之后，将z作为sigmoid的输入使得z的值转化为在0-1之间的值，然后计算概率，最后根据sigmod函数的特点也就是当输入为零的时候，函数值为0.5，以0.5为分界线来划分数据的类型。
1.1梯度上升法

梯度上升算法用来求函数的最大值，而梯度下降算法用来求函数的最小值。
求一个函数的最大值，在数学中，是不是通过对函数求导，然后算出导数等于0，或者导数不存在的位置作为极值，然后如果极值只有一个开区间内是不是极值就是最大值。但是在实际应用中却不是这么简单的去求最大值，而是通过迭代的方式一步一步向最值点靠近，最后得到最值。这也就是梯度上升法的思想

其中，m为样本的总数，y(i)表示第i个样本的类别，x(i)表示第i个样本，需要注意的是θ是多维向量，x(i)也是多维向量。
梯度上升迭代公式为：

代码实现：

1. import matplotlib.pyplot as plt
2. import numpy as np
3. import warnings
4. warnings.filterwarnings('ignore')
6. def loadDataSet():
7.     dataMat = []   #创建数据列表
8.     labelMat = []   #创建标签列表
9.     fr = open('testSet.txt')    #打开文件
10.     for line in fr.readlines():     #逐行读取
11.         lineArr = line.strip().split()  #去回车，放入列表
12.         dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])]) #添加数据
13.         labelMat.append(int(lineArr[2]))   #添加标签
14.     fr.close()   #关闭文件
15.     return dataMat, labelMat    #返回
17. def sigmoid(inX):
18.     return  1.0/(1+np.exp(-inX))
20. #dataMatIn，它是一个2维NumPy数组，每列分别代表每个不同的特征,每行则代表每个训练样本
21. def gradAscent(dataMatIn,classLabels):
22.     dataMatrix=np.mat(dataMatIn)#转换成numpy的mat
23.     labelMat=np.mat(classLabels).transpose()#为了便于矩阵运算，需要将该行向量转换为列向量，做法是将原向量转置，再将它赋值给labelMat
24.     m,n=np.shape(dataMatrix) #返回dataMatrix的大小。m为行数,n为列数。
25.     alpha=0.001 #向目标移动的步长
26.     maxCycles=500 #maxCycles是迭代次数
27.     weights=np.ones((n,1))#权重初始化都为1
28.     for k in range(maxCycles):
29.         h=sigmoid(dataMatrix*weights)#梯度上升矢量化公式
30.         error=(labelMat-h)#相当于公式中的y(i)-h(x(i))
31.         weights=weights+alpha*dataMatrix.transpose()*error#公式里面的累加在这里使用矩阵相乘来实现(矩阵相乘的过程中先乘再加)
32.     return weights.getA()  #返回权重数组
34. def plotBestFit(weights):
35.     dataMat, labelMat = loadDataSet() #加载数据集
36.     dataArr = np.array(dataMat) #转换成numpy的array数组
37.     n = np.shape(dataMat)[0]  #数据个数
38.     xcord1 = []; ycord1 = []   #正样本
39.     xcord2 = []; ycord2 = []   #负样本
40.     for i in range(n):    #根据数据集标签进行分类
41.         if int(labelMat[i]) == 1:
42.             xcord1.append(dataArr[i,1]); ycord1.append(dataArr[i,2]) #1为正样本
43.         else:
44.             xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2])#0为负样本
45.     fig = plt.figure()
46.     ax = fig.add_subplot(111) #添加subplot
47.     ax.scatter(xcord1, ycord1, s = 20, c = 'red', marker = 's',alpha=.5)#绘制正样本
48.     ax.scatter(xcord2, ycord2, s = 20, c = 'green',alpha=.5)  #绘制负样本
49.     x = np.arange(-3.0, 3.0, 0.1)
50.     y = (-weights[0] - weights[1] * x) / weights[2]
51.     ax.plot(x, y)
52.     plt.title('BestFit')  #绘制title
53.     plt.xlabel('X1'); plt.ylabel('X2') #绘制label
54.     plt.show()
56. if __name__ == '__main__':
57.     dataMat, labelMat = loadDataSet()
58.     weights = gradAscent(dataMat, labelMat)
59.     plotBestFit(weights)
60.     #plotDataSet()

复制代码

正常运行之前报了如下错误AttributeError: partially initialized module ‘matplotlib.cbook’ has no attribute ‘deprecated’ (most likely due to a circular import)
然后就是你安装的matplotlib版本太高了，可以将以前的版本给卸载了，然后安装一个版本较低的matplotlib，然后就可以解决
测试结果

总结：代码运行的大致思路如下：

• 拿到数据之后通过loadDataSet函数将数据的坐标和类别分别存放在两个数组中
  
• 然后通过拿到的坐标数组和分类数组，首先先将数组转换为矩阵，方便后面矩阵的运算，然后根据sigmoid函数将所有数据转化为0-1之间的数据，也即公式中的h(xi)的值，然后用每个样本的类标签减去梯度上升的矢量值，最后带入公式算出当前的权重值
  
• 绘图
  
  
  
  • 然后绘图的时候定义两种数据类型的x坐标的数组和对应的y坐标的数组，然后通过数据集，将每种类别对应的x坐标和y坐标分开存入对应的数组中
    
  • 然后绘制一张空的面板，添加坐标系，然后将每个同一类别的点画在面板上，并且同一类别的点的颜色相同
    
  • 最后画拟合的直线，横坐标的范围已知，然后取sigmoid 函数为0。0是两个分类(类别1和类别0)的分界处。因此，我们设定 0 = w0x0 + w1x1 + w2x2，然后解出X2和X1的关系式(即分隔线的方程，注意X0=1)。
    
  • 最后根据函数方程和x的值将直线画在面板上即可
    
  
  

1.2改进的梯度上升算法

改进的第一点也就是改变每次向目标点靠近的步长的值，最初的时候也能稍微大点，然后随着迭代次数的增加，也就是离目标点越来越近，此时每次向前的步长也越来越小；第二点就是样本的选取也是随机的；第三点就是原来计算出来的h是一个100乘以1的矩阵，而现在算出来的是一个0-1之间的数值；第四点就是原来计算回归系数的时候使用一个3*100的矩阵乘以一个100*1的一个矩阵，现在是三个数值的乘积

1. #改进的梯度上升算法
2. #迭代次数150是根据上面的代码测试出来的
3. def stocGradAscent1(dataMatrix, classLabels, numIter=150):
4.     m,n = np.shape(dataMatrix)  #返回dataMatrix的大小。m为行数,n为列数。
5.     weights = np.ones(n)   #参数初始化
6.     for j in range(numIter):
7.         dataIndex = list(range(m))
8.         for i in range(m):
9.             alpha = 4/(1.0+j+i)+0.01  #降低alpha的大小，每次减小1/(j+i)。
10.             randIndex = int(random.uniform(0,len(dataIndex)))  #随机选取样本
11.             h = sigmoid(sum(dataMatrix[randIndex]*weights)) #选择随机选取的一个样本，计算h
12.             error = classLabels[randIndex] - h   #计算误差
13.             weights = weights + alpha * error * dataMatrix[randIndex]  #更新回归系数
14.             del(dataIndex[randIndex])  #删除已经使用的样本
15.     return weights

复制代码

1.3回归系数与迭代次数的关系

1. #2.为改进之前查看回归系数与迭代次数的关系
2. def gradAscent1(dataMatIn, classLabels):
3.     dataMatrix = np.mat(dataMatIn)   #转换成numpy的mat
4.     labelMat = np.mat(classLabels).transpose()  #转换成numpy的mat,并进行转置
5.     m, n = np.shape(dataMatrix)    #返回dataMatrix的大小。m为行数,n为列数。
6.     alpha = 0.01  #移动步长,也就是学习速率,控制更新的幅度。
7.     maxCycles = 500    #最大迭代次数
8.     weights = np.ones((n,1))
9.     weights_array = np.array([])
10.     for k in range(maxCycles):
11.         h = sigmoid(dataMatrix * weights)  #梯度上升矢量化公式
12.         error = labelMat - h
13.         weights = weights + alpha * dataMatrix.transpose() * error
14.         weights_array = np.append(weights_array,weights)
15.     weights_array = weights_array.reshape(maxCycles,n)
16.     return weights.getA(),weights_array #将矩阵转换为数组，并返回
18. #改进之后的
19. def stocGradAscent(dataMatrix, classLabels, numIter=150):
20.     m,n = np.shape(dataMatrix)                                                #返回dataMatrix的大小。m为行数,n为列数。
21.     weights = np.ones(n)                                                       #参数初始化
22.     weights_array = np.array([])                                            #存储每次更新的回归系数
23.     for j in range(numIter):
24.         dataIndex = list(range(m))
25.         for i in range(m):
26.             alpha = 4/(1.0+j+i)+0.01                                            #降低alpha的大小，每次减小1/(j+i)。
27.             randIndex = int(random.uniform(0,len(dataIndex)))                #随机选取样本
28.             h = sigmoid(sum(dataMatrix[randIndex]*weights))                    #选择随机选取的一个样本，计算h
29.             error = classLabels[randIndex] - h                                 #计算误差
30.             weights = weights + alpha * error * dataMatrix[randIndex]       #更新回归系数
31.             weights_array = np.append(weights_array,weights,axis=0)         #添加回归系数到数组中
32.             del(dataIndex[randIndex])                                         #删除已经使用的样本
33.     weights_array = weights_array.reshape(numIter*m,n)                         #改变维度
34.     return weights,weights_array

复制代码

测试结果

让我们分析一下。我们一共有100个样本点，改进的随机梯度上升算法迭代次数为150。而上图显示15000次迭代次数的原因是，使用一次样本就更新一下回归系数。因此，迭代150次，相当于更新回归系数150*100=15000次。简而言之，迭代150次，更新1.5万次回归参数。从上图左侧的改进随机梯度上升算法回归效果中可以看出，其实在更新2000次回归系数的时候，已经收敛了。相当于遍历整个数据集20次的时候，回归系数已收敛。训练已完成。
上图右侧的梯度上升算法回归效果，梯度上升算法每次更新回归系数都要遍历整个数据集。从图中可以看出，当迭代次数为300多次的时候，回归系数才收敛。凑个整，就当它在遍历整个数据集300次的时候已经收敛好了。
2.从疝气病症状预测病马的死亡率

代码实现

1. def classifyVector(inX, weights):
2.     prob = sigmoid(sum(inX*weights))
3.     if prob > 0.5: return 1.0
4.     else: return 0.0
6. def colicSklearn():
7.     frTrain = open('horseColicTraining.txt')  #打开训练集
8.     frTest = open('horseColicTest.txt') #打开测试集
9.     trainingSet = []; trainingLabels = []
10.     testSet = []; testLabels = []
11.     for line in frTrain.readlines():#取出训练集中的每一行的所有数据
12.         currLine = line.strip().split('\t')#每一行里数据的划分以空格划分
13.         lineArr = []
14.         for i in range(len(currLine)-1):#遍历每一行的每个元素
15.             lineArr.append(float(currLine[i]))#lineArr里存放的就是每一行中所有的数据
16.         trainingSet.append(lineArr)#将每一行的数据存放在训练集中
17.         trainingLabels.append(float(currLine[-1]))#拿到每一行的最后一列也即数据类别
18.     for line in frTest.readlines():#取出测试集中每一行的所有数据
19.         currLine = line.strip().split('\t')
20.         lineArr =[]
21.         for i in range(len(currLine)-1):
22.             lineArr.append(float(currLine[i]))
23.         testSet.append(lineArr)
24.         testLabels.append(float(currLine[-1]))
25.     classifier = LogisticRegression(solver='liblinear',max_iter=10).fit(trainingSet, trainingLabels)
26.     test_accurcy = classifier.score(testSet, testLabels) * 100
27.     print('正确率:%f%%' % test_accurcy)

复制代码

测试结果