机器学习-logistic回归原理与实现

一、什么是logistic回归？

logistic回归又叫对数几率回归，适合数值型的二值型输出的拟合，它其实是一个分类模型，比如根据患者的医疗数据判断它是否能被治愈。

二、logistic回归数学原理与算法实现

我们考虑1个输入的$n$维数据$x=(x_1,x_2,\dots ,x_n)$，我们对输入数据进行线性加权得到
$$g(x)=w_0+w_1{x}_1+\dots +w_n{x}_n=w^{T}x$$
前面说到，logistic回归用于而分类，假设得到的类别为0或者1，那么可以使用sigmoid函数处理输入，这个函数类似于阶跃函数但是又是连续型函数，看下这个函数长什么样

将$g(x)$作为sigmoid函数的输入，得到
$$f(x)=\frac{1}{1+e^{-g(x)}}$$
上式做个简单的变换
$$\ln \frac{f(x)}{1-f(x)}=w^{T}x$$
将$f(x)$视为类后验概率估计$P(y=1|x)$，则上式可以重写为
$$\ln \frac{P(y=1|x)}{P(y=0|x)}=w^{T}x$$
那么从而可以得到
$$\begin{gathered} P(y=1|x)=f(x) \\ P(y=0|x)=1-f(x) \end{gathered}$$
所以得到一个观测值的概率为
$$P(y|x,w)=[f(x){]}^y\cdot [1-f(x){]}^{1-y}$$
设输入数据为 X = [ x 11 x 12 … x 1 n x 21 x 22 … x 2 n ⋮ ⋮ ⋮ … ⋮ x m 1 x m 2 … x m n ] = { x 1 , x 2 , … , x m } X=\left[\begin{aligned}&{x_{11}} & {x_{12}} & {\dots} & {x_{1 n}} \\ &{x_{21}} &{x_{22}} & {\dots} & {x_{2 n}} \\ {\vdots} & {\vdots} & {\vdots} & {\dots} & {\vdots} \\ &{x_{m 1}} & {x_{m 2}} & {\dots} & {x_{m n}}\end{aligned}\right]=\{x_1,x_2,\ldots,x_m\} X=⎣⎢⎢⎢⎢⎢⎡​⋮​x11​x21​⋮xm1​​x12​x22​⋮xm2​​…………​x1n​x2n​⋮xmn​​⎦⎥⎥⎥⎥⎥⎤​={x1​,x2​,…,xm​}，$y=\left[\begin{array}{r}c{y}_1\\ y_2\\ ⋮\\ y_m\end{array}\right]$，$x_i$表示第$i$个输入数据，上式的似然函数为
$$L(w)=\prod _{i=1}^m{\left[f\left(x_i\right)\right]}^{y_i}{\left[1-f\left(x_i\right)\right]}^{1-y_i}$$
然后我们的目标是求出使这一似然函数的值最大的参数估计，最大似然估计就是求出参数$w_0,w_1,\dots ,w_n$，使得上式取得最大值，对上式两边取对数得到
$$\ln L(w)=\sum _{i=1}^m\left(y_i\ln \left[f\left(x_i\right)\right]+\left(1-y_i\right)\ln \left[1-f\left(x_i\right)\right]\right)$$

2.1 梯度上升法估计参数

我们考虑$\ln L(w)$中间的一部分，对$w_k$求导得到
$$\begin{array}{rl}& {\left(y_i\ln \left[f\left(x_i\right)\right]+\left(1-y_i\right)\ln \left[1-f\left(x_i\right)\right]\right)}^{\prime }\\ & =\frac{y_i}{f\left(x_i\right)}\cdot {\left[f\left(x_i\right)\right]}^{\prime }+\left(1-y_i\right)\cdot \frac{-{\left[f\left(x_i\right)\right]}^{\prime }}{1-f\left(x_i\right)}\\ & =\left[\frac{y_i}{f\left(x_i\right)}-\frac{1-y_i}{1-f\left(x_i\right)}\right]\cdot {\left[f\left(x_i\right)\right]}^{\prime }\\ & =\left(f\left(x_i\right)-y_i\right)g^{\prime }(x)\\ & =x_{ik}\left[f\left(x_i\right)-y_i\right]\end{array}$$
那么
$$\frac{\partial \ln L\left(w_k\right)}{\partial w_k}=\sum _{i=1}^m{x}_{ik}\left[f\left(x_i\right)-y_i\right]=0$$
我们使用**梯度上升法(Gradient ascent method)**求解参数$w$，其迭代公式为
$$w=w+\alpha \nabla \ln L(w)$$
梯度已经在上面计算过了，即
$$\nabla \ln L(w)=\frac{\partial \ln L\left(w\right)}{\partial w}$$
定义损失误差为
$$E=f(x)-y$$
所以我们现在可以得到
$$w=w+\alpha \sum _{i=1}^m{x}_{ik}\left[f\left(x_i\right)-y_i\right]=w+\alpha X^{T}E$$

代码实现

from numpy import *
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties


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(dataMatrix, classLabels):
    dataMatrix = mat(dataMatrix)  # convert to NumPy matrix
    labelMat = mat(classLabels).transpose()  # convert to NumPy matrix
    m, n = shape(dataMatrix)
    alpha = 0.01
    maxCycles = 500
    weights = ones((n, 1))
    weights_list = list()
    for k in range(maxCycles):  # heavy on matrix operations
        h = sigmoid(dataMatrix*weights)  # matrix mult
        error = (labelMat - h)  # vector subtraction
        weights = weights + alpha * dataMatrix.transpose() * error  # matrix mult
        weights_list.append(weights)
    return weights, weights_list                                            

    
