R语言画单因素和多因素逻辑回归森林图 多因素logistic回归森林图

一、Logistic回归

介绍Logistic回归模型之前，首先看一个并不常见的概率分布，即Logistic分布。设X是连续随机变量，X服从Logistic分布是指X具有下列的分布函数和密度函数：

式中，$µ$为位置参数，$γ$> 0为形状参数。Logistic的分布的密度函数f(x)和分布函数F(x)的图形如下图所示。其中分布函数属于Logistic函数，其图形为一条S形曲线。该曲线以点($µ$,1/2)为中心对称，即满足：

曲线在中心附近增长较快，在两端增长速度较慢。形状参数$γ$的值越小，曲线在中心附近增长得越快。

二、Logistic回归模型

线性回归的应用场合大多是回归分析，一般不用在分类问题上。原因可以概括为以下两个：

• 1）回归模型是连续型模型，即预测出的值都是连续值（实数值），非离散值；
• 2）预测结果受样本噪声的影响比较大。

2.1 LR模型表达式

LR模型表达式为参数化的Logistic函数（默认参数$µ$= 0,$γ$= 1）

其中h_θ(x)作为事件结果y=1的概率取值。这里x ∈ R_n+1, y ∈ {1, 0}, θ ∈ R_n+1是权值向量。其中权值向量w中包含偏置项，即w = (w₀, w₁ , · · · , w_n) ， x = (1, x₁, x₂, · · · , x_n)

2.2 理解LR模型

2.2.1 对数几率

一个事件发生的几率（odds）是指改事件发生的概率与该事件不发生的概率的比值。如果事件发生的概率是$p$，那么该事件的几率为$p/1-p$，该事件的对数几率（log odds）或者logit函数是：

对LR而言，根据模型表达式可以得到：

即在LR模型中，输出y=1的对数几率是输入x的线性函数。或者说输出y=1的对数几率是由输入x的线性函数表示的模型，即LR模型。

2.2.1 函数映射

除了从对数几率的角度理解LR外，从函数映射也可以理理解LR模型。

考虑对输入实例x进行分类的线性表达式θ^T ，其值域为实数域。通过LR模型表达式可以将线性函数θ^Tx的结果映射到(0,1)区间，取值表示为结果为1的概率（在二分类场景中）。

线性函数的越接近于无穷大，概率值就越接近 1；反之，其值越接近于负无穷，概率就越接近于0，这样的模型就是LR模型。

LR本质上还是线性回归，知识特征到结果的映射过程中加一层函数映射（即sigmoid函数），即先把特征线性求和，然后使用sigmoid函数将线性和约束至（0,1）之间，结果值用于二分或者回归预测。

2.2.3 概率解释

LR模型多用于解决二分类问题，如广告是否被点击（是/否）、商品是否被购买（是/否）等互联网领域中常见的应用场景。但是实际场景中，我们又不把它处理理成“绝对的”分类问题，而是用其预测值作为事件发生的概率。

这里从事件、变量以及结果的角度给予解释。

我们所能拿到的训练数据统称为观测样本。问题：样本是如何生成的？

一个样本可以理解为发生一次事件、样本生成的过程即为事件发生的过程。对于0/1分类问题来讲，产生的结果有两种可能，符合伯努利试验的概率假设。因此、我们可以说样本的生成过程即为伯努利试验过程、产生的结果符合伯努利分布、这里我们假设结果为1的概率为h_θ(x)，那么对于第i个样本概率公式表示如下：

将上面两个公式合并在一起，可得到第i个样本正确预测的概率：

上式是对一个样本进行建模的数据表达。对于所有的样本，假设每条样本生成过程独立，在整个样本空间中（N个样本）的概率分布（即似然函数）为：

通过极大似然估计（Maximum Likelihood Evaluation，简称MLE）方法求概率参数。具体地，第三节给出了通过随机梯度下降法（SGD）求参数。

三、模型参数估计

3.1 Sigmoid函数

上图所示即为sigmoid函数，它的输入范围为$−∞→ +∞$，而值域刚好为(0, 1)，正好满足概率分布为(0, 1)的要求。用概率去描述分类器，自然要比阈值要来的方便。而且它是一个单调上升的函数，具有良好的连续性，不存在不连续点。此外非常重要的，sigmoid函数求导后为：

