基本原理

找个合适的决策函数即$h$函数（hypothesis），用来预测对输入数据的判断结果。这个过程较为关键，你要对数据有一定的了解或者知道预测函数的“大概”形式，比如线性函数还是非线性函数。
构造损失函数，该函数表示预估的结果（$h$）与实际（$y$）的差异。同时综合考虑所有样本点的“损失”，将Cost求平均，记为$J(\theta)$函数。
通过梯度下降法求$J(\theta)$函数的最小值，使得是预测与真实值之间的偏差达到最小。

求解过程

为什么使用Logistic函数

现实生活中，二分类问题最常遇到：判断一笔交易是欺诈还是非欺诈；一封邮件是正常邮件还是垃圾邮件等，此时目标变量即类别标签为离散值，通常记作$\{1,0\}$，分别表示正类和负类。而最简单的分类器就是线性分类器，通过找到一个超平面(Hyperplan)，将两个不同的类区分开。对于这个超平面，可以用如下的线性函数表示：

$w^T+b=0$

当 $w^T+b<0$ 时，即 $\sum\limits_{i=1}^{n} w_ix_i<-b$ ，此时的样本类别 $y=0$

当 $w^T+b>0$ 时，即 $\sum\limits_{i=1}^{n} w_ix_i>-b$ ，此时的样本类别 $y=1$

这不就是我们熟悉的感知机么，即划定一个阈值（-b），通过比较样本与阈值的大小来进行分类。

虽然这个模型简单直观，但存在几个问题。一，假设阈值为10，现在有一个样本进来，最后计算得出10.01，是归为1好还是0好呢？二，这个函数在阈值这个点上不连续，求导很不方便。

因此，逻辑回归引入了Logistic函数，其数学函数是：

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

对应的函数曲线如下图：

通过Logistic函数的图像，发现优点有：

它将线性回归输出的很大范围的数，例如负无穷到正无穷，压缩到(0,1)之间，且远离0的地方函数值会很快接近0/1，这个性质使我们能够以概率的方式来解释
它是一个单调上升函数，具有良好的连续性，不存在不连续点

决策函数

我们通过Logistic函数，构造了一个预测函数$h$为：

$h_\theta (x)=g(\theta^{T}x)=\frac{1}{1+e^{-\theta^{T}x}}$

$h_\theta (x)$函数的值有特殊含义，它表示结果取1的概率。因此对于输入X分类结果为类别1与类别0的概率分别是

$P(y=1|x;\theta )=h_\theta(x)$ $P(y=0|x;\theta )=1-h_\theta(x)$

我们可以把上面两个式子合并成一个，如下：

$P(y|x;\theta )=h_\theta(x)^{y}(1-h_\theta(x))^{(1-y)}$

对于Logistic Regression算法来讲，我们需要求解分割超平面中的参数，即权重矩阵$W$与偏置向量$b$。

损失函数

我们用极大似然法对参数进行估计，即找到一组参数，使得在这组参数下，数据的似然度越大。假设训练集有$n$个独立的训练样本

$\{(x^{(1)},y^{(1)}) ,(x^{(2)},y^{(2)}),...,(x^{(n)},y^{(n)})\}$

每一个样本被观测到的概率为

$P(y=1|x)^{y}(1-P(y=1|x))^{1-y}$

那么我们的整个样本集，也就是$n$个独立的样本出现的似然函数为

$L(\theta)=\prod\limits P(y=1|x)^{y}(1-P(y=1|x))^{1-y}$

取对数可得到对数似然度

$l(\theta)=\sum\limits ylog(h_\theta(x))+(1-y)log(1-h_\theta(x))$

另一方面，我们更常遇到的是损失函数的概念，主要有0-1损失，log损失等。其中log损失在单个点上的定义为

$-log(P(Y|X))=-ylog(p(y|x))-(1-y)(1-log(p(y|x)))$

若取在整个数据集上的平均log损失，得到

$J(\theta)=-\frac{1}{N} l(\theta)$

即在Logistic Regression中，最大化似然函数与最小化log损失函数是等价的。

损失函数的优化方法