def plot_weights(weights_list):
    font = FontProperties(fname=r"/System/Library/Fonts/PingFang.ttc", size=12)
    fig = plt.figure(figsize=(8, 8))
    x = range(len(weights_list))
    w0 = [item[0, 0] for item in weights_list]
    w1 = [item[1, 0] for item in weights_list]
    w2 = [item[2, 0] for item in weights_list]
    plt.subplot(311)
    plt.plot(x, w0, 'r-', label="w0")
    plt.ylabel("w0")
    plt.subplot(312)
    plt.plot(x, w1, 'g-', label="w1")
    plt.ylabel("w1")
    plt.subplot(313)
    plt.plot(x, w2, 'b-', label="w2")
    plt.ylabel("w2")
    plt.xlabel("迭代次数", FontProperties=font)
    plt.show()


def plotBestFit(weights):
    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(figsize=(6, 4))
    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.transpose())
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()


def plot_sigmoid():
    x = arange(-60.0, 60.0, 1)
    y = sigmoid(x)
    fig = plt.figure(figsize=(8, 4))
    plt.xlabel('x')
    plt.ylabel('sigmoid(x)')
    plt.plot(x, y.transpose())
    plt.show()


if __name__ == "__main__":
    data_mat, label_mat = loadDataSet()
    weights, weights_list = gradAscent(data_mat, label_mat)
    #plot_weights(weights_list)
    #plotBestFit(weights)
    #plot_sigmoid()

注意：在上程序的第45行使用的是$w_0+w_1{x}_1+w_2{x}_2=0$等式，而不是$w_0+w_1{x}_1+w_2{x}_2=1$，为什么呢？我们观察sigmoid函数可以发现0值是函数决策边界，而$g(x)=w_0+w_1{x}_1+w_2{x}_2$又是sigmoid函数的输入，所以另$g(x)=0$便可以得到分类边界线。

上述程序的运行结果如下图所示。

迭代次数与最终参数的变化如下图所示。

2.2 改进的随机梯度上升法

梯度上升算法在每次更新回归系数(最优参数)时，都需要遍历整个数据集。假设，我们使用的数据集一共有100个样本，回归参数有 3 个，那么dataMatrix就是一个100×3的矩阵。每次计算h的时候，都要计算dataMatrix×weights这个矩阵乘法运算，要进行100×3次乘法运算和100×2次加法运算。同理，更新回归系数(最优参数)weights时，也需要用到整个数据集，要进行矩阵乘法运算。总而言之，该方法处理100个左右的数据集时尚可，但如果有数十亿样本和成千上万的特征，那么该方法的计算复杂度就太高了。因此，需要对算法进行改进，我们每次更新回归系数(最优参数)的时候，能不能不用所有样本呢？一次只用一个样本点去更新回归系数(最优参数)？这样就可以有效减少计算量了，这种方法就叫做随机梯度上升算法(Stochastic gradient ascent)。

def stocGradAscent(dataMatrix, classLabels, numIter=150):
    dataMatrix = mat(dataMatrix)
    labelMat = mat(classLabels).transpose()
    m,n = np.shape(dataMatrix)                                                
    weights = np.ones((n,1))  
    weights_list = 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(random.uniform(0,len(dataIndex)))              
            h = sigmoid(sum(dataMatrix[randIndex]*weights))                   
            error = classLabels[randIndex] - h                               
            weights = weights + alpha * dataMatrix[randIndex].transpose() * error 
            del(dataIndex[randIndex])  
            weights_list.append(weights)                                      
    return weights, weights_list 

该算法第一个改进之处在于，alpha在每次迭代的时候都会调整，并且，虽然alpha会随着迭代次数不断减小，但永远不会减小到0，因为这里还存在一个常数项。必须这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。如果需要处理的问题是动态变化的，那么可以适当加大上述常数项，来确保新的值获得更大的回归系数。另一点值得注意的是，在降低alpha的函数中，alpha每次减少$1/(j+i)$，其中j是迭代次数，i是样本点的下标。第二个改进的地方在于跟新回归系数(最优参数)时，只使用一个样本点，并且选择的样本点是随机的，每次迭代不使用已经用过的样本点。这样的方法，就有效地减少了计算量，并保证了回归效果。

使用随机梯度上升法获得的分类结果如下图所示

迭代次数与最终参数的变化如下图所示，可以看到，在 8000 次以后，各参数基本趋于稳定，这个过程大约迭代了整个矩阵 80次，相比较于原先的 300 多次，大幅减小了迭代周期。

三、sklearn算法实现

sklearn实现 logistic 回归使用的是linear_model模型，还是使用上述的数据，代码如下

def logistic_lib(dataMatrix, classLabels):
    from sklearn import linear_model
    model_logistic_regression = linear_model.LogisticRegression(solver='liblinear',max_iter=10)
    classifier = model_logistic_regression.fit(dataMatrix, classLabels)
    # 回归系数
    print(model_logistic_regression.coef_)
    # 截距
    print(model_logistic_regression.intercept_)
    # 准确率
    accurcy = classifier.score(dataMatrix, classLabels) * 100
    print(accurcy)

输出结果如下

[[ 2.45317293  0.51690909 -0.71377635]]
[2.45317293]
97.0

参考

[1] CSDN博客-logistic回归原理解析–一步步理解
[2] CSDN博客-logistic回归原理及公式推导
[3] 李航. 统计学习方法, 清华大学出版社