Logistic回归

假设现在有一些数据点，我们用一条直线对这些点进行拟合（该直线称为最佳拟合直线），这个拟合的过程就称为回归。

利用Logistic（逻辑斯蒂）回归是一个分类模型而不回归模型。其进行分类的主要思想是：根据现有数据对分类边界线建立回归公式，以此进行分类。这里的“回归”一词源于最佳拟合，表示要找到最佳拟合参数。而最佳拟合参数就是在训练分类器时，通过最优化算法获得。

首先，逻辑斯蒂回归是一种线性分类器，针对的是线性可分问题！

先看两个函数logit和sigmoid函数.

几率比：指特定事件发生的几率。如下式所示。其中p为事件1发生的概率:$\frac{p}{1-p}$.

进一步地，我们可以logit函数，为几率比的对数函数$logit(p)=log\frac{p}{1-p}$,在logistic回归中,$q=p(y=1|x)$.

而对于某一样本属于特定类别的概率，为logit函数的反函数，称为logistic函数，由于它的图像呈S形，有时也称为sigmoid函数：$\pi (z)=\frac{1}{1+e^{-z}}$,在logistic回归中,$z=w^{T}x+b=w_1{x}_1+w_2{x}_2+\cdots +w_n{x}_n+b$,这里要知道${\pi }^{{}^{\prime }}(z)=\pi (z)(1-\pi (z))$.

1 Logistic回归模型

下面正式地来讲Logistic回归模型。

考虑具有个独立变量的向量$x=(x_1,x_2,\cdots ,x_n),$,设条件慨率$p(y=1|x)=p$为根据观测量相对于某事件x发生的概率.

那么Logistic回归模型可以表示为$p(y=1|x)=\pi (z)=\frac{1}{1+e^{-(w^{T}x+b)}}$,

其中,$z=w^{T}x+b=w_1{x}_1+w_2{x}_2+\cdots +w_n{x}_n+b$.

那么在x条件下y不发生的概率为

$p(y=0|x)=1-p(y=1|x)=1-\frac{1}{1+e^{-(w^{T}x+b)}}=\frac{e^{-(w^{T}x+b)}}{1+e^{-(w^{T}x+b)}}$.

所以事件发生与不发生的概率之比为$\frac{p(y=1|x)}{p(y=0|x)}=\frac{p}{1-p}=e^{w^{T}x+b}$,这个比值称为事件的发生比（the odds of experiencing an event），简记为odds。

对odds取对数得到$\ln (\frac{p}{1-p})=wx+b$.

先给出求解参数w的思路，下面我们便按照这个思路进行求解参数w：

1、为求解参数w，我们需要定义一个准则函数 J(w)，利用准则函数求解参数w.
2、我们通过最大似然估计法定义准则函数J(w).
3、接下来通过最大化准则函数J(w)便可求出参数w的迭代表达式.
4、为了更好地使用数据求出参数w，我们将第三步得到的w的迭代时向量化。

步骤一、求解准则函数J(w).

先将上面p(y|x)的两个式子合并: $p(y|x,w)=(\pi (z){)}^y(1-\pi (z){)}^{1-y}$.

上式是个二分类问题，所以y只是取两个值0或是1。

设有m组样本如下:

$(x_1^{(1)},x_2^{(1)},\cdots ,x_n^{(1)}),(x_1^{(2)},x_2^{(2)},\cdots ,x_n^{(2)}),\cdots ,(x_1^{(m)},x_2^{(m)},\cdots ,x_n^{(m)})$,

因为各个观测样本之间相互独立，那么它们的联合分布为各边缘分布的乘积,得到似然函数为

$L(w)=\prod _{i=1}^m(\pi (z^{(i)}){)}^{y^{(i)}}(1-\pi (z^{(i)}){)}^{1-y^{(i)}}$,其中$z^{(i)}=w_0+w_1{x}_1^{(i)}+w_2{x}_2^{(i)}+\cdots +w_n{x}_n^{(i)}$.

