逻辑回归原理

1、逻辑回归1. 基本原理Logistic Regression和Linear Regression的原理是相似的，可以简单的描述为这样的过程：（1）找一个合适的预测函数，一般表示为h函数，该函数就是我们需要找的分类函数，它用来预测输入数据的判断结果。这个过程时非常关键的，需要对数据有一定的了解或分析，知道或者猜测预测函数的“大概”形式，比如是线性函数还是非线性函数。（2）构造一个Cost函数（损失函数），该函数表示预测的输出（h）与训练数据类别（y）之间的偏差，可以是二者之间的差（h-y）或者是其他的形式。综合考虑所有训练数据的“损失”，将Cost求和或者求平均，记为J()函数，表示所有训练数据

2、预测值与实际类别的偏差。（3）显然，J()函数的值越小表示预测函数越准确（即h函数越准确），所以这一步需要做的是找到J()函数的最小值。找函数的最小值有不同的方法，Logistic Regression实现时有的是梯度下降法（Gradient Descent）。2. 具体过程1.1 构造预测函数Logistic Regression虽然名字里带“回归”，但是它实际上是一种分类方法，用于两分类问题（即输出只有两种）。根据步骤，需要先找到一个预测函数（h），显然，该函数的输出必须是两个值（分别代表两个类别），所以利用了Logistic函数（或称为Sigmoid函数），函数形式为： 对应的函数图像是

3、一个取值在0和1之间的S型曲线（图1）。图1接下来需要确定数据划分的边界类型，对于图2和图3中的两种数据分布，显然图2需要一个线性的边界，而图3需要一个非线性的边界。接下来我们只讨论线性边界的情况。图2图3对于线性边界的情况，边界形式如下：构造预测函数为：h(x)函数的值有特殊的含义，它表示结果取1的概率，因此对于输入x分类结果为类别1和类别0的概率分别为：3.2 构造Cost函数Andrew Ng在课程中直接给出了Cost函数及J()函数如式（5）和（6），但是并没有给出具体的解释，只是说明了这个函数来衡量h函数预测的好坏是合理的。实际上这里的Cost函数和J()函数是基于最大似然估计推导得

4、到的。下面详细说明推导的过程。（4）式综合起来可以写成：取似然函数为：对数似然函数为：最大似然估计就是要求得使l()取最大值时的，其实这里可以使用梯度上升法求解，求得的就是要求的最佳参数。但是，在Andrew Ng的课程中将J()取为（6）式，即：因为乘了一个负的系数-1/m，所以J()取最小值时的为要求的最佳参数。3.3 梯度下降法求J()的最小值求J()的最小值可以使用梯度下降法，根据梯度下降法可得的更新过程：式中为学习步长，下面来求偏导：上式求解过程中用到如下的公式：因此，（11）式的更新过程可以写成：因为式中本来为一常量，所以1/m一般将省略，所以最终的更新过程为：另外，补充一下，3.

5、2节中提到求得l()取最大值时的也是一样的，用梯度上升法求（9）式的最大值，可得：观察上式发现跟（14）是一样的，所以，采用梯度上升发和梯度下降法是完全一样的，这也是机器学习实战中采用梯度上升法的原因。3.4 梯度下降过程向量化关于更新过程的vectorization，Andrew Ng的课程中只是一带而过，没有具体的讲解。机器学习实战连Cost函数及求梯度等都没有说明，所以更不可能说明vectorization了。但是，其中给出的实现代码确是实现了vectorization的，图4所示代码的32行中weights（也就是）的更新只用了一行代码，直接通过矩阵或者向量计算更新，没有用for循环，

6、说明确实实现了vectorization，具体代码下一章分析。文献3中也提到了vectorization，但是也是比较粗略，很简单的给出vectorization的结果为：且不论该更新公式正确与否，这里的(.)是一个求和的过程，显然需要一个for语句循环m次，所以根本没有完全的实现vectorization，不像机器学习实战的代码中一条语句就可以完成的更新。下面说明一下我理解机器学习实战中代码实现的vectorization过程。约定训练数据的矩阵形式如下，x的每一行为一条训练样本，而每一列为不同的特称取值：约定待求的参数的矩阵形式为： 先求x.并记为A：求h(x)-y并记为E：g(A)的参数

