Logistic Regression 原理及推导 python实现

一、问题引入

首先，Logistic回归是一种广义的线性回归模型，主要用于解决二分类问题。比如，现在我们有N个样本点，每个样本点有两维特征x1和x2，在直角坐标系中画出这N个样本的散点图如下图所示，

蓝色和红色分别代表两类样本。现在我们的目标是，根据这N个样本所表现出的特征以及他们各自对应的标签，拟合出一条直线对两类样本进行分类，直线的上侧属于第一类，直线的下侧属于第二类。那么我们如何寻找这条直线呢？我们知道，一条直线可以用简单的公式表示

$$y=\theta _{0}+\theta_{1}x_{1}+\theta_{2}x_{2}+... = \theta ^{T}x$$ 参数$\theta^{T}$的选择决定了直线的位置，如果我们选取了一组参数$\theta$
导致直线的位置是这样的

那肯定不合理，因为两类样本完全没有被分开，而如果我们得到了这样一条直线

两类样本虽然存在一些小错误，但是基本上被分开了。由此，我们可以看到，Logistic Regression问题最终变成了求解最优参数$\theta$的问题。

二、原理

样本的每一维特征的取值在经过参数$\theta$线性组合之后取值范围是实数集（-inf, inf），而要想对实数进行二分类就要通过一个函数将实数投影到某个有限区间上，在有限区间内找到一个阈值，大于这个阈值分为第一类，小于等于这个阈值分为第二类。LR找到的这个投影函数就是sigmoid函数

$$g(z)=\frac{1}{1+e^{-z}}$$

值域为（0，1），当x=0时，函数值为0.5。在实际分类中，由于假设样本均匀分布，所以阈值通常选取为0.5。

现在我们有N个样本$x_{1},x_{2},x_{3}...$，每个样本有一个类别标签y(y=0 or y=1)与之对应，我们知道它们的对应关系可以由参数$\theta$来估计，因此用极大似然估计法来求解我们我们的目标$\theta$。

首先，把问题变成一个概率问题：
在某个$x$和$\theta$的取值下，$y=1$的概率为$h_{\theta }(x)$，

$$P(y=1|x;\theta )=h_{\theta }(x)$$
在某个$x$和$\theta$的取值下，$y=0$的概率为$1-h_{\theta }(x)$， $$P(y=0|x;\theta )=1-h_{\theta }(x)$$
由于$y$只有两种取值：0和1，因此综合两种情况，对于每一个样本点来说， $$P(y|x;\theta )=(h_{\theta }(x))^{y}(1-h_{\theta }(x))^{1-y}$$
考虑样本集中的所有样本点，由于每个样本之间相互独立，因此它们的联合分布等于各自边际分布之积，
$$L(\theta )=\coprod_{i=1}^{m}P(y_{i}|x_{i};\theta )=(h_{\theta }(x_{i}))^{y_{i}}(1-h_{\theta }(x_{i}))^{1-y_{i}}$$

这就是我们求解$\theta$需要的似然函数，我们通过他来求解在$\theta$为何值时，$x$取某个值出现某个$y$的概率最大。

对$J(\theta)取对数$，因为$ln(x)$和$x$单调性相同

$$l(\theta )=lnL(\theta )=\sum_{i=1}^{m}(y_{i}lnh_{\theta }(x_{i})+(1-y_{i})ln(1-h_{\theta }(x_{i})))$$

给出损失函数$J(\theta)=-\frac{1}{m}l(\theta)$，对$J(\theta)$求偏导，

理应令求偏导后的表达式等于零求极值，但是无法解析求解，因此用梯度下降法逐渐迭代，找到局部最优解。为什么梯度下降法能够做到呢？

可以看到$\theta$的取值和$J(\theta)$存在着一一对应的关系，让$\theta$沿着$J(\theta)$梯度的方向减小，可以最快速的逼近$J(\theta)$的最小值，但其实往往找到的是极小值，局部最优。

由此可以得到$\theta$的更新公式

三、Logistic回归算法步骤

1. 初始化回归系数$\theta$为1
2. 重复下面步骤直至收敛
   {
   计算整个数据集的梯度
   使用$\alpha$ x gradient更新回归系数$\theta$
   }
3. 返回回归系数$\theta$

四、python实现
现在，我们来解决一个实际问题：分类0/1数字。

在工程目录中有train和test两个文件夹，里面分别存有若干txt文件，每个txt文件用32*32的0/1矩阵表示一个数字0或者1，我们的目标是用train中的文件训练出一个logistic回归模型，也就是得到一组$\theta$，能够用来分类test中的数据。
假设train中有m个文件，将每个文件中的数据拼接成一个1*1024（32*32=1024）的向量，相当于m个样本每个样本有1024维特征，这样训练出的$\theta$也有1024维。同样，test中的n个文件也得到n*1024的矩阵，用之前训练好的$\theta$进行分类。

"""
learned on Wed Feb 15 22:11:00 2017
@author maggie

description:
loadData function
load all the data from the given file, and return as mat

sigmoid function
sigmoid functon in logistic regression

gradAscent function
get the theta vector according to the gradient ascent method

classify function
using the theta vector to classify the test dataset

digitRecognition function
initialize all 

"""

#/usr/bin/python
from numpy import *
from os import listdir

def loadData(dir):
    trainfileList = listdir(dir)
    m = len(trainfileList)
    dataArray = zeros((m, 1024)) #store the data
    labelArray = zeros((m, 1))   #store the label
    for i in range(m):
        tempArray = zeros((1, 1024))
        filename = trainfileList[i]
        fr = open(‘%s/%s‘ %(dir, filename))
        for j in range(32):
            linestr = fr.readline()
            for k in range(32):
                tempArray[0, 32*j+k] = int(linestr[k])
        dataArray[i,:] = tempArray
        filename0 = filename.split(‘.‘)[0]
        label = filename0.split(‘_‘)[0]
        labelArray[i] = int(label)
    return dataArray, labelArray

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

def gradAscent(dataArray, labelArray, alpha, maxCycles):
    dataMat = mat(dataArray)  #size : m x n
    labelMat = mat(labelArray)  #size : m x 1
    m, n = shape(dataMat)
    weigh = ones((n, 1))  #initialize the theta vector
    for i in range(maxCycles):
        h = sigmoid(dataMat * weigh)
        error = labelMat - h #size : m x 1
        weigh = weigh + alpha * dataMat.transpose() * error  #update the theta vector
    return weigh

def classify(testDir, weigh):
    dataArray, labelArray = loadData(testDir)
    dataMat = mat(dataArray)
    labelMat = mat(labelArray)
    h = sigmoid(dataMat * weigh)
    m = len(h)
    error = 0.0
    for i in range(m):
        if int(h[i]) > 0.5:
            print int(labelMat[i]), ‘is classified as : 1‘
            if int(labelMat[i]) != 1:
                error += 1
                print ‘error‘
        else:
            print int(labelMat[i]), ‘is classified as : 0‘
            if int(labelMat[i]) != 0:
                error += 1
                print ‘error‘
    print ‘error rate is ‘, ‘%.4f‘ %(error/m)

def digitRecognition(trainDir, testDir, alpha=0.07, maxCycles=10):
    data, label = loadData(trainDir)
    weigh = gradAscent(data, label, alpha, maxCycles)
    classify(testDir, weigh)
    print weigh


if __name__ == ‘__main__‘:
    digitRecognition(‘train‘,‘test‘)

得到的结果，可以看到错误率为0.018.

(train和test数据以及代码可以点击此处下载）(https://github.com/zjsghww/MachineLearning)