Logistic回归模型的Python及C++实现

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

优点：计算代价不高，易于理解和实现
缺点：容易欠拟合，分类精度可能不高
适用数据类型：数值型和标称型数据

1.1逻辑斯谛分布

分布函数为：

$$F(x)=P(X\leq x)= \frac{1}{1+e^{-(x-u)/\gamma}}$$
密度函数为：
$$f(x)=\frac{e^{-(x-u)/\gamma}}{\gamma(1+e^{-(x-u)/\gamma})^2}$$
其中， $u$为位置参数，$\gamma\gt{0}$为形状参数
逻辑斯谛分布的密度函数 $f(x)$和分布函数 $F(x)$的图形如图所示

逻辑斯谛函数的分布函数是一条S形曲线(sigmoid curve),该曲线以点 $(u,\frac{1}{2})$为中心对称，即满足 $$F(-x+u)-\frac{1}{2}=-F(x+u)+\frac{1}{2}$$
曲线在中心附近增长速度较快，在两端增长速度较慢，形状参数 $\gamma$的值越小，曲线在中心附近增长得越快。

1.2基于最优化方法的最佳回归系数确定

sigmoid函数的输入记为z,则sigmoid函数为：

$$\sigma(z)=\frac{1}{1+e^{-z}}$$
其中 $z=w_0x_0+w_1x_1+...+w_nx_n$ ，采用向量的形式可以写成 $z=w^Tx$，其中的 $x$是分类器的输入数据，向量$w$也就是我们要找的最佳系数

1.3二项逻辑斯谛回归模型

二项逻辑斯谛回归模型(binomial logistic regression model)是一种分类模型，由条件概率$P(Y|X)$表示，形式为逻辑斯谛分布。这里随机变量X取值为实数，随机变量Y为1或0。我们通过监督学习的方法来估计模型参数。
定义(逻辑斯谛回归模型)二项逻辑斯谛回归模型是如下的条件概率分布：

$$P(Y=1|x)=\frac{1}{1+e^{-z}}=\frac{e^z}{1+e^z}\\ =\frac{e^{w^Tx}}{1+e^{w^Tx}}$$
$$P(Y=0|x)=1-P(Y=1|x)\\ =\frac{e^{-z}}{1+e^{-z}}=\frac{1}{1+e^z}\\ =\frac{1}{1+e^{w^Tx}}$$
这里， $x\in \bf{R}^n$是输入， $Y \in${0,1}是输出， $w \in \bf{R}$是参数
一个事件的几率(odds)是指该事件发生的概率与该事件不发生的概率的比值。如果 $P(Y=1|x)=p$ ,则 $odds=\frac{p}{1-p}$,则对数几率为：
$$ln{\frac{P(Y=1|x)}{1-P(Y=1|x)}}=w^Tx$$
也就是说，在逻辑斯谛回归模型中，输出 $Y=1$的对数几率是输入 $x$线性函数，由$P(Y=1|x)=\frac{w^Tx}{1+w^Tx}$可以看出，线性函数的值越接近正无穷，概率值就越接近1；线性函数的值越接近负无穷，概率值就越接近0

1.4模型参数估计

逻辑斯谛回归模型学习时，对于给定的训练数据集

$$T=\{(x_1,y_1),(x_2,y_2),...,(x_M,y_M)\}$$ 其中， $x_i \in \bf{R}^n, \it y_i \in \{0,1\}$可以应用极大似然估计模型参数，从而得到逻辑斯谛回归模型，设 $P(Y=1|x)=\pi(x), P(Y=0|x)=1-\pi(x)$
似然函数为:
$$\prod_{i=1}^{M}[\pi(x_i)]^{y_i}[1-\pi(x_i)]^{1-y_i}$$ 对数似然函数为：
$$L(w)=\sum_{i=1}^{M}[y_i\ln\pi(x_i)+(1-y_i)\ln(1-\pi(x_i))]\\ =\sum_{i=1}^{M}[y_i\ln{\frac{\pi(x_i)}{1-\pi(x_i)}}+\ln(1-\pi(x_i))]\\ =\sum_{i=1}^{M}[y_i(w^Tx_i)-\ln(1+e^{w^Tx})]$$ 对 $L(w)$求极大值，得到 $w$的估计值，这样问题就变成了以对数似然函数为目标函数的最优化问题，逻辑斯谛回归中通常采用的方法是梯度下降法及拟牛顿法,其中
$$w^Tx=w_0+w_1x^1+w_2x^2+...,+w_Nx^N$$
$x_i^{j}表示第i个观测量x的第j个变量$