Logistic Regression主要使用梯度下降法（Gradient Descent）对损失函数进行优化。它有许多优点，例如只要求解一阶导数，计算成本小，使得能在大规模数据集上应用。它是一种迭代求解的方法，通过在每一步沿着梯度方向，不断调整参数的值来逼近最优值。步骤如下：

选择下降方向即梯度方向 $\nabla {J(\theta)}$
选择步长，更新参数 $\theta_j = \theta_j - \alpha \frac{\partial}{\partial \theta_j}J(\theta)$
重复上述步骤，直至满足终止条件

其中损失函数$J(\theta)$在方向 $\theta_j$上的导数为：

$\frac{\partial{J}}{\partial{\theta_j}} = \frac{1}{n}\sum_{i=1}^{n} (g(\theta^T x^{(i)}) - y^{(i)})x_j^{(i)}$

或者写成矩阵的形式：

$\frac{\partial}{\partial\theta}J(\theta) = X^T(h_{\theta}(X) - Y )$

【矩阵法推导梯度】
由于代数法推导过于繁琐，这里我们用矩阵法推导梯度：

$J(\theta) = -Y\bullet logh_{\theta}(X) - (E-Y)\bullet log(E-h_{\theta}(X))$

用$J(\theta)$对$\theta$向量求导，得：

$\frac{\partial}{\partial\theta}J(\theta) = -Y \bullet X^T\frac{1}{h_{\theta}(X)}h_{\theta}(X)(1-h_{\theta}(X)) + (E-Y)\bullet X^T\frac{1}{1-h_{\theta}(X)}h_{\theta}(X)(1-h_{\theta}(X))$ $\frac{\partial}{\partial\theta}J(\theta) = -Y \bullet X^T(1-h_{\theta}(X)) + (E-Y)\bullet X^T h_{\theta}(X)$ $\frac{\partial}{\partial\theta}J(\theta) = X^T (h_{\theta}(X)-Y )$

其中，除了矩阵求导的链式法则，我们还用到了下面三个公式：

$\frac{\partial}{\partial X}logX = 1/X$ $\frac{\partial}{\partial z}g(z) = g(z)(1-g(z))$ $\frac{\partial}{\partial\theta}X\theta = X^T$

算法实现

这里主要参考的是《机器学习实战》中的代码。

训练样本100条，每一列分别表示$X0、X1、Y$，形式如下：

-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

加载代码如下：

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])]) #X0设为1.0，构成拓充后的输入向量
        labelMat.append(int(lineArr[2]))
    return dataMat,labelMat

可视化的代码：

def plotBestFit(weights):
    """
    画出数据集和逻辑斯谛最佳回归直线
    """
    import matplotlib.pyplot as plt
    dataMat,labelMat=loadDataSet()
    dataArr = array(dataMat)
    n = shape(dataArr)[0] 	#行数
    xcord1 = []; ycord1 = []
    xcord2 = []; ycord2 = []
    for i in range(n):
        if int(labelMat[i])== 1:	#标签为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')
    if weights is not None:
        x = arange(-3.0, 3.0, 0.1)
        y = (-weights[0]-weights[1]*x)/weights[2]   #决策边界：令w0*x0 + w1*x1 + w2*x2 = 0，其中x0=1，解出x1和x2的关系
        ax.plot(x, y)                               
    plt.xlabel('X1'); plt.ylabel('X2');
    plt.show()

梯度上升算法

def sigmoid(inX):
    return 1.0/(1+exp(-inX))
 
def gradAscent(dataMatIn, classLabels):
    """
    逻辑斯谛回归梯度上升优化算法
    """
    dataMatrix = mat(dataMatIn)             #转换为 NumPy 矩阵数据类型
    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)     # 100*3 3*1 --> 100*1
        error = (labelMat - h)       
        # sigma (样本点的输出与真实的差异)*该样本点在方向j上的取值       
        weights += alpha * dataMatrix.transpose() * error  #3*100 100*1 --> 3*1
    return weights

调用方法：

1
2
3

dataArr, labelMat = loadDataSet()
weights = gradAscent(dataArr, labelMat)
plotBestFit(weights)

效果示意：

做成动画展现出来

