机器学习系列(1)_逻辑回归初步

一、总述

逻辑回归是应用很是普遍的一个分类机器学习算法，它将数据拟合到一个logit函数(或者叫作logistic函数)中，从而可以完成对事件发生的几率进行预测。

二、由来

    要说逻辑回归，咱们得追溯到线性回归，想必你们对线性回归都有必定的了解，即对于多维空间中存在的样本点，咱们用特征的线性组合去拟合空间中点的分布和轨迹。以下图所示：

    线性回归能对连续值结果进行预测，而现实生活中常见的另一类问题是，分类问题。最简单的状况是是与否的二分类问题。好比说医生须要判断病人是否生病，银行要判断一我的的信用程度是否达到能够给他发信用卡的程度，邮件收件箱要自动对邮件分类为正常邮件和垃圾邮件等等。

    固然，咱们最直接的想法是，既然可以用线性回归预测出连续值结果，那根据结果设定一个阈值是否是就能够解决这个问题了呢？事实是，对于很标准的状况，确实能够的，这里咱们套用Andrew Ng老师的课件中的例子，下图中X为数据点肿瘤的大小，Y为观测结果是不是恶性肿瘤。经过构建线性回归模型，如hθ(x)所示，构建线性回归模型后，咱们设定一个阈值0.5，预测hθ(x)≥0.5的这些点为恶性肿瘤，而hθ(x)<0.5为良性肿瘤。

    但不少实际的状况下，咱们须要学习的分类数据并无这么精准，好比说上述例子中忽然有一个不按套路出牌的数据点出现，以下图所示：

    你看，如今你再设定0.5，这个断定阈值就失效了，而现实生活的分类问题的数据，会比例子中这个更为复杂，而这个时候咱们借助于线性回归+阈值的方式，已经很难完成一个鲁棒性很好的分类器了。

    在这样的场景下，逻辑回归就诞生了。它的核心思想是，若是线性回归的结果输出是一个连续值，而值的范围是没法限定的，那咱们有没有办法把这个结果值映射为能够帮助咱们判断的结果呢。而若是输出结果是 (0,1) 的一个几率值，这个问题就很清楚了。咱们在数学上找了一圈，还真就找着这样一个简单的函数了，就是很神奇的sigmoid函数(以下)：

 

    若是把sigmoid函数图像画出来，是以下的样子：

Sigmoid Logistic Function

    从函数图上能够看出，函数y=g(z)在z=0的时候取值为1/2，而随着z逐渐变小，函数值趋于0，z逐渐变大的同时函数值逐渐趋于1，而这正是一个几率的范围。

    因此咱们定义线性回归的预测函数为Y=WTX，那么逻辑回归的输出Y= g(WTX)，其中y=g(z)函数正是上述sigmoid函数(或者简单叫作S形函数)。

 

三、断定边界

    咱们如今再来看看，为何逻辑回归可以解决分类问题。这里引入一个概念，叫作断定边界，能够理解为是用以对不一样类别的数据分割的边界，边界的两旁应该是不一样类别的数据。

    从二维直角坐标系中，举几个例子，大概是以下这个样子：

    有时候是这个样子：

    

    甚至多是这个样子：

    上述三幅图中的红绿样本点为不一样类别的样本，而咱们划出的线，不论是直线、圆或者是曲线，都能比较好地将图中的两类样本分割开来。这就是咱们的断定边界，下面咱们来看看，逻辑回归是如何根据样本点得到这些断定边界的。

    咱们依旧借用Andrew Ng教授的课程中部分例子来说述这个问题。

    回到sigmoid函数，咱们发现：  

        当g(z)≥0.5时, z≥0;

        对于hθ(x)=g(θTX)≥0.5, 则θTX≥0, 此时意味着预估y=1;

        反之，当预测y = 0时，θTX<0;

    因此咱们认为θTX =0是一个决策边界，当它大于0或小于0时，逻辑回归模型分别预测不一样的分类结果。

    先看第一个例子hθ(x)=g(θ0+θ1X1+θ2X2)，其中θ0 ,θ1 ,θ2分别取-3, 1, 1。则当−3+X1+X2≥0时, y = 1; 则X1+X2=3是一个决策边界，图形表示以下，恰好把图上的两类点区分开来：

    例1只是一个线性的决策边界，当hθ(x)更复杂的时候，咱们能够获得非线性的决策边界，例如：

    这时当x12+x22≥1时，咱们断定y=1，这时的决策边界是一个圆形，以下图所示：

    因此咱们发现，理论上说，只要咱们的hθ(x)设计足够合理，准确的说是g(θTx)中θTx足够复杂，咱们能在不一样的情形下，拟合出不一样的断定边界，从而把不一样的样本点分隔开来。