1.3梯度算法

利用梯度算法的迭代公式

$$w:=w+\alpha\nabla_wf(w)$$
利用偏微分公式给系数 $w$的每个分量迭代求值
$$\frac{\partial{\ln L(w)}}{\partial{w_k}}=\sum_{i=1}^{M}x_i^{k}(y_i-\pi(x_i))$$

1.4 训练算法：使用梯度算法求最佳参数

给出100个样本点，每个点包含两个数值型特征，每个回归系数初始化为1
重复R次：
计算整个数据集的梯度
使用$\alpha\nabla_wf(w)$更新回归系数的向量
返回回归系数

1.5 Python算法

创建logRegres.py文件

from numpy import *

def loadDataSet():
    dataMat=[];labelMat=[]
    fr=open('testSet.txt')
    for line in fr.readlines():
        lineArr=line.strip().split()
        dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])
        labelMat.append(int(lineArr[2]))
    return dataMat,labelMat


def sigmoid(inX):
    return 1.0/(1+exp(-inX))


def gradAscent(dataMatIn,classLabels):
    dataMatrix=mat(dataMatIn)               
    labelMat=mat(classLabels).transpose()   
    m,n=shape(dataMatrix)
    alpha=0.001
    maxCycles=500
    weights=ones((n,1))
    for k in range(maxCycles):
        h=sigmoid(dataMatrix*weights)
        error=(labelMat-h)
        weights=weights+alpha*dataMatrix.transpose()*error
    return weights

在python提示符下，敲入下面的代码：

>>>import logRegres
>>>dataAr,labelMat=logRegres.loadDataSet()
>>>logRegres.gradAscent(dataArr,labelMat)

结果为
matrix([[4.12414349],
[0.48007329],
[-0.6168482]])

def plotBestFit(wei):
    import matplotlib.pyplot as plt
    weights=wei.getA()
    dataMat,labelMat=loadDataSet()
    dataArr=array(dataMat)
    n=shape(dataArr)[0]
    xcord1=[];ycord1=[]
    xcord2=[];ycord2=[]
    for i in range(n):
        if int(labelMat[i]==1):
            xcord1.append(dataArr[i,1])
            ycord1.append(dataArr[i,2])
        else:
            xcord2.append(dataArr[i,1])
            ycord2.append(dataArr[i,2])
    fig=plt.figure()
    ax=fig.add_subplot(111)
    ax.scatter(xcord1,ycord1,s=30,c='red',marker='s')
    ax.scatter(xcord2,ycord2,s=30,c='green')
    x=arange(-3.0,3.0,0.1)
    y=(-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x,y)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

1.6 C++算法

#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include <cstring>
#include <stdio.h>
#include <algorithm>
#include <cmath>

using namespace std;

void loadDataset(vector<vector<double>> &dataMat,vector<int> &labelMat,const string &filename)
{

    ifstream file(filename);
    string line;
    while(getline(file,line))
    {
        vector<double> data;
        double x1,x2;
        int label;
        sscanf(line.c_str(),"%lf  %lf  %d",&x1,&x2,&label);
        data.push_back(1);
        data.push_back(x1);
        data.push_back(x2);
        dataMat.push_back(data);
        labelMat.push_back(label);
    }
}

