05 Logistic Regression

1 Logistic 原理

逻辑回归从字面上看是一个回归问题，很容易和回归分析联系到一起。这是字面意思给人的错觉，其实它不是回归问题！它是最典型的分类问题，经典的二分类，比如检查邮件是否是垃圾邮件，检查一个零件是否是合格件，判断一个西瓜是否成熟等等，这些日常生活中经常用于判断是非黑白的问题，都能用到这个经典的二分类之逻辑回归。如下图：

这样一个线性可分的问题，需要这么一个超平面将其分开。而逻辑回归，使用阈值函数，将样本映射到不同类别中。常用的映射函数为Sigmoid函数形式如下：
$$f(x)=\frac{1}{1+e^{-x}}$$它的图像如下：

从图像可以看出，其函数值值域在（0,1），在 0 附近的变化较为明显。其导数为：
$$f^{\prime }(x)=\frac{e^{-x}}{(1+e^{-x}{)}^2}=f(x)[1-f(x)]$$
对于线性回归模型：$$h_{\theta }(x)={\theta }^{T}x$$
将线性模型带入到g(x)中，得到最终的逻辑回归模型：
$$h_{\theta }(x)=g({\theta }^{T}x)=\frac{1}{1+e^{-{\theta }^{T}x}}$$
对于输入向量X，其属于正例的概率：
$$P(y=1|X;\theta )=h_{\theta }(x)$$
对于输入向量X,其属于负例的概率为：
$$P(y=0|X;\theta )=1-h_{\theta }(x)$$

整合起来，就是其属于类别c的概率为：
$$P(y=c|X;\theta )=(h_{\theta }(x){)}^y(1-h_{\theta }(x){)}^y$$

要求参数，方法与线性回归一样，先求出由所有样本组成的似然函数，转化为对数似然，转化为梯度下降，根据梯度下降的方向更新迭代权重参数，选取代价函数本次与上次的差小于某个阈值的迭代思路求解。
其中，似然函数为，m表示样本个数
$$L(\theta )=\prod _{i=1}^m(h_{\theta }(x){)}^y(1-h_{\theta }(x){)}^{1-y}$$
对数似然函数：
$$log(L(\theta ))=\sum _{i=1}^{m}y(h_{\theta }(x))+(1-y)(1-h_{\theta }(x))$$
然后，利用对数极大似然估计，即求上式的极大值，引入因子 -1/m，转化为求下式的极小值：
$$J(\theta )=-\frac{1}{m}log(L(\theta ))$$
对上式求偏导后的结果如下，在此不一 一详细写出推导过程了，但是说明下式的参数含义，其中${\theta }_j$表示第i个样本的第j个特征的权重参数，注意如果有100个特征，自然地就有100个权重参数，等号左边是对${\theta }_j$求偏导的，因此这个式子相当于对某个样本的所有特征参数都适应于下面的公式。
$$\frac{\partial }{\partial {\theta }_j}J(\theta )=\frac{1}{m}\sum _{i=1}^m(h_{\theta }(x)-y_i){\theta }_i^j$$

2 代码

先生成第一张图的数据，也就是实验数据。然后使用梯度下降法求出俩个特征对应的权重参数。

设定一个学习率迭代参数，当与前一时步的代价函数与当前的代价函数的差小于阈值时，计算结束，我们将得到3个权重参数，其中包括两个特征的权重参数，和偏置项的权重参数。
梯度下降求逻辑回归模型的权重参数的基本思路如下：

1. ‘model’ 建立的逻辑回归模型：包括Sigmoid映射, 也就是上面的 $h_{\theta }(x)$
2. ‘cost’ 代价函数,也就是上诉的$J(\theta )$
3. ‘gradient’ 梯度公式,也就是上诉的$\frac{\partial }{\partial {\theta }_j}J(\theta )$
4. ‘theta update’ 参数更新公式
5. ‘stop stratege’ 迭代停止策略：代价函数小于阈值法

#!/user/bin/env python
#-*- coding:utf-8 -*-
#by: LiuWei
"""逻辑斯提克回归"""
import numpy as np
import matplotlib.pyplot as plt
import time

#按照一定规律均匀分布含有两个特征的数据点
def createData(samplecnt,coef=1.0,intercept=0.05):
    x1 = np.random.uniform(0,1,samplecnt)
    x2 = np.random.uniform(0,1,samplecnt)
    index = (x2-intercept)/x1 <coef
    x1_pos = x1[index]
    x2_pos = x2[index]
    index = (x2-intercept)/x1 >=coef
    x1_neg = x1[index]
    x2_neg = x2[index]
    plt.xlabel("w1")
    plt.ylabel("w2")
    plt.scatter(x1_pos,x2_pos)
    plt.scatter(x1_neg,x2_neg)
    regx = np.linspace(0,1,samplecnt)
    regy = coef*regx+intercept
    #plt.plot(regx,regy,color='g')
    plt.show()
    return x1_pos,x1_neg,x2_pos,x2_neg