def gradAscent(dataMatIn, classLabels, history_weight):
    dataMatrix = mat(dataMatIn)             #转换为 NumPy 矩阵数据类型
    labelMat = mat(classLabels).transpose() #转换为 NumPy 矩阵数据类型
    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 += alpha * dataMatrix.transpose() * error  #矩阵内积
        history_weight.append(copy(weights))
    return weights
def draw_line(weights):
    x = arange(-5.0, 5.0, 0.1)
    y = (-weights[0]-weights[1]*x)/weights[2]   #令w0*x0 + w1*x1 + w2*x2 = 0，其中x0=1，解出x1和x2的关系
    line.set_data(x, y)
    return line,
 
def init():
    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])
    ax = fig.add_subplot(111)
    ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
    ax.scatter(xcord2, ycord2, s=30, c='green')
    plt.xlabel('X1'); plt.ylabel('X2');
 
    return draw_line(ones((n,1)))
 
def animate(i):
    return draw_line(history_weight[i])
# 调用
history_weight = []
dataMat,labelMat=loadDataSet()
gradAscent(dataMat, labelMat, history_weight)
fig = plt.figure()
currentAxis = plt.gca()
ax = fig.add_subplot(111)
line, = ax.plot([], [], 'b', lw=2)
anim = animation.FuncAnimation(fig, animate, init_func=init, frames=len(history_weight), interval=30, repeat=False,blit=True)
plt.show()

可视化

随机梯度上升算法

梯度下降法每次更新权值向量$W$是需要遍历所有样本点，计算代价高。而随机梯度下降采用的方法是一次仅用一个样本点数据来计算梯度，从而更新权值向量$W$。

算法：

def stocGradAscent0(dataMatrix, classLabels, history_weight):
    """
    随机梯度上升算法
    """
    dataMatrix = array(dataMatrix)
    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 + dataMatrix[i] * alpha * error
        history_weight.append(copy(weights))
    return weights

其中$h$与$error$均为数值，没有矩阵运算。

可视化

虽然每次迭代的复杂度小的多（仅使用一个样本点），但同时使得最后的分类效果不如梯度上升算法，于是我们对随机梯度算法进行改进。

改进随机梯度上升

每次迭代式减小步长$\alpha$ ，使算法参数尽可能收敛

def stocGradAscent1(dataMat, classLabels, history_weight, numIter=150):
    """
    改进的随机梯度上升算法
    """
    dataMatrix = array(dataMat)
    m,n = shape(dataMatrix)
    weights = ones(n)                                           #初始化为单位矩阵
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.0001                          #步长递减，但是由于常数存在，所以不会变成0
            randIndex = int(random.uniform(0,len(dataIndex)))   #总算是随机了
            #print(j,i,randIndex,dataIndex,dataIndex[randIndex])
            h = sigmoid(sum(dataMatrix[randIndex]*weights))
            error = classLabels[randIndex] - h
            weights = weights + alpha * error * dataMatrix[randIndex]
            history_weight.append(copy(weights))
            del(dataIndex[randIndex])                           #删除这个样本，以后就不会选到了
    return weights

可视化

总结

优点

形式简单，模型的可解释性非常好。从特征的权重可以看到不同的特征对最后结果的影响。
训练速度较快，资源占用小。分类的时候，计算量仅仅只和特征的数目相关，且内存只需要存储各个维度的特征值。
方便输出结果调整。逻辑回归可以很方便的得到最后的分类结果，因为输出的是每个样本的概率分数，我们可以很容易的对这些概率分数进行cutoff，也就是划分阈值

缺点

形式非常的简单（非常类似线性模型），很难去拟合数据的真实分布。
LR对样本分布敏感，很难处理数据不平衡的问题，所以要注意样本的平衡性（通常需要过采样）。
对于非线性切分，处理比较麻烦。必须在原始数据上做一些非线性变换，比如把X做个平方项，X1*X2等。
逻辑回归本身无法筛选特征。有时候，我们会用gbdt来筛选特征，然后再上逻辑回归。

参考资料

寒小阳 - 机器学习系列(1)_逻辑回归初步
 美团点评 - Logistic Regression 模型简介
 logistic回归详解(二）：损失函数（cost function）详解
 逻辑回归的常见面试点总结