7、A为一列向量，所以实现g函数时要支持列向量作为参数，并返回列向量。由上式可知h(x)-y可以由g(A)-y一次计算求得。再来看一下（15）式的更新过程，当j=0时：同样的可以写出j，综合起来就是：综上所述，vectorization后更新的步骤如下：（1）求A=x.；（2）求E=g(A)-y；（3）求:=-.x.E,x表示矩阵x的转置。也可以综合起来写成：前面已经提到过：1/m是可以省略的。4. 代码分析图4中是机器学习实战中给出的部分实现代码。图4sigmoid函数就是前文中的g(z)函数，参数inX可以是向量，因为程序中使用了Python的numpy。gradAscent函数是梯度上升的实

8、现函数，参数dataMatin和classLabels为训练数据，23和24行对训练数据做了处理，转换成numpy的矩阵类型，同时将横向量的classlabels转换成列向量labelMat，此时的dataMatrix和labelMat就是（18）式中的x和y。alpha为学习步长，maxCycles为迭代次数。weights为n维（等于x的列数）列向量，就是（19）式中的。29行的for循环将更新的过程迭代maxCycles次，每循环一次更新一次。对比3.4节最后总结的向量化的更新步骤，30行相当于求了A=x.和g(A)，31行相当于求了E=g(A)-y，32行相当于求:=-.x.E。所以这

9、三行代码实际上与向量化的更新步骤是完全一致的。4.1代码示例使用鸢尾花数据集进行代码运行from sklearn.datasets import load_irisimport matplotlib.pyplot as pltimport numpy as npiris = load_iris()data = iris.datatarget = iris.target#print data:10#print target10:X = data0:100,0,2.y = target0:100print X:5print y-5:label = np.array(y)index_0 = np.w

10、here(label=0)plt.scatter(Xindex_0,0,Xindex_0,1,marker=x,color = b,label = 0,s = 15)index_1 =np.where(label=1)plt.scatter(Xindex_1,0,Xindex_1,1,marker=o,color = r,label = 1,s = 15)plt.xlabel(X1)plt.ylabel(X2)plt.legend(loc = upper left)plt.show()逻辑回归模型的类：import numpy as npclass logistic(object):def _

11、init_(self):self.W = Nonedef train(self,X,y,learn_rate = 0.01,num_iters = 5000):num_train,num_feature = X.shape#init the weightself.W = 0.001*np.random.randn(num_feature,1).reshape(-1,1)loss = for i in range(num_iters):error,dW = pute_loss(X,y)self.W += -learn_rate*dWloss.append(error)if i%2