四、代价函数与梯度降低

    咱们经过对断定边界的说明，知道会有合适的参数θ使得θTx=0成为很好的分类断定边界，那么问题就来了，咱们如何断定咱们的参数θ是否合适，有多合适呢？更进一步，咱们有没有办法去求得这样的合适参数θ呢？

    这就是咱们要提到的代价函数与梯度降低了。

 

    所谓的代价函数Cost Function，实际上是一种衡量咱们在这组参数下预估的结果和实际结果差距的函数，好比说线性回归的代价函数定义为:

固然咱们能够和线性回归类比获得一个代价函数，实际就是上述公式中hθ(x)取为逻辑回归中的g(θTx)，可是这会引起代价函数为“非凸”函数的问题，简单一点说就是这个函数有不少个局部最低点，以下图所示：

而咱们但愿咱们的代价函数是一个以下图所示，碗状结构的凸函数，这样咱们算法求解到局部最低点，就必定是全局最小值点。

      所以，上述的Cost Function对于逻辑回归是不可行的，咱们须要其余形式的Cost Function来保证逻辑回归的成本函数是凸函数。

      咱们跳过大量的数学推导，直接出结论了，咱们找到了一个适合逻辑回归的代价函数:

      Andrew Ng老师解释了一下这个代价函数的合理性，咱们首先看当y=1的状况：

       若是咱们的类别y = 1, 而断定的hθ(x)=1，则Cost = 0，此时预测的值和真实的值彻底相等，代价本该为0；而若是判断hθ(x)→0，代价->∞，这很好地惩罚了最后的结果。

       而对于y=0的状况，以下图所示，也一样合理：

       下面咱们说说梯度降低，梯度降低算法是调整参数θ使得代价函数J(θ)取得最小值的最基本方法之一。从直观上理解，就是咱们在碗状结构的凸函数上取一个初始值，而后挪动这个值一步步靠近最低点的过程，以下图所示：

    

 

      咱们先简化一下逻辑回归的代价函数：

       从数学上理解，咱们为了找到最小值点，就应该朝着降低速度最快的方向(导函数/偏导方向)迈进，每次迈进一小步，再看看此时的降低最快方向是哪，再朝着这个方向迈进，直至最低点。

       用迭代公式表示出来的最小化J(θ)的梯度降低算法以下：

五、代码与实现

    咱们来一块儿看两个具体数据上作逻辑回归分类的例子，其中一份数据为线性断定边界，另外一份为非线性。

    示例1。

    第一份数据为data1.txt，部份内容以下：

    咱们先来看看数据在空间的分布，代码以下。

 

[python] view plain copy

1. from numpy import loadtxt, where  
2. from pylab import scatter, show, legend, xlabel, ylabel  
3.   
4. #load the dataset  
5. data = loadtxt('/home/HanXiaoyang/data/data1.txt', delimiter=',')  
6.   
7. X = data[:, 0:2]  
8. y = data[:, 2]  
9.   
10. pos = where(y == 1)  
11. neg = where(y == 0)  
12. scatter(X[pos, 0], X[pos, 1], marker='o', c='b')  
13. scatter(X[neg, 0], X[neg, 1], marker='x', c='r')  
14. xlabel('Feature1/Exam 1 score')  
15. ylabel('Feature2/Exam 2 score')  
16. legend(['Fail', 'Pass'])  
17. show()  

 

 

    获得的结果以下：

    下面咱们写好计算sigmoid函数、代价函数、和梯度降低的程序：

 