以下的推导中会用到，带来了很大的便利。

3.2 参数估计推导

上一节的公式不仅可以理解为在已观测的样本空间中的概率分布表达式。如果从统计学的角度可以理解为参数$θ$似然性的函数表达式（即似然函数表达式）。就是利用已知的样本分布，找到最有可能（即最大概率）导致这种分布的参数值；或者说什么样的参数才能使我们观测到目前这组数据的概率最大。参数在整个样本空间的似然函数可表示为：

为了方便参数求解，对这个公式取对数，可得对数似然函数：

最大化对数似然函数其实就是最小化交叉熵误差（Cross Entropy Error）。先不考虑累加和，我们针对每一个参数w_j求偏导：

最后，通过扫描样本，迭代下述公式可求得参数：

其中$a$表示学习率，又称学习步长。此外还有Batch GD，共轭梯度，拟牛顿法（LBFGS），

ADMM分布学习算法等都可以用来求解参数。另作优化算法一章进行补充。以上的推导是LR模型的核心部分，在机器学习相关面试中，LR模型公式推导可能是考察频次最高的一个点。要将其熟练推导。

3.3 分类边界

知道如何求解参数后，我们看一下模型得到的最后结果是什么样的。

假如我们的决策函数为：

选择0.5作为阙值是一个一般的做法，实际应用时特定的情况可以选择不同的阙值，如果对正例的判别准确率要求高，可以选择阙值大一些，对正例的召回要求高，则可以选择阙值小一些。很容易看出来$θ$^Tx>0时，y=1，否者y=0。$θ$^Tx=0是模型隐含的分类平面（在高维空间中，我们说的是超平面）。所以说逻辑回归本质上是一个线性模型，但是这不意味着只有线性可分的数据能通过LR模型求解，实际上，我们可以通过特征变换的方式把低维空间转换到高维空间（kernel trick），而在低维空间不可分的数据，到高维空间中线性可分的几率会高一些。下面两个图的对比说明了线性可分类曲线和非线性曲线（通过特征映射）。

左图是一个线性可分的数据集，右图在原始空间中线性不可分，但是在特征转换[x₁, x₂] => [x₁,x₂,x²₁ , x²₂ , x₁x₂]后的空间是线性可分的，对应的原始空间中分类边界为一条椭圆曲线。不过，通常使用的kernel都是隐式的，也就是找不到显式地把数据从低维映射到高维的函数，而只能计算高维空间中的数据点的内积。在这种情况下，logistic regression模型就不能再表示成w^Tx+b的形式（原始形式primal form），而只能表示成$∑$_ia_i<x_i, x> + b的形式（对偶形式dualform）。忽略b，则原始形式的模型参数只有w，只需要一个数据点那么多的存储量；而对偶形式的模型不仅需要存储各个a_i，还要存储训练数据x_i本身，这个存储量就大了。

SVM也具有原始形式和对偶形式，相比之下SVM的对偶形式是稀疏的，即只有支持向量的a_i非零，才需要存储相应的x_i，所以，在非线性可分的情况下，SVM用的更多。

3.4 示例代码

• 几种参数推导代码比较

def loadDataSet():
    dataMat=[]
    labels=[]
    with open('testSet.txt','r') as f:
        for line in f.readlines():
            linArr=line.strip().split()
            dataMat.append([1.0,float(linArr[0]),float(linArr[1])])
            labels.append(int(linArr[2]))
    return dataMat,labels
def sigmoid(X): # Sigmoid函数
    return 1.0/(1+np.exp(-X))
# 梯度上升
def gradAscent(datamain,classLabels):
    dataMatrix=np.mat(datamain)
    labelMat=np.mat(classLabels).transpose() # 矩阵转换
    m,n=np.shape(dataMatrix)
    alpha=0.001  # 学习率
    maxCycles=500
    weights=np.ones((n,1))  # 初始化权重
    for k in range(maxCycles):
        h=sigmoid(dataMatrix*weights)
        error=labelMat-h
        weights=weights+alpha*(dataMatrix.transpose()*error) #更新权重
    return weights
datasets,labels=loadDataSet()
weights=gradAscent(datasets,labels)

