logistic回归详解

逻辑斯谛回归（logistic regression）是统计学习中的经典分类方法，虽然带有回归的字眼，但是该模型是一种分类算法，逻辑斯谛回归是一种线性分类器，针对的是线性可分问题。利用logistic回归进行分类的主要思想是：根据现有的数据对分类边界线建立回归公式，以此进行分类。这里的“回归”一词源于最佳拟合，表示要找到最佳拟合参数集，因此，logistic训练分类器时的做法就是寻找最佳拟合参数，使用的是最优化方法，接下来我们都会讲到。

一、Logistic回归函数：

我们先说一个概念，事件的几率（odds），是指该事件发生的概率与该事件不发生的概率的比值。如果事件发生的概率是p，那么该事件的几率是p/(1-p)。取该事件发生几率的对数，定义为该事件的对数几率（log odds）或logit函数：

事件发生的概率p的取值范围为[0,1]，对于这样的输入，计算出来的几率只能是非负的，而通过取对数，便可以将输出转换到整个实数范围内。
那我们将输出转换到整个实数范围内的目的是什么呢？因为这样，我们就可以将对数几率记为输入特征值的线性表达式 ：

其中，p(y =1|x)是条件概率分布，表示当输入为x时，实例被分为1类的概率，依据此概率我们能得到事件发生的对数几率。但是，我们的初衷是做分类器，简单点说就是通过输入特征来判定该实例属于哪一类别或者属于某一类别的概率。所以我们取logit函数的反函数，令的线性组合为输入，p为输出，经如下推导

公式1就是logistic函数。大家应该对Φ(x)很熟悉，是一个sigmoid函数，类似于阶跃函数的S型生长曲线。

上图给出了sigmoid函数的曲线图。当x为0时，sigmoid函数值为0.5。随着x的增大，对应的sigmoid函数的值将逼近于1；而随着x的减小，sigmoid函数的值将逼近于0。

       那么对于公式1，我们可以这样解释：为了实现logistic回归分类器，我们可以在每个特征上都乘以一个回归系数，然后把所有的结果值相加，将这个总和带入sigmoid函数中。进而得到一个范围在0-1之间的数值。最后设定一个阈值，在大于阈值时判定为1，否则判定为0。以上便是逻辑斯谛回归算法是思想，公式就是分类器的函数形式。

二、最佳回归系数w的确定（极大似然估计+最优化方法）

       上一节我们已经确定了logistic分类函数，有了分类函数，我们输入特征向量就可以得出实例属于某个类别的概率。但这里有个问题，权重w（回归系数）我们是不确定的。正如我们想的那样，我们需要求得最佳的回归系数，从而使得分类器尽可能的精确。

       如何才能获得最佳的回归系数呢？这里就要用到最优化方法。在很多分类器中，都会将预测值与实际值的误差的平方和作为损失函数（代价函数），通过梯度下降算法求得函数的最小值来确定最佳系数

from numpy import *
import matplotlib.pyplot as plt


######################################################读取数据############################################################
def loadDataSet():
    '''数据集的前两个值分别为X1和X2,第三个值是数据对应的类别标签，为了方便计算，把X0的值设置成了1.0'''
    dataMat = []
    labelMat = []
    fr = open('testSet.txt')
    for line in fr.readlines():
        lineArr = line.strip().split()  #按行读取并按空格拆分，每一行为一个列表， ['-0.017612', '14.053064', '0']
        # 为了方便计算，我们将 X0 的值设为 1.0 ，也就是在每一行的开头添加一个 1.0 作为 X0
        # 因为线性回归化式为 H(x) = W0 + W1*X1 + W2*X2，即为 (W0, W1, W2)*(1.0, X1, X2),其中 (W0, W1, W2) 即为所求回归系数 W，所以模拟了一个数据列全为1.0
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])]) #[[1.0, -0.017612, 14.053064], [1.0, -1.395634, 4.662541],........]含有100个元素的列表，每一个元素都为一个列表，包含3个元素
        labelMat.append(int(lineArr[2])) #100个元素（类别标签）的列表
    return dataMat, labelMat