[python] view plain copy

1. def sigmoid(X):  
2.     '''''Compute sigmoid function '''  
3.     den =1.0+ e **(-1.0* X)  
4.     gz =1.0/ den  
5.     return gz  
6. def compute_cost(theta,X,y):  
7.     '''''computes cost given predicted and actual values'''  
8.     m = X.shape[0]#number of training examples  
9.     theta = reshape(theta,(len(theta),1))  
10.       
11.     J =(1./m)*(-transpose(y).dot(log(sigmoid(X.dot(theta))))- transpose(1-y).dot(log(1-sigmoid(X.dot(theta)))))  
12.       
13.     grad = transpose((1./m)*transpose(sigmoid(X.dot(theta))- y).dot(X))  
14.     #optimize.fmin expects a single value, so cannot return grad  
15.     return J[0][0]#,grad  
16. def compute_grad(theta, X, y):  
17.     '''''compute gradient'''  
18.     theta.shape =(1,3)  
19.     grad = zeros(3)  
20.     h = sigmoid(X.dot(theta.T))  
21.     delta = h - y  
22.     l = grad.size  
23.     for i in range(l):  
24.         sumdelta = delta.T.dot(X[:, i])  
25.         grad[i]=(1.0/ m)* sumdelta *-1  
26.     theta.shape =(3,)  
27.     return  grad  

 

 

    咱们用梯度降低算法获得的结果断定边界是以下的样子：

    最后咱们使用咱们的断定边界对training data作一个预测，而后比对一下准确率：

 

[python] view plain copy

1. def predict(theta, X):  
2.     '''''Predict label using learned logistic regression parameters'''  
3.     m, n = X.shape  
4.     p = zeros(shape=(m,1))  
5.     h = sigmoid(X.dot(theta.T))  
6.     for it in range(0, h.shape[0]):  
7.         if h[it]>0.5:  
8.             p[it,0]=1  
9.         else:  
10.             p[it,0]=0  
11.     return p  
12. #Compute accuracy on our training set  
13. p = predict(array(theta), it)  
14. print'Train Accuracy: %f'%((y[where(p == y)].size / float(y.size))*100.0)  

 

    计算出来的结果是89.2%

 

    示例2.

    第二份数据为data2.txt，部份内容以下：

    咱们一样把数据的分布画出来，以下：

    咱们发如今这个例子中，咱们没有办法再用一条直线把两类样本点近似分开了，因此咱们打算试试多项式的断定边界，那么咱们先要对给定的两个feature作一个多项式特征的映射。好比说，咱们作了以下的一个映射：

    代码以下：

 

[python] view plain copy

1. def map_feature(x1, x2):  
2.     ''''' 
3.     Maps the two input features to polonomial features. 
4.     Returns a new feature array with more features of 
5.     X1, X2, X1 ** 2, X2 ** 2, X1*X2, X1*X2 ** 2, etc... 
6.     '''  
7.     x1.shape =(x1.size,1)  
8.     x2.shape =(x2.size,1)  
9.     degree =6  
10.     mapped_fea = ones(shape=(x1[:,0].size,1))  
11.     m, n = mapped_fea.shape  
12.     for i in range(1, degree +1):  
13.         for j in range(i +1):  
14.             r =(x1 **(i - j))*(x2 ** j)  
15.             mapped_fea = append(<span style="font-family: Arial, Helvetica, sans-serif;">mapped_fea</span><span style="font-family: Arial, Helvetica, sans-serif;">, r, axis=1)</span>  
16.     return mapped_fea  
17. mapped_fea = map_feature(X[:,0], X[:,1])  

 

      接着作梯度降低：

 

[python] view plain copy