然后我们的目标是求出使这一似然函数的值最大的参数估计，最大似然估计就是求出参数w,b,使得L(w)取得最大值，对函数L(w)取对数得到

$J(w)=\frac{1}{m}\ln L(w)=\frac{1}{m}\sum _{i=1}^m(y^{(i)}\ln (\pi (z^{(i)}))+(1-y^{(i)})\ln (1-\pi (z^{(i)})))$.

现在开始L(w)分别对$w_1,w_2,\cdots ,w_n$求偏导,下面只计算对$w_j$的偏导.

$\frac{\partial }{\partial w_j}\ln L(w)=\sum _{i=1}^m(y^{(i)}\frac{1}{\pi (z^{(i)})}-(1-y^{(i)})\frac{1}{1-\pi (z^{(i)})})\frac{\partial }{\partial w_j}\pi (z^{(i)})$,

其中,$\frac{1}{\partial w_j}\pi (z^{(i)})=\frac{\partial \pi (z^{(i)})}{\partial z^{(i)}}\frac{\partial z^{(i)}}{\partial w_j}$.

由sigmoid函数的导数公式得${\pi }^{{}^{\prime }}(z^{(i)})=\pi (z^{(i)})(1-\pi (z^{(i)}))$,从而,$\frac{1}{\partial w_j}\pi (z^{(i)})=\pi (z^{(i)})(1-\pi (z^{(i)}))x_j$,带入上式得,

$\frac{\partial }{\partial w_j}\ln L(w)=\sum _{i=1}^m(y_i\frac{1}{\pi (z^{(i)})}-(1-y_i)\frac{1}{1-\pi (z^{(i)})})\pi (z^{(i)})(1-\pi (z^{(i)}))x_j^{(i)}=\sum _{i=1}^m(y_i-\pi (z^{(i)}))x_j^{(i)}$.

步骤二 梯度下降算法求解参数w.

这里我们使用梯度下降算法求解参数w，因此参数w的迭代式为：

$w:=w+\Delta w,\Delta W=-\alpha \nabla J(w)$,这里w是一个n行1列的向量,$w=(w_0,w_1,\cdots ,w_n{)}^T$,

$w_j$的更新公式如下:
$w_j:=w_j+\Delta w_j,\Delta w_j=-\alpha \frac{\partial J(w)}{\partial w_j}$,
其中,$\frac{\partial J(w)}{w_j}=\sum _{i=1}^m(y_i-\pi (z^{(i)}))x_j^{(i)}$.

因此，更新权重的公式如下:

$w_j:=w_j-\alpha \sum _{i=1}^m(y_j-\pi (z))x_j,j=0,1,2,\cdots ,n$.

矩阵版本的更新公式为:

$w:=w-\alpha \nabla J(w)=w-\alpha X^{T}E$,这里X和E的值如下所示:

$
X=\begin{bmatrix}
1 & x_1^{(1)} & x_2^{(1)} & \cdots & x_n^{(1)} \
1 & x_1^{(2)} & x_2^{(2)} & \cdots & x_n^{(2)} \
\vdots & \vdots & \vdots & \ddots & \vdots \
1 & x_1^{(m)} & x_2^{(m)} & \cdots & x_n^{(m)}
\end{bmatrix}
$,
$y=\left[\begin{array}{c}{y}^{(1)}\\ y^{(2)}\\ ⋮\\ y^{(m)}\end{array}\right]$
$,
w=\begin{bmatrix}
w_0 \ w_1 \ w_2 \ \vdots \ w_n
\end{bmatrix}
$,
$E=y-sigmoid(X\ast w)$.

步骤三 梯度上升算法的伪代码如下:

每个回归系数初始化为1
重复R次:
\qquad 计算整个数据集的梯度
\qquad 使用$\eta \times gradient$更新回归系数的向量
返回回归系数

代码如下:

def sigmoid(x):
    return 1.0/(1+np.exp(-x))