#组合成原始数据
def combine_data(x1_pos,x1_neg,x2_pos,x2_neg):
    x1_pos_1 = x1_pos.reshape(-1,1)
    x2_pos_1 = x2_pos.reshape(-1,1)
    x_pos = np.concatenate((x1_pos_1,x2_pos_1),axis=1)
    x_pos_shape =  np.shape(x_pos)
    y_pos = np.ones(x_pos_shape[0])
    y_pos = y_pos.reshape(-1,1)
    data_pos = np.concatenate((x_pos,y_pos),axis=1)
    x1_neg_1 = x1_neg.reshape(-1,1)
    x2_neg_1 = x2_neg.reshape(-1,1)
    x_neg = np.concatenate((x1_neg_1,x2_neg_1),axis=1)
    x_neg_shape = np.shape(x_neg)
    y_neg = np.zeros(x_neg_shape[0])
    y_neg = y_neg.reshape(-1,1)
    data_neg = np.concatenate((x_neg,y_neg),axis=1)
    data = np.vstack((data_pos,data_neg))
    data = np.random.permutation(data)
    return data

#生成数据
x1_pos,x1_neg,x2_pos,x2_neg=createData(200)
data=combine_data(x1_pos,x1_neg,x2_pos,x2_neg)
# print(data[:10,:])
X = data[:,:2]
y = data[:,-1]

#偏移量 b shape=(200,1)
b = np.ones(200)
#将偏移量与2个特征值组合 shape = (200,3)
X = np.column_stack((b,X))

# #model 
def sigmoid(x):
    return 1/(1+ np.exp(-x))
def model(theta,X):
    theta = np.array(theta)
    return sigmoid( X.dot(theta) )

#cost
def cost(m,theta,X,y):
    ele = y*np.log(model(theta,X)) + (1-y)*np.log(1-model(theta,X))
    item_sum = np.sum(ele)
    return -item_sum/m

#gradient
def gradient(m,theta,X,y,cols):
    grad_theta = []
    for j in range(cols):
        grad = (model(theta,X) - y).dot(X[:,j])
        grad_sum = np.sum(grad)
        grad_theta.append(grad_sum/m)
    return np.array(grad_theta)

#theta update
def theta_update(grad_theta,theta,sigma):
    return theta - sigma * grad_theta

#stop stratege
def stop_stratege(cost,cost_update,threshold):
    return cost-cost_update < threshold

#逻辑回归算法
def LogicRegression(X,y,threshold,m,xcols):
    start = time.clock()
    #设置权重参数的初始值
    theta = np.zeros(xcols)
    #迭代步数
    iters = 0
    #记录代价函数的值
    cost_record=[]
    #学习率
    sigma = 0.01
    cost_val = cost(m,theta,X,y)
    cost_record.append(cost_val)
    while True:
        grad = gradient(m,theta,X,y,xcols)
        #参数更新
        theta = theta_update(grad,theta,sigma)
        cost_update = cost(m,theta,X,y)
        if stop_stratege(cost_val,cost_update,threshold):
            break
        iters=iters+1
        cost_val = cost_update
##        print("cost_val:%f" %cost_val)
        cost_record.append(cost_val)
    end = time.clock()
    print("LogicRegressionconvergence duration: %f s" % (end - start))
    return cost_record, iters,theta

cost_record, iters,theta= LogicRegression(X,y,1e-7,200,3)#调用逻辑回归函数
plt.scatter(x1_pos,x2_pos)
plt.scatter(x1_neg,x2_neg)
wp = np.linspace(0.0,1.0,200)
plt.plot(wp,-(theta[0]+theta[1]*wp)/theta[2],color='g')
plt.show()

调用函数，收敛时得到的权重为：

array([  1.4441927 ,  17.07777938, -17.26372013])

参数的含义：第一个权重参数为偏置项，第二、三个权重参数相当，只不过贡献方向相反而已。

画出决策面如下：

效果还不错！

总结

1. 如果代价函数的最后稳定的值，确认比较大，比如0.5，说明模型中一定存在某些bug，比如在调试过程中，将标签值错误地被赋值了第三列，实际应该为第四列，所以导致最后迭代终止时的成本值为0.50。
2. 学习率直接关系到迭代速度，如果学习率太小，迭代下降的会很慢，相反会比较快。
3. 迭代终止的策略选取，一般会根据迭代次数，成本函数前后两次的差小于某个阈值等，如果选取的终止策略不当，会导致看似收敛，实际成本值还很大的情况。