# 画出决策边界
import matplotlib.pyplot as plt
import matplotlib
def ploBestFit(wei):
    weights=wei
    dataMat,labelMat=loadDataSet()
    dataArr=np.array(dataMat)
    n=np.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=np.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()
ploBestFit(weights.getA())

在全部数据集上梯度上升得到的分类情况：

# 随机梯度上升，占用更少的计算资源
def stocGradAscent(dataMatrix,classLabels):
    m,n=np.shape(dataMatrix)
    learning_rate=0.01
    weights=np.ones(n)
    for i in range(m):
        h=sigmoid(np.sum(dataMatrix[i]*weights)) # 选取一个数据而非选取全部数据
        error=classLabels[i]-h
        weights=weights+learning_rate*error*dataMatrix[i]
    return weights
weights=stocGradAscent(np.array(datasets),labels)
ploBestFit(weights)

随机梯度上升算法效果并不是很好，对随机梯度上升算发法进行优化

def stocGradeAscent1(dataMatrix,classLabels,num=150):
    m,n=np.shape(dataMatrix)
    weights=np.ones(n)
    for j in range(num): # 迭代次数
        dataIndex=list(range(m))
        for i in range(m):
            lr=4/(1.0+j+i)+0.01
            randIndex=int(np.random.uniform(0,len(dataIndex))) # 随机选取数据集数据
            h=sigmoid(np.sum(dataMatrix[randIndex]*weights))
            error=classLabels[randIndex]-h
            weights=weights+lr*error*dataMatrix[randIndex]
            del(dataIndex[randIndex]) # 删除已选数据
    return weights
weights=stocGradeAscent1(np.array(datasets),labels,1000)
ploBestFit(weights)

可见，优化后的随机梯度上升算法和梯度上升算法效果相当

• “从疝气病症预测病马的死亡率” 示例代码

def classifyVector(X,weights):
    prob=sigmoid(np.sum(X*weights))
    if prob>0.5:
        return 1.0
    return 0.0
def colicTest():
    ftr=open("horseColicTraining.txt",'r')
    fte=open("horseColicTest.txt",'r')
    train=[];train_labels=[]
    for line in ftr.readlines():
        currline=line.strip().split('\t')
        linArr=[]
        for i in range(21):
            linArr.append(float(currline[i]))
        train.append(linArr)
        train_labels.append(float(currline[21]))
    trainWeights=stocGradeAscent1(np.array(train),train_labels)
    errorCount=0;numTestVecx=0.0
    for line in fte.readlines():
        currline=line.strip().split('\t')
        linArr=[]
        for i in range(21):
            linArr.append(float(currline[i]))
        if classifyVector(np.array(linArr),np.array(trainWeights))!=int(currline[21]):
            errorCount+=1
        numTestVecx+=1
    errorRate=errorCount/numTestVecx
    print("the error rate of this test is %f"%(errorRate))
    fte.close()
    ftr.close()
    return errorRate
def mulTest():
    errorSum=0.0
    for i in range(10):
        errorSum+=colicTest()
    print("the avg error rate is %f"%(errorSum/10.0))     
mulTest()   
输出：
the error rate of this test is 0.358209
the error rate of this test is 0.373134
the error rate of this test is 0.402985
the error rate of this test is 0.656716
the error rate of this test is 0.298507
the error rate of this test is 0.373134
the error rate of this test is 0.343284
the error rate of this test is 0.313433
the error rate of this test is 0.253731
the error rate of this test is 0.388060
the avg error rate is 0.376119

四、延伸

4.1 生成模型与判别模型

逻辑回归是一种判别模型，表现为直接对条件概率P(y|x)建模，而不不关心背后的数据分布P(x, y)

。而高斯贝叶斯（Gaussian Naive Bayes）是⼀一种生成模型，先对数据的联合分布建模，再通过

贝叶斯公式来计算属于各个类别的后验概率，即

通常假设P(x|y)是高斯分布，P(y)是多项式分布，相应的参数可以通过最大似然估计得到。如果

我们考虑二分类问题，通过简单的变化可以得到：

如果σ₁= σ₀，二次项会抵消，我们得到一个简单的线性关系：

上式进一步可以得到：