1. def cost_function_reg(theta, X, y, l):  
2.     '''''Compute the cost and partial derivatives as grads 
3.     '''  
4.     h = sigmoid(X.dot(theta))  
5.     thetaR = theta[1:,0]  
6.     J =(1.0/ m)*((-y.T.dot(log(h)))-((1- y.T).dot(log(1.0- h)))) \  
7.             +(l /(2.0* m))*(thetaR.T.dot(thetaR))  
8.     delta = h - y  
9.     sum_delta = delta.T.dot(X[:,1])  
10.     grad1 =(1.0/ m)* sumdelta  
11.     XR = X[:,1:X.shape[1]]  
12.     sum_delta = delta.T.dot(XR)  
13.     grad =(1.0/ m)*(sum_delta + l * thetaR)  
14.     out = zeros(shape=(grad.shape[0], grad.shape[1]+1))  
15.     out[:,0]= grad1  
16.     out[:,1:]= grad  
17.     return J.flatten(), out.T.flatten()  
18. m, n = X.shape  
19. y.shape =(m,1)  
20. it = map_feature(X[:,0], X[:,1])  
21. #Initialize theta parameters  
22. initial_theta = zeros(shape=(it.shape[1],1))  
23. #Use regularization and set parameter lambda to 1  
24. l =1  
25. # Compute and display initial cost and gradient for regularized logistic  
26. # regression  
27. cost, grad = cost_function_reg(initial_theta, it, y, l)  
28. def decorated_cost(theta):  
29.     return cost_function_reg(theta, it, y, l)  
30. print fmin_bfgs(decorated_cost, initial_theta, maxfun=500)  

 

      接着在数据点上画出断定边界：

 

[python] view plain copy

1. #Plot Boundary  
2. u = linspace(-1,1.5,50)  
3. v = linspace(-1,1.5,50)  
4. z = zeros(shape=(len(u), len(v)))  
5. for i in range(len(u)):  
6.     for j in range(len(v)):  
7.         z[i, j]=(map_feature(array(u[i]), array(v[j])).dot(array(theta)))  
8. z = z.T  
9. contour(u, v, z)  
10. title('lambda = %f'% l)  
11. xlabel('Microchip Test 1')  
12. ylabel('Microchip Test 2')  
13. legend(['y = 1','y = 0','Decision boundary'])  
14. show()  
15. def predict(theta, X):  
16.     '''''Predict whether the label 
17.     is 0 or 1 using learned logistic 
18.     regression parameters '''  
19.     m, n = X.shape  
20.     p = zeros(shape=(m,1))  
21.     h = sigmoid(X.dot(theta.T))  
22.     for it in range(0, h.shape[0]):  
23.         if h[it]>0.5:  
24.             p[it,0]=1  
25.         else:  
26.             p[it,0]=0  
27.     return p  
28. #% Compute accuracy on our training set  
29. p = predict(array(theta), it)  
30. print'Train Accuracy: %f'%((y[where(p == y)].size / float(y.size))*100.0)  

 

       获得的结果以下图所示：

       咱们发现咱们获得的这条曲线确实将两类点区分开来了。

六、总结

       最后咱们总结一下逻辑回归。它始于输出结果为有实际意义的连续值的线性回归，可是线性回归对于分类的问题没有办法准确而又具有鲁棒性地分割，所以咱们设计出了逻辑回归这样一个算法，它的输出结果表征了某个样本属于某类别的几率。

       逻辑回归的成功之处在于，将本来输出结果范围能够很是大的θTX 经过sigmoid函数映射到(0,1)，从而完成几率的估测。

        而直观地在二维空间理解逻辑回归，是sigmoid函数的特性，使得断定的阈值可以映射为平面的一条断定边界，固然随着特征的复杂化，断定边界多是多种多样的样貌，可是它可以较好地把两类样本点分隔开，解决分类问题。

       求解逻辑回归参数的传统方法是梯度降低，构造为凸函数的代价函数后，每次沿着偏导方向(降低速度最快方向)迈进一小部分，直至N次迭代后到达最低点。

 

 

七、补充

 

      本文的2份数据可在http://pan.baidu.com/s/1pKxJl1p上下载到，分别为data1.txt和data2.txt，欢迎你们本身动手尝试。

      关于逻辑回归的完整ipython notebook示例代码能够在个人github上(https://github.com/HanXiaoyang/ML_examples/tree/master/logistic_regression)下载到，欢迎指正。