12、00=0:print i=%d,error=%f %(i,error)return lossdef compute_loss(self,X,y):num_train = X.shape0h = self.output(X)loss = -np.sum(y*np.log(h) + (1-y)*np.log(1-h)loss = loss / num_traindW = X.T.dot(h-y) / num_trainreturn loss,dWdef output(self,X):g = np.dot(X,self.W)return self.sigmod(g)def sigmod(self,X

13、):return 1/(1+np.exp(-X)def predict(self,X_test):h = self.output(X_test)y_pred = np.where(h=0.5,1,0)return y_pred训练测试，并且可视化跟踪的损失lossimport matplotlib.pyplot as plty = y.reshape(-1,1)#add the x0=1one = np.ones(X.shape0,1)X_train = np.hstack(one,X)classify = logistic()loss = classify.train(X_train,y)p

14、rint classify.Wplt.plot(loss)plt.xlabel(Iteration number)plt.ylabel(Loss value)plt.show()可视化决策边界label = np.array(y)index_0 = np.where(label=0)plt.scatter(Xindex_0,0,Xindex_0,1,marker=x,color = b,label = 0,s = 15)index_1 =np.where(label=1)plt.scatter(Xindex_1,0,Xindex_1,1,marker=o,color = r,label = 1

15、,s = 15)#show the decision boundaryx1 = np.arange(4,7.5,0.5)x2 = (- classify.W0 - classify.W1*x1) / classify.W2plt.plot(x1,x2,color = black)plt.xlabel(X1)plt.ylabel(X2)plt.legend(loc = upper left)plt.show()4.1代码示例第二部分1. 数据准备下面的数据来自机器学习实战中的示例：-0.14.0-1.4.1-0.6.0-1.7.00.11.02.0-.0-0.10

16、.0上面的数据一共是3列10行，其中前两列为x1和x2的值，第3列表示y的值；10行表示取了10个样本点。我们可以将这些数据当做训练模型参数的训练样本。见到训练样本就可以比较直观的理解算法的输入，以及我们如何利用这些数据来训练逻辑回归分类器，进而用训练好的模型来预测新的样本（检测样本）。从逻辑回归的参数形式，式子（1）我们可以看到逻辑回归模型中有两个待定参数a（x的系数）和b（常数项），我们现在给出来的数据有两个特征x1, x2，因此整个模型就增加了一项：ax1+ cx2+ b。为了形式上的统一，我们使用带下标的a表示不同的参数（a0表示常数项b并作x0的参数，a1、a2分别表示x1和x2的参

17、数），就可以得到：$ a_0x_0+ a_1x_1+ a_2x_2 $这样统一起来后，就可以使用矩阵表示了（比起前面展开的线性表示方式，用矩阵表示模型和参数更加简便，而且矩阵运算的速度也更快）：$ beginbmatrix a_0 & a_1 & a_2 endbmatrix beginbmatrix x_0 x_1 x_2 endbmatrix = a mathrm T X$将上面的式子带入到（1）式，我们就可以得到逻辑回归的另一种表示形式了：$p(x; a) = frac11 + e-a mathrm T X （2）$此时，可以很清楚的看到，我们后面的行动都是为了确定一个合适的a（一个参数

18、向量），使得对于一个新来的X（也是一个向量），我们可以尽可能准确的给出一个y值，0或者1.注：数据是二维的，也就是说这组观察样本中有两个自变量，即两个特征（feature）。2. 训练分类器就像上面说的，训练分类器的过程，就是根据已经知道的数据（训练样本）确定一个使得代价函数的值最小的a（参数向量/回归系数）的过程。逻辑回归模型属于有监督的学习方法，上面示例数据中的第3列其实是训练样本提供的标准答案。也就是说，这些数据是已经分好类的（两类，0或者1）。在训练阶段，我们要做的就是利用训练样本和（2）式中的模型，估计一个比较合适的参数a，使得仅通过前面两列数据（观察值/测量值）就可以估计一个值h(

19、a)，这个值越接近标准答案y，说明我们的模型预测的越准确。下面是估计回归系数a的值的过程，还是借鉴了机器学习实战中的代码，做了少量修改：其中计算参数梯度，即代价函数对每个参数的偏导数（下面代码中的第36-38行），的详细推导过程可以参考这里1 2 Created on Oct 27, 20103 Logistic Regression Working Module 4 author: Peter 5 6 from numpy import * 7 import os 8 9 path = D:MechineLearningMLiA_SourceCodemachinelearninginacti

20、onCh0510 training_sample = trainingSample.txt11 testing_sample = testingSample.txt12 13 # 从文件中读入训练样本的数据，同上面给出的示例数据14 # 下面第20行代码中的1.0表示x0 = 115 def loadDataSet(p, file_n):16 dataMat = ; labelMat = 17 fr = open(os.path.join(p, file_n)18 for line in fr.readlines():19 lineArr = line.strip().split()20 da

21、taMat.append(1.0, float(lineArr0), float(lineArr1) # 三个特征x0, x1, x221 labelMat.append(int(lineArr2) # 标准答案y22 return dataMat,labelMat23 24 def sigmoid(inX):25 return 1.0/(1+exp(-inX)26 27 # 梯度下降法求回归系数a，由于样本量少，我将迭代次数改成了1000次28 def gradAscent(dataMatIn, classLabels):29 dataMatrix = mat(dataMatIn) #con

22、vert to NumPy matrix30 labelMat = mat(classLabels).transpose() #convert to NumPy matrix31 m,n = shape(dataMatrix)32 alpha = 0.001 # 学习率33 maxCycles = 100034 weights = ones(n,1)35 for k in range(maxCycles): # heavy on matrix operations36 h = sigmoid(dataMatrix*weights) # 模型预测值, 90 x 137 error = h - l

23、abelMat # 真实值与预测值之间的误差, 90 x 138 temp = dataMatrix.transpose()* error # 交叉熵代价函数对所有参数的偏导数, 3 x 139 weights = weights - alpha * temp # 更新权重40 return weights41 42 # 下面是我自己写的测试函数43 def test_logistic_regression():44 dataArr, labelMat = loadDataSet(path, training_sample) # 读入训练样本中的原始数据45 A = gradAscent(da

24、taArr, labelMat) # 回归系数a的值46 h = sigmoid(mat(dataArr)*A) #预测结果h(a)的值47 print(dataArr, labelMat)48 print(A)49 print(h)50 # plotBestFit(A.getA()51 52 test_logistic_regression()上面代码的输出如下： 一个元组，包含两个数组：第一个数组是所有的训练样本中的观察值，也就是X，包括x0, x1, x2；第二个数组是每组观察值对应的标准答案y。(1.0, -0., 14., 1.0, -1., 4., 1.0, -0., 6.5386

25、2, 1.0, -1., 7., 1.0, 0., 11., 1.0, 0., 7., 1.0, 0., 12., 1.0, -2.46015, 6., 1.0, 0., 9., 1.0, -0., 10., 0, 1, 0, 0, 0, 1, 0, 1, 0, 0) 本次预测出来的回归系数a，包括a0, a1, a2 1.-0. -0. 根据回归系数a和（2）式中的模型预测出来的h(a)。这里预测得到的结果都是区间(0, 1)上的实数。 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 标准答案是0, 1，如何将预测到的结果与标准答案y进行比较呢？取0.5作为阈值，大于该值的样本就

26、划分到1这一组，小于等于该值的样本就划分到0这一组，这样就可以将数据分为两类。检查一下结果可以看到，我们现在分出来的1这一类中包括原来y=1的两个样本，另一类包括原来y=0的所有样本和一个y=1的样本（分错了）。鉴于我们选择取的样本比较少（只有10个），这样的效果其实还算非常不错的！3. 结果展示上面已经求出了一组回归系数，它确定了不同类别数据之间的分割线。可以利用X内部（x1与x2之间的关系）的关系画出该分割线，从而更直观的感受到分类的效果。添加下面一段代码：1 # 分类效果展示，参数weights就是回归系数 2 def plotBestFit(weights): 3 import mat

27、plotlib.pyplot as plt 4 dataMat,labelMat=loadDataSet()5 dataArr = array(dataMat)6 n = shape(dataArr)0 7 xcord1 = ; ycord1 = 8 xcord2 = ; ycord2 = 9 for i in range(n):10 if int(labelMati)= 1:11 xcord1.append(dataArri,1); ycord1.append(dataArri,2)12 else:13 xcord2.append(dataArri,1); ycord2.append(dat

28、aArri,2)14 fig = plt.figure()15 ax = fig.add_subplot(111)16 ax.scatter(xcord1, ycord1, s=30, c=red, marker=s)17 ax.scatter(xcord2, ycord2, s=30, c=green)18 x = arange(-3.0, 3.0, 0.1)19 y = (-weights0-weights1*x)/weights2 # x2 = f(x1)20 ax.plot(x, y)21 plt.xlabel(X1); plt.ylabel(X2);22 plt.show()将上面的test_logistic_regression()函数中的最后一句注释去掉，调用p