################################################分类器的分类（转换）函数#####################################################
def sigmoid(inX):
    return 1.0 / (1 + exp(-inX))



################################### 梯度上升算法，用来计算出最佳回归系数w(weights)#############################################
#输入数据特征与数据的类别标签，dataMatIn存放的是100*3的矩阵
    """
    第一个参数(dataMatIn)是3维数组，每列代表每个不同特征，每行代表每个训练样本
    第二个参数(labelMat)是类别标签，1*100的行向量，为便于计算，将其转换为列向量，即进行转置，并赋值给labelMat
    转化矩阵[[0,1,0,1,0,1.....]]为[[0],[1],[0].....]
    """
def gradAscent(dataMatIn, classLabels):
    dataMatrix = mat(dataMatIn) #转为100*3的矩阵 [ 1.0000000e+00 -1.7612000e-02  1.4053064e+01]
    labelMat = mat(classLabels).transpose()  #转为100*1的矩阵，transpose() 行列转置函数,将行向量转化为列向量   =>  矩阵的转置
    # m->数据量、样本数100    n->特征数3
    m, n = shape(dataMatrix)# shape()函数是numpy.core.fromnumeric中的函数，它的功能是查看矩阵或者数组的维数
    alpha = 0.001  # 步长，向函数增长最快的方向的移动量，即学习率
    maxCycles = 500# 迭代次数
    weights = ones((n, 1))  # 转为3*1的矩阵，元素为1的矩阵赋给weihts，即回归系数初始化为1


    # 循环 maxCycles次, 每次都沿梯度向真实值 labelMat 靠拢
    for k in range(maxCycles):
        # 求当前sigmoid函数的值
        h = sigmoid(dataMatrix * weights)  # 矩阵相乘 包含了300次的乘积 [100*3] * [3*1] = [100*1]
        error = (labelMat - h)  # 向量减法，计算真实类别与预测类别的差值，h是一个列向量，列向量的元素个数等于样本数，即为100
        weights = weights + alpha * dataMatrix.transpose() * error  # 矩阵相乘，dataMatrix.transpose()* error 就是梯度f(w)，按照该差值的方向调整回归系数
    return weights




###################################################分析数据：画出决策边界####################################################
# 画出数据集和Logistic回归最佳拟合直线
def plotBestFit(weights):
    dataMat, labelMat = loadDataSet()
    dataArr = array(dataMat)
    n = shape(dataArr)[0]   # n->数据量、样本数
    # xcord1,ycord1代表正例特征1,正例特征2
    # xcord2,ycord2代表负例特征1,负例特征2
    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()

    # “111”表示“1×1网格，第一子图”，“234”表示“2×3网格，第四子图”。
    ax = fig.add_subplot(111)
    ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
    ax.scatter(xcord2, ycord2, s=30, c='green')

    # 画线
    # 设定边界直线x和y的值,并且以0.1为步长从-3.0到3.0切分。
    x = arange(-3.0, 3.0, 0.1)#[-3.000000000 -2.90000000 -2.80000000 -2.70000000
    """
    y的由来？
    首先理论上是这个样子的:
    dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
    w0*x0+w1*x1+w2*x2=f(x)
    x0最开始就设置为1， x2就是我们画图的y值。0是两个类别的分界处
    所以： w0+w1*x+w2*y=0 => y = (-w0-w1*x)/w2   
    """
    y = (-weights[0] - weights[1] * x) / weights[2]  #最佳拟合直线

    ax.plot(x, y)
    plt.xlabel('X1');
    plt.ylabel('X2');
    plt.show()



# 输出运用梯度上升优化算法后得到的最理想的回归系数的值
if __name__ =='__main__':
    dataArr,labelMat=loadDataSet()
    weights=gradAscent(dataArr,labelMat)
    plotBestFit(weights.getA())