double scalarProduct(vector<double> &w,vector<double> &x)
{
    double ret=0;
    for(int i=0;i<w.size();i++)
        ret+=w[i]*x[i];
    return ret;
}

double sigmoid(double z)
{
    return exp(z)/(1+exp(z));
}

vector<vector<double>> matTranspose(vector<vector<double>> &dataMat)
{
    vector<vector<double>> ret(dataMat[0].size(),vector<double>(dataMat.size(),0));
    for(int i=0;i<ret.size();i++)
        for(int j=0;j<ret[0].size();j++)
            ret[i][j]=dataMat[j][i];
    return ret;
}

void  gradAscent(vector<double> &weight,
                vector<vector<double>> &dataMat,vector<int> &labelMat)
{
    int maxCycles=500;
    double alpha=0.001;
    vector<vector<double>> dataMatT=matTranspose(dataMat);
    while(maxCycles>0)
    {
        vector<double> h;
        vector<double> error;
        for(auto &data:dataMat)
            h.push_back(sigmoid(scalarProduct(data,weight)));
        for(int i=0;i<labelMat.size();i++)
            error.push_back(labelMat[i]-h[i]);
        for(int i=0;i<weight.size();i++)
            weight[i]+=alpha*scalarProduct(dataMatT[i],error);
        maxCycles--;
    }

}

int main()
{
    vector<vector<double>> dataMat;
    vector<int> labelMat;
    string filename("testSet.txt");
    loadDataset(dataMat,labelMat,filename);
    vector<double> weight(dataMat[0].size(),1);
    gradAscent(weight,dataMat,labelMat);
    for(auto v:weight)
        cout<<v<<endl;  
}

1.7 随机梯度

梯度算法在每次更新回归系数时需要遍历整个数据集，如果样本数和特征数太多，则该方法的计算复杂度就太高了，一种改进方法是一次仅用一个样本点来更新回归系数，该方法称为随机梯度上升算法，由于可以在新样本到来时对分类器进行增量式更新，因而随机梯度上升算法是一个在线学习算法，与“在线学习相对应”，一次处理所有数据被称作是“批处理”
随机梯度上升算法：
所有回归系数初始化为1
对数据集中每个样本
计算该样本的梯度
使用$\alpha \times \nabla_wf(w)$更新回归系数
返回回归系数

1.8 随机梯度上升算法的Python实现

1.8.1 未优化的随机梯度算法

def stocGradAscent0(dataMatrix,classLabels):
    m,n=shape(dataMatrix)
    alpha=0.01
    weights=ones(n)
    for i in range(m):
        h=sigmoid(sum(dataMatrix[i]*weights))
        error=classLabels[i]-h
        weights=weights+alpha*error*dataMatrix[i]
    return weights

分类效果：

可以看到，分类的效果并不像前面那么完美，但前者是在整个数据集上迭代了500次才得到的，而此算法只是遍历了一遍训练数据。
一个判断优化算法优劣的可靠方法是看它是否收敛，我们对上述的梯度上升算法进行修改，使其在整个数据集上运行50次即可发现分类效果大幅度改进

1.8.2 一种改进措施

def stocGradAscent1(dataMatrix,classLabels,numIter=150):
    m,n=shape(dataMatrix)
    weights=ones(n)
    for j in range(numIter):
        dataIndex=range(m)
        for i in range(m):
            alpha=4/(1.0+j+i)+0.01
            randIndex=int(random.uniform(0,len(dataIndex)))
            h=sigmoid(sum(dataMatrix[randIndex]*weights))
            error=classLabels[randIndex]-h
            weights=weights+alpha*error*dataMatrix[randIndex]
            del(dataIndex[randIndex])
    return weights

第一处改进是alpha在每次迭代的时候都会调整，可以避免数据波动。
第二处改进是随机选取样本来更新回归系数，这种方法可以减少周期性的波动

1.9 随机梯度算法的C++实现