def gradAscent(data,classLabels):
    """
    Input:
        m,n=data.shape
        data的第一列是全1,即将w_0+w_1*x_1+w_2*x_2+...+w_{n-1}*x_{n-1}写作w_0*x_0+w_1*x_1+w_2*x_2+...+w_{n-1}*x_{n-1},其中x_0=1
        m表示样本数,比如100组样本
        n表示变量的个数,即x_0,x_1,x_{n-1},其中x_0恒等于1

        classLabels为m行1列
    Output:
        weights为n行1列,其中w_0为常数项,w_i为x_i的系数
    """
    data=np.matrix(data)
    labelMat=np.matrix(classLabels)
    m,n=data.shape
    alpha=0.001
    maxCycles=500
    weights=np.asmatrix(np.ones((n,1)))
    for k in range(maxCycles):
        h=sigmoid(data*weights)
        error=labelMat-h
        weights=weights+alpha*data.T*error
    return weights

加载数据,运行函数如下:

import pdb
#pdb.set_trace()
import pandas as pd
dataset=pd.read_csv('d:/testSet.txt',sep='\t',header=None)
labels_pd=dataset.loc[:,2]
labels=np.matrix(labels_pd).T
data_pd=dataset.loc[:,0:1]
data=pd.concat([pd.DataFrame(np.ones(100)),data_pd],axis=1)
w=gradAscent(data,labels)
#结果如下:
matrix([[ 4.12414349],
        [ 0.48007329],
        [-0.6168482 ]])

2 分析结果

def plotBestFit(data,labels,weights):
    import matplotlib.pyplot as plt
    import numpy as np
    data=np.matrix(data)
    labels=np.matrix(labels)
    m,n=data.shape
    x1=[];x2=[];y1=[];y2=[]
    for i in range(m):
        if int(labels[i,0])==1:
            x1.append(data[i,1]);y1.append(data[i,2])
        else:
            x2.append(data[i,1]);y2.append(data[i,2])
    fig=plt.figure()
    ax=fig.add_subplot(111)
    ax.scatter(x1,y1,s=30,c='red',marker='s')
    ax.scatter(x2,y2,s=30,c='green')
    x=np.arange(-3.0,3.0,0.1)
    y=(-weights[0,0]-weights[1,0]*x)/weights[2,0]
    ax.plot(x,y)
    plt.xlabel('x1');plt.ylabel('x2')
    plt.show()
plotBestFit(data,labels,w)

结果如下:

3 随机梯度上升

梯度上升算法在每次更新回归系数时都需要遍历整个数据集，该方法在处理小数据时还尚可，但如果有数十亿样本和成千上万的特征，那么该方法的计算复杂度太高了，改进方法便是一次仅用一个数据点来更新回归系数，此方法便称为随机梯度上升算法！由于可以在更新样本到来时对分类器进行增量式更新，因而随机梯度上升算法是一个“在线学习算法”。而梯度上升算法便是“批处理算法”！

3.1 代码和结果

def stocGradAscent(data,classLabels):
    data=np.array(data)
    labelMat=np.array(classLabels)
    m,n=data.shape
    alpha=0.01
    weights=np.ones(n,)

    for i in range(m):
        h=sigmoid(sum(data[i]*weights))
        error=labelMat[i,0]-h
        weights=weights+alpha*data[i]*error

    return weights

导入数据并画图如下:

import pandas as pd
dataset=pd.read_csv('d:/testSet.txt',sep='\t',header=None)
labels_pd=dataset.loc[:,2]
labels=np.matrix(labels_pd).T
data_pd=dataset.loc[:,0:1]
data=pd.concat([pd.DataFrame(np.ones(100)),data_pd],axis=1)
w=stocGradAscent(data,labels)

w=np.matrix(w).T
plotBestFit(data,labels,w)

可以看到,分类器错分了大概三分之一的样本.

3.2 结果分析