可以看到，这个概率和逻辑回归中的形式是一样的，这种情况下高斯贝叶斯和LR会学习到同一个模型。实际上，在更一般的假设（P(x|y)的分布属于指数分布族）下，我们都可以得到类似的结 论。

4.2 多分类

如果y不是在[0, 1]中取值，而是在K个类别中取值，这时问题就变为一个多分类问题。有两种方式可以出处理该类问题：一种是我们对每个类别训练一个二元分类器（One-vs-all），当K个类别不是互斥的时候，比如用户会购买哪种品类，这种方法是合适的。如果K个类别是互斥的，即 y = i的时候意味着 y 不能取其他的值，比如用户的年龄段，这种情况下 Softmax 回归更合适⼀一些。Softmax 回归是直接对逻辑回归在多分类的推广，相应的模型也可以叫做多元逻辑回归

（Multinomial Logistic Regression）。模型通过 softmax 函数来对概率建模，具体形式如下：

而决策函数为

对于的损失函数为

类似的，我们也可以通过梯度下降或其他高阶方法来求解该问题，这里不再赘述。

4.3 应用

这里以预测用户对品类的购买偏好为例，介绍一下美团是如何用逻辑回归解决工作中问题的。该问题可以转换为预测用户在未来某个时间段是否会购买某个品类，如果把会购买标记为1，不会购买标记为0，就转换为一个二分类问题。我们用到的特征包括用户在美团的浏览，购买等历史信息，见下表

其中提取的特征的时间跨度为30天，标签为2天。生成的训练数据大约在7000万量级（美团一个月有过行为的用户），我们人工把相似的小品类聚合起来，最后有18个较为典型的品类集合。如果用户在给定的时间内购买某一品类集合，就作为正例。有了训练数据后，使用Spark版的LR算法对每个品类训练一个⼆二分类模型，迭代次数设为100次的话模型训练需要40分钟左右，平均每个模型2分钟，测试集上的AUC也大多在0.8以上。训练好的模型会保存下来，用于预测在各个品类上的购买概率。预测的结果则会用于推荐等场景。

由于不同品类之间正负例分布不不同，有些品类正负例分布很不均衡，我们还尝试了不同的采样方法，最终目标是提高下单率等线上指标。经过一些参数调优，品类偏好特征为推荐和排序带来了超过1%的下单率提升。

此外，由于LR模型的简单高效，易于实现，可以为后续模型优化提供一个不错的baseline，我们在排序等服务中也使用了了LR模型。

逻辑回归的数学模型和求解都相对比较简洁，实现相对简单。通过对特征做离散化和其他映射，逻辑回归也可以处理非线性问题，是一个非常强大的分类器。因此在实际应用中，当我们能够拿到许多低层次的特征时，可以考虑使用逻辑回归来解决我们的问题。

4.4 LR与SVM

两种方法都是常见的分类算法，从目标函数来看，区别在于逻辑回归采用的是logistical loss，svm采用的是hinge loss。这两个损失函数的目的都是增加对分类影响较大的数据点的权重，减少与分类关系较小的数据点的权重。SVM的处理方法是只考虑support vectors，也就是和分类最相关的少数点，去学习分类器。而逻辑回归通过非线性映射，大大减小了离分类平面较远的点的权重，相对提升了与分类最相关的数据点的权重。两者的根本目的都是一样的。此外，根据需要，两个方法都可以增加不同的正则化项，如l1,l2等等。所以在很多实验中，两种算法的结果是很接近的。但是逻辑回归相对来说模型更简单，好理解，实现起来，特别是大规模线性分类时比较方便便。而SVM的理解和优化相对来说复杂一些。但是SVM的理论基础更加牢固，有一套结构化风险最小化的理论基础，虽然一般使用的人不太会去关注。还有很重要的⼀一点，SVM转化为对偶问题后，分类只需要计算与少数几个支持向量的距离，这个在进行复杂核函数计算时优势很明显，能够大大简化模型和计算量。

两者对异常的敏感度也不一样。同样的线性分类情况下，如果异常点较多的话，无法剔除，首先LR中每个样本都是有贡献的，最大似然后会自动压制异常的贡献，SVM+软间隔对异常还是比较敏感，因为其训练只需要支持向量，有效样本来就不高，一旦被干扰，预测结果难以预料。