void stocGradAscent(vector<double> &weight,
                    vector<vector<double>> &dataMat,vector<int> &labelMat,int numIter=150)
{
    double alpha=0.01;
    double h=0.0;
    int i=0;
    int j=0;
    double error=0.0;
    vector<int> randIndex;
    for(i=0;i<dataMat.size();i++)
        randIndex.push_back(i);

    for(int k=0;k<numIter;k++)
    {
        random_shuffle(randIndex.begin(),randIndex.end());

        for(i=0;i<dataMat.size();i++)
        {
            alpha=4/(1+k+i)+0.01;
            h=sigmoid(scalarProduct(dataMat[randIndex[i]],weight));
            error=labelMat[randIndex[i]]-h;
            for(j=0;j<weight.size();j++)
            {
                weight[j]+=alpha*error*dataMat[randIndex[i]][j];
            }
        }
    }
}

2.0训练数据

-0.017612   14.053064   0
-1.395634   4.662541    1
-0.752157   6.538620    0
-1.322371   7.152853    0
0.423363    11.054677   0
0.406704    7.067335    1
0.667394    12.741452   0
-2.460150   6.866805    1
0.569411    9.548755    0
-0.026632   10.427743   0
0.850433    6.920334    1
1.347183    13.175500   0
1.176813    3.167020    1
-1.781871   9.097953    0
-0.566606   5.749003    1
0.931635    1.589505    1
-0.024205   6.151823    1
-0.036453   2.690988    1
-0.196949   0.444165    1
1.014459    5.754399    1
1.985298    3.230619    1
-1.693453   -0.557540   1
-0.576525   11.778922   0
-0.346811   -1.678730   1
-2.124484   2.672471    1
1.217916    9.597015    0
-0.733928   9.098687    0
-3.642001   -1.618087   1
0.315985    3.523953    1
1.416614    9.619232    0
-0.386323   3.989286    1
0.556921    8.294984    1
1.224863    11.587360   0
-1.347803   -2.406051   1
1.196604    4.951851    1
0.275221    9.543647    0
0.470575    9.332488    0
-1.889567   9.542662    0
-1.527893   12.150579   0
-1.185247   11.309318   0
-0.445678   3.297303    1
1.042222    6.105155    1
-0.618787   10.320986   0
1.152083    0.548467    1
0.828534    2.676045    1
-1.237728   10.549033   0
-0.683565   -2.166125   1
0.229456    5.921938    1
-0.959885   11.555336   0
0.492911    10.993324   0
0.184992    8.721488    0
-0.355715   10.325976   0
-0.397822   8.058397    0
0.824839    13.730343   0
1.507278    5.027866    1
0.099671    6.835839    1
-0.344008   10.717485   0
1.785928    7.718645    1
-0.918801   11.560217   0
-0.364009   4.747300    1
-0.841722   4.119083    1
0.490426    1.960539    1
-0.007194   9.075792    0
0.356107    12.447863   0
0.342578    12.281162   0
-0.810823   -1.466018   1
2.530777    6.476801    1
1.296683    11.607559   0
0.475487    12.040035   0
-0.783277   11.009725   0
0.074798    11.023650   0
-1.337472   0.468339    1
-0.102781   13.763651   0
-0.147324   2.874846    1
0.518389    9.887035    0
1.015399    7.571882    0
-1.658086   -0.027255   1
1.319944    2.171228    1
2.056216    5.019981    1
-0.851633   4.375691    1
-1.510047   6.061992    0
-1.076637   -3.181888   1
1.821096    10.283990   0
3.010150    8.401766    1
-1.099458   1.688274    1
-0.834872   -1.733869   1
-0.846637   3.849075    1
1.400102    12.628781   0
1.752842    5.468166    1
0.078557    0.059736    1
0.089392    -0.715300   1
1.825662    12.693808   0
0.197445    9.744638    0
0.126117    0.922311    1
-0.679797   1.220530    1
0.677983    2.556666    1
0.761349    10.693862   0
-2.168791   0.143632    1
1.388610    9.341997    0
0.317029    14.739025   0