将stocGradAscent修改一下,将上面的代码执行num次,保存每次迭代时的weights的值,最后做出weights的三个量随迭代次数的变化曲线图.

def stocGradAscent(data,classLabels):
    data=np.array(data)
    labelMat=np.array(classLabels)
    m,n=data.shape
    alpha=0.5
    weights=np.ones(n,)
    num=200
    weights_list=[]
    for j in range(num):
        for i in range(m):
            h=sigmoid(sum(data[i]*weights))
            error=labelMat[i,0]-h
            weights=weights+alpha*data[i]*error
            weights_list.append(list(weights))
    return weights,weights_list

加载数据,运行程序:

import pandas as pd
dataset=pd.read_csv('d:/testSet.txt',sep='\t',header=None)
labels_pd=dataset.loc[:,2]
labels=np.matrix(labels_pd).T
data_pd=dataset.loc[:,0:1]
data=pd.concat([pd.DataFrame(np.ones(100)),data_pd],axis=1)
w,w_list=stocGradAscent(data,labels)

最后,画出weights的三个量随迭代次数的变化曲线图如下:

import matplotlib.pyplot as plt
w_list_mat=np.array(w_list)
fig=plt.figure()
ax1=fig.add_subplot(311)
ax1.plot(w_list_mat[:,0])
plt.ylabel('x0')
ax2=fig.add_subplot(312)
ax2.plot(w_list_mat[:,1])
plt.ylabel('x1')
ax3=fig.add_subplot(313)
ax3.plot(w_list_mat[:,2])
plt.ylabel('x2')
plt.show()

可以看出,回归系数经过大量迭代才能达到稳定值,并且仍然有局部的波动现象.

3.3 改进的随机梯度上升算法

随机梯度上升算法虽然大大减少了计算复杂度，但是同时正确率也下降了！所以可以对随机梯度上升算法进行改进！改进分为两个方面：

改进一:对于学习率alpha采用非线性下降的方式使得每次都不一样.

改进二:每次使用一个数据，但是每次随机的选取数据，选过的不在进行选择.

def stocGradAscent1(data,classLabels,numIter=200):
    data=np.array(data)
    labelMat=np.array(classLabels)
    m,n=data.shape
    alpha=0.5
    weights=np.ones(n,)
    weights_list=[]
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            alpha=4/(1.0+j+i)+0.01
            randIndex = int(np.random.uniform(0,len(dataIndex)))
            h=sigmoid(sum(data[randIndex]*weights))
            error=labelMat[randIndex,0]-h
            weights=weights+alpha*data[randIndex]*error
            del(dataIndex[randIndex])
            weights_list.append(list(weights))
    return weights,weights_list

import pandas as pd
dataset=pd.read_csv('d:/testSet.txt',sep='\t',header=None)
labels_pd=dataset.loc[:,2]
labels=np.matrix(labels_pd).T
data_pd=dataset.loc[:,0:1]
data=pd.concat([pd.DataFrame(np.ones(100)),data_pd],axis=1)
w,w_list=stocGradAscent1(data,labels)

w=np.matrix(w).T
plotBestFit(data,labels,w)

画出weights的三个值随迭代次数的变化曲线图如下:

import matplotlib.pyplot as plt
w_list_mat=np.array(w_list)
fig=plt.figure()
ax1=fig.add_subplot(311)
ax1.plot(np.linspace(0,50,150),w_list_mat[:,0])
plt.ylabel('x0')
ax2=fig.add_subplot(312)
ax2.plot(np.linspace(0,50,150),w_list_mat[:,1])
plt.ylabel('x1')
ax3=fig.add_subplot(313)
ax3.plot(np.linspace(0,50,150),w_list_mat[:,2])
plt.ylabel('x2')
plt.show()

可以看到,本图没有出现之前随机梯度算法参数weights的周期性波动,这归功于stocGradAscent1()里的样本随机选择机制,其次,可以看到,除了x1的微小波动外,x0,x2的很快就收敛,即改进的随机梯度算法收敛